**Research article**
02 Jan 2019

**Research article** | 02 Jan 2019

# Stochastic modeling of flow and conservative transport in three-dimensional discrete fracture networks

I-Hsien Lee Chuen-Fa Ni Fang-Pang Lin Chi-Ping Lin and Chien-Chung Ke

^{1,2},

^{1,2},

^{3},

^{1},

^{4}

**I-Hsien Lee et al.**I-Hsien Lee Chuen-Fa Ni Fang-Pang Lin Chi-Ping Lin and Chien-Chung Ke

^{1,2},

^{1,2},

^{3},

^{1},

^{4}

^{1}Graduate Institute of Applied Geology, National Central University, Taoyuan city, 32001, Taiwan^{2}Center for Environmental Studies, National Central University, Taoyuan City, 32001, Taiwan^{3}Data Computing Division, National Center for High Performance Computing, Hsinchu City, 30076, Taiwan^{4}Sinotech Engineering Consultants Inc., Taipei City, 11494, Taiwan

^{1}Graduate Institute of Applied Geology, National Central University, Taoyuan city, 32001, Taiwan^{2}Center for Environmental Studies, National Central University, Taoyuan City, 32001, Taiwan^{3}Data Computing Division, National Center for High Performance Computing, Hsinchu City, 30076, Taiwan^{4}Sinotech Engineering Consultants Inc., Taipei City, 11494, Taiwan

**Correspondence**: Chuen-Fa Ni (nichuenfa@geo.ncu.edu.tw)

**Correspondence**: Chuen-Fa Ni (nichuenfa@geo.ncu.edu.tw)

Received: 24 Jul 2018 – Discussion started: 21 Aug 2018 – Revised: 05 Dec 2018 – Accepted: 06 Dec 2018 – Published: 02 Jan 2019

This study presents the stochastic Monte Carlo
simulation (MCS) to assess the uncertainty of flow and conservative transport
in 3-D discrete fracture networks (DFNs). The MCS modeling workflow involves
a number of developed modules, including a DFN generator, a DFN mesh
generator, and a finite element model for solving steady-state flow and
conservative transport in 3-D DFN realizations. The verification of the
transport model relies on the comparison of transport solutions obtained from
HYDROGEOCHEM model and an analytical model. Based on 500 DFN realizations in
the MCS, the study assesses the effects of fracture intensities on the
variation of equivalent hydraulic conductivity and the exhibited behaviors of
concentration breakthrough curves (BTCs) in fractured networks. Results of
the MCS show high variations in head and Darcy velocity near the specified
head boundaries. There is no clear stationary region obtained for the head
variance. However, the transition zones of nonstationarity for *x*-direction
Darcy velocity is obvious, and the length of the transition zone is found to
be close to the value of the mean fracture diameter for the DFN realizations.
The MCS for DFN transport indicates that a small sampling volume in DFNs can
lead to relatively high values of mean BTCs and BTC variations.

Successful characterizations of flow and contaminant transport in fractured geologic formations depend on adequate descriptions of complex geometrical structures, which comprise a wide variety of fractures and their connections (Ahmed et al., 2015; Pichot et al., 2012; Weng et al., 2014). The fracture characteristics can be quantified by using various statistical parameters, including the fracture orientation, length, shape, and permeability alongside the fracture intensity and connectivity (Bonnet et al., 2001; Botros et al., 2008; Bour et al., 2002; Koike et al., 2015; Stephens et al., 2015). These commonly used parameters represent fracture networks at sites of interest and bridge gaps between limited field observations and flow and transport implementations for site-specific issues.

Using the discrete fracture network (DFN) approach to characterize the flow and transport in fractured media is a challenging task for practical applications. Intensive research over the past 3 decades has led to the development of numerous models that are based on the DFN approach to model the flow or transport in fractured formations (Cacas et al., 1990a, b; de Dreuzy et al., 2013; Hyman et al., 2015a; Liu and Neretnieks, 2006; Long et al., 1985; Pichot et al., 2012; Xu and Dowd, 2010). Advanced 3-D DFN approaches typically include procedures of fracture generation, DFN meshing, flow and transport, or particle tracking (de Dreuzy et al., 2013; Erhel et al., 2009; Hyman et al., 2014; Pichot et al., 2012; Xu and Dowd, 2010; Zhang, 2015; Trinchero et al., 2016; Fourno et al., 2019). Particle tracking algorithms are usually preferred to simulate DFN transport and have recently been widely implemented to evaluate the time resistance of contaminants for fractured formations (Hyman et al., 2015a, b; Makedonska et al., 2015; Painter et al., 2008; Stalgorova and Babadagli, 2015; Wang and Cardenas, 2015). The objective of the Lagrangian approach is to avoid numerical difficulties in solving the advection dispersion equation (ADE) in complex DFN domains. Such DFN transport models use the released particles to represent the contaminant with a specified mass or concentration. Many previous studies have discussed the issues in treating particle movement in fracture networks (Hyman et al., 2015a, b; Johnson et al., 2006; Makedonska et al., 2015; Painter et al., 2008; Park et al., 2003; Wang and Cardenas, 2015; Zafarani and Detwiler, 2013).

Over the years, many studies have focused on developing flow and transport models and integrating DFN simulation workflows for 3-D fracture networks (Hyman et al., 2014, 2015a; Lee and Ni, 2015). Specifically, the DFN transport was mainly modeled based on Lagrangian approaches such as particle tracking and random walk algorithms (e.g., Makedonska et al., 2015; Painter et al., 2008; Stalgorova and Babadagli, 2015; Wang and Cardenas, 2015). Numerical solutions to the ADE based on the Eulerian approach have not been widely implemented because of computational issues, such as numerical dispersion and convergence in the model for complex fracture connections (Odling, 1997; Berrone et al., 2018).

