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

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


Introduction
Characterizing the flow of fluids and the associated transport of dissolved chemicals through subsurface fractured porous media is a critical step for many industrial and civil engineering endeavors, including the geological sequestration CO 2 [29,33], drinking water aquifer management [39,55], the environmental restoration of contaminated fractured media [45,46,54], hydrocarbon extraction [26,44], and the permanent disposal of spent nuclear fuel and high-level radioactive waste [20,34].A cornerstone of simulating flow, transport, or heat transfer in fractured porous media is the adopted system representation, a computational mesh and solving the governing equations of interest on said mesh.Recent advances in discrete representations of fractures, such as discrete fracture network (DFN) and discrete fracture matrix (DFM) models, have vastly improved our capabilities to model physical and chemical processes on realistic fracture networks, but are still limited in several regards.DFN models only represent the fractures, neglecting important contributions from the surrounding rock (i.e., matrix), which can be critical in a variety of applications, including enhanced geothermal and unconventional reservoirs.DFM models are unequivocally the most accurate representation of the true fractured system [6], but are still limited by difficulties in meshing algorithms and solving the governing equations, which is complicated by the coupling between fracture and matrix domains.While DFN models have been used to model tens of thousands of fractures, DFM models have only been used to model hundreds [31].Consequently, we frequently employ continuum models, also known as equivalent continuum/continuous porous medium models (ECPM), where the fracture network properties are upscaled into equivalent properties of a porous medium when the fracture number is large and the matrix processes cannot be neglected.They are typically used due to their ease of implementation and wide compatibility with existing multi-physics codes.However, they are known to smear out important length scales and it is not fully understood whether or how their accuracy depends on the underlying fracture network properties and the resolution of the computational mesh in relation.Previous work by Jackson et al. [32] showed that ECPMs are able to represent the underlying DFNs, in terms of reproducing the effective permeability.However, they only considered two different well-connected DFNs of relatively high density and no flow contribution from the matrix.More recently, Kottwitz et al. [37] systematically analyzed the relationship between the resolution of the ECPM mesh and the effective permeability in steady state flow simulations.They found monotonic convergence of the effective permeability as the mesh size was decreased and suggest minimum and maximum cell sizes related to fracture sizes needed in continuum meshes to preserve the underlying network characteristics.The networks considered in their study were networks with higher fracture intensity/density.Sweeney et al. [52] also showed that upscaling for a single fracture is accurate regardless of mesh size.Taken together with prior studies, this suggests that changing the topology of the effective network is responsible for differences when considering different mesh resolutions for the same fracture network.For instance, at a coarse resolution, properties of two fractures may be upscaled to the same cell in the continuum mesh, even if they do not intersect in the underlying network, in effect creating a "false connection" in the network.These types of changes in topology have never been systematically investigated, but are critical to understanding whether we can actually use continuum models with confidence for a spectrum of DFNs, not just ones with high intensity/density.In principle, if two fractures are close to one another and do not intersect, then the meshing generation parameters need to be selected so that the background element sizes are smaller than the distance between the fractures to eliminate a false connection.Within a stochastically generated fracture network, attempts to resolve these distances, which will quickly become below the size resolvable with finite precision arithmetic, are futile regardless of the heroic efforts put forth.Thus, false connections are unavoidable when upscaling a stochastically generated network.Moreover, multiple false connections can occur within a single control volume.Deviation of connectivity/topology of the continuum mesh from the fracture network could influence the local and global flow and transport within the network.Specifically, global changes in topology could occur when false connections result in a continuum mesh with a connected pathway of fracture cells between flow boundaries, i.e., percolates, when the underlying DFN does not percolate.However, characterizing the impact of these false connections has yet to be systematically performed, and and general guidelines have yet to be established.
In this work, we have assessed the impacts of so-called false connections on the accuracy and convergence of flow and transport in fractured systems represented by an ECPM model.We found that selecting a mesh resolution such that the global topology of the upscaled mesh matches the fracture network is essential.If the upscaled mesh has a connected pathway of fracture (higher permeability) cells but the fracture network does not, then the estimates for effective permeability and solute breakthrough will be incorrect.Local false connections between fractures due to a coarse mesh result in more solute dispersion, but to a smaller degree than if there is mismatch in global connectivity, which is consistent with prior work.To this end, we provide some principle-based guidelines.Foremost, macro-scale connectivity is an essential property of the network to be represented in an upscaled continuum.Likewise, while false connections cannot be eliminated entirely, they can be managed by choosing appropriate mesh resolution and refinement for a given network to obtain more accuracte estimations of effective permeability and breakthrough curve behavior.These general guidelines are refined both in terms of the difference between matrix and fracture permeabilities and network densities.

