Assessment of uncertainty in discrete fracture network modeling using probabilistic distribution method

Many engineering activities such as the geological disposal of nuclear waste and shale gas development are considered in the research of fractured media (Moreno et al. 2006). In the case of the geological disposal of high-level radioactive waste, fractured rocks are the main contamination transmission media, and the fracture characteristics determine the areas that become contaminated (Spycher et al. 2003; Schwartz 2012; Follin et al. 2013). The determination of the contaminated areas is important for further development of the site. Sensitive objects that require protection should be moved away as soon as a contaminated area is identified. Some field tracer tests can then be conducted in the contaminated areas as a basis for locating a radioactive waste disposal reservoir. In view of the foregoing researches, solute transport in geological disposal of radioactive waste has increasingly attracted the attention of researchers (Bredehoeft 1997; Browning et al. 2003; Ophori 2004; Kelkar et al. 2013).

The groundwater flow characteristics induced by complex rock structures in fractured media are described by discontinuity, heterogeneity, and anisotropy. Numerical simulation of particle tracking is an important research method used to solve complex solute transport problems in fractured media, and is increasingly being applied to the determination of contaminated areas (Frampton & Cvetkovic 2007, 2009; Reeves et al. 2008a, 2008b, 2008c; Rahmatian et al. 2015). Based on previous relevant studies, there are three main numerical models used for the simulation of fractured media, namely, the equivalent continuum model, dual-continuum model, and discrete fracture network (DFN) model (Bodin et al. 2007; Blessent et al. 2011). The equivalent continuum model is used to derive relevant equations based on the porous media flow theory (Reeves et al. 2008a, 2008b, 2008c), and it is employed in the study of nuclear waste disposal in the USA (Browning et al. 2003; Kelkar et al. 2013), Canada (Ophori 2004), and Sweden (Selroos et al. 2013). The dual-continuum model has been used to model radionuclide migration in nuclear waste disposal sites in Sweden (Schwartz 2012) and the Yucca Mountain (Spycher et al. 2003).

In the gradual development of the DFN model simulation, the effect of rock matrix has been ignored, and the fractures network is treated as the only flow area to describe the true status of the fluid flow through the fracture network (Berkowitz 2008; Follin et al. 2013). A DFN model considers the flow of groundwater in fractures, and does not consider the effect of matrix. It effectively describes the preferential channel of flow in the fractured medium. Thus, a DFN model adequately considers the heterogeneity and discontinuity characteristics of a fractured rock. When the matrix consists of dense rock, the permeability can be ignored and the model would produce highly precise simulation results (Blessent et al. 2011). To develop a DFN model, clear fracture geometry and spatial position information are required (Selroos & Painter 2012). The required detailed geometrical description of the fracture network geometry makes DFN models computationally demanding (Dverstorp et al. 1992). For geological disposal site-scale investigations, DFN models can be used to adequately represent the anisotropic and discontinuity features of fractured media (Frampton & Cvetkovic 2011).

Groundwater modeling commonly employs a single conceptual model (Refsgaard et al. 2006; Ye et al. 2010). Uncertainty in data, model structure, and model parameters can propagate throughout model runs, causing the model output to substantially deviate from the expected response of the model (Guzman et al. 2015). There are two main different uncertainties when the three main types of numerical models aforementioned are employed. The first type of uncertainty is the parametric uncertainty, which is due to uncertainty in each of the actual model, cannot be ignored. This uncertainty is commonly addressed by generating random values for the input parameters of the flow and contaminant transport models (Huang et al. 2011). Researchers have also used the Monte Carlo stochastic method and Bayesian approach to reduce parametric uncertainty (Feyen et al. 2003; Singh et al. 2010). Furthermore, a model-averaging method for assessing the uncertainly of a groundwater conceptual model has been developed (Ye et al. 2010) and applied (Woldt et al. 1992; Feyen & Caers 2006; Dong et al. 2013). MODPATH for particle tracking using the flow simulation from MODFLOW has likewise been used to apply the particle tracking method to the case of an equivalent continuum model (Reeves et al. 2008a, 2008b, 2008c), wherein 100 simulated particle trajectories were integrated and depicted in a figure (Langevin 2003). The goal of all the foregoing work was to reduce or assess the parametric uncertainty. The second type of uncertainty is engendered by the distinct difference between a DFN model and an equivalent continuum model. When a stochastic simulation is performed using equivalent continuum models, the flow areas of the models are the same because each model forms a continuum. Conversely, when DFN models are employed, the model structures are strongly dependent on the distribution of the fractures, which are generated by the stochastic method. Each DFN model thus has a different fracture network, and this leads to uncertainty of the model structures.