With the advantages of computational technologies, the stochastic modeling of flow and Eulerian-based transport in 3-D DFNs has become a feasible task. It is an important issue to quantify flow and transport uncertainties based on available DFN properties. The objectives of this study are to develop and implement numerical models for stochastic modeling of flow and conservative transport in 3-D DFNs. The stochastic Monte Carlo simulation (MCS) is conducted to assess the flow and transport uncertainty induced by the 3-D DFNs. In this study, we first assess the developed ADE model by comparing the solutions of simple porous fractures with those from the HYDROGEOCHEM finite element model (Yeh et al., 2004) and the analytical model developed by Wexler (1992). Then, we use the MCS to evaluate the equivalent hydraulic conductivity for specified DFN statistical parameters. The collected flow and transport realizations enable the analyses of flow and transport uncertainties in the fractured simulation domain. The simulation results are expected to provide general insight into the evaluations of flow and transport uncertainty based on the available DFN geometrical properties.

In this study, the fractures in a DFN are considered to be porous media with impermeable surfaces that are connected to the formation matrix. The two impermeable surfaces of a fracture are considered to be two rough parallel plates that enable fluids to pass through the fracture at a relatively high velocity (e.g., Kwicklis and Healy, 1993; Lee and Ni, 2015; Pruess and Tsang, 1990). The following presents the mathematical formulas, a brief description of the mesh generation, and the finite element models for simulating the 3-D DFN flow and transport.

## 2.1 Flow and transport equations

The mathematical formulation for the DFN consists of flow and transport in a set of 2-D porous fracture plates connected in a 3-D domain. The coupling of flow and transport in porous media has been widely investigated in fields that are related to groundwater hydrology (Dagan, 1989; Hartley and Joyce, 2013; Yeh et al., 2004; Zheng and Bennett, 2002). Based on the concept of mass conservation and Darcy's law, the equations for solving the steady-state and depth-averaged hydraulic head for 2-D porous fractures can be expressed as

subject to the boundary conditions

and

where *h*(*x*) is the hydraulic head, *K*(*x*) is the hydraulic conductivity,
and *b*(*x*) is the aperture for fractures. For saturated flow, the locations
of the fracture has been taken into account in Eq. (1). We assumed that the
flows in fractures were parallel to the fracture and the relatively small
fracture apertures insignificantly influenced the vertical flow in the 2-D
fractures. The notation *Q*(*x*) represents the sources or sinks applied to
the 2-D porous fractures. The Cartesian coordinate *x* (*x*=1 and 2)
represents the *x* and *y* directions in the horizontal modeling domain of
the fractures. Moreover, *h*_{D} represents the prescribed head at
the Dirichlet boundary Γ_{D}, and *q*_{N} is the
specific flux at the Neumann boundary Γ_{N}. The notation
** n** is a unit vector that is normal to the boundary Γ

_{N}. In this study, we consider the DFNs to be fully saturated. The

*z*coordinates for fractures represent the elevation heads, which have been considered in the calculations of hydraulic heads. Therefore, the solutions of hydraulic heads in 2-D fractures in the 3-D domain can be obtained from Eqs. (1) to (3).

Similar to the flow simulation, the depth-averaged conservative solute transport equation for saturated fractured porous media is governed by the ADE and can be written as (e.g., Dagan, 1989; Ni et al., 2009; Zheng and Bennett, 2002)

subject to the initial and boundary conditions

and

where *c*(*x*,*t*) is the volumetric solute concentration that is measured in
the liquid phase, and *Q*_{c}(*x*,*t*) represents the rate where the
volumetric solute concentration is injected (source) or extracted (sink) from
the DFN. The notation $v\left(x\right)=-K\left(x\right)\mathrm{\nabla}h\left(x\right)/n\left(x\right)$ is the seepage
velocity, and *n*(*x*) is the effective porosity in the porous fractures.
Calculating the seepage velocities at nodes relies on the obtained hydraulic
heads at element centers in the DFN (i.e., the solution of Eqs. 1 to 3). This
study used an improved approach proposed by Yeh (1981) to obtain the seepage
velocity at nodes to consider the global mass conservation of flow in the
simulation domain. In Eqs. (5) to (7), *c*_{0} represents the initial
concentration in the entire modeling domain Ω, *c*_{D} is the
specified concentration at the Dirichlet boundary, and *q*_{c} is the
dispersive flux at the Neumann boundary. Moreover, *D*_{e}(*x*) in
Eq. (4) is considered the macro-dispersion coefficient, which is evaluated
based on the seepage velocity (Zheng and Bennett, 2002)

where *δ*_{ij} is the Kronecker delta and *α*_{L}=*α*_{L}(*x*) is the longitudinal dispersivity in the principal
flow direction. *α*_{T}=*α*_{T}(*x*) represents the
transverse dispersivity, which is perpendicular to the longitudinal
dispersivity. Notations *v*_{i} and *v*_{j} are the seepage velocities in
different directions in the porous fractures, and $\stackrel{\mathrm{\u203e}}{v}$ represents
the magnitude of the seepage velocity. In Eq. (8), *D*_{0}(*x*) is the
effective molecular diffusion coefficient. In this study, *D*_{0}(*x*) in
our model is assumed to be negligible, which implies that the dispersion is
dominated by advective transport and mechanical dispersion (Zheng and
Bennett, 2002; Bear and Cheng, 2010).

## 2.2 DFN connections and the unstructured mesh generation

The information of fracture orientations enables the direct simulation of flow and transport in a series of 2-D fractures and the reproduction of the flow and transport behaviors for a 3-D DFN system. This study defines a DFN without isolated fractures as an effective discrete fracture network (DFNe). Figure 1 shows the definitions of the individual fracture in a 3-D DFNe. Based on the long axis of an elliptical porous fracture, the positive trend and plunge angles are defined as clockwise from the north and downward from the horizontal plane, respectively. In this study, the intersections of a fracture and the simulation boundaries have to be identified (Fig. 1) before the mesh generation is implemented. We generate the fracture length for the long and short axes of each fracture in the 3-D DFN based on the uniform distribution. The larger value of the two generated radii is used to obtain the long axis of the elliptical fractures. In addition, isolated fractures that are not connected to other fractures and simulation boundaries are removed for computational efficiency.