Numerical Simulations
Fracture network generation, upscaling, flow and transport simulations are performed using the dfnWorks computational suite [24].dfnWorks combines the feature rejection algorithm for meshing fram [23,38], the LaGriT meshing toolbox [40], and the parallelized subsurface flow and reactive transport code pflotran [42].fram is used to generate three-dimensional fracture networks.LaGriT is used to create a computational mesh representation of the DFN in parallel.The networks are upscaled to an equivalent continuum porous media using the upscaled discrete fracture matrix (UDFM) model of Sweeney et al. [52].pflotran is used to numerically integrate the governing equations to obtain steady state flow field fields within the domains, and then simulate transport therein.

Network Generation
We generate a set of generic three-dimensional fracture networks at various densities.The networks are generic in that they do not represent a partic-ular field site, but their parameters are based on field site characteristics of fractured crystalline rock [9].The networks are generated in a cubic domain with sides of length L = 50 meters.All networks are generated using the same parameters, except the number of fractures, which we use to control the network density.The networks consist of a single fracture family.Each fracture is a planar disc with an aspect ratio of one, i.e., each fracture is a circle.The radius r of the fractures are sampled from a truncated power-law distribution with a probability density function of and parameters: decay exponent α = 1.8, lower cut off r 0 = 1, and upper cut off r u = 10m.This style DFN generation is referred to as a Poissonian generation opposed to a kinematic generation where fractures nucleate and grow [13,14,41,53].Fracture centers are uniformly distributed throughout the domain.During the generation of the network, we expand the domain slightly (5 meters) in each direction to avoid low-density issues near domain boundaries.
The fracture orientations are uniformly distributed across the unit sphere, which mimics disordered fractures networks observed in the field [22,35].The orientations are sampled from a three dimensional von Mises Fisher distribution, In ( 2), µ is the mean direction vector of the fracture family, T denotes transpose, and κ ≥ 0 is the concentration parameter that determines the degree of clustering around the mean direction.Values of κ close to zero lead to a uniform distribution of points on the sphere while larger values create points with a small deviation from mean direction.To obtain uniformly random orientations, we set κ = 0.1 and mean normal vector of (0, 0, 1).The distribution is sampled using the method detailed in [57].
Variations in the hydraulic properties of the fractures are created by positively correlating each fracture's hydraulic aperture to its radius via The parameters of the relationship are based on field data [51], again in a fractured crystalline rock.This perfectly correlated relationship between fracture size and hydraulic aperture is a common choice for DFN models [8,15,18,25,34,56].Under this model setup, the hydraulic aperture is constant within each fracture, but it varies between fractures with larger fractures having wider hydraulic apertures than smaller fractures.We do not include hydraulic aperture variability into these simulations because our goal is to upscale the networks to an equivalent continuum, and such properties are integrated out.
The set of DFNs is ordered by network density.We adopt a definition of network density based on the percolation parameter p provided by de Dreuzy et al. [16], which is defined as dr min(r, αL)p r (r) .
We define the critical percolation density value (denoted p c ) to be the minimum number of fractures required so there is almost surely a connected cluster of fractures that spans the whole domain, i.e., the value of N such that p in ( 4) is equal to unity [5,11,12,48].For the domain size and truncated power-law distribution parameters considered here, p c ≈ 1000 fractures.This definition of density provides a normalization factor so that networks can be placed in a relative context using a dimensionless percolation parameter where p is the number of fractures in a particular network during generation.We consider values of p = 0.5, 0.75, 1.0, 1.5 and 2.0 (Figure 1).The fractures are colored by their radius to demonstrate the wide range of length scales in the simulations.We consider two sub-cases for each density value.In the first, we consider the network with all of the generated fractures.In the second, we only consider isolated fractures and clusters relative to the flow conditions.i.e., we retain only clusters of fractures that connect the inflow and outflow boundaries.We refer to these cases as isolated fractures retained and isolate fractures removed, respectively.The second case, where isolated fractures are removed, is a common pre-processing step in DFN modeling used to reduce the computational cost of the simulation.These isolated sub-networks do not contribute to network flow, i.e., simulating flow within only the fracture network domain, and can therefore be removed without influencing the simulation.However, the impact of removing/including isolated fractures within the upscaled system is relatively un-characterized, and performing this characterization is a goal of this study.Detecting clusters that connect the domain and isolated clusters is performed using a graph-based algorithm within dfnWorks, cf.[27] and [28] for details of the graph-based techniques.If a network connects the inflow and outflow boundaries, then we say the network is percolating, else the network is non-percolating.The p = 0.5 and 0.75 networks are non-percolating, while all other networks p ≥ 1 are percolating.
Table 1 reports the dimensionless network density (p ), initial number of fractures (N ), number of non-isolated fractures ( N ), percentage of the total number of fractures connecting the inflow/outflow boundaries ( N /N ), fracture intensity with isolated fractures retained (P 32 ), and fracture intensity with isolated fractures removed ( P 32 ).The fracture intensity is defined as where S f is fracture surface area and the sum is taken over all fractures in the network, and V is the total volume of the domain.P 32 is reported in units of m −1 .Since all network generation parameters are fixed, the initial number of fractures placed in the domain increases in proportion to the selected density, cf.(5).In the case of p = 0.5 and 0.75, there are no fractures remaining in the domain once isolated fractures have been removed and P 32 = 0. Therefore, we only consider the cases where all fractures are retained for those two densities, as the upscaled continuum media would otherwise be uniformly matrix cells.
For the percolating networks, p ≥ 1, the percentage of the initial network retained increases slightly with density, but not in the same linear fashion as the total number of fractures initially placed into the domain.Likewise, the fracture intensity increases once the network percolates.