Few researchers have considered the uncertainty problem of model structures (Refsgaard et al. 2006). The traditional approach to model uncertainty analysis, which considers only a single conceptual model, may fail to adequately sample the relevant space of plausible conceptual models. Therefore, multiple conceptual models are often developed to represent the groundwater systems (Reeves et al. 2010, 2013; Pham & Tsai 2016), and multiple realizations of DFN models have been developed to assess the uncertainty of the model structures (Robertson 1990; Reeves et al. 2013). In addition, by means of the particle tracking method, two-dimensional (2D) DFN realizations have been used to determine the typical distributions of chemical solutes (Berkowitz & Scher 1998). To investigate the effects of different injection modes on the transport properties of a kilometer-scale site, three-dimensional (3D) DFN models with multi-realizations were applied to fractured granite at the Forsmark site in Sweden (Hyman et al. 2015). However, in the limited works on the uncertainty of model structures, the results of the particle tracking modeling are presented by the simply superimposed display method. This is because of the difficulty of assessing the uncertainty in DFN modeling and quantifying the potentially contaminated areas.

In view of the foregoing work, this study focuses on model structure uncertainty in DFN modeling using the probabilistic distribution method (PDM). The employed method enables illustration of the potentially contaminated area using probabilistic distribution maps, as well as quantification of the uncertainty of the occurrence of contamination. The particle tracking is modeled using the random walk method. The PDM was further applied to a potential high-level radioactive waste disposal reservoir site in Beishan, China, and the high-possibility contaminated areas were determined. This paper includes a description of the procedure that was used to build the employed DFN model.

As the study area, the Beishan area is located in the northwestern Gansu Province in China (Figure 1). The area has been selected as a potential site for the disposal of high-level radioactive waste (Guo et al. 2001; Ju et al. 2006; Dong et al. 2009; Wang 2010). The crust in the area is block structure, with the crust thickness of 47 km–50 km. No earthquakes with Ms > 4.75 were recorded in history. Since the Tertiary's period, it is a slowly uplifting area without obvious differential movement. Those geological characteristics show that the crust in the Beishan area is stable. In Beishan area, the dominant rocks are porphyritic monzonitic granite and tonalite. The deep intact rock of Beishan area is of high density, low porosity, high strength, low strain and high brittleness (Wang 2014).

The Beishan area is an arid Gobi desert area, with an average annual precipitation of 70 mm, with annual evaporation of 3,000 mm. More importantly, there is no perennial stream and other surface water bodies in the area. Atmospheric precipitation is the main source of the groundwater recharge in this area. Therefore, the Beishan area is poor in groundwater resources. Pumping tests carried out by local geological teams in the 1980s in the area have shown that, for most of the wells, the outflow rates are less than 50 m³/d (Ju et al. 2006). The present water table in the potential site area is 28–46 m below ground surface. Chemical analysis shows that the groundwater in Beishan area is of Cl-SO₄-Na type, with a pH value of 7–9, while the total dissolved solid is larger than 2 g/L (Wang 2014). Evaporation is the main discharge in this area. According to the regional hydrogeological data, the groundwater surrounding the leakage source flows from west to east.

Xinchang, as one of the eight granite sections selected as potential sites for disposal of high-level radioactive waste in Beishan area (Ju et al. 2006), is located in the middle-southern part of Beishan area. Granite is the most common rock type in the area (Sun et al. 2015), which has a low permeability. Surface geological mapping, surface geophysical survey and borehole investigation have shown the good integrity of the granitic rock mass in Xinchang section. In such low-permeability rocks, the groundwater flow rate and flow direction are mainly controlled by small-scale fractures, the locations, and the hydraulic properties of the rocks. BS18 is one of the deep boreholes drilled in the area (Zhao et al. 2013). The surrounding rock near BS18 was the specific study area. BS18 is 603 m deep and located in the Xinchang section, northeast in the Beishan area (Chen et al. 2007; Wang 2014). Considering the meter-scale of the study, a DFN model was useful for investigating the particle transport.