Figure 2 shows an example of a DFNe connection for two intersected fractures. Figure 2a shows the generated mesh in the 3-D domain for two intersecting fractures. Mesh generation begins with the generation of initial fracture meshes for each fracture plate (i.e., Fractures 1 and 2 in Fig. 2c and e, respectively). Figure 2b displays two intersecting fractures that were individually rotated back to the 2-D horizontal plane. In Fig. 2b the plunge and trend values differ for each fracture in a 3-D DFNe, so the fractures in a 2-D coordinate system might not overlap. The intersections for each fracture are also located in different areas (see Fig. 2b). However, the length of the intersection should be identical to that of the intersecting fractures. The fracture intersections and simulation boundaries are recorded for our unstructured mesh generation model.

The mesh generation starts with generating initial mesh for each fracture. The mesh generator allows users to define intervals of mesh boundary nodes. In this study, the Delaunay triangulation algorithm is used for generating the initial meshes. The special treatment of fracture intersections rely on the boundary recovery algorithm. In the process of boundary recovery for fracture intersections, we allow the node interval be reduced with a defined ratio according to the smallest node interval along the edges of the connected fractures. The ratio can be 1∕10 or a smaller value, depending on the problem. The details of the mesh generation algorithm for connected fractures are available in the study of Lee and Ni (2015).

## 2.3 Numerical finite element solution to the DFNe

To solve the governing equations of flow and transport for the DFNe framework, we employ the Galerkin finite element method and the biconjugate gradient matrix solver to solve Eqs. (1) and (4). The linear function for hydraulic heads and the concentration at the nodes surrounding an element of the triangular mesh system can be represented with

and

respectively, where the notation *N*_{α}(*x*) is the shape function determined by
the mesh in fractures in the 2-D porous fracture plates in a 3-D DFNe, and
the shape function has the following formula:

where *a*_{0} and *a*_{i} are coefficients of the shape functions. In this
study, the coefficients were determined based on the following formulas as
applied to nodes in a triangular element:

where $J=\mathrm{2}\times \phantom{\rule{0.125em}{0ex}}A$, and *A* represents the area of an element. Notation
*x*_{ij} represents the coordinate value in the *i* direction at the *j*th
node in a triangular element. In this study, the solution to the DFNe
followed similar processes to what are used in classical finite element
methods. However, the system of equations to be solved in the DFNe is
different and relatively complex compared to a typical 2-D problem, because
the nodes along the fracture intersections in the solution processes could
introduce additional terms in the coefficient matrices for the system of
equations. The roles of the additional terms in the coefficient matrices are
to build connections between fractures through the element nodes along the
lines of intersection. More intersections in the DFNe would yield more complex
coefficient matrices for the system of equations. Figure 2c–e also
demonstrate an example for the connection of nodes and elements along a
fracture intersection. Suppose that the Fractures 1 and 2 in Fig. 2c and
e have a line of intersection (dashed line). Elements 1 to 5 in
Fracture 1 and elements 6 to 10 in Fracture 2 share the same nodes (nodes 48, 49,
and 50) along the line of intersection. Let us focus on the node 49. For all the
elements in Fracture 1 and Fracture 2 that are connected to the node 49, the
calculation of the coefficient matrix for node 49 must rely on integrating
weightings from shape functions of the connected elements and the associated
nodes. This inclusion of all nodal information in the matrix system for heads
could resolve the detailed behavior of flow in elements that are connected to
the fracture intersections. The mass flux of concentration near the
intersections follow a similar procedure to build the coefficient matrices
for the ADE solutions.

The features of the HYDROGEOCHEM model are not for DFN flow and transport modeling. For simple cases such as a single horizontal fracture plate or cross-shaped porous fracture, it is possible to simulate a fracture and matrix system using the HYDROGEOCHEM model if one can use small mesh sizes to resolve the fracture apertures and matrix system. In addition, the differences of the hydraulic conductivity between the fracture and matrix need to be large enough to minimize the influence of the matrix. Because the HYDROGEOCHEM was developed based on the finite element method, the numerical dispersion might be similar to the developed model in this study. This study further uses a two-dimensional analytical solution proposed by Wexler (1992) to conduct verification of the developed model. The comparison is limited to the case with advection and dispersion in a horizontal porous fracture plate.

Based on the verified DFN flow and transport model, we then conducted 500 MCS
realizations to assess the upscaled flow behaviors with various fracture
intensities for 2-D profiles (i.e., *P*_{21}) and 3-D rock volumes (i.e.,
*P*_{32}). *P*_{21} and *P*_{32} are defined as the length of fracture traces per
unit area and the area of fractures per unit volume, respectively. The MCS
flow realizations are further used to assess the effects of DFN properties on
the flow and transport uncertainties. In this study, statistical structures
that are relevant to the distributions of the fracture properties included
Poisson and uniform distributions (Lee and Ni, 2015; Xu and Dowd, 2010).

## 3.1 Transport model verification by using HYDROGEOCHEM model

This study employs two cases in a 2 m × 2 m × 2 m
fractured rock domain (Fig. 3), including a horizontal fracture plate
(Case 1) and a cross-shaped fracture network (Case 2), to verify the
developed transport model by using the HYDROGEOCHEM model. The test cases
represent fracture sizes that are considerably larger than the controlled
modeling domain. Local mesh refinements in the HYDROGEOCHEM model are
required to resolve fractures in a rock volume to obtain fracture structures
in HYDROGEOCHEM (see Fig. 3a and c). However, the hydraulic conductivity
values are different to represent a horizontal fracture plate (Fig. 3a) and a
cross-shaped fracture network (Fig. 3c). For the two test cases, we assume a
uniform fracture aperture of 0.001 m. The fracture hydraulic conductivity is
1.0 m d^{−1} for the fracture plates in
the developed model. This fracture hydraulic conductivity value in
HYDROGEOCHEM model is also applied to the elements that represent the
fracture locations (Fig. 3a and c). The matrix of hydraulic conductivity in the
HYDROGEOCHEM model is 10^{−5} m d^{−1} to create a high variety of
hydraulic conductivity values for the fractures and rock matrix. In the test
cases, the effective porosity for the fractures is at a relatively large
constant value of 0.43. A relatively small isotropic dispersivity of 0.001 m
is used to evaluate the advection-dominated transport. The boundary
conditions along the boundaries parallel to the flow direction are specified
to be no-flow boundary conditions, except for the cross-shaped fracture
network case, where a slightly upward flow along the vertical fracture is
introduced (Fig. 3c and d). In Fig. 3c and d, the boundary conditions are
the same as in Fig. 3a and b, but a constant head of 9.01 m is specified on the top
side of the vertical fracture (*x*=1.0 m and *z*=1.0 m). Such a constant
head of 9.01 on the top side of the vertical fracture can produce stress for
upward flow in the vertical fracture. The Neumann boundary condition is
assigned in these cases as the transport boundary condition. In the test
cases, we release an initial Gaussian distribution plume in the horizontal
fracture plate. The time step for the transport simulation is 0.1 days
throughout the simulation time (5.0 days). Similar to the transport solution
in the developed model, the ADE solution technique in the HYDROGEOCHEM model
is the Eulerian-based approach for comparison purposes. Table 1 lists the
flow and transport parameters for the test cases.