UDFM
Once the networks are generated, we upscale them into a single continuum model using the UDFM method presented by Sweeney et al. [52].There are two main steps to the UDFM method, which are (1) creating the variable resolution continuum mesh and (2) combining the hydraulic attributes of the fractures, e.g., apertures; surface areas; and volumes, with the background matrix permeability and porosity values.We highlight the main aspects of each step below.Complete details of the model along with verification examples are provided in [52].

Meshing
Given a DFN, the UDFM method uses octree refinement to produce a spatially variable Delaunay tetrahedral mesh where the mesh resolution is a function of the distance to the fractures in the DFN, i.e., the mesh is refined more in closer proximity to fractures.The higher resolution close to the fractures allows for the representation of small length scales, which, in turn, allows for the resolution of local gradients in flow and transport simulations, the highest of which occur close to the fractures.The control parameters for the mesh generation are the edge length of the initial mesh element size (l) and the number of refinement levels desired in the final continuum mesh (orl ).
The initial mesh size is also the largest in the domain.The meshing algorithm initializes a uniform hexahedral mesh with edge length l.The cells in the initial mesh that intersect fractures in the DFN are tagged as fracture cells, while the cells that are not intersected by a fracture are tagged as matrix cells.All of the fracture cells and their neighbors (cells that share a face with a fracture cell) are refined into eight sub-hexagons with edge length l/2 to create a child mesh.This process is performed recursively until the desired orl is obtained.The final refined hexahedral mesh is then converted to a Delaunay tetrahedral mesh, whose dual mesh is a Voronoi tesselation that can be used in two-point flux simulation codes.
For these simulations we fix l = 5 m and solely consider the impact of varying orl values.It should be noted that different choices of l and orl can give the same resolution near the fractures.For example, a mesh generated with values of l = x and orl = y will have the same finest mesh resolution as a mesh generated with values of l = 2x and orl = y + 1, or, l = x/2 and orl = y − 1, but not necessarily the same total number of cells.Generally, it is best to balance the computational load of meshing with that of simulation, which means the user must determine a priori the appropriate contribution of each parameter to the resolution and refinement, which itself might change depending on the network.From prior work, we know going beyond orl = 4 creates a bottleneck in the meshing algorithm, but starting with an l too small will defeat the purpose of the refinement to begin with, creating a mesh with too many cells for typical simulation tools [52].Consequently, we select octree refinement levels of orl = 1, 2, 3, and 4. Figure 2 shows the mesh for the p = 1.5 network.As the octree refinement level increases, the fractures are better resolved in the continuum mesh.While larger values of orl are more computationally demanding, they allow for a better continuum representation of the DFN.

