Abstract
In this study, we analyze the impact of transforming a discrete fracture network (DFN) model to an equivalent continuum model (ECM) on flow and solute transport characteristics. The analysis was conducted in a setting derived from – though simplified relative to – in-situ fracturing conditions at the Forsmark site. The geometrical structure of the DFN model considered is combined successively with a highly simplified transmissivity model (single constant value) in order to isolate spatial and structural effects, and with a more realistic model in which transmissivities depend on both fracture size and orientation (indirectly reflecting mechanical stress conditions). Two upscaling methods are considered: a simplified geometry-based method and a numerical flow-based method, in which local ECM cell properties are directly informed by local flow and transport simulations. For both approaches, we quantify the impact of ECM resolution. Specifically, we assess the sensitivity of local properties as well as local and global flow and transport indicators, to both the upscaling method and the ECM grid resolution. The results demonstrate that the transformation from DFN to ECM overestimates hydraulic conductivity, underestimates the geometric porosity but to a lesser extent and overestimates the transport porosity. This also artificially reduces flow path variability and tortuosity, except at very high resolutions. Consequently, average transfer times are shorter in ECMs than in DFNs, with discrepancies increasing as grid resolution coarsens. Similar trends are observed for first arrival times and mode (the peak of the distribution), but to a lesser extent. DFNs are also more likely to have very long transport times. Finally, we show that ECMs derived from geometry-based upscaling are highly sensitive to grid resolution, whereas flow-based upscaling exhibits significantly lower sensitivity.
1 Introduction
Discrete fracture network (DFN) modelling is often considered the most appropriate approach for quantitatively describing geometry, groundwater flow, and solute transport in fractured rock formations (; ; ; ). While DFN modelling may be an appropriate approach, it also presents significant challenges. One challenge concerns data requirements, as DFN approaches require more extensive input data than continuum approaches based on representative elementary volumes. Another challenge relates to computational demands: for a given model domain size, DFN models are typically more computationally expensive than continuum representations. The latter issue is the topic of the present paper, where we specifically examine whether the computational demands of DFN models can be relaxed by upscaling the DFN to an equivalent continuum representation. For such upscaling to be considered successful, the resulting continuum models must reproduce the flow and transport characteristics of the underlying DFN.
Modelling flow and transport processes in a heterogeneous porous medium requires the specification of several physical properties, including hydraulic conductivity, porosity, and storativity. It should be noted that upscaling flow and upscaling transport may require different approaches and may not yield equivalent results (; ). Upscaling hydraulic conductivity values from the measurement scale to the numerical grid scale has been extensively investigated in the context of heterogeneous porous media (e.g., ; ; ). Upscaling in poorly fractured media, where fractures are embedded in a quasi-impermeable matrix, introduces additional challenges. Fractures are essentially two-dimensional features that must be represented within three-dimensional volumes; they span multiple length scales and orientations, and the occurrence of connected fracture clusters between hydraulically active boundaries is essential for defining an equivalent permeability.
Broadly speaking, upscaling methods for fractured and porous-fractured media can be divided into two categories: analytical approaches (e.g., ; ; ; ) and numerical approaches (; ; ; ; ; ; ). In their review, emphasize anisotropy and the representation of equivalent properties in tensorial or non-tensorial forms. introduced an early flow-based upscaling approach for fractured rocks, assuming flow concentrated within fractures surrounded by an impermeable rock matrix. Their method applies linear head boundary conditions and several numerical experiments to derive equivalent permeability tensors. proposed two variants of the flow-based approach for two-dimensional fractured and porous rocks, using single- and multi-boundary integration and accounting for a permeable rock matrix in addition to the fracture network. also assumed a permeable matrix around the fractures and developed a pressure-flux volume-averaging method to derive a full equivalent block permeability tensor without a priori assumptions regarding the orientation of the eigenvectors, requiring numerical experiments under classical permeameter conditions. extended the multi-boundary method to three-dimensional porous fractured rocks and compared several flow-based and geometry-based methods, but for DFN models containing only a few tens of fractures. We note that, for most of the flow-based approaches cited above, especially for porous fractured rocks, the focus lies on issues raised by numerical schemes under consideration for the ECM calculation, and fractures larger than the ECM meshes. Consequently, the DFN complexity is constrained by fracture–matrix meshing considerations. As a result, very few studies address the large size range and high DFN complexity investigated in the present paper.
In this study, we assess under which conditions upscaled ECMs provide reliable representations of flow and advective transport, and we highlight differences between alternative upscaling methods. Specifically, we compare a geometry-based upscaling proposed by , coupled with a finite-volume numerical scheme for ECM calculations, with a single-boundary integration flow-based upscaling method. The underlying hypothesis is that flow-based upscaling more effectively preserves key characteristics such as permeability, anisotropy, and transport regimes. The performance of the different upscaling techniques is evaluated by direct comparison of the ECM results with those obtained from the underlying DFN models.
The DFN model used in this study reflects conditions at the Forsmark site in Sweden (; ; ), which has been selected for a planned spent nuclear fuel repository. The site is characterized by sparsely fractured crystalline rock, for which a DFN representation is considered appropriate. Such sparsely fractured formations pose challenges for upscaling from DFN to ECM, due to limited connectivity and likely anisotropic flow behavior with bottlenecks in the fracture network ().
The paper is structured as follows. Section 2 describes the upscaling methods considered. Section 3 presents the numerical setup of the test cases. The results are discussed in Section 4, followed by concluding remarks in Section 5.
2 The DFN-to-ECM transformation
2.1 Parameterization of an ECM (element properties)
An ECM represents a rock mass in three dimensions using properties defined per volumetric element. For the simulation of groundwater flow and advective and non-reactive solute transport, the key parameters are hydraulic conductivity and porosity. Hydraulic conductivity is defined in tensor form to account for flow anisotropy, while porosity is defined as a scalar quantity. The hydraulic conductivity is hence represented by the following matrix:
The tensorial representation imposes additional constraints on its coefficients (; ; ), as summarized in Equation 2. The conductivity matrix must be invertible and symmetric; furthermore, the block determinant for indices and where must be strictly positive. The same constraint applies to the determinant of the matrix itself.
The six components of the ECM hydraulic conductivity property can be reduced to three if the orientations of the principal axes are known in advance. Under isotropic conditions, the tensor further reduces to a single scalar value.
Spatial variations in ECM properties are defined at the resolution scale of a regular grid composed of cubic cells. This resolution can be uniform (one value, e.g., Figures 1a,b) or variable, and is informed by the DFN itself (Figures 1c,d).
FIGURE 1
The resolution of the ECM grid must be considered in relation to the intrinsic characteristic length scales of the DFN. These include the minimum and maximum values of the fracture length distribution, the average spacing between fractures (2/P32), and the connectivity scale. For DFNs characterized by power-law length distributions, as is the case here, these characteristic scales are well defined and widely documented.
2.2 Calculation of cell properties
When representing a DFN by equivalent properties within an ECM cell, DFN connectivity, both within the cells and with the cell faces, is a key issue. Discretisation inherently introduces artificial connectivity, either by merging initially distinct flow paths that are closer than the cell size or by adding new connected paths that do not exist in the original DFN. To mitigate these effects, the transformation from DFN to ECM begins with a filtering step. In this step, only fracture clusters connected to the external boundaries of the full domain are retained, while isolated clusters are removed, as they do not contribute to flow and transport.
After the DFN connectivity filtering, we consider two upscaling methods for transforming the DFN into ECM properties. The first is a simplified geometry-based method (hereafter referred to as geom-based), and the second is a numerical method based on direct flow simulations (hereafter referred to as flow-based). These approaches are described below.
In the geom-based approach, a DFN is transformed into grid-cell properties using the volume–intersection method proposed by . This process is briefly summarized here. Each conductive fracture in the DFN is represented as a tabular volume with a prescribed aperture. All fractures intersecting a given cell are identified, and the intersecting volume between each fracture and the cell is computed. Each selected fracture contributes to the grid cell value of a variable by an amount equal to the intersecting volume times the value of the variable in question. The contributions from all the intersecting fractures are summed and normalized by the total cell volume to obtain the upscaled cell value (see, e. g., Equation 3 for the porosity variable). The finite-volume numerical scheme implemented in DarcyTools (; ) employs staggered grids. Scalar variables like pressure and salinity are defined at cell centers, while vector quantities such as velocity are defined at cell faces. Porosity is hence defined at cell centers while hydraulic conductivity is defined at cell faces. For hydraulic conductivity, three offset grids are used – one for each direction – each shifted by a half-cell. This leads to the definition of three directional conductivity components (), corresponding to cell faces whose normal vectors align with the Cartesian direction . The principal directions of the hydraulic conductivity tensors are by definition the ones of the ECM Cartesian cubic grid (off-diagonal terms , with and the Cartesian direction, are therefore equal to zero by definition).
To summarize, the cell porosity of the geom-based upscaling is:where is the cell volume, is the volume of the entire fracture, and the fraction of the fracture volume within the cell, and the summation index over all the fractures that intersect the cell. The equivalent hydraulic conductivities (in the direction ) are computed following the same principle.
For the flow-based approach, grid cell components are estimated using direct numerical simulations, replacing the geom-based calculations. In the general formulation expressed below, three flow simulations are required to derive the six components of the equivalent permeability tensor (Equation 1).
Each flow simulation is based on an imposed linear head boundary over all the cell faces (hereafter denoted GRD, for imposed head GRaDient). Successive GRD conditions are applied with the head gradient being collinear with each of the Cartesian directions of the ECM grid, and the total flow rate per cell face being computed (see Figure 2) and leading to a total of 18 measurements. For a flow simulation and head gradient in the direction , the cell face flow rates and ( and subscripts refer to the two opposite faces of the cell in the direction , where is a cartesian direction) are computed. Note that the GRD boundary condition does not force flows through opposite sides to be equal (i.e., ). The hydraulic conductivity components , and (line in Equation 1) are next deduced by using the following expression (Equation 4):
FIGURE 2
By averaging the flow rates across opposite cell faces, the total of 18 flow-rate measurements is reduced to 9 components (Equation 1). A further, commonly applied simplification to satisfy the symmetry requirement of the equivalent conductivity tensor is to reduce the number of independent coefficients from 9 to 6 by averaging opposite off-diagonal terms. In the present study, the procedure is simplified further by fixing the tensor principal directions a priori to coincide with the ECM grid axes (for , components are neglected). As a result, the number of parameters is reduced from 6 to 3. If a cell has no connectivity between any of the three pairs of opposite faces, it is assigned a null hydraulic conductivity (and porosity). If at least one pair of opposite faces is connected, the cell is set as conductive, and any unconnected diagonal components are assigned a residual (minimum) value. Note that flow simulations under GRD boundary conditions, as described here at the cell scale, are also performed at the full DFN and ECM grid scales to enable direct comparison between DFN and ECM results derived using the two upscaling options. In this case, the 18 large scale (entire domain sides) flow-rate measurements are reduced to six main values following the same simplification procedure.
The flow-based definition of cell porosity is adapted slightly from . It is defined as a scalar quantity obtained by averaging the results of three transport simulations performed along the three Cartesian directions of the cell. For each direction, permeameter (PRM) conditions are set. Unlike GRD conditions (Figure 2), the lateral cell faces (i.e., faces not orthogonal to the head gradient) are assigned no-flow boundary conditions. After computing the head field, a transport simulation is performed (see Annex 7.3), involving an instantaneous, flow-weighted injection of particles at the inflow face and recording of particle arrival times at the outflow face. For each tested direction, the flow-based porosity is calculated using Equation 5 (where, denotes the mean arrival time, the flow rate through the cell, the distance between entry and exit faces, and the area of the face. Only one scalar porosity value is then obtained as the arithmetic average of the three directional values as defined from Equation 5 (unconnected directions are not considered).
The flow-based porosity is the volume explored by the particles in the three transport tests, which is to be smaller than, or equal to, total (or geometric) porosity. We verified that the number of particles was sufficient to explore areas with very low flow rates; however, dead ends associated with boundary conditions are not counted. This underestimation may be artificial if local cell dead ends are not dead ends at the global level. The advantage of this method is that it respects average transport times.
2.3 Numerical methods
Laminar flow conditions with Newtonian fluids (Darcy conditions) are assumed. DFN numerical flow and transport simulations are performed with DFN.lab (https://fractorylab.org/dfnlab-software/). ECM numerical flow and transport are performed in DarcyTools, based on a finite volume scheme (; ); and in DFN.lab by using a 3D first-order Hybrid Discontinuous Galerkin finite element method, noted HDG (; ). The HDG method ensures mass conservation, whilst simultaneously providing constant element-wise hydraulic head on both faces and within the elements. Tetrahedral and hexahedral mesh can be used, noted HDG-tet and HDG-hex (the later provided by the library NGsolve (). Both HDG-tet and HDG-hex schemes use Raviart-Thomas (RT) basis functions to discretize the flow field. DFN transport simulations are based on a Lagrangian approach and advective particle tracking without diffusion and sorption or reactivity components. At the scale of a DFN, dispersion emerges from the heterogeneity of local velocities.
3 Case study and numerical experiments
The study is being conducted in the context of the Forsmark site in Sweden (; ), using fracturing conditions and a DFN model recipe representative of the site. Flow and transport are modelled within a simple cubic domain with an extent of 500 m (Figure 3).
FIGURE 3
The DFN model recipe has five fracture orientation sets (see Table 1). Fracture sizes follow a power-law distribution, with a common range (minimum and maximum sizes) for all the sets, while the scaling exponent and density term are defined separately for each set (Table 1). The recipe is assumed to describe all transmissive (i.e., open) fractures; non-transmissive fractures are excluded from the outset. Open fracture intensities, per set and per fracture size range, are defined. Over the full range [; ], the total fracture intensity () is equal to 1.02 . The evolution of fracture intensity with increasing truncation length (fracture cutoff diameter ) can be readily computed. As increases beyond the minimum value (equal to , the fracture intensity decreases significantly, from 1.02 to 0.053 and to 0.037 for equal to , 3 m, and 12 m, respectively. Simultaneously, the DFN percolation parameter – a statistical indicator of DFN connectivity – decreases only slightly, from 3.41 to 3.25 over the entire domain. Values above the percolation threshold (i.e., above ∼2.5) indicate that the DFN remains connected at the domain scale. The minor decrease in percolation parameter, as the DFN cutoff length increases, reflects that large-scale connectivity is primarily influenced by the largest fractures. These variations with fracture cutoff size, in fracture intensity and percolation parameter, are directly linked to the scaling exponent of the DFN recipe.
TABLE 1
| Set name | Fisher distributions for orientations | Power-law size distribution | Intensity term | ||||
|---|---|---|---|---|---|---|---|
| | Trend (°) | Plunge (°) | Dispersion | Radius min () | Radius max () | Scaling exponent () | |
| NS | 87 | 2 | 21.7 | 0.038 | 564 | 2.50 | 0.142 |
| NE | 135 | 3 | 21.5 | 0.038 | 564 | 2.70 | 0.345 |
| NW | 41 | 2 | 23.9 | 0.038 | 564 | 3.10 | 0.133 |
| EW | 190 | 1 | 30.6 | 0.038 | 564 | 3.10 | 0.081 |
| SH | 343 | 80 | 8.2 | 0.038 | 564 | 2.38 | 0.316 |
| Total | | | | | | | 1.02 |
DFN model recipe () of transmissive fractures with the parameters from .
The DFN recipe is completed by two transmissivity models (Table 2): one with constant transmissivity, labelled , and one with distributed transmissivities, labelled . In the latter model, transmissivity is correlated to the fracture size (size to the power of ) and fracture orientation reflecting the in-situ stresses () and Table 3. Under the thrust-fault stress regime assumed for the Forsmark site, fractures experiencing lower normal stress exhibit higher transmissivity. As a result, the sub-horizontal fractures tend to be more transmissive than the vertical ones. The distributed transmissivity model is considered more realistic than the constant-transmissivity model. The parameters of the transmissivity laws are not further calibrated to in-situ hydraulic data, as absolute flow rates and arrival times are only compared across modelling cases. For the sake of simplicity, hydraulic, mechanical, and transport apertures () are assumed equal. Fracture properties for flow and transport are hence fully defined by a single parameter, the transmissivity, from which a unique aperture is derived using the classical cubic law.
TABLE 2
| Label | Reference | Definition |
|---|---|---|
| Constant | ||
| Stress and size | , With the fracture diameter and the normal stress on the fracture, the hydraulic aperture, the fluid density, the gravity and the fluid viscosity. |
Fracture transmissivity models.
TABLE 3
| (MPa) | (MPa) | (MPa) | Trend (°) |
|---|---|---|---|
| 5.3 | 13.6 | 23.9 | 145 |
Thrust fault regime stress field representative of depths of a few hundred meters.
An arbitrary number of ten DFN realizations of the DFN geometrical recipe are generated and analyzed in this study. Each realization is first examined independently; ECM performance is evaluated by direct comparison with the corresponding DFN realization. Boxplot representations (e.g., Figure 8) are subsequently used to illustrate the variability in results arising from the DFN stochasticity and to identify qualitative trends.
Each DFN realization is combined successively with the two transmissivity laws, resulting in a total of 20 DFN cases. For each DFN realization, several ECM transformations are studied. These include single-resolution grids (S5, S10 and S50, applied to all ten DFN realizations) and multiresolution grids (DG1 and DG2, applied only to DFN realization 1). Both grid types are combined with the flow-based and geom-based upscaling approaches. The full set of combinations is summarized in Table 6.
Two types of hydraulic boundary conditions are considered. The GRD boundary conditions (Figure 2) are used both at the cell scale for flow-based upscaling (Section 2.2) and at the domain scale. In addition, classical permeameter (PRM) boundary conditions are applied in the transport simulations.
4 Results–performance of ECM methods
4.1 ECM cell properties
ECM grid properties are computed for all DFN realizations, upscaling approaches, and grid resolution combinations (Tables 4–6). The sensitivity of local property distributions – hydraulic conductivities (, , ) and porosity () at the cell scale – to grid resolution is illustrated for the flow-based upscaling using DFN realization 1. Results are shown for both single-resolution and multi-resolution grids (Figures 4a,b) for hydraulic conductivities, and c) and d) for porosities). First, at a fixed ECM grid resolution (i.e., cell size), the distributions of local permeabilities in the x, y and z directions are broadly similar for the constant transmissivity model (see, e.g., the light-blue square symbols for S50 in Figure 4a). When the transmissivity model is used, the distribution of local (vertical direction) exhibits a slightly higher shift compared to and (S50 in light-blue square symbols for transmissivity law in Figure 4b). This reflects the increased transmissivity of sub-horizontal fractures. However, these differences remain modest compared to the much stronger influence of ECM grid resolution. For the coarsest grid, the conductivity distribution is the narrowest, spanning two orders of magnitude and exhibiting the lowest range of values. Increasing the resolution shifts the distribution towards higher conductivities and a wider range (up to almost four orders of magnitude). Despite their structural differences (single- versus multi-resolution), the DG1 and S5 grids have a similar maximum resolution of 3.9–5 m and produce distributions that are very close to each other. When the constant transmissivity model is replaced by the distributed TSL model, which correlates with lengths and orientations, the resulting ECM distributions widen further, by up to five orders of magnitude for the same DFN realization. Although the anisotropy component increases (z differs from x and y), the dominant trend remains the first-order dependence on grid resolution. At the domain scale, this apparent resolution scale effect is counterbalanced by the relative increase of the grid connectivity. As cell size decreases, the number of cells increases and the proportion of empty cells grows (as illustrated qualitatively in Figure 1). In addition, the nature of the DFN representation within a cell evolves: cells increasingly contain single fractures rather than clusters of fractures as resolution increases.
FIGURE 4
TABLE 4
| Label | Resolution type | Domain () | Cell size | Applied to |
|---|---|---|---|---|
| S50 | Single | 50 | 10 DFN realizations (*2) | |
| S10 | Single | 10 | 10 DFN realizations (*2) | |
| S5 | Single | 5 | 10 DFN realizations (*2) | |
| DG1 | Multi3 | [3.9 to 15.6] | Only for DFN realization 1 | |
| DG2 | Multi4 | [1.95 to 15.6] | Only for DFN realization 1 |
ECM grids summary.
TABLE 5
| Geometric recipe parameters | Nb realizations | Transmissivity law | ECM grids |
|---|---|---|---|
| See Table 1 | 10 | Realization 1: S5, S10, S50, DG1, DG2 Other realizations: S5, S10, S50 | |
| | 10 | Realization 1: S5, S10, S50, DG1, DG2 Other realizations: S5, S10, S50 |
DFN and realizations.
TABLE 6
| Grid Label | Connectivity filter | Grid type | Cell size(s) | ∼ Cells | Upscaling | Hydraulic boundary conditions | Flow and transport solver | Transport |
|---|---|---|---|---|---|---|---|---|
| DFN | 6 | | | | | GRD | DFN.lab | x |
| DFN | 2 | | | | | PRM | DFN.lab | x |
| S50 | 6 | Single | 50 | Geom-based | GRD | DarcyTools | | |
| S10 | 6 | Single | 10 | Geom-based | GRD | DarcyTools | | |
| S5 | 6 | Single | 5 | Geom-based | GRD | DarcyTools | | |
| 1.25 | 6 | Single | 1.25 | tbf | Geom-based | GRD | DarcyTools | |
| DG1 | 6 | Multi 3 | [3.9–15.6] | Geom-based | GRD | DarcyTools | | |
| DG2 | 6 | Multi 4 | [1.95–5.6] | Geom-based | GRD | DarcyTools | | |
| S50 | 6 | Single | 50 | Flow-based | GRD | DarcyTools | | |
| S10 | 6 | Single | 10 | Flow-based | GRD | DarcyTools | | |
| S5 | 6 | Single | 5 | Flow-based | GRD | DarcyTools | | |
| DG1 | 6 | Multi 3 | [3.9–15.6] | Flow-based | GRD | DarcyTools | | |
| DG2 | 6 | Multi 4 | [1.95–15.6] | Flow-based | GRD | DarcyTools | | |
| S50 | 6 | Single | 50 | Flow-based | GRD | HDG-tet | | |
| S10 | 6 | Single | 10 | Flow-based | GRD | HDG-tet | | |
| S5 | 6 | Single | 5 | Flow-based | GRD | HDG-tet | | |
| S50 | 6 | Single | 50 | Flow-based | PRM | HDG-hex | x | |
| S10 | 6 | Single | 10 | Flow-based | PRM | HDG-hex | x | |
| S5 | 6 | Single | 5 | Flow-based | PRM | HDG-hex | x | |
| DG1 | 6 | Multi 3 | [3.9–15.6] | Flow-based | PRM | HDG-hex | x | |
| DG2 | 6 | Multi 4 | [1.95–15.6] | Flow-based | PRM | HDG-hex | x | |
| S50 | 2 | Single | 50 | Flow-based | PRM | HDG-hex | x | |
| S10 | 2 | Single | 10 | Flow-based | PRM | HDG-hex | x | |
| S5 | 2 | Single | 5 | Flow-based | PRM | HDG-hex | x | |
| S50 | 6 | Single | 50 | Flow-based | GRD | HDG-hex | | |
| S10 | 6 | Single | 10 | Flow-based | GRD | HDG-hex | | |
| S5 | 6 | Single | 5 | Flow-based | GRD | HDG-hex | |
Cases for the DFN realization 1.
The distributions of local porosities for single-resolution grids (Figures 4c,d) show an increase in cell porosity values when increasing the resolution. As the resolution increases, the number of cells increases, the proportion of empty cells increases, and the number of fractures per non-empty cell decreases, as the average spacing between fractures and the size of each fracture become much larger than the size of the cell itself. At the extreme, cells may be small enough to contain no more than one fracture (or two at fracture intersections), allowing the definition of a maximum end-member porosity per cell (corresponding to the ratio of a single through-going fracture volume to the cell volume). Conversely, at coarser resolutions, cell sizes exceed both fracture spacing and fracture size. For these conditions, the cells are active only when they contain connected fracture clusters. As for hydraulic conductivity, replacing the constant transmissivity model (Figure 4c) with the distributed transmissivity model (Figure 4d) leads to a broader porosity distribution.
4.2 Definition of indicators from large-domain conditions
Several flow and transport indicators are defined to evaluate the performance of the ECM transformation under different cases. All indicators are derived from large-scale simulations, after the computation of domain-scale flow fields under GRD or PRM hydraulic boundary conditions. Indicators are evaluated either locally, at the cell scale, or globally, across the faces of the entire domain.
At the domain scale, equivalent permeability components are defined from the total flow rates crossing the faces of the large domain after computing the flow field under GRD hydraulic boundary conditions. This definition is similar to the one used at the ECM cell scale (see Section 2.2), but applied at the domain level and consistently for both DFN and ECM representations. Six domain-scale equivalent permeability components are hence defined (see Section 4.3.2).
For the local-scale evaluation, a set of measurement probes is distributed throughout the entire domain. Each probe consists of a delimited planar surface (e.g., a square) through which local flow () is measured as the volumetric flow rate crossing the probe surface, divided by the probe area and head gradient (yielding a quantity with units in ). Identical probes are used in both DFN and ECM simulations. They are superimposed on the grid such that, in each cell, three probes are centered at the cell center and oriented successively along the , and directions, corresponding to the cell-face orientations and GRD head gradients.
In addition to flow-based indicators, transport performance is evaluated using domain-scale porosity and several key components of the particle arrival time distributions. These include the first arrival time, mode, mean arrival time, and tailing exponent characterizing late-time arrivals ().
4.3 Flow numerical experiments under GRD conditions
We first evaluate the ability of the ECM transformation to reproduce the flow-rate structure, prior to analyzing transport behavior. This comparison is based on domain-scale GRD flow simulations and analyses performed at local and domain scales. In total, nearly 130 different configurations are defined (each tested under three to five large-scale hydraulic boundary conditions), combining ten DFN geometric realizations, two transmissivity models, five ECM grids, two upscaling methods (geom- and flow-based), and three directions for the imposed GRD head gradient. First, we focus on local flow-rate indicators (Section 4.3.1), before considering their large-scale implications using domain-scale indicators (Section 4.3.2).
4.3.1 Distribution of local flow rates
A systematic deviation is observed between the local flow rate distributions obtained from the DFN and those obtained from the ECM. In general, ECM-derived distributions are narrower and shifted toward larger values. This is illustrated in Figure 5a for a representative case (DFN realization 1, transmissivity model T1), using the single-resolution S5 grid with flow-based upscaling. The deviation is clearly visible, but the overlap between the two distributions remains substantial. A one-to-one comparison of local flow rates (Figure 5b) reveals the same systematic trend, but with considerable scatter. Notably, the smallest flow rates present in the DFN are not reproduced in the ECM. Similar deviations were observed across all cases, with magnitudes that vary depending on the ECM parameters and the DFN realization.
FIGURE 5
A one-to-one comparison of local flow rates obtained using flow-based and geom-based upscaling is shown in Figure 6 for the same DFN realization (realization 1) and the DG1 multiresolution ECM grid. The geom-based ECM flow rates tend to be larger than the flow-based ones.
FIGURE 6
For a given resolution, the ECM transformation tends to overestimate local flow rates, with a stronger overestimation associated with the geom-based upscaling. This is summarized in Equation 6:
The effect becomes even more pronounced as the resolution coarsens, a point to which we return when discussing domain-scale behavior.
4.3.2 Domain-scale values
The large-scale flow indicators compiled here are summarized in Figures 7–9. For each flow simulation, six equivalent domain scale permeability components are derived from the flow rates across the domain faces (Figure 2): the diagonal components, denoted (i.e., , , ), and the off-diagonal components, denoted (i.e., , , ).
FIGURE 7
Figure 7 compares the DFN and ECM (S50 grid) equivalent conductivities for the ten DFN realizations. Variability between realizations can reach a factor of 2.5 for the DFN and the ECM. The vertical conductivity () is also consistently lower than the horizontal components ( and ). Overall, the ECM conductivities overestimate the DFN ones by a factor of roughly 2–3 in all directions and for all ten realizations.
The DFN and ECM values are averaged over the ten DFN realizations and represented as boxplots in Figures 8, 9, sorted by increasing characteristic ECM cell size. For single-resolution grids, this characteristic size is constant, while for multiresolution grids it is taken as the minimum cell size. Figure 8 shows that permeability is increasingly overestimated as cell size grows. Whereas results are comparable for the finest resolution tested, they diverge markedly towards coarser ECM resolutions. For flow-based upscaling, the overestimation is only weakly sensitive to cell size and stays within a factor of 1.8–2.5. For geom-based upscaling, however, the deviation significantly increases with cell size and reaches up to an order of magnitude. When the transmissivity model is switched from constant values (T1) to a size- and orientation-dependent formulation, the same trends persist, although they are slightly attenuated. The geom-based approach implicitly assumes that all fractures fully penetrate grid cells and contribute to connectivity and fluid flow. This assumption is valid only for quasi-infinite fractures, and not for fractures of finite size distributions or for resolutions coarser than fracture size and average spacing. As cell size increases, this assumption is increasingly violated, leading to progressively larger overestimations of upscaled hydraulic conductivity. No such assumption is required for the flow-based approach.
FIGURE 8
FIGURE 9
Finally, we compare domain-scale directional permeabilities in the , and directions to assess whether the simplified geom-based upscaling affects the anisotropy structure of the permeability field and the prediction of large-scale flow rates. For this purpose, an anisotropy indicator is defined as the ratio of the maximum to the minimum equivalent permeability (among the three values and ). The results are plotted in Figure 9. The DFN anisotropy ratio (left column in Figure 9) reflects the anisotropy of the DFN geometrical structure, with several orientation sets of various weight (see section 7.1) that will also interact with the fracture transmissivity law (orientation and size dependency of the fracture transmissivities). The DFN anisotropy ratio lies between about 1 and 2.2 for T1, and 2–10 for TSL. For the ECM configurations the anisotropy ratio is similar to the DFN ratio only for the flow-based upscaling. It is resolution dependent, with an increasing underestimation toward coarser resolution, for the geom-based upscaling.
4.4 Results for flow experiments under changing hydraulic boundary conditions
Before moving on to the transport analysis, in which conventional PRM conditions are chosen to facilitate the interpretation of arrival-time distributions, we evaluate whether changes in the hydraulic boundary conditions (GRD versus PRM) affect macroscopic parameters relevant for transport properties. Under PRM conditions, a single direction is chosen for the external head gradient (e.g., the direction), with fixed heads applied on the two domain faces normal to this direction. In contrast to GRD conditions, no-flow boundaries are imposed on the remaining faces (those normal to the and directions), such that only two faces are hydraulically active and used to compute connectivity clusters.
One DFN realization (realization 1) and one transmissivity model () are selected for this comparison. First, the previously computed flow-based ECM grids (S5, S10, and S50) are reused for flow simulations under PRM conditions in the direction. The results are shown in Figure 10 (empty disc symbols labelled ECM from DFN6 in the legend) and summarized in Table 7. As previously observed under GRD conditions, the ECM overestimates flow rates, and the magnitude of this overestimation increases with cell size. The observed deviation is, however, larger under PRM conditions; for example, for grid S5 the estimation reaches a factor of approximately 4, compared with about 2 under GRD conditions (see Figure 8).
TABLE 7
| Model | Connected cluster | HBC | ||||
|---|---|---|---|---|---|---|
| | | | ||||
| DFN | 2 | PRM | 294 | 126 | ||
| S5 | 6 | PRM | 1,164 | 41.5 | ||
| S10 | 6 | PRM | 2,121 | 22 | ||
| S50 | 6 | PRM | 3,670 | 9.4 | ||
| S5 | 2 | PRM | 897 | 41.9 | ||
| S10 | 2 | PRM | 1,318 | 27.3 | ||
| S50 | 2 | PRM | 2,423 | 11.5 |
Summary of domain scale key indicators for DFN realization 1 and PRM flow conditions.
FIGURE 10
This difference can be explained by the underlying connectivity structure of the DFN. Due to corner fractures (see schematic in Figure 10), the GRD DFN flow rate must always be greater than or equal to the PRM flow rate under similar head gradient conditions. Indeed, corner fractures, or corner-clusters of fractures, bring additional paths between adjacent faces of the domain which potentially contribute to increasing all the flow rate components (, see Figure 2). In the DFN example shown in Figure 11, the GRD flow rate is approximately 1,300 (dash-dotted line), whereas the corresponding PRM flow rate is close to 300 (dashed line). If such corner-connected paths are not filtered out during ECM construction, they artificially increase hydraulic conductivity in certain cells and may remain active flow paths in the ECM, even though they are inactive in the DFN under PRM conditions.
FIGURE 11
To investigate this effect, a second set of ECMs is generated from the same DFN, but constructed using the DFN connectivity obtained under PRM conditions. In this case, corner-connected paths are removed, while fracture clusters connecting the two opposite active faces of the domain are retained. This ECM configuration is labelled ECM from DFN2 in Figure 10. Under these conditions, the observed deviation between DFN and the ECM flow rates is reduced from a factor of approximately 4 to about 3 (still larger than the factor of ∼2 under GRD conditions; see Figure 8).
The question of whether a unique ECM can be defined independently of boundary conditions is discussed further below. For the transport simulations, we have adopted PRM conditions, as they allow a consistent definition of entry and exit faces for all particles.
4.5 Results for transport experiments under PRM conditions
Next, we analyze the results of transport simulations to evaluate how the simplifications of the ECM approach – namely, the cell size of single-resolution grids – impact solute transport. PRM hydraulic boundary conditions are applied in the x direction, and the flow is simulated using DFN. lab/HDG-hex. For each DFN or ECM configuration, travel-time distributions are calculated successively under PRM conditions along the x, y, and z directions. For each test, hundreds of thousands of particles are injected at the inlet boundary, and the travel time and distance are recorded when a particle reaches the opposite boundary. Particle injection at the inlet is flux-weighted to avoid non-stationary effects (). Travel-time distributions (breakthrough curves) and transport indicators (Section 4.2) are then derived and analyzed. Five different realizations of the DFN recipe (realizations 1–5), combined with the two transmissivity distributions ( and ), are analyzed.
Figure 12 shows the results for a DFN realization (realization 1) with constant transmissivity (). The resulting PDF of arrival times (black line in Figures 12a,b) fits well with the Stretched Inverse Gamma (SIG) model defined by . The tail of the PDF at long travel times follows a power law, with an exponent (the parameter of the SIG) equal to 2.7. The DFN transport porosity (column in Table 7) and geometrical porosity (column in Table 7) are then compared. The transport porosity is defined as the ratio between the measured total flow rate (column in Table 7) and the mean arrival time (column in Table 7) as indicated in Equation 5. The geometrical porosity is defined as the open volume of the connected cluster (i.e., for each fracture, surface times aperture) divided by the total volume. For DFN realization 1, the transport porosity reaches only two-thirds of the available (geometrical) porosity, indicating that one-third of the open volume is not visited by particles due to very low flow rates, mostly in dead-end structures. The DFN is next transformed into an ECM using both the connectivity-2 and connectivity-6 filters (prior to deriving ECM properties) and single-resolution grids S50 (Figure 12a) and S10 (Figure 12b). As expected with flow-based methods, porosity is lower than that of the DFN (see Section 2.2). It decreases as the grid cell size increases, approaching the transport porosity of the DFN (Table 7; Figure 13b). For all ECM simulations, the flow-based porosity is very similar to the transport porosity. This shows that the flow-based method accurately identifies the volume with flow. Whatever the grid resolution (i.e., 10 or 50), the travel-time distributions are similar for both filtered DFNs along the fastest paths, up to times slightly larger than the mean travel time. However, the PDF slope at long travel times is strikingly different. The ECM built from the connectivity-2 DFN displays more late arrival times and, consequently, a larger mean arrival time. This can be attributed to the fact that the connectivity-6 ECM contains more flow paths and is therefore more prone to local overestimations, leading to fewer slow paths. In contrast, the fastest paths likely occur away from no-flow boundaries and are therefore similar for both connectivity filtering configurations. Beyond the connectivity-2 versus connectivity-6 comparison, all ECM breakthrough curves show lower early arrival, peak, mean, and late arrival times than the corresponding DFN results.
FIGURE 12
FIGURE 13
The analysis is extended to a total of five DFN realizations (realizations 1–5), each successively transformed into ECMs using the three grids S5, S10, and S50, with only connectivity-6 retained for the pre-filtering step. This results in a total of 20 transport simulations. DFN and ECM equivalent permeabilities, expressed as large-scale flow rates, are shown in Figure 13a. Variations between DFN realizations can reach a factor of three. As discussed in Section 4.3.2, the ECM transformation results in an overestimation of the flow rate, with the magnitude increasing as grid resolution decreases and reaching up to a factor of ten for the coarsest grid. ECM porosity, which is calculated as local transport porosity, lies between geometric porosity of the network and transport porosity of the DFN (Figure 13b). For the largest cell size (S50), ECM porosity is approximately equal to the transport porosity of the DFN, reflecting the calculation method. For smaller cell sizes (S10 and S5), however, it is higher, indicating that the method identifies a larger volume with non-zero flow than actually exists in the DFN.
Details of the arrival-time distributions are shown in Figure 14. All travel-time distributions exhibit a peak (i.e., the PDF mode) at travel times 1.5–4 times greater than the first arrival time, as well as a long power-law tail at larger travel times (Figure 14c). This shape is well described by the Stretched Inverse Gamma (SIG) distribution, whose exponents and characterize the travel-time distribution (). The fitting parameters of the SIG model are indicators of the distribution dispersion. In particular, , the exponent of the long-tail power law, controls the contribution of the longest times to the mean time: the smaller is, the larger the contribution. Although the ECM arrival-time distributions are all systematically shifted towards shorter times, important nuances appear when comparing first arrival times, peak times, and mean times (Figures 14a,b). The ECM transformation reduces travel times for all transmissivity models and grid sizes, which is consistent with the overestimation of flow rates and the underestimation of porosity shown in Figure 13. Among the three characteristic times (first arrival, mode, and mean), the first arrival time in ECMs is close to that of the DFNs. In contrast, the mean arrival time shows the largest deviation from the DFN reference and can differ by up to a factor of ten for some models (Figure 14b). Simulations using the TSL transmissivity model exhibit greater variability in characteristic times and larger deviations from DFN results than those using the constant transmissivity model T1. The tail indicator is consistently lower for DFN simulations than for any ECM grid, although differences diminish as grid resolution increases (Figure 14c). Low values indicate a higher probability of long travel times and lead to a strong separation between the mean and the PDF mode, with mean-to-mode ratios of up to 7 for the lowest values (Figure 14c). As increases, this ratio approaches 1, indicating that the contribution of the long tail to the mean becomes negligible. This is quantitatively well captured by the SIG model (Figure 14d). For values smaller than 3, transport becomes anomalous, with dispersion that varies continuously with scale, unlike normal transport (; ). This regime applies to the DFN simulations and to some simulations performed on the finest ECM grid. All other cases exhibit normal transport behavior.
FIGURE 14
5 Discussion and conclusions
In this study, we analyzed the impact of transforming a DFN into an equivalent continuum model (ECM) for modelling flow and solute transport. The analysis was conducted in a specific–though still simplified–context inspired by fracture conditions at the Forsmark site, where fluid flow and solute transport occur exclusively within the fracture network and the surrounding matrix can be considered impervious. The geometrical structure of the DFN was successively combined with (i) a simple transmissivity model using a single constant value, to highlight spatial and structural effects, and (ii) a more realistic model in which transmissivity is a function of both fracture size and orientation, indirectly reflecting mechanical stresses. Two upscaling approaches were considered: a simplified geom-based method and a numerical flow-based method, in which local cell properties are directly informed by local flow and transport simulations. For both approaches, we quantified the impact of varying the ECM resolution by analyzing the sensitivity of local properties as well as local and global flow and transport indicators to the upscaling method and ECM grid size. Our results demonstrate that the transition from DFN to ECM systematically overestimates hydraulic conductivity, underestimates the geometric porosity to a lesser extent, and overestimates the transport porosity. This also artificially reduces the variability and tortuosity of flow paths, except at very high resolutions. Furthermore, the ECM approach increases velocities at the domain scale and shortens arrival times, while significantly narrowing the distribution of late arrivals, thereby reducing the power-law tail exponent. First arrival times are relatively insensitive to ECM parameters, whereas mode and mean arrival times show strong sensitivity to both grid resolution and upscaling method.
Both geom-based and flow-based upscaling methods have pros and cons that depend on DFN complexity, on the upscaling resolution scale relative to the DFN itself, and whether the resolution can be imposed or tuned for a given application. The geom-based approach is much quicker than the flow-based approach when generating an ECM grid at a given resolution, because the analytical calculations are simpler and meshing is not required. However, this apparent efficiency can come at the cost of accuracy if the chosen resolution is too coarse relative to the DFN complexity. Geom-based upscaling yields results comparable to the DFN only when cells are small enough, such that most non-void cells contain a maximum of one crossing fracture. In this case, using a multi-resolution grid–with resolution locally adapted to fracture presence and with empty cells removed–can limit computational capacity. Nevertheless, this configuration closely resembles the original DFN, as the entire fracture network is effectively reproduced within the grid. Although the flow-based upscaling can be more time-consuming during the grid definition phase–especially at very fine resolutions–it provides more accurate results across all ECM resolutions. This advantage is particularly important when grid cells contain clusters of connected fractures, a situation that commonly arises when the resolution is constrained by computational limitations and upscaling is required to reduce DFN complexity.
The representation of connectivity and flow paths emerges as a key issue in the transformation of DFN to ECM. At the cell scale, distinct pathways are effectively merged, artificially increasing connectivity and, consequently, the inferred effective hydraulic conductivity. This effect diminishes only in cells containing a single connected pathway, but persists as soon as multiple fractures and pathways are present. Except for the single-fracture-cell configuration, geom-upscaling is therefore more strongly affected by connectivity issues. In addition, ECM sometimes creates paths that do not exist in the original DFN, for example, through artificial connections between adjacent or opposite cell faces, resulting in an overestimation of total flow. It is possible to reduce this overconnectivity by identifying isolated DFN clusters associated with the boundary conditions applied to the model. This is particularly relevant when using permeameter-type conditions with impermeable lateral boundaries; however, the resulting ECM grid derived from a filtered DFN is then only suitable for the specific boundary conditions at hand.
Finally, several aspects warrant further investigation. When it is not possible to construct an ECM with one fracture per cell, the impact of the DFN characteristic scales (fracture size distributions, minimum and maximum fracture sizes, connectivity length scales, and average spacing) on permeability and porosity upscaling should be investigated further. In addition, under these broader conditions, the assumption adopted here of representing hydraulic conductivity using three predefined coefficients aligned with the Cartesian axes should be revisited, particularly with respect to the interplay between fracture connectivity, orientation, and preferential flow directions.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
CD: Conceptualization, Validation, Supervision, Methodology, Writing – original draft, Writing – review and editing. QC: Methodology, Conceptualization, Visualization, Investigation, Validation, Writing – review and editing. US: Conceptualization, Software, Investigation, Writing – review and editing. RLG: Writing – review and editing, Methodology, Validation, Investigation, Software. BP: Methodology, Investigation, Software, Validation, Writing – review and editing. J-OS: Supervision, Funding acquisition, Writing – review and editing, Writing – original draft, Conceptualization, Resources, Methodology. PD: Methodology, Writing – review and editing, Validation, Supervision, Conceptualization, Writing – original draft.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This research was funded by the Swedish Nuclear Fuel and Waste Management Company (SKB). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.
Conflict of interest
Authors CD, QC, RLG and BP were employed by Itasca, Factory. Author US was employed by Computer-Aided Fluid Engineering AB.
Author J-OS was employed by Swedish Nuclear Fuel and Waste Management Company (SKB).
The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Abbreviations
DFN, discrete fracture network; ECM, equivalent continuum model; HBC, hydraulic boundary conditions; HDG, hybrid discontinuous galerkin method; PDF, probability density function; PRM, PeRMeameter hydraulic boundary conditions; GRD, head GRaDient hydraulic boundary conditions; SIG, Stretched Inverse Gamma; DFN2, DFN connected cluster with PRM conditions; DFN6, DFN connected cluster with GRD conditions; S5, regular grid of 100x100x100 blocks with block size of 5m; S10, regular grid of 50x50x50 blocks with block size of 10m; S50, regular grid of 10x10x10 blocks with block size of 50m; DG1, multi-resolution grid of 1.1 million blocks with minimum size of 3.9m; DG2, multi-resolution grid of 6 million blocks with minimum size of 1.95m; T1, constant transmissivity (T) model with T = 1m2/s; TSL, variable transmissivity (T) model with T = f(σ,l), l the fracture length and σ the applied normal stress.
References
1
BaxterS.HartleyL.HoekJ.MyersS.TsitsopoulosV.WilliamsT. (2019). Upscaling of brittle deformation zone flow and transport properties, SKB Report R-19-01. Svensk Kärnbränslehantering AB.
2
BearJ. (1988). Dynamics of fluids in porous media. Courier Corporation.
3
BerkowitzB.ScherH. (1997). Anomalous transport in random fracture networks. Phys. Rev. Lett.79 (20), 4038–4041. 10.1103/PhysRevLett.79.4038
4
ChenM.BaiM.RoegiersJ.-C. (1999). Permeability tensors of anisotropic fracture networks. Math. Geol.31 (4), 335–373. 10.1023/A:1007534523363
5
ChenT.ClauserC.MarquartG.WillbrandK.MottaghyD. (2015). A new upscaling method for fractured porous media. Adv. Water Resources80, 60–68. 10.1016/j.advwatres.2015.03.009
6
ChenT.ClauserC.MarquartG.WillbrandK.HillerT. (2018). Upscaling permeability for three-dimensional fractured porous rocks with the multiple boundary method. Hydrogeol. J.26 (6), 1903–1916. 10.1007/s10040-018-1744-z
7
CockburnB.DongB.GuzmánJ.RestelliM.SaccoR. (2009). A hybridizable discontinuous galerkin method for steady-state convection-diffusion-reaction problems. SIAM J. Sci. Comput.31 (5), 3827–3846. 10.1137/080728810
8
DarcelC.Le GocR.LavoineE.DavyP.Mas IvarsD.SykesE.et al (2024). Coupling stress and transmissivity to define equivalent directional hydraulic conductivity of fractured rocks. Eng. Geol.342, 107739. 10.1016/j.enggeo.2024.107739
9
DavyP.DarcelC.Le GocR.MunierR.SelroosI.-O.Mas IvarsD. (2018). DFN, why, how and what for, concepts, theories and issues. Proceedings of 2nd international discrete fracture network engineering conference, Seattle, USA, June 20-22 2018, 2018. American Rock Mechanics Association.
10
DavyP.Le GocR.DarcelC.SelroosJ. O. (2023). Scaling of fractured rock flow. Propos. Indicators Selection DFN Based Flow Models Comptes Rendus Géoscience355, 667–690. 10.5802/crgeos.174
11
DavyP.Le GocR.DarcelC.PinierB.SelroosJ.-O.Le BorgneT. (2024). Structural and hydrodynamic controls on fluid travel time distributions across fracture networks. Proc. Natl. Acad. Sci.121 (47), e2414901121. 10.1073/pnas.2414901121
12
DentzM.CortisA.ScherH.BerkowitzB. (2004). Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport. Adv. Water Resour.27 (2), 155–173. 10.1016/j.advwatres.2003.11.002
13
FollinS. (2008). Bedrock hydrogeology forsmark, site descriptive modelling, SDM-site forsmark, SKB report R-08-95. Svensk Kärnbränslehantering AB.
14
FollinS.LevénJ.HartleyL.JacksonP.JoyceS.RobertsD.et al (2007). Hydrogeological characterisation and modelling of deformation zones and fracture domains, forsmark modelling stage 2.2, R-07-48. Svensk Kärnbränslehantering AB.
15
HoteitH.FiroozabadiA. (2008). Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures. Adv. Water Resour.31 (1), 56–73. 10.1016/j.advwatres.2007.06.006
16
HymanJ. D.KarraS.MakedonskaN.GableC. W.PainterS. L.ViswanathanH. S. (2015). dfnWorks: a discrete fracture network framework for modeling subsurface flow and transport. Comput. Geosciences84, 10–19. 10.1016/j.cageo.2015.08.001
17
JacksonC. P.HochA. R.TodmanS. (2000). Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium. Water Resour. Res.36 (1), 189–202. 10.1029/1999WR900249
18
LangP.PalusznyA.ZimmermanR. (2014). Permeability tensor of three‐dimensional fractured porous rock and a comparison to trace map predictions. J. Geophys. Res. Solid Earth119 (8), 6288–6307. 10.1002/2014JB011027
19
LeiQ.LathamJ.-P.TsangC.-F. (2017). The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks. Comput. Geotechnics85, 151–176. 10.1016/j.compgeo.2016.12.024
20
LiZ.NguyenS. (2025). Modelling flow and transport in fractured crystalline rocks by an upscaled equivalent continuous porous media method. Geomechanics Energy Environ.41, 100625. 10.1016/j.gete.2024.100625
21
LongJ. C. S.RemerJ. S.WilsonC. R.WitherspoonP. A. (1982). Porous media equivalents for networks of discontinuous fractures. Water Resour. Res.18 (3), 645–658. 10.1029/WR018i003p00645
22
NeumanS. P. (2005). Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeol. J.13 (1), 124–147. 10.1007/s10040-004-0397-2
23
OdaM. (1985). Permeability tensor for discontinuous rock masses. Geotechnique35 (4), 483–495. 10.1680/geot.1985.35.4.483
24
PainterS.CvetkovicV. (2005). Upscaling discrete fracture network simulations: an alternative to continuum transport models. Water Resour. Res.41, W02002. 10.1029/2004WR003682
25
PouyaA.CourtoisA. (2002). Définition de la perméabilité équivalente des massifs fracturés par des méthodes d'homogénéisation. Comptes Rendus Geosci.334 (13), 975–979. 10.1089/neu.2016.4789
26
PouyaA. (2005). Tenseurs de perméabilité équivalente d'un domaine hétérogène fini. Comptes Rendus Geosci.337 (6), 581–588. 10.1016/j.crte.2005.02.002
27
RenardP.AbabouR. (2022). Equivalent permeability tensor of heterogeneous media: upscaling methods and criteria (review and analyses). Geosciences12 (7), 269. 10.3390/geosciences12070269
28
RenardP.de MarsilyG. (1997). Calculating equivalent permeability: a review. Adv. Water Resour.20 (5), 253–278. 10.1016/S0309-1708(96)00050-4
29
Sanchez-VilaX.GuadagniniA.CarreraJ. (2006). Representative hydraulic conductivities in saturated groundwater flow. Rev. Geophys.44 (3). 10.1029/2005RG000169
30
Sánchez‐VilaX.GirardiJ. P.CarreraJ. (1995). A synthesis of approaches to upscaling of hydraulic conductivities. Water Resour. Res.31 (4), 867–882. 10.1029/94WR02754
31
SchöberlJ. (2014). C++ 11 implementation of finite elements in NGSolve. Institute for Analysis and Scientific Computing, Vienna University of Technology.
32
SelroosJ.-O.Mas IvarsD.MunierR.HartleyL.LibbyS.DarcelC.et al (2022). Methodology for discrete fracture network modelling of the forsmark site. Volume I: concepts, data and interpretation methods, SKB R-20-11. Svensk Kärnbränslehantering AB.
33
SKB (2008). Site description of forsmark at completion of the site investigation phase SDM-site forsmark. SKB TR-08-05. Svensk Kärnbränslehantering AB.
34
SvenssonU. (2001). A continuum representation of fracture networks. Part I: method and basic test cases. J. Hydrology250, 170–186. 10.1016/S0022-1694(01)00435-8
35
SvenssonU.FerryM. (2014). DarcyTools: a computer code for hydrogeological analysis of nuclear waste repositories in fractured rock. J. Appl. Math. Phys.2, 365–383. 10.4236/jamp.2014.26044
36
SvenssonU.KuylenstiernaH.-O.FerryM. (2010). DarcyTools version 3.4-Concepts, methods and equations, SKB R-07-38. Svensk Kärnbränslehantering AB.
37
WenX.-H.Gómez-HernándezJ. J. (1996). Upscaling hydraulic conductivities in heterogeneous media: an overview. J. Hydrology183 (1), ix–xxxii. 10.1016/S0022-1694(96)80030-8
38
WillmannM.CarreraJ.Sánchez‐VilaX. (2008). Transport upscaling in heterogeneous aquifers: what physical parameters control memory functions?Water Resour. Res.44, W12437. 10.1029/2007WR006531
39
ZhangX.SandersonD. J. (2002). Numerical modelling and analysis of fluid flow and deformation of fractured rock masses. Elsevier.
Summary
Keywords
DFN, ECM, flow, fracture, transport, upscaling
Citation
Darcel C, Courtois Q, Svensson U, Le Goc R, Pinier B, Selroos J-O and Davy P (2026) Transformation of discrete fracture networks into equivalent continuum models for sparsely fractured rocks: comparing flow-based and geometry-based upscaling for flow and transport. Front. Nucl. Eng. 5:1760157. doi: 10.3389/fnuen.2026.1760157
Received
03 December 2025
Revised
22 January 2026
Accepted
09 February 2026
Published
22 April 2026
Volume
5 - 2026
Edited by
Yuankai Yang, Forschungszentrum Juelich, Germany
Reviewed by
Tao Chen, Shandong University of Science and Technology, China
Zhenze Li, Canadian Nuclear Safety Commission (CNSC), Canada
Updates
Copyright
© 2026 Darcel, Courtois, Svensson, Le Goc, Pinier, Selroos and Davy.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Caroline Darcel, c.darcel@itasca.fr
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.