^{a} The specified boundary conditions for HYDROGEOCHEM and the developed
model are applied to the intersection between fractures and the western (or
eastern) boundary of the simulation domain.
^{b} The specified boundary conditions for HYDROGEOCHEM and the developed
model are applied to the intersection between the vertical fracture and the
top boundary of the simulation domain.
^{c} The seepage velocity at each node is evaluated based on the Darcy
flux obtained at the node. n/a: not applicable.

## 3.2 Transport model verification by using analytical solution

The study of Wexler (1992) considers a horizontal two-dimensional domain. The
simulation domain is similar to the first case in Sect. 3.1. However, the
study of Wexler (1992) considers an *x*-direction uniform flow and constant
values of longitudinal and transverse dispersion coefficients in the
simulation domain. The transport process in the solution of Wexler (1992)
also involves the first-order decay. The transport equation has the following
formula:

with initial and boundary conditions

and

where the *c*(*x*,*t*) is the concentration, and *D*_{x} and *D*_{y} are the
longitudinal and transverse dispersion coefficients. The *v*_{x} is the
uniform seepage velocity in the *x* direction. Notation *n**λ* in Eq. (15)
is the first-order decay coefficient for the model. In the test example, the
specified concentration *c*_{D} is applied along the inlet boundary
(i.e., *x*=0) in the interval of *Y*_{1} and *Y*_{2}. To fit the condition
of the model in this study, we have neglected the decay process for
comparison purpose. This approximation yields the closed-form solution:

The calculation of Eq. (21) requires numerical approximations. The study of
Wexler (1992) suggested the Gauss–Legendre iteration algorithm to obtain the
solution. However, their results indicated that a small *x* value might lead
to numerical errors for the iterations.

In the test example (Case 3), the horizontal porous fracture plate has the
size of 2 m × 2 m. We follow the assumptions applied to the
solution in the study of Wexler (1992). The zero concentration is used as
the initial condition, and the concentration of 1.0 is specified in the
interval between *Y*_{1}=0.75 and *Y*_{2}=1.25 m along the inlet boundary
(i.e., *x*=0). The uniform seepage velocity in the *x* direction is
0.1 m d^{−1}. The longitudinal and transverse dispersivities for the case
are 0.1 m. Note that this dispersivity value was determined based on the
study of scale-dependent dispersivity proposed by Gelhar et al. (1992). The
interested scale for the study is on the order of 1 m.
Therefore, the isotropic dispersivity of 0.1 m was specified for the
transport simulation.

## 3.3 Flow upscaling behaviors

The equivalent hydraulic conductivity for a specified representative
elementary volume (REV) is the basis for conducting flow upscaling for
practical problems that cover simulation domains on the order of hundreds of
meters to several kilometers. Similar to the fractured rock volume in the
test case, we generate 500 DFN realizations to assess the flow and upscaling
behaviors for various fracture intensities and the associated fracture
properties (e.g., fracture locations, plunges and trends, and sizes). Table 2
lists the parameters for generating the DFN realizations. In this study, the
*P*_{21} values for different DFN realizations are calculated based on the
downstream boundary profile at *x*=2.0 m (Fig. 4). In the numerical
example, the criteria of the nodal spaces along the fracture boundaries are
fixed at 0.1 m, and the nodal spaces for fracture intersections have a
reduction rate of 1∕10 for all the DFNe realizations. Figure 4 also
shows the conceptual model alongside one DFNe realization and the associated
head field as a numerical example. Table 3 lists the hydrogeological
parameters for the flow and transport simulations.

The calculation of the equivalent hydraulic conductivity for a rock system considers the concept of mass conservation applied to a REV. The flow passing through the 3-D control volume can be represented by the following formula:

where *Q*_{frac} (*L*^{3}∕*T*) is the net flow rate for a specified
direction of a DFNe in a defined rock control volume. The hydraulic head
gradient of a control volume is defined as the ratio of the head difference
Δ*h* (*L*) to the flow length Δ*L* (*L*) between two boundaries
of the control volume. The notation *A* (*L*^{2}) represents the flow area
for the control volume. The equivalent hydraulic conductivity *K*_{eq} can be determined if the flow rate of a fractured rock is obtained from the
DFNe flow model.

## 3.4 Sampling volumes and flow and transport uncertainties

In this numerical example, we investigate the effects of averaging strategies
(i.e., along vertical lines or on profiles perpendicular to flows) on the
observations of breakthrough curves (BTCs) in fractured formations. This numerical example involves
using the random DFNe flow realizations from the example in Sect. 3.3.
Table 3 lists the parameters applied to the transport simulations. Figure 5
shows the boundary conditions for the transport simulations. Constant
concentration values of 1.0 and 0.0 (kg m^{−3}) are specified on
boundaries at *x*=0.0 and 2.0 (m), respectively. The constant concentration
values of 0.0 (kg m^{−3}) at the outflow can represent the scenario that a
large water body was connected to the fractured rock at the outflow boundary.
The remaining boundaries are zero-concentration flux boundary conditions.
Each time step in the simulation produces 500 concentration values at the
computational nodes based on the 500 DFNe realizations. Then, we collect the
concentration values and present the BTCs with means and standard deviations
(SDs) to provide the variation bandwidths.