This section describes the methods used to generate and mesh the 3D DFN model of Beishan area. The modeling of the particle tracking and implementation of the PDM are also described. DFN model simulation using particle tracking generally involves the following three steps: (1) generating a realization of the network based on a fracture system model, (2) calculating the flow field for the imposed boundary conditions, and (3) moving particles by advection through the network (Parney & Smith 1995). The present study particularly involved a fourth step, namely, postprocessing of the particle tracking results by the PDM. The probability of the contaminant particle occurrence (P) is a binary variable (one for occurrence, zero for the opposite) in a DFN model. However, with the increase of the realizations, the P is calculated using the total times of particle occurrence by dividing the total number of realizations. Actually, the P values are related to the uncertainty in DFN modeling, and a high P value represents a high probability of the occurrence of contamination particle.

Based on the field data (borehole logging), fracture network properties are approximated by the best-fit theoretical statistical distributions (Starzec & Andersson 2002). To construct a model with complex fractures network, the Monte Carlo random generation method has been widely used (Prickett et al. 1981; Berkowitz & Scher 1998; Outters 2003). The Monte Carlo method was also employed in this study. For numerical problems involving high-dimension (Hastings 1970) and complex fracture network, a high-performance workstation is required for the application of the Monte Carlo method, which is included in the FracMan software (Dershowitz et al. 1993). The employed 3D DFN model contained a large number of two-dimensional intersecting fractures, which could be considered as plates with different apertures and sizes. Each fracture was partitioned into 2D triangular finite elements by Meshmaster (a FracMan subroutine), and the meshes could be refined and edited by Edmesh (another FracMan subroutine). A Galerkin finite element method was used to calculate the head distributions on the joint nodes. Twenty realizations of the DFN models were generated separately, and each model was used to conduct the particle tracking modeling, which was described below in detail.

The Lagrangian approach has been used in recent years to solve the solute transport problem (Park et al. 2001; Kim et al. 2004; Trinchero et al. 2016). It focuses on the changes in the number of particles at a mesh node and the solute particle trajectories, which properly reflect the heterogeneity of the fractured media. The contaminant movement in geological formations can be investigated by examining the various transitions of the particles, which represent, e.g. the dissolved solutes (Berkowitz et al. 2006). The particle tracking method has been used to determine the transport characteristics of groundwater flow in a fracture network (Parney & Smith 1995; Tsang et al. 1996; Painter 2002; Berkowitz 2008; Willmann et al. 2013). The random-walk theory is mainly presently employed for particle tracking in a DFN model (Prickett et al. 1981; Berkowitz et al. 2006; Bodin et al. 2007; Neuman & Tartakovsky 2009), as was done in this study. Each fracture was divided into several 2D triangular meshes. Each particle represented a fraction of the total solute mass. They were introduced into the model at specified solute sources and removed at encountered sink nodes in the cases of steady-state problems. In each time step, the particle movement is in accordance with the convective transport theory, but also includes random dispersive movement.

It should be noted that, if the particle remains within or exits the mesh, the new position of the particle would be calculated using the remaining or intercept time, i.e. ⁠, respectively. The relevant time is used to determine the boundary exit point coordinates. All the above procedures were implemented in this study in the MAFIC module of FracMan (Outters 2003). The module was developed for the simulation of transient flow and solute transport through 3D rock masses using DFN models, and has been applied to Oskarshamn in south-east Sweden (Starzec & Andersson 2002).

Here, we utilized the PDM to assess the uncertainty of the model structures, especially in randomly generated DFN models. The presently employed method enables the determination of the probability of contamination occurring in each zone of a DFN model. The particle tracking method is used to examine the solute transport process in a DFN model, and the PDM was utilized based on the particle position data for each time step. Considering the extreme case, the solute has a very strict concentration standard – one particle of the solute in place will make it a contamination area.

Considering there is bias when using a single realization of DFN to make an inference, we employed the multi-realization method to assess the uncertainty of the model structures. Multi-realization means that many DFN models are generated based on the same statistical data of the existing fractures. A mass of particles is placed at the same source point in these models. For each realization, the particle tracking is simulated using the same boundary and initial model conditions. The particle position data for each time step is outputted from MAFIC for the following process, which describes the MATLAB work flow of the PDM:

• (a)

  Define a rectangular region. The size of the region should be the same as the size of the DFN model. Divide the region into a large number of small rectangular grid cells, and size of the grid cell was chosen according to the requirement of the study and computer capacity. The number of stochastic realizations was set to be N.