False Connections
When mapping a DFN to an ECPM, it is possible to create false connections between fractures.We define a false connection in the continuum mesh as when two fractures that do not intersect in the DFN are upscaled into the same control volume cell in the continuum mesh.Figure 3 presents a twodimensional pictorial explanation of a false connection.The original network of fractures, grey lines, is shown in the upper left sub-figure.There is only one pathway in the network between the left and right domain boundaries.The upper right sub-figure shows a continuum mesh where fracture cells in the mesh are grey and matrix cells are white.Note that in this cartoon, a uniform quad mesh is used for simplicity.In this case, the underlying topology of the DFN is preserved by the upscaled mesh.The lower right sub-figure shows another continuum mesh with coarser resolution.Here, a false connection is created in the bottom of the domain, indicated by a red star, which allows for fracture flow along two pathways rather than one.
A primary design goal of the UDFM method is to minimize the occurrence of false connections, through the use of spatially variable mesh resolution, while also minimizing the number of mesh elements, and in turn computational cost of the simulation.These two desires compete with one another.Selecting a value of orl allows for some control of the number of false connections, but Fig. 2: Mesh refinement levels applied to a fracture network with density p = 1.5.As the octree refinement level increases, we can see the individual fractures getting better resolved.This makes the meshing more computationally demanding, but reduces the number of false connections and allows for more accurate representation of the flow and transport properties of the system.
cannot remove them entirely.The advantage of the UDFM method compared to a uniform resolution mesh is that to achieve the same number of false connections in the domain requires drastically fewer elements, due to the variable refinement.Table 2 reports false connections and mesh information for the different refinement levels across all density values.The total number of false connection (#f), the number of cells with false connections (fc), the total number of Voronoi cells (vc), the percentage of Voronoi cells with false connections (fc/vc × 100).For comparison, the number of cells in an equivalent uniform resolution hexahedral mesh (n) are also included along with a computational cost comparison (vc/n × 100) in terms of reduction of number of mesh cells, which in turn reduces the degrees of freedom in the computational physics simulation and over all cost.By equivalent uniform resolution hexahedral we mean meshing the domain with hexahedrons of size equal to the smallest hexahedral produced using the octree refinement level.The equivalent number of hexahedral cells (n) is given by where L is the domain size length, L = 50.The plus one in the first expression is for end points of the domain edges.
Fig. 3: Example of how false connections can change connectivity in ECPM mesh.In the original network, there is only one percolating pathway connecting the left and right boundaries.When a fine ECPM mesh is used to represent the network, the correct connectivity is preserved, whereas when a coarser ECPM mesh is used to represent the network, an artificial connection appears in the lower half of the domain.
As orl increases, the number of false connections decreases monotonically, as expected.The difference in the number of false connections between subsequent orl values decreases with increasing orl values, i.e., the difference between orl = i and orl = i + 1 is larger than the difference between orl = i + 1 and orl = i + 2. This is true for all observed densities, when isolated fractures are both retained and removed.Removal of isolated fractures also impacts the number of false connections.In the case of retaining the isolated fractures, there are higher percentages of Voronoi cells that contain a false connection when compared to the cases where isolated fractures are removed.However, the difference between the two cases reduces as orl increases.By orl = 3, the two values are quite similar.Likewise, the UDFM mesh size relative to the equivalent hexahedral mesh is comparable when isolated fractures are retained at low orl values, but significantly smaller when isolated fractures are removed.Even for orl = 1, the UDFM mesh has 50% fewer cells than the equivalent hexahedral mesh.As orl increases the UDFM meshes reduce in size relative to the equivalent hexahedral mesh even more, for both cases.
Network density also influences the number of false connections along with the orl values.As the density increases, the number of false connections increases, which is expected as there are more fractures in the domain, and they can be closer to one another.Table 2 provides measurements of local/cell-based false connections.However, false connections can also have an influence on a global/network-scale connectivity.It is possible that a series of local-false connections can combine and create a pathway between the inflow and outflow boundaries of fracturecells through the UDFM domain where the underlying network does not connect, cf., the bottom pathway created in the coarse mesh example provided in Figure 3. Once the UDFM mesh is created, we query the connectivity of the mesh and whether the fracture cells percolate through the domain, i.e., if there is a set of fracture cells in the UDFM mesh that connects the inflow and outflow boundaries.A comparison of the UDFM mesh and DFN percolation is provided in Table 3.For the percolating networks, p ≥ 1, the topological properties match for all values of orl.However, for the non-percolating networks, p =0.5, 0.75, the topology of the UDFM only matches the DFN for the highest level of refinement orl = 4.
Table 3: Comparison of the UDFM mesh and DFN percolation.Legend: "+" indicates the UDFM percolates, "-" indicates the UDFM mesh does not percolate.Green -UDFM mesh topology coincides with the DFN topology for percolating networks.Yellow -UDFM mesh topology coincides with the DFN topology for non-percolating networks.Red -UDFM mesh topology does not coincide with the DFN topology.

Upscaling/Effective Hydraulic properties
Once the domain is meshed, the upscaled/effective hydraulic properties, permeability and porosity, of each cell need to be determined.In the UDFM model, we take the spectral radius of a full rank dyadic permeability tensor constructed using the linear superposition of the fracture apertures, porosities, and coordinate transformation matrix to be the equivalent permeability of the Voronoi cell.Permeability and porosity at each fracture cell is based on the volume of fractures in the control volume, which is calculated using the intersection area of each fracture within each control volume, as well as the fracture orientations, defined by the components of their normal vectors.The aperture of the fractures (b f , [m]) are given by equation ( 3), the surface area of a fracture within a control volume (a f , [m 2 ]) is directly measured using computational geometry, and the volume of a single fracture (v f , [m 3 ]) is given by their product Dividing v f by the total cell volume (v c , [m 3 ]) provides the cell-based fracture porosity (φ f , [-]) The components of the normal vector of the fracture n f = {n 1 , n 2 , n 3 } are used to define a coordinate transformation tensor for that particular fracture Then equivalent fracture permeability tensor (K F , [m 2 ]) is defined at each control volume in the continuum mesh by linear superposition of these terms taken over all fractures intersecting with that cell Notice that K F is a 3x3 full rank second-order (dyadic) tensor and the formulation is loosely based on the cubic law.Given the inability of most flow and transport codes to handle an anisotropic permeability tensor on unstructured meshes, we take the maximum eigenvalue (λ) of K F to be a single equivalent fracture permeability for the cell, where ρ(K F ) denotes the spectral radius of K F defined as Likewise, the total volume v F of N intersecting fractures in the control volume is the summation of individual volumes which is used to determine the equivalent porosity of a fracture cell (φ F , [-]) Finally, the equivalent permeability of the fracture cell is the weighted arithmetic average of the matrix permeability k m fracture and where the equivalent porosity is given by φ F .Equations ( 15) and ( 16) ensure we model the correct volumetric flux out of each cell, while limiting the transport to the fractures.Sweeney et al. [52] verified this approach is conceptually consistent with analytical solutions.We consider a matrix porosity value of φ m = 0.01.These values fall within observed ranges of crystalline rock, of which our fracture network generation parameters are also loosely based [30,58].The matrix cells in the mesh are directly assigned values of k m and φ m .

Flow and transport simulations
Once the single continuum mesh with hydraulic parameters is generated, flow and transport are then simulated.We use the massively parallel flow and reactive code pflotran to perform the simulations [42].We simulate solute transport through steady isothermal flow within the domain.We drive flow through the domain using a pressure difference of 0.001 MPa across the x axis.Neumann, no-flow, boundary conditions are applied along lateral boundaries and gravity is not included in these simulations.The distribution of pressure is modeled using Darcy's Law where The values on the right hand side of k eff are directly observable from the simulations.We model solute transport using an Eulerian formulation.We consider conservative, decaying, and adsorbing tracers.The movement of a conservative (non-reactive) tracer with concentration C(t, x) is governed by where D is the dispersion coefficient that is defined here solely by the molecular diffusion coefficient and is set to 10 −9 m 2 /s.The movement of the decaying tracer (ϕC * ) is governed by where λ is the decay rate constant [1/s], which is used to impose a half life of 100 years.The movement of the adsorbing tracer (ϕC ) with retardation coefficient (R, [-]) is governed by The retardation coefficient R is expressed as where s l is liquid saturation [-], ρ w is water density [kg/m 3 ], and K D is the distribution coefficient [kg/m 3 ].The porosity will vary across the domain based on the upscaled porosity values.The distribution coefficient represents the ration of sorbed to aqueous concentrations.Here, we set the distribution coefficient to R = 4 • 10 3 .All tracers are subject to the following initial and boundary conditions.Initially, there is no tracer in the domain (Ω) Tracer is introduced uniformly into the domain across the entire inlet face (Γ I ) using a pulse injection where δ(0, x) is the Dirac delta function applied to the inlet face at time t = 0.The inflow boundary is assigned a zero-gradient boundary condition to prevent solute diffusing out of the domain in that direction.The outflow boundary (Γ O ) is assigned an adsorbing boundary condition.All other boundaries are reflective.Transport simulations are run for 10 8 years.We measure the tracer mass flow (mol/yr) through the outflow boundary as a function of time, which we refer to as the breakthrough curve (BTC).We normalize these measurements by the initial tracer mass (mol) injected into the system.Time is normalized by the peak arrival time of the BTC for orl = 1 for the conservative tracer.These normalizations facilitate easier comparisons between tracer types (conservative, decaying, and sorbing), orl values, fracture densities, and whether isolated fractures are retained or removed.

Results
In this section, we begin by reporting the measurement of effective permeability.Next, we discuss the breakthrough curves obtained in our ensemble of networks.

Effective permeability
We populate the right hand side of (18) using the numerical solution to obtain an estimate for the effective permeability k eff of each realization.For these samples, we considered a matrix permeability of 1 • 10 −16 m 2 .The values are plotted in the top row of Figure 4 with "isolated fractures retained" realizations on the left and "isolated fractures removed" on the right.Each line corresponds to a density value p .Values of k eff increase with increasing p values due to larger numbers of fractures with higher permeability being placed into the domain.The values of k eff monotonically decrease with increasing orl values, in all cases.Recall that the background matrix permeability is 10 −16 m 2 , which provides a lower bound to k eff .The difference between k eff at subsequent orl levels decreases as orl increases, which indicates convergence to a relatively stable value should occur with a sufficient number of refinement levels.Due to computational limitations, expanding to higher orl to identify this value was not possible for the higher density values p .Nonetheless, we examine the convergence rate by computing the relative difference/error (e i ) of orl =1,2,3 compared to the value obtained with orl =4 using the definition where k i eff denotes the effective permeability at orl = i.These values are plotted in the lower row of Figure 4.In the case of "isolated fractures retained", the error values seem to decrease exponentially with orl values.The reported error for orl =1 ranges from ≈ 15 − 80 times that of orl = 4, with the highest values observed for the non-percolating networks, p < 1.As we increase the number of refinement levels, the error decreases, and by orl = 3, the estimates are all less than 2, which given the small absolute values, O(10 −14 − 10 −15 ), are likely acceptable.In the case of "isolated fractures removed", the error values also declines exponentially with orl.ere, the estimates provided at the lower orl levels are much closer to the orl = 4 value when compared to the "isolated fractures retained" counter parts.For orl = 3, the values are less than 1.
Figure 5 reports the effective permeability for additional matrix permeability values of 1•10 −12 m 2 , 1•10 −14 m 2 , 1•10 −16 m 2 , and 1•10 −18 m 2 for p = 0.50 and p = 0.75, but with isolated fractures retained.We only consider the samples because they are below the percolation threshold and the topology of the UDFM mesh and DFN do not match.Similar to the previous samples, the effective permeabilities decrease with orl value.However, the rate of decreasing, depends on the background matrix permeability.For the highest values, there is relatively little change with orl value when compared to the lower values.A typical fracture permeability value in our model setup is around 1 • 10 −9 m 2 .So, the case with k m = 1 • 10 −12 m 2 , the domain is relatively homogeneous.However, for the lowest matrix permeability values, 1 • 10 −18 m 2 , the domain is highly heterogeneous, with permeability values ranging close to 10 orders of magnitude in the domain.In turn, if a macro-scale false connection forms, as in these cases, there is highly permeability channel through the domain.This results in effective permeability errors on the order of 4 orders of magnitude.Most notably, the effective permeability for the 1 • 10 −18 m 2 case does not stabilize as the values for the other ones do.

Breakthrough curves
The breakthrough curves (BTC) of the tracers through the outlet face of the domain as a function of time are discussed in this section.We begin with the two densities with networks that do not percolate (p = 0.5 & 0.75).Then we discuss the three densities that do percolate (p = 1.0, 1.5 & 2.0).

p < 1
Figure 6 (a) shows the breakthrough curves in the fracture network with density p = 0.5.We begin the descriptions with the conservative tracer, which is shown as solid lines.For octree refinement level orl =1, there is an initial peak and then a plateau of tracer slowly exiting the domain.For octree refinement levels orl =2 and 3, there are two distinct peaks.There is an initial peak corresponding to transport through the fractures and then there is a second peak corresponding to transport through the matrix.At orl =4, there is only a single peak arriving at the same time as the secondary peaks in orl =2 and 3. Recall that for this density, there was not only a significant number of false connections for orl =1, 2, and 3, but the global topology of the UDFM model did not match the DFN, i.e., fracture cells in the UDFM mesh percolate through the domain while the DFN does not percolate, cf.Table 3.Thus, there is a connected pathway of fracture cells in the UDFM model for orl =1, 2, and 3, which leads to the early peak.However, physically speaking, the tracer must travel through the matrix, which is the case for orl =4.In the case of the decaying tracer, orl =1, 2, and 3 also show an initial peak in tracer concentration due to transport through the fracture cells but the tracer decays before we observe the second peak, which corresponds to travel through the matrix.
For orl =4, there is a peak of tracer passing through the matrix, but is significantly delayed.The sorbing tracer breakthrough curves roughly resemble the conservative tracer curves, being significantly delayed due to adsorption, though they are more spread out than the conservative and decaying tracers.At the exit through the matrix, the peak height of tracer concentration has reduced significantly due to delay.At orl =1, there a wide spread in arrival times for all tracers.As we increase the octree refinement level from orl =2 to orl =4 for the fracture network with removed isolated fractures, the actual behavior of the conservative tracer is  revealed -the initial peak in tracer concentration disappears because there is no fracture connection from the inflow to the outflow boundary.
The breakthrough curves for a fracture network with density p = 0.75 is shown in Figure 6 (b).Like p = 0.5, the DFN does not percolate, but orl =1, 2 and orl =3 UDFM meshes do.In general, the behavior is quite similar to what was observed for the p = 0.5 network.Once the percolation structure of the UDFM model matches the DFN, then matrix dominated transport is observed.The falsely percolating simulations show a dual peak structure that indicates an initial peak through a connected fracture network and a second peak due to transport through the matrix at the exit.The profile of the BTC for the sorbing tracers are delayed and somewhat different from conservative tracer, i.e., they have less pronounced peaks.

p ≥ 1
Figure 7 shows the BTCs for p = 1.0:(a) isolated fractures retained and (b) isolated fractures removed.Colors and line styles are the same as in Figure 6.Rather than seeing a stark change in the BTC behavior due to a change in global connectivity, we observe the influence of changes in local connectivity and mesh refinement.Although the global connectivity of the UDFM model and DFN correspond for all orl values, the number of local false connections rapidly declines with increasing orl value, cf.Table 2.For the case where isolated fractures are retained, 59% of control volumes have a false connection at orl =1, compared to less than 1% at orl =4.In the case of isolated fractures are retained and orl =1, we see a single peak arrival with a long tail for all tracer types.As the orl increases, we see the formation of a second, matrix-  transport-dominated peak.This latter peak is likely the result of better defined pathways through the system due to mesh refinement.
For the case where isolated fractures are removed, the convergence with rising orl is more pronounced.This is true especially true for transport through the matrix where the tails of the distribution collapse onto one another.Recall that the number of local false connections is much smaller compared to the case where isolated fractures are retained.
Figure 8 depicts the BTCs for p = 1.5.The BTC show more or less the same behavior as in Figure 7 with respect to how the BTC is influenced by the orl value and whether or not isolated fractures are retained or removed.A small difference is that the latter peak, which corresponds to matrix transport, does not show the same collapse, nor the same valley between the two peaks.
Figure 9 shows the BTCs for p = 1.5.Here, convergence of the BTC with increasing orl is much clearer than in the other cases, especially for the case where isolated fractures are removed.Although the case where isolated fractures are retained has more cells overall, the local false connections and dispersion influences convergence of the breakthrough curves.

Matrix Permeability Variations
Figure 10 and 11 show the breakthrough curves for p = 0.5 and p = 0.75, respectively, for the four matrix permeability values that we considered.Recall that the domain for the k m = 10 −12 m 2 values are relative homogeneous, and that is reflected in the BTC.The BTC here, appear more like a pulse traveling through a low-variance porous media than and fractured system.For orl = 1, the pulse is wide, and then narrows with rising orl.The permeability is suffi- Fig. 9: Breakthrough curves of Eulerian tracer in a fracture network with density p = 2.0.c is the concentration of the tracer at time t, and c 0 is the total amount of tracer in the pulse.We examine the solute transport for fracture networks with (a) isolated fractures retained, and (b) isolated fractures removed.
ciently high that we do not see the effects the decay rate in the decaying tracer.
As the matrix permeability rises, the influence of the embedded fractures can be seen in the form of the false early breakthroughs due to the network-scale false percolation.These results suggest that not only is the topology of the UDFM mesh important, but there is an interplay between the heterogeneity of the domain and topology.We have provided a set of numerical simulations to characterize the influence of mesh refinement level on flow and transport properties in continuum representations of fractured porous media.We generated a set of generic threedimensional fracture networks at various densities, where the radii of the fractures were sampled from a truncated power-law distribution, with parameters  loosely based on field site characterizations.We considered five network densities, defined using a dimensionless version of density based on percolation theory.Two values were below the critical percolation value, one at the value, and two above.We considered two sub-cases for each density value.In the first, we considered the network with all of the generated fractures.In the second, we only considered non-isolated fractures and clusters relative to the flow conditions.i.e., we retained only clusters of fractures that connected the inflow and outflow boundaries.Once the networks were generated, we upscaled them into a single continuum model using the UDFM method presented by Sweeney et al. [52].We considered four octree refinement levels of orl = 1,2,3, and 4. We considered steady isothermal pressure-driven flow through each domain and then simulated conservative, decaying, and adsorbing tracers using a pulse injection into the domain.
Keeping the isolated fractures resulted in more false connections when compared to the cases where isolated fractures were removed, and thus more refinement levels were needed to resolve the true structure of the fracture network.This was most apparent for the low densities (p = 0.5 and p = 0.75), where we found that even a refinement level of orl =3 was not sufficient to eliminate the false connections in the network and the UDFM mesh percolated, when in fact, the DFN did not.These global false connections dramatically affected the values of effective permeability.The estimates provided by orl =1 and 2, were multiple orders of magnitude higher than that provided by orl =4, which did not percolate.In the cases where the topology matched between the UDFM model and DFN, p ≥ 1, the differences in the estimates were smaller.Likewise, when we considered a range of effective permeability values, we saw that not only is the topology of the UDFM mesh important, but there is an interplay between the heterogeneity of the domain and topology.For cases with lower matrix permeability, the error in effective permeability estimations was much higher than in the case with higher matrix permeability.
In terms of transport, the conservative tracer breakthrough curves for orl =1,2,3 for p = 0.5, and p = 0.75 showed erroneous peaks caused by artificially connected pathways in the system.Once the topology of the UDFM mesh and DFN match, then only later time breakthrough is observed, which is consistent with the physical system.When increasing the density of the fracture network, the number of possible percolation paths increased and thus the effect of the false connections to the shape of the breakthrough curves decreased.While the global topology matched, there were still changes in local connectivity of the mesh, and more dispersion in the breakthrough curves was observed.We observed a dual peak structure, where the first peak was slightly higher than the second one, that is shortly before most of the solute had left the domain.The first peak was due to the flow and transport through the fracture cells, while the second peak was due to transport through the matrix.The breakthrough curves for all refinement levels and densities decayed sharply after the second peak and overlapped for all fracture networks above the percolation threshold.For the decaying tracer, the breakthrough curves are narrower and in most cases, disappeared from the system prior to the expected matrix dominated breakthrough.In the case of the non-percolating network, the decaying tracer never broke through the matrix at the expected time due to the decay rate.Otherwise, the decaying tracer behaved like the conservative tracer.The sorbing tracer showed delayed first arrival times (approx.10 4 ) in comparison to the conservative and decaying tracers for the percolating networks.When looking at the non-percolating networks, we observed a drastic increase of the first arrival time depending on the refinement level.For example, the first arrival time was at about 10 4 for refinement level orl =1, and up to 10 6 for refinement levels orl =3 and orl =4.
In the case of higher matrix permeability, the BTC appear more like a pulse traveling through a low-variance porous media than a fractured system.For orl = 1, the pulse is wide, and then narrows with rising orl.As the matrix permeability rises, the influence of the embedded fractures can be seen in the form of the false early breakthroughs due to the network-scale false percolation.These results again emphasize that not only is the topology of the UDFM mesh important, but there is an interplay between the heterogeneity of the domain and topology.
The results of the simulations point to a few key considerations when upscaling fracture networks.Foremost, selecting a mesh resolution so that the global topology of the upscaled mesh matches the fracture network is essential.If the upscaled mesh has a connected pathway of fracture (higher permeability) cells but the fracture network does not, then the estimates for effective permeability and solute breakthrough are incorrect.Previous work that studied the consistency of ECPMs relative to DFNs did not consider these global topological changes in their analyses [32,37].To prevent such situations, identifying percolation in the fracture network and upscaled mesh is an essential part of the modeling process.Such queries can be easily performed using graph theory software such as networkX [21], among others.Local false connections also impact transport behavior, but to a smaller degree than if the global connectivity is incorrect, which is consistent with previous work [32,37].These small short circuits in the system, which occur more frequently if isolated fractures are retained, increase dispersion of the solute plume.We did not discriminate between dispersion caused by false connections themselves and numerical dispersion caused by mesh size and numerical scheme, but this would be a problem worth investigating in future work.False connections cannot be eliminated entirely, but they can be managed with mesh resolution.While using a uniform mesh resolution is simple, adopting that methodology can lead to intractable numbers of small cells.On the other hand, adopting octree meshing to obtain sufficient levels of refinement leads to fewer computational cells (up to a 90% reduction in overall cell count), but it is more laborious to create the mesh.
The computational representation of fractured porous media remains a complex part of subsurface flow and transport modeling.While numerical models are rapidly advancing with the introduction of multi-dimensional DFM models, kinematic DFN, embedded fracture matrix models, and coupled thermohydro-mechanical-chemical models [1,2,3,4,7,10,13,14,17,19,31,36,41,43,47,49,50,53,55], mesh generation remains a key step for all modeling endeavors.Therefore, continued characterization of how meshing techniques impact physical phenomena beyond the standard convergence of the numerical scheme remains requisite for modeling endeavors.Extensions of this research should focus on how the influence of meshing methods interplay with various upscaling techniques as well as the aforementioned modeling methodologies.
Campaign of the U.S. Department of Energy Office of Nuclear Energy.Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525.The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.Assigned LA-UR-23-20307.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Fig. 1 :
Fig.1: Network realizations from each density considered in this publication (p = 0.5, 0.75, 1.0 and 2.0).The fracture radius is color-coded to highlight the range of length scales in the simulations.When the network's density is below the percolation threshold (p < 1), the fracture network often does not connect between the inflow and outflow boundaries.For networks above the percolation threshold (p > 1), we can clearly see at least one path connecting the inflow and outflow boundaries.As the density increases, the number of fractures increases and thus the number of possible paths though the fracture network.
q is the Darcy flux [m/s], k is the permeability [m 2 ], µ is the fluid viscosity [Pa s], and P is the fluid pressure [Pa].Simulations are performed with the fluid temperature at 20 • C, which corresponds to a viscosity of µ = 8.9e−4 Pa s for the pressure values considered.Once steady conditions are obtained, we can invert (17) and obtain a numerical estimation for the effective permeability of the block

Fig. 4 :
Fig. 4: Effective permeability of samples at different orl values with a matrix permeability of 1 • 10 −16 m 2 .(a) Isolated Fractures Retained -Effective Permeability (b) Isolated Fractures Removed -Effective Permeability (c) Isolated Fractures Retained -Error factor of effective permeability.(d) Isolated Fractures Removed -Relative error of effective permeability.

Fig. 7 :
Fig. 7: Breakthrough curves of Eulerian tracer in a fracture network with density p = 1.0.(a) isolated fractures retained and (b) isolated fractures removed.

Table 1 :
Network Characterization: dimensionless network density (p ), initial number of fractures (N ), number of non-isolated fractures ( N ), percentage of the total number of fractures connecting the inflow/outflow boundaries ( N /N ), fracture intensity with isolated fractures retained (P 32 ), fracture intensity with isolated fractures removed ( P 32 ).

Table 2 :
UDFM mesh information for network densities p at different orl values where n is the equivalent number of hexahedral cells.Number of false connections (#f), number of cells with false connections (fc), total number of Voronoi cells (v), percentage of Voronoi cells with false connections (fc/vc[%]), and relative mesh size compared to equivalent hexahedral mesh (vc/n[%]).
Fig.8: Breakthrough curves of Eulerian tracer in a fracture network with density p = 1.5.c is the concentration of the tracer at time t, and c 0 is the total amount of tracer in the pulse.We examine the solute transport for fracture networks with (a) isolated fractures retained, and (b) isolated fractures removed.