The *x* and *y* coordinates for the vertical sampling lines of concentration
are (0.5, 1.0), (1.0, 1.0), and (1.5, 1.0). We allow the vertical lines to
have a fixed diameter to represent wells in practical problems, and the wells
are assumed to be open in the rock volume (Fig. 5). The BTCs along
the wells are calculated by using the flux weighted concentration (i.e.,
$v(x,t)c(x,t)$) at nodes involved in the wells. In this study, the well
diameter is fixed to 4 in (0.102 m) to calculate the mean BTCs and BTC uncertainties
along the wells.

We further define four profiles to assess the flow and transport
uncertainties. The profiles can be considered as a series of wells installed
along the profiles. The vertical profile (*y*=1.0 m) along the flow
direction (see Fig. 5) is used for assessing spatial variability in flow and
transport. This *y*=1.0 m profile has the width of 4 in. along the
profile. To obtain the head and transport uncertainties along the vertical
profile at each *x* location, we average the solutions of flow and transport
along the *z* direction. The sampling *x* step is the same as the profile
width. Solutions of MCS at each *x* location are collected based on the nodal
heads and concentration fluxes in the cuboid.

Three profiles perpendicular to the flow direction have the same *x*
locations as the sampling wells, i.e., *x*=0.5, 1.0, and 1.5 (m). In this
study, we also consider the profiles to have the same widths as the well
diameter (4 in.) for comparing vertical line (i.e., the well) and profile
(i.e., a series of wells) sampling strategies. The representative heads and
BTCs of the profiles are calculated by averaging the fracture nodal heads and
concentrations in the profile volumes.

This study focused on a relatively small fractured rock volume that was 2 m × 2 m × 2 m in size for the model verification and stochastic flow and transport modeling. The modeling domain might not be limited to this size when implementing this model. The following sections provide the results of the transport model verification and the uncertainties obtained from realizations of MCS.

## 4.1 Transport model verification

Figures 6 and 7 show a comparison of the concentration distributions for the horizontal fracture plate and cross-shaped fractures. The results in Figs. 6 and 7 show that identical solutions were obtained from both the developed DFN model and the HYDROGEOCHEM model. All the temporal and spatial variations in the plume were determined, and the solutions from the developed and HYDROGEOCHEM models were found to be identical. Figure 7 shows the concentration distributions after 3.0 days (Fig. 7a and b) and 5.0 days (Fig. 7c and d) when using the HYDROGEOCHEM and developed models for a cross-shaped fracture system. With a specified small upward flow applied to the vertical fracture, portions of the concentrations moved upward near the fracture line of intersection (Fig. 7b and d). This slightly upward flow relied on the constant head of 9.01 m that was applied on the top side of the vertical fracture. Again, the developed and HYDROGEOCHEM models were found to be identical for the cross-shaped fracture network.

Figure 8 shows the comparison of the solute concentration obtained from
analytical (dashed lines) and the developed (solid lines) models. Figure 8a–d show the concentration distribution at time 5.0, 7.5, 10.0, and
15.0 days, respectively. The concentration for the contour lines are for the
relative concentration (*c*∕*c*_{D}) of 0.1, 0.5, and 0.7. The
continuous release of the concentration along the central portion of the *x*=0
boundary leads to a high concentration area near the release location. In the
test example the uniform flow was assigned for the fracture plate. We set the
uniform seepage velocities in the *x* direction and *y* direction to be 0.1 and
0 m d^{−1} in all simulation domains. The uniform flow in the *x* direction
then leads to an isotropic dispersion for the simulation case. In summary, the
results of the developed model agree with those of the analytical model.

## 4.2 Flow upscaling behaviors

Figure 9 shows the results of the flow simulations when applied to the
500 DFN realizations. Figure 9a shows that the number of fractures increased
with the effective fracture intensity ${P}_{\mathrm{32}}^{e}$, where ${P}_{\mathrm{32}}^{e}$ is
the effective fracture intensity for total fracture intensity ${P}_{\mathrm{32}}^{t}$.
The effective fracture intensity represents the fracture intensity that is
evaluated without the isolated fractures in a rock volume. A higher number of
fractures could create greater variations in the 3-D fracture intensity.
However, such behavior was not valid for the 2-D effective fracture intensity
${P}_{\mathrm{21}}^{e}$ (see Fig. 9b). Similar numbers of fractures might vary in a
wide range of ${P}_{\mathrm{21}}^{e}$ values. This is because of the ${P}_{\mathrm{21}}^{e}$
calculation that relies on counting the trace length (i.e., fracture
intersections between fractures and the rock volume boundary) on the
downstream boundary profile (i.e., *x*=2.0 m) of the simulation domain. We
found that different sampling profiles at fixed *x* locations exhibited
similar patterns for the comparison of variables ${P}_{\mathrm{21}}^{e}$ and
${P}_{\mathrm{32}}^{e}$. Figure 9c summarizes the variations in the equivalent
hydraulic conductivity with different fracture intensities. The results
revealed that the high fracture intensity generally created a high equivalent
hydraulic conductivity. The results in Fig. 9c shows that the equivalent
hydraulic conductivity values were 2 to 3 orders of magnitude lower than the
specified fracture hydraulic conductivity value. Figure 9c shows that the
${P}_{\mathrm{32}}^{e}$ and equivalent hydraulic conductivity is highly correlated,
i.e., high values of ${P}_{\mathrm{32}}^{e}$ can yield high equivalent hydraulic
conductivity. Figure 9d further shows the comparison of ${P}_{\mathrm{32}}^{e}$ and
${P}_{\mathrm{32}}^{t}$. The results indicated that small numbers of fractures can
result in large variations of ${P}_{\mathrm{32}}^{t}$ when compared to a known
${P}_{\mathrm{32}}^{e}$. Such results implied that the fracture diameter used in the
DFN generations might be relatively small so that many isolated fractures
were removed for flow and transport simulations. The behavior could also lead
to high variations of flow and transport simulations.

## 4.3 Sampling volumes and flow and transport uncertainties