• (b)

  Set a loop for reading the data for each realization. A variable (cnt) is used to count the frequency of the particle occurrence. For one-realization data, the position of each particle is projected to the corresponding zone grid (cnt = 1). If many particles appear in the same grid within the one-realization data, 1 should still be assigned to cnt. When the next realization data is being read and particles appear in the grids, execute cnt=cnt+1.

• (c)

  The probability of a particle occurring in a grid cell can be expressed as P=cnt/N. The probabilistic distribution figure enables the determination of the area most likely to be contaminated in the fractured media.

It should be noted that the output data comprises only the particle position data for a given time step. When the cell size is smaller than the distance between adjacent particles (Figure 3), the passed grids would be ignored if only the particle positions were considered. We then calculated the intersection of the particle trajectory and the grid line and determined the midpoint of the segment that connected the adjacent intersection. The passed grids are identified based on the location of the midpoints.

The employed method required that we first identified the distribution characteristics of the fractures to develop the numerical model. The modeling results for one realization of the DFN model contained uncertainty. The PDM was thus used to assess the uncertainty of the model structure and calculate the high-possibility contaminated areas.

Well imaging is a technology used to obtain internal images of a borehole and it enables researchers to determine the distribution of deep underground fractures and their occurrences (Dickens 1994; Nolen-Hoeksema et al. 2012). The method can detect all the components in the borehole through the intuitive observation. The whole borehole can be imaged and expanded into a floor plan and 3D column picture. Through high definition images, the apertures of fractures can be obtained. At last, the structure of the hole can be quantitatively and visually. The well imaging data of BS18 were used to develop the DFN model.

The fractures intersected by the well were logged by well imaging, and the fracture depths were also recorded. Figure 4 shows the variation of the cumulative fracture percentage with increasing depth. At the point where the cumulative fracture curve appears to twist, the number of fractures per unit length of the scanline (linear density: P10) begins to change. Based on the P10, the fractures are divided into five zones, which are distinguished by different colors in Figure 4. The well imaging data provide information about the fractures that are intersected by BS18. As can be observed from Figure 4, the actual data is in good agreement with the simulation results for one sample realization. This indicates that the developed DFN model properly represents the actual fractured media. It can also be seen from Figure 5(a) that the fracture intersections concentrate at several focus points, while Figure 5(b) shows that the fractures in Zone 5 are scattered. Because only eight fractures were logged in Zone 5, the statistical results for the zone may also be decentralized.

Beishan is a potential high-level radioactive waste disposal reservoir site. Such reservoirs are usually located at depths of 300–1,000 m (Wang 2010). In this study, it was assumed that the underground repository would be located at a depth of 550 m, which corresponds to the domain of Zone 5. Because of the stochastic distribution of the simulated fractures, the radioactive waste leakage source was assumed to be located at a depth of 540–560 m, which is also within Zone 5. The fractures in Zone 5 were thus considered for further investigation. The location of each zone and the fracture statistical data of each zone used for generating the model were presented in Table 1. The size of the employed model was 600 m (x) × 600 m (y) × 93 m (z).

The study area lies in an arid region and the location of the repository is quite deep. The surface and bottom of the model were thus assigned no-flux boundaries. According to the regional hydrogeological data, the groundwater surrounding the leakage source flows from west to east. The south and north boundaries were consequently also all considered no-flux boundaries, while the east and west boundaries were considered constant-head boundaries. The hydraulic gradient in the area is 0.493% and, based on field data, the heads of the west and east boundaries were respectively set to 1,681.5 and 1,678.5 m. The mass of each particle was 1E-6 kg, and 1,000 particles were released at t = 0 day instantly. The simulation time was 8,000,000 days.

Fractures in each DFN models are divided into plenty of triangular meshes. Meshes in a sample DFN model is shown in Figure 7. Particles transport pathways in the triangular meshes were calculated by the random walk method. The sample model in Figure 7 contains 183,728 triangular elements, which implies a high-performance workstation is necessary for this study.

The detailed parameters for generating DFN model, the meshing and solver parameters are shown in Table 2.