Figure 10 shows a realization of the transport simulation based on the DFNe flow field in Fig. 4. These results revealed that the continuous concentration released on the left boundary gradually migrated along the connected fractures. The spatial distributions of the concentration on the fracture plates were highly variable. This finding validated the concept that was proposed by Park et al. (2003), who stated that local flow cells contribute less to flow and contaminant transport in fracture formations.

Figure 11 shows the mean concentration BTCs (solid lines) and the associated
SD intervals (i.e., ±SD from the means) at sampling wells (Fig. 11a to
c) and along profiles (Fig. 11d to f) at the locations *x*=0.5, 1.0, and 1.5 m,
respectively. Comparisons of the means and ±SD BTCs at the wells and
along the profiles are shown in Fig. 11g–i, where the solid lines
represent the mean concentrations along the wells and the dashed lines are
the average concentrations along the profiles. We collected 500 BTC
realizations (gray lines) from the flow simulations and recursively estimated
the means and SD of the concentrations at different time steps. The results
in Fig. 11 clearly show that the maximum mean concentration BTCs for the
different averaging strategies might not reach the maximal concentration that
is released on the boundaries. The maximum concentrations of the mean BTCs
were 60 % to 80 % of the concentration at the maximum released
source. A longer simulation time might be required to obtain the maximum
concentration. Similar results were obtained with the same averaging strategy
at different monitoring locations. In general, small sampling volumes (i.e.,
the wells) obtain relatively large values of concentration means for
specified times. In addition, such small sampling volumes can also lead to
higher variations of BTCs (i.e., SD). In this study, the differences between
two sampling strategies were insignificant.

The acceptable number of the MCS realizations was decided based on the
comparison of statistical moments for different numbers of MC realizations.
With the realization numbers up to 500, the overall trends of head and
velocity variances were obvious, except for some variations along the
selected profiles. Figure 12 shows the distribution flow uncertainties along
the flow direction. Please note that the results were based on collections of
vertical averaged heads and Darcy velocities at each *x* location along *y*=1.0 m. The results in Fig. 12 showed the high nonstationarity of the head and
*x*-direction Darcy velocity. The distribution of the head variance exhibited
high variation at the inlet boundary, and the head variance gradually
decreased to a small value at the outlet boundary (Fig. 12a and b). This was
an interesting result, because the behavior was different from cases in
heterogeneous porous media, which showed similar head variations near inlet
and outlet boundaries (Li et al., 2004; Ni and Li, 2005, 2006). We believed
that the extremely high head variations near the DFN inlet boundary could be
induced by the generated DFN realizations. In the DFN flow simulation, the
constant head of the inlet boundary was applied to the fractures that
connect to the inlet boundary. The number of intersections at boundary
controlled the inlet flow for a generated DFN. Therefore, the uncertainty of
the head must involve the uncertainty of fractures connected to the inlet
boundary.

Figure 12c and d show the distribution of velocity variance in
*x* direction. In general, the highest value of the *x*-direction velocity
variance was 1 order of magnitude smaller than the highest value of the
head variance. The *x*-direction Darcy velocity also showed high variations
near inlet and outlet boundaries. The variances of *x*-direction Darcy
velocity at boundaries relies on calculations of head gradients and fracture
connections at boundaries. These two parameters are not deterministic for
MCS. These integrated uncertainties from head gradients and fracture
connections therefore lead to the increase of velocity uncertainties near
the boundaries. Away from the inlet and outlet boundaries, we found that the
velocity variance along the *x* direction showed a relatively stationary
zone. This behavior had been observed in groundwater flow in porous media (Ni
et al., 2010, 2011). The distance of the variance transition zone was close
to 0.5 m from the boundaries. Such a value was close to the mean fracture
diameter used for generating DFN realizations.

Figure 12e–h show the means and variances of the Darcy velocities in the
*y* and *z* directions. The results indicated that the velocity variations
in the *y* and *z* directions were relatively stationary. We found that the
boundary-induced nonstationarity was not obvious for velocity variations in
*y* and *z* directions. The stationary variances of the velocities in *y*
and *z* directions were similar to that of the *x*-direction velocity. In
Fig. 12 the fluctuations of head and velocity variances indicated that more
DFN realizations can improve the distributions of head and velocity
variances.

Figure 13 shows the distributions of transport uncertainty at different times
along the centerline profile (*y*=1.0 m). This study specified a constant
concentration on the inlet boundary. Similar to the situation of the constant
head condition in flow simulations, transport results also showed high
nonstationarity near the inlet boundary. However, the variances can propagate
to the downstream area with time. In the transport simulations, we specified
a constant concentration of 0 (kg m^{−3}) on the outlet boundary. The
uncertainty of concentration at the outlet boundary increased with time
because of the arrivals of concentration fronts. Similar to the condition on
the inlet boundary, the number of fracture intersections might influence the
increase of the concentration variation at the outlet boundary.

Numerical solutions to the ADE based on the Eulerian approach have not been widely implemented because of heavy computational costs. This study developed and implemented numerical models for stochastic modeling of flow and conservative transport in 3-D DFNs. The developed ADE model was validated by comparing the solutions of simple porous fractures with those from the HYDROGEOCHEM finite element model and an analytical model. When testing the transport model, identical temporal and spatial solutions were obtained from the developed model and the HYDROGEOCHEM model based on a Gaussian-type initial plume that was released in the porous fracture plate. For a simplified case, an analytical 2-D transport solution exists. The developed model accurately produced the results of concentration distributions in a horizontal fracture plate.

The MCS flow simulations showed that different fracture intensities can result in variations in the equivalent hydraulic conductivity that were 2 to 3 orders of magnitude lower than the fracture hydraulic conductivity values.

Simulations of transport in 3-D DFNs revealed that the maximum concentration of mean BTCs for different averaging strategies might not have reached the concentration. The sampling strategies along the wells and profiles yielded similar BTCs patterns. Based on the MCS, the means and SD for the two sampling strategies were observed at different sampling locations. MCS results showed that a smaller sampling volume can lead to relatively large values of mean concentrations and concentrations of SD for specified times.

MCS flow and transport showed that the distribution of the head variance
exhibited high variations at the inlet boundary, and the head variance
gradually decreased to a small value at the outlet boundary. The extremely
high head variation near the DFN inlet boundary could be induced by the generated
DFN realizations. No stationary zone for the head variance was obtained based
on collected MCS realizations. In the study, the value of the highest
*x*-direction velocity variance was 1 order of magnitude smaller than the
value of the highest head variance. The velocity variance along the
*x* direction showed a relatively stationary zone away from the inlet and
outlet boundaries. The distance of the variance transition zone was close to
0.5 m, which was the value of the mean fracture diameter for DFN generations.
Results of the velocity variations in the *y* and *z* directions appeared
relatively stationary along the flow direction. The boundary-induced
nonstationarity was not obvious for velocity variations in *y* and
*z* directions. The MCS of transport modeling showed high nonstationarity of
concentration variation near the inlet boundary, and the variance can
propagate downstream with time. The simulation results can provide
general insight into the evaluations of flow and transport uncertainty based
on the available DFN geometrical properties. One can employ the fracture
statistics to evaluate the possible flow paths or contaminant transport in
fractured rocks. The insufficient data from sites were represented by the
uncertainty and can be used for risk assessments or engineering designs.

The data that support the models in this study are available from the corresponding author upon request.

CFN and IHL conceived and designed the main idea of paper. IHL, CFN and FPL developed the model and analyzed the results. CPL and CCK contributed the data preparation and figures. IHL and CFN wrote the paper. All Authors read and approved the final draft.

The authors declare that they have no conflict of interest.

This research was partially supported by the Ministry of Science and
Technology of the Republic of China under grants MOST 103-2221-E-008-049-MY3, NSC
102-2116-M-008-010, and NSC 100-2625-M-008-005-MY3 as well as by the Institute of
Nuclear Energy Research, the Atomic Energy Council of the Executive Yuan, under grant
NL1030099, and by the Soil and Groundwater Pollution Remediation Fund in
2017.

Edited by: Philippe Ackerer

Reviewed by: two anonymous referees

Ahmed, R., Edwards, M. G., Lamine, S., Huisman, B. A. H., and Pal, M.: Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, J. Comput. Phys., 284, 462–489, 2015.

Bear, J. and Cheng, A. H. D.: Modeling groundwater flow and contaminant transport, Springer, New York, 2010.

Berrone, S., Fidelibus, C., Pieraccini, S., Scialò, S., and Vicini, F.: Unsteady advection-diffusion simulations in complex Discrete Fracture Networks with an optimization approach, J. Hydrol., 566, 332–345, 2018.

Bonnet, E., Bour, O., Odling, N. E., Davy, P., Main, I., Cowie, P., and Berkowitz, B.: Scaling of fracture systems in geological media, Rev. Geophys., 39, 347–383, 2001.

Botros, F. E., Hassan, A. E., Reeves, D. M., and Pohll, G.: On mapping fracture networks onto continuum, Water Resour. Res., 44, W08435, https://doi.org/10.1029/2007WR006092, 2008.

Bour, O., Davy, P., Darcel, C., and Odling, N.: A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (Hornelen Basin, Norway), J. Geophys. Res.-Sol. Ea., 107, ETG 4-1–ETG 4-12, https://doi.org/10.1029/2001JB000176, 2002.

Cacas, M. C., Ledoux, E., de Marsily, G., Barbreau, A., Calmels, P., Gaillard, B., and Margritta, R.: Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model, Water Resour. Res., 26, 491–500, 1990a.

Cacas, M. C., Ledoux, E., de Marsily, G., Tillie, B., Barbreau, A., Durand, E., Feuga, B., and Peaudecerf, P.: Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model, Water Resour. Res., 26, 479–489, 1990b.

Dagan, G.: Flow and transport in porous formations. Springer-Verlag GmbH & Co. KG, 1989.

de Dreuzy, J.-R., Pichot, G., Poirriez, B., and Erhel, J.: Synthetic benchmark for modeling flow in 3-D fractured media, Comput. Geosci., 50, 59–71, 2013.

Erhel, J., de Dreuzy, J. R., Beaudoin, A., Bresciani, E., and Tromeur-Dervout, D.: A parallel scientific software for heterogeneous hydrogeoloy, Parallel Computational Fluid Dynamics 2007, Springer Berlin Heidelberg, 39–48, 2009.

Fourno, A., Ngo, T.-D., Noetinger, B., and La Borderie, C.: FraC: A new conforming mesh method for discrete fracture networks, J. Comput. Phys., 376, 713–732, 2019.

Gelhar, L. W., Wetly, C., and Rehfeldt, K. R.: A critical review of data on field-scale dispersion in aquifers, Water Resour. Res., 28, 1955–1974, 1992.

Hartley, L. and Joyce, S.: Approaches and algorithms for groundwater flow modeling in support of site investigations and safety assessment of the Forsmark site, Sweden, J. Hydrol., 500, 200–216, 2013.

Hyman, J. D., Gable, C. W., Painter, S. L., and Makedonska, N.: Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: A feature rejection algorithm for meshing strategy, SIAM J. Sci. Comput., 36, A1871–A1894, 2014.

Hyman, J. D., Karra, S., Makedonska, N., Gable, C. W., Painter, S. L., and Viswanathan, H. S.: dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport, Comput. Geosci., 84, 10–19, 2015a.

Hyman, J. D., Painter, S. L., Viswanathan, H., Makedonska, N., and Karra, S.: Influence of injection mode on transport properties in kilometer-scale three-dimensional discrete fracture networks, Water Resour. Res., 51, 1–20, 2015b.

Johnson, J., Brown, S., and Stockman, H.: Fluid flow and mixing in rough-walled fracture intersections, J. Geophys. Res.-Sol. Ea., 111, B12206, https://doi.org/10.1029/2005JB004087, 2006.

Koike, K., Kubo, T., Liu, C., Masoud, A., Amano, K., Kurihara, A., Matsuoka, T., and Lanyon, B.: 3-D geostatistical modeling of fracture system in a granitic massif to characterize hydraulic properties and fracture distribution, Tectonophysics, 660, 1–16, 2015.