The particle tracking method and PDM were applied to the solute transport modeling of the Beishan area. Twenty realizations of the DFN model were performed separately, with 1,000 particles placed at the leakage source in each realization. The 3D particles' trajectories calculated in one sample model are shown in Figure 8. It can be found that the trajectories extend along the flow direction. To implement the PDM, a displayed side was first chosen and divided into a large number of rectangular grids, with each grid assigned a probability based on the number of particle occurrences. The employed grid cell size was 5 m × 5 m. Sample particle trajectories in the DFN model at 8E6 days are shown in Figure 9(a) and 9(c). As shown in this figure, the trajectories diverged. In addition, the fracture trends in this region were not concentrated, and this resulted in the solute being significantly transported in the direction of the regional groundwater, and diffused in the south-north direction. If the generated DFN model is used for simulation, the areas passed by the particles (red broken lines) would be determined to be certainly contaminated areas. A comparison of the results of the contaminated areas obtained by a single simulation with those obtained by the PDM method clearly reveals the inadequacy of the former method.

Twenty realizations of the DFN model were implemented for the investigation described here. The probabilistic distributions at 5E6 days for different numbers of realizations are shown in Figure 10, and the relative red zones in the figure indicate the very high probability of particle occurrences. We were generally more interested in the distance over which the pollutants might be transported in a region. It was therefore assumed that, when the P > 0.5, the domain of the grid had a high possibility to be contaminated. We also assumed that areas with P > 0.8 were certainly contaminated. With an increasing number of realizations N, the potentially contaminated areas and the high-probability areas gradually increase, although the potentially contaminated areas are gradually refined internally. The foregoing shows that employed PDM can be used to describe the status of the contaminant particle transport in DFN models.

The probability of particle occurrence was also quantified for each grid cell which was divided by the PDM. The obtained probabilistic distributions of 20 realizations are shown in Figure 11(a)–11(c), where the trajectories are colored according to the probabilities. The potentially contaminated areas (P > 0) at 5E6 days were observed to be larger than those calculated using one realization (Figure 9(a) and 9(c)). The certainly contaminated areas were quite small in Figure 11(a)–11(c), while the high-possibility contaminated areas were found to be discrete in Figure 11(a)–11(c). Different types of contaminated areas at 5E6 days were more identifiable than the areas passed by the particles for one realization because the uncertainty of the model structure was considered.

The PDM was also used to quantify the contamination probability changes over time, as shown in Figure 11. As can be seen from the figure, the potentially contaminated areas (blue areas) and high-possibility contaminated areas gradually increase with time. This is particularly obvious from the xz side view, with the high-possibility contaminated areas occupying most of the flow domain after 5E6 days. Particularly noteworthy is that the high-possibility contaminated areas near the release point are more concentrated than the farther areas after 5E6 days. The potentially contaminated areas and high-possibility contaminated areas are thus increasingly more dispersed with increasing distance from the release point.

The variations of the different types of contaminated areas over time from the xy and xz views are calculated. After 5E5 days of particle transport, the potentially contaminated areas covered 46,400 m² in the xy plan view, while the high-possibility contaminated areas covered 1,950 m². The total high-possibility contaminated area in the xz view was larger than that in the xy view at the monitoring time, but there were hardly any certainly contaminated areas owing to the additional discrete fractures of the model. However, certainly contaminated areas appeared in the xz view after 5E6 days, implying that a longer time allows for the possibility of definite contamination.

This paper utilized the PDM to assess the uncertainty of Beishan DFN models. In the application to solute transport, the particle tracking modeling was done using the random walk method and the PDM was implemented with multi-realization of fracture network. The employed method not only enables the determination of areas that are potentially contaminated by the solute using probabilistic distribution maps, but also affords quantification of the uncertainty of the contamination.

The PDM was applied to a potential high-level radioactive waste disposal reservoir site in Beishan, China. The procedure for developing the employed DFN model was described. The differences between the solute transport results for one realization and multiple realizations were analyzed, as well as the sensitivity to the number of realizations. Furthermore, the variations of contaminated areas given the different contamination levels over time were established. The results showed that the employed PDM can be used to describe the status of the contaminant particle transport in DFN models. The high-possibility contaminated areas near the contaminant release point were more concentrated than the farther areas after 5E6 days. In addition, the potentially contaminated and high-possibility contaminated areas were increasingly more dispersed with increasing distance from the release point. After 5E5 days of particle transport, the potentially contaminated areas covered a total of 46,400 m² in the xy plan view, while the high-possibility contaminated areas covered 1,950 m². A longer time was found to increase the certainly contaminated areas. The PDM combined with DFN modeling can provide a reference for the design of the potential geological disposal of high-level radioactive waste.

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB10030601). Partial support was also provided by the National Natural Science Foundation of China (Grant No. 41202182) and the Youth Innovation Promotion Association CAS (Grant No. 2016063).