Kwicklis, E. M. and Healy, R. W.: Numerical investigation of steady liquid water flow in a variably saturated fracture network, Water Resour. Res., 29, 4091–4102, 1993.

Lee, I. H. and Ni, C.-F.: Fracture-based modeling of complex flow and
CO_{2} migration in three-dimensional fractured rocks, Comput.
Geosci., 81, 64–77, 2015.

Li, S. G., Liao, H. S., and Ni, C. F.: Stochastic modeling of complex nonstationary groundwater systems, Adv. Water Resour., 27, 1087–1104, 2004.

Liu, L. and Neretnieks, I.: Analysis of fluid flow and solute transport through a single fracture with variable apertures intersecting a canister: Comparison between fractal and Gaussian fractures, Phys. Chem. Earth, 31, 634–639, 2006.

Long, J. C. S., Gilmour, P., and Witherspoon, P. A.: A model for steady fluid flow in random three-dimensional networks of disc-shaped fractures, Water Resour. Res., 21, 1105–1115, 1985.

Makedonska, N., Painter, S., Bui, Q., Gable, C., and Karra, S.: Particle tracking approach for transport in three-dimensional discrete fracture networks, Computat. Geosci., 19, 1123–1137, 2015.

Ni, C. F. and Li, S. G.: Simple closed form formulas for predicting groundwater flow model uncertainty in complex, heterogeneous trending media, Water Resour. Res., 41, W11503, https://doi.org/10.1029/2005WR004143, 2005.

Ni, C. F. and Li, S. G.: Modeling groundwater velocity uncertainty in complex composite media, Adv. Water Resour., 29, 1866–1875, 2006.

Ni, C. F., Yeh, T. C. J., and Chen, J. S.: Cost-effective hydraulic tomography surveys for predicting flow and transport in heterogeneous aquifers, Environ. Sci. Technol., 43, 3720–3727, 2009.

Ni, C. F., Li, S. G., Liu, C. J., and Hsu, S. M.: Efficient conceptual framework to quantify flow uncertainty in large-scale, highly nonstationary groundwater systems, J. Hydrol., 384, 297–307, 2010.

Ni, C. F., Lin, C. P., Li, S. G., and Chen, J. S.: Quantifying flow and remediation zone uncertainties for partially opened wells in heterogeneous aquifers, Hydrol. Earth Syst. Sci., 15, 2291–2301, https://doi.org/10.5194/hess-15-2291-2011, 2011.

Odling, N. E.: Scaling and connectivity of joint systems in sandstones from western Norway, J. Struct. Geol., 19, 1257–1271, 1997.

Painter, S., Cvetkovic, V., Mancillas, J., and Pensado, O.: Time domain particle tracking methods for simulating transport with retention and first-order transformation, Water Resour. Res., 44, W01406, https://doi.org/10.1029/2007WR005944, 2008.

Park, Y.-J., Lee, K.-K., Kosakowski, G., and Berkowitz, B.: Transport behavior in three-dimensional fracture intersections, Water Resour. Res., 39, 1215, https://doi.org/10.1029/2002WR001801, 2003.

Pichot, G., Erhel, J., and de Dreuzy, J.: A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks, SIAM J. Sci. Comput., 34, B86–B105, 2012.

Pruess, K. and Tsang, Y. W.: On two-phase relative permeability and capillary pressure of rough-walled rock fractures, Water Resour. Res., 26, 1915–1926, 1990.

Stalgorova, E. and Babadagli, T.: Modified Random Walk–Particle Tracking
method to model early time behavior of EOR and sequestration of CO_{2}
in naturally fractured oil reservoirs, J. Petrol. Sci. Eng., 127, 65–81,
2015.

Stephens, M. B., Follin, S., Petersson, J., Isaksson, H., Juhlin, C., and Simeonov, A.: Review of the deterministic modelling of deformation zones and fracture domains at the site proposed for a spent nuclear fuel repository, Sweden, and consequences of structural anisotropy, Tectonophysics, 653, 68–94, 2015.

Trinchero, P., Painter, S., Ebrahimi, H., Koskinen, L., Molinero, J., and Selroos, J.-O.: Modelling radionuclide transport in fractured media with a dynamic update of Kd values, Comput. Geosci., 86, 55–63, 2016.

Wang, L. and Cardenas, M. B.: An efficient quasi-3-D particle tracking-based approach for transport through fractures with application to dynamic dispersion calculation, J. Contam. Hydrol., 179, 47–54, 2015.

Weng, X., Kresse, O., Chuprakov, D., Cohen, C.-E., Prioul, R., and Ganguly, U.: Applying complex fracture model and integrated workflow in unconventional reservoirs, J. Petrol. Sci. Eng., 124, 468–483, 2014.

Wexler, E. J.: Analytical solutions for one-, two-, and three-dimensional solute transport in ground-water systems with uniform flow, USGS Open-File Report, 89–56, 1992.

Xu, C. and Dowd, P.: A new computer code for discrete fracture network modelling, Comput. Geosci., 36, 292–301, 2010.

Yeh, G., Sun, J., Jardine, P., Burger, W., Fang, Y., Li, M., and Siegel, M.: HYDROGEOCHEM 4.0: A Coupled Model of Fluid Flow, Thermal Transport and HYDROGEOCHEMical Transport through Saturated-Unsaturated Media: Version 4.0. ORNL/TM-2004/103, Oak Ridge National Laboratory, Oak Ridge, 2004.

Yeh, G.-T.: On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow, Water Resour. Res., 17, 1529–1534, 1981.

Zafarani, A. and Detwiler, R. L.: An efficient time-domain approach for simulating Pe-dependent transport through fracture intersections, Adv. Water Resour., 53, 198–207, 2013.

Zhang, Q.-H.: Finite element generation of arbitrary 3-D fracture networks for flow analysis in complicated discrete fracture networks, J. Hydrol., 529, Part 3, 890–908, 2015.

Zheng, C. and Bennett, G. D.: Applied contaminant transport modeling, Wiley-Interscience New York, 2002.