Implementation of a Time-Domain Random-Walk Method into a Discrete Element Method to Simulate Nuclide Transport in Fractured Rock Masses

Wu, Yuexiu; Liu, Quansheng; Chan, Andrew H. C.; Liu, Hongyuan

doi:https://doi.org/10.1155/2017/5940380

Geofluids

On this page

Abstract Introduction Conclusion Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2017 | Article ID 5940380 | https://doi.org/10.1155/2017/5940380

Implementation of a Time-Domain Random-Walk Method into a Discrete Element Method to Simulate Nuclide Transport in Fractured Rock Masses

Yuexiu Wu,¹Quansheng Liu,¹Andrew H. C. Chan,²and Hongyuan Liu²

Academic Editor: Paolo Fulignati

Received12 May 2017

Accepted20 Jun 2017

Published15 Aug 2017

Abstract

It is essential to study nuclide transport with underground water in fractured rock masses in order to evaluate potential radionuclide leakage in nuclear waste disposal. A time-domain random-walk (TDRW) method was firstly implemented into a discrete element method (DEM), that is, UDEC, in this paper to address the pressing challenges of modelling the nuclide transport in fractured rock masses such as massive fractures and coupled hydromechanical effect. The implementation was then validated against analytical solutions for nuclide transport in a single fracture and a simple fracture network. After that, the proposed implementation was applied to model the nuclide transport in a complex fracture network investigated in the DECOVALEX 2011 project to analyze the effect of matrix diffusion and stress on the nuclide transport in the fractured rock masses. It was concluded that the implementation of the TDRW method into UDEC provided a valuable tool to study the nuclide transport in the fractured rock masses. Moreover, it was found that the total travel time of the nuclide particles in the fractured rock masses with the matrix diffusion and external stress modelled was much longer than that without the matrix diffusion and external stress modelled.

1. Introduction

1.1. Significance

With the development of nuclear waste disposal, nuclide transport with underground water in fractured rock masses has been an important issue in recent decades [1]. Especially after the nuclear leakage accident in Fukushima in Japan [2], this becomes more and more pressing. Therefore, it is essential to understand the principal transport mechanisms of radionuclides with underground water in fractured rock masses and evaluate the potential of radionuclide leakage from the repositories to the biosphere.

1.2. Literature Review about Nuclide Transport Simulation

Once the radionuclide is released from the container, it will flow with underground water in fractured rock masses to the biosphere. Therefore, water flow process in fractured rock masses should be studied firstly and foremost. Till now, five models have been found in literatures to simulate the water flow in fractured rock masses, that is, continuous model [3–5], dual-continuum model [6–8], equivalent continuum model [9–11], discrete fracture network (DFN) [12–14], and hybrid model [15–17]. In continuous model, rock masses are treated as continuum porous media with porosity and permeability and continuum-mechanics formulations can be used to describe flow and transport in media. It is a simple way but does not reflect the real flow path and velocity distribution in rock masses. In dual-continuum model, the rock matrix and fractures are represented as two overlapping and interactive continuum media with different mechanical and flow characteristics. However, this model cannot be applied to disconnected fractures in fractured media and the evaluation of fluid exchange between the matrix and the fractured domain is a complex task. In equivalent continuum model, fractured rock masses are treated as equivalent continua for large-scale analyses with equivalent mechanical and hydraulic properties that are measured or calculated in Representative Element Volume (REV). Nevertheless, REV may not exist for fractured rock masses in general. Moreover, flow rate in fractures is much larger than that in rock matrix. Thus, the DFN seems to be more suitable to simulate the water flow in fractured rock masses, irrespective of its high demand of computer resource. As a matter of fact, DFN has already been implemented into several numerical methods to simulate water flow in fractured rock masses, such as distinct deformation analysis (DDA) [18, 18–22], numerical manifold method (NMM) [23, 24], extended finite element method (XFEM) [25–27], and discrete element method (DEM) [28, 29]. Compared with other numerical methods, DEM coupled with DFN shows the advantages of simulating water flow through fractured rock masses, since it is capable of simulating coupled hydromechanical effect and, especially, massive fractures.

Secondly, radionuclide transport process associated with water flow through fractured rock masses should be analyzed. According to former researches [30, 31], the difficulty of discovering nuclide transport mechanism lies in solving the partial differential transport equation. As mentioned by Delay et al. [32], Lagrangian schemes (particle tracking methods) are more accurate and less computationally intensive than Eulerian schemes. The main Lagrangian schemes include random-walk particle tracking (RWPT) method [33, 34], convolution-based particle tracking (CBPT) method [35, 36], and continuous-time random-walk (CTRW) method [37, 38]. For very heterogeneous media such as fractured rock mass, flow velocities may span several orders of magnitude. With the RWPT method, a large number of jumps are required in low-velocity areas to move the particles significantly, which may cause the calculations to become inefficient in terms of computational costs. The CBPT method is valid for mainly steady-state flow field and linear transport process. Some discrepancies are caused when the CTRW method is employed to the advection-dispersion equation. Aiming at removing these constraints, Delay and Bodin [39] firstly proposed a new method called time-domain random-walk (TDRW) method to simulate solute transport in heterogeneous media, which was further extended to simulate solute transport in fractured rock masses [40–42].

As can be seen from the literature review above, DEM and the TDRW method have the advantages of modelling water flow and solute transport in fractured rock masses. Thus, the TDRW method was firstly implemented into the DEM software UDEC in this paper.

1.3. Objectives and Organization of Paper

The primary objective of this paper is to propose a new numerical tool to simulate nuclide transport with underground water in fractured rock mass. In order to achieve this objective, the paper is organized as follows.

Theoretical introduction about hydromechanical coupling in UDEC is given in Section 2; Section 3 discusses nuclide transport equation and theory about time-domain random-walk method; the implementation of TDRW into UDEC is given in Section 4; verification of implementation program with analytical solution is given in Section 5; Section 6 discusses application of implementation program to DECOVALEX (DEmonstration of COupled models and their VALidation against EXperiment) 2011 project to analyze the effect of matrix diffusion and stress on nuclide transport process; and conclusion is given in Section 7.

Throughout this study, it was concluded that, due to the effects of both matrix diffusion and stress, it took much more time for particles to outflow the fractured rock masses. Moreover, the proposed numerical method is a valuable numerical tool to study nuclide transport process through fractured rock masses, which has also wide applications in the field of underground water pollution, oil reservoir engineering, and so on.

2. Principle of Flow in UDEC

In UDEC [43], fractured rock masses are divided into many deformable blocks whose boundaries are fractures. Each deformable block is meshed into triangle elements, as shown in Figure 1. Motion is first calculated at the grid points of triangle elements within blocks and the stresses within the elements are then obtained according to the block material constitutive relations. Moreover, a fully coupled mechanical-hydraulic analysis can also be performed in UDEC, in which matrix is regarded as impermeable and fluid flows only through fractures.

2.1. Mechanical Behavior

For any grid point in a triangle element, the total force is obtained as a sum of four terms:where are the external applied loads, are the internal stresses contributed by zones adjacent to grid point , are the contact forces from all the adjacent blocks, and are gravity forces. The internal stress is calculated using the following equation:where is element stress tensor and is the unit vector outward normal to the integral path , which follows the closed polygonal line around the grid point . The gravity force is calculated using the following equation:where is the lumped gravitational mass, defined as the sum of one-third of the masses of triangles connected to the grid point, and is the gravitational acceleration.

The contact force exists only when grid points are on block boundaries (like grid point in Figure 1) and is calculated according to the contact type and joint constitutive model. Taking the Coulomb slide model as an example, the contact force at the grid point from an adjacent block is divided into two components: (normal contact force) and (shear contact force):where is the normal displacement increment for each contact, is joint normal stiffness, and is the contact area, andwhere is the shear displacement increment and is joint shear stiffness. If the shear contact force exceeds the shear strength, it is modified according to the following equation:where is joint cohesion and is joint frictional strength.

The motion of the grid point is obtained according to Newton’s second law:where is the velocity of grid point and is time.

After the calculation of displacement of each grid point, the strain and stress of each triangle element are obtained according to the block constitutive model:where is time step.where and are Lame’s constant and is the bulk strain.

2.2. Fluid Flow in Fractures

Depending on the different contact type, the flow rate in contact domain is calculated differently. For a point contact (from a domain with pressure to a domain with pressure ),where is the flow rate, is a point contact permeability factor, and is the pressure difference calculated using the following equation:where is the fluid density; is the acceleration of gravity; and and are -coordinates of the domain centers which are perpendicular to the joint surface. For edge-edge contact (an edge of one element contacts with an edge of another element), the flow rate in contact domain is calculated using the following equation:where is a joint permeability factor (whose theoretical value is one-twelfth dynamic viscosity of fluid), is the contact hydraulic aperture, and is the length assigned to the contact between the domains. Fluid velocity can be obtained using (13) according to the flow rate calculated according to (10) and (12):where is the fluid crossing area.

2.3. Coupled Mechanical-Hydraulic Analysis

When the fluid flows in the contact domain, the fluid pressure is regarded as an external force applied on the contact domain, which is added into . The external force is calculated using the following equation:where is the fluid pressure in a contact domain, is the normal direct of contact, and is the contact length. With the motion of each block, the contact area will change following (15), which causes flow rate variation: where is the initial fracture aperture and is the normal displacement of contact area.

3. Nuclide Transport

3.1. Nuclide Transport Equation

It is known that nuclide transport is the consequence of several physical mechanisms such as advection, hydrodynamic dispersion, matrix diffusion, decay, and adsorption. As mentioned by Bodin et al. [42], the nuclide transport can be described in three media: () fracture, () matrix, and () stagnant zone. The stagnant zones in the fracture plane act as an additional “nonflowing” pore space available for solute diffusion. Thus, transport in the stagnant zone is ignored in this study. The nuclide transport equation in 2D fracture can be written as follows [42]:where is the concentration of nuclide in fracture, is the fluid velocity in fracture, is half of the aperture, is the concentration of nuclide in rock matrix, is the decay constant of nuclide, the term represents the decay effect, is the hydrodynamic dispersion coefficient in fracture, is a retardation factor accounting for the sorption of nuclide onto fracture surface, the last term on the right of the equation represents the matrix diffusion effect, and and are the space coordinates along and perpendicular to the fracture, respectively, as shown in Figure 2. Not considering the effect of decay and matrix diffusion, the equation becomes the advection-dispersion-sorption (ADS) equation.

The retardation factor is defined using the following equation [42]:where is the surface sorption coefficient of nuclide onto the fracture surface. The nuclide transport equation in rock matrix can be written as follows [42]:where is the apparent-diffusion coefficient of the nuclide in the matrix, which is calculated as ; is effective diffusion coefficient; is matrix porosity; is matrix density; and is volumetric-sorption coefficient of the nuclide in the matrix.

3.2. Time-Domain Random-Walk Method

3.2.1. Advection-Dispersion-Sorption Equation

In TDRW method, the solution of (16) can be divided into three steps: () advection-dispersion-sorption (ADS) transport equation, () decay effect, and () matrix diffusion effect. The ADS equation can be rewritten as (19) [42], considering the change of variable ():By applying an equivalent change of variable, (19) can be written with Fokker-Planck formalism:According to the solution of random-walk approach, the mean and variance of nuclide transport time distribution for a displacement of are obtained using (21a) and (21b), respectively [42]:For Peclet number larger than 10, it is shown that the travel time distribution is lognormal for a displacement of [30]:where is the travel time of the nuclide for a displacement and is a random number that follows a standard normal distribution.

For very low flow velocities in short fractures (i.e., the Peclet number is less than 10), the assumption of a lognormal distribution is not accurate. Bodin et al. [30] suggested that (22b) should be modified:The total travel time in a fracture of length can be written asWhen Peclet number is less than 10, the variable is replaced by .

3.2.2. Matrix Diffusion

To account for the interaction between ADS and matrix diffusion in a bond of length , the fluid velocity and dispersion coefficient in (22a) and (22b) were replaced by an apparent fluid velocity and an apparent dispersion coefficient [42]:where is a dispersivity constant and is a retardation factor accounting for matrix diffusion:where , , is the matrix porosity, is the bulk density of matrix, is the volumetric-sorption coefficient of nuclide in matrix, and is the effective diffusion coefficient of nuclide in matrix.

The modified (22a) and (22b) with and in place of and , respectively, were used to calculate the mean and variance of the log transform (or ) and according to (22b) and (22c), respectively. The total travel time for the coupled ADS and matrix diffusion is calculated using (27) for [42]:where is the particle residence time in the fracture for pure advection delayed by sorption reactions onto the fracture walls, expressed as and is a random number drawn from a uniform distribution between 0 and 1. For , the variable in (27) is replaced by .

3.2.3. Decay Effect

If the linear decay model is used, the nuclide mass decreases according to the following equation [42]:where and represent the masses of the nuclide at node and , respectively.

For simplicity, the decay effect is not considered here, which accords with the situation in DECOVALEX 2011 project.

3.2.4. Mass Partition Model at Fracture Intersection

Because the nuclide was tracked in fractures from the inflow boundary to the outflow boundary in TDRW method, the mass partition model at fracture intersection is vital in the calculation. There are three models: the perfect-mixing model [44], the stream-tube model [45], and the diffusional-mixing model [46]. The focus of this paper is not the mass partition law of solute at fracture intersection. For simplicity, the perfect-mixing model was adopted here:where is the probability of particle transition from an inlet bond to an outlet bond ; is the flow rate in the outlet bond ; and is the sum of the flow rates over all the outlet bonds connected to the junction of interest. It can be seen from (29) that mass partition is only related to the outlet fluxes.

4. Implementation of TDRW into UDEC

In TDRW method, the fractures are regarded as the sum of many bonds and a continuous connection of the bonds from inlet to outlet boundaries is considered as a flow path (e.g., the dashed line shown in Figure 3). In a bond, the fluid flow is a 1D flow with a constant fluid velocity. Therefore, the first thing, before we do the particle tracking calculation, is to get the flow paths of particle transport.

In UDEC, jointed rock masses are divided into blocks by fractures. The geometric information of fracture and block, such as corner coordinates and their connection, is stored in list in files (e.g., BLOCK.FIN, DOMAIN.FIN, and CONTACT.FIN). According to the above files, the coordinates of fracture intersections can be obtained easily, as well as their connective information. The TDRW method is implemented into UDEC using the FISH language, as shown in Figure 4.

After the hydromechanical calculation with UDEC, the final flow velocity at each fracture intersection can be obtained and this is the input data for TDRW simulation. At the beginning, nuclide particles are injected into the fractured rock masses through the inlet boundary. According to the inlet flow rate, the injection location of each particle is chosen. The particle transport time in each bond is calculated with ((22a), (22b), (22c))–(28). At the fracture intersection, the perfect-mixing model is adopted to decide the next step of particles until they arrive at the outlet boundary. If the injected mass is and the total number of particles is , mass of each particle is expressed as

If the collected particle number on the outlet boundary is in a period of time , the outlet concentration is written aswhere is the total outlet fluid flow rate.

The following steps are adopted to obtain the curve of the relative concentration versus time:(1)Calculate the maximum and minimum travel times of all particles.(2)Divide the travel time period (from minimum travel time to maximum travel time) into parts, in which the time interval is .(3)Collect particles (number is ) in each time interval and calculate their concentration with (31).

5. Verification

Two numerical particle transport tests were performed to verify the implementation of the TDRW method into UDEC, that is, particle transports in a single fracture and in a simple facture network.

5.1. Analytical Solution of the Nuclide Transport Equation

Analytical solution for nuclide transport in a single fracture has been deduced in literatures [47–49]. The analytical solution is to be used later to verify the numerical results obtained with combination of DEM and TDRW. For simplicity, the analytical solution for advection-dispersion transport equation with constant concentration injection is adopted, which can be written as follows:where is the constant concentration injection; and is a complementary error function.

5.2. Calibration Test 1: Particle Transport in a Single Fracture

In the calibration test 1 of particle transport in a single fracture, the injection mass at the inlet boundary of the fracture was 0.001 g, and the total number of nuclide particles was 30000. The trace of the fracture was 10 m, fluid velocity was 4.0 × 10⁻⁴ m/s, hydrodynamic dispersion coefficient in the fracture was 2.0 × 10⁻⁵ m²/s, and the aperture of the fracture was 2.5 × 10⁻⁴ m. Thus, the Peclet number () was 200. The analytical solution can be obtained using (32). The corresponding numerical solution was obtained using the TDRW method implemented into UDEC, which is compared in Figure 5 with the analytical solution. From Figure 5, the maximum and minimum differences between numerical and analytical solutions are 11% and 0.03%, respectively. The reason for this difference is that the numerical results are iterated. In general, the results from UDEC with the TDRW method implemented agreed well with the analytical solution.

5.3. Calibration Test 2: Particle Transport in a Simple Fracture Network

The simple fracture network was shown in Figure 6. A hydromechanical calculation was performed with UDEC at first. The water heads at top and bottom boundaries were 100 m and 99 m, respectively, and two lateral sides were impermeable. The aperture of the fractures was 2.5 × 10⁻⁶m. The other parameters were the same as those in test 1. The analytical solution was the summarizing solution of (32) for each fracture. Figure 7 was the comparison between the numerical solution and analytical solution for the particle transport in this simple fracture network. From Figure 7, the maximum and minimum differences between numerical and analytical solutions are 10.1% and 0, respectively.

From the comparisons shown in Figures 5 and 7, it is concluded that the numerical solution obtained using the TDRW method implemented into UDEC was close to the analytical solution. Thus, it is proven that the TDRW method implemented into UDEC can well simulate nuclide particle transport in fracture networks.

6. Application to Particle Transport in Discrete Fracture Network

6.1. Calculation Models and Material

The proposed method was used to evaluate the effect of matrix diffusion and stress ratio on particle transport process in the fracture network studied in the DECOVALEX 2011 project. The calculation model was shown in Figure 8, which contained 771 fractures. The hydromechanical boundary condition was depicted in Figure 9. The hydraulic heads at the top and bottom boundaries were 30 m and 10 m, respectively. On the two lateral sides, a linear distribution of hydraulic head was applied with a maximum value of 30 m and a minimum value of 10 m; in order to analyze the effect of matrix diffusion, the same hydromechanical boundaries were applied for the case with matrix diffusion and the case without matrix diffusion. Accounting for the effect of stress ratio, simulations without stress () and with different stress ratios were performed. For the cases with stress, a constant vertical stress of 5 MPa was specified on the top and bottom boundaries. And a horizontal stress, with ratios of vertical/horizontal stress varying from 1 to 2, 3, and 5, was applied, respectively (). For each calculation, 10000 particles were injected at a time along the inflow boundary. At each injected point, which is the intersection of fracture and inflow boundary, the injecting particle number is proportional to the local flow rate. This is equivalent to injecting a constant concentration at the inflow boundary.

The block material property was GPa, and GPa.

The normal behavior of the fractures in the calculation model was nonlinear, as shown in Figure 10 [28]. Other parameters of the fractures were listed in Table 1 [28]. The aperture of the fractures followed a lognormal distribution with a mean value of 65 and a second moment of 1.0. The minimum and maximum apertures were 1 and 200 , respectively. The residual aperture was assumed to be 20% of the initial aperture.

The hydrodynamic dispersion coefficient in fractures was 1.0 × 10⁻⁵ m²/s. The hydromechanical calculation was done firstly and the fluid velocity could then be obtained. After that, the TDRW simulation was performed.

6.2. Results and Discussion

6.2.1. Fluid Flow

Flow rates in fractures which intersect the bottom, right, and left outlet boundaries with different applied stress were shown in Figures 11–13. The location in Figures 11–13 was the coordinates of intersections of fractures along the corresponding boundary. On the bottom outlet boundary, the maximum flow rate and average flow rate without stress applied were 1.16 × 10⁻⁸ m³/(s·m) and 4.06 × 10⁻⁹ m³/(s·m), respectively, which were larger than those with stress applied. From Figure 11, it was shown that the main outlet positions (intersections with flow rate larger than the average flowrate) changed comparing the results not considering stress with results considering stress effect. When the stress was applied, the main outlet positions were same for different stress ratio. However, the maximum and average values of flow rate decreased with the increasing of stress ratio. On the right boundary, the maximum and average values of flow rate in intersections without stress were 8.59 × 10⁻⁸ m³/(s·m) and 3.13 × 10⁻⁸ m³/(s·m), respectively, which were obviously larger than those considering stress ratio. On the left boundary, the maximum and average values of flow rate in intersections without stress were 5.38 × 10⁻⁷ m³/(s·m) and 8.02 × 10⁻⁸ m³/(s·m), respectively, and the same effect of stress ratio on fluid flow was observed. This was because the fractures were compressed to be narrow with the increasing of external stress and the fluid flow became smaller.

6.2.2. Nuclide Transport

(1) Without Matrix Diffusion. After the nuclide transport calculations with different stress ratio, nuclides which outflow the outlet boundary at a certain time were collected. The ratio of number of leaking nuclides from the outlet boundary to the total number of nuclides was plotted versus traveling time (i.e., breakthrough curve) in Figure 14. From Figure 14, it was shown that the total traveling times with different stress boundaries were 1.97 × 10⁴ s, 2.58 × 10⁴ s, 3.12 × 10⁴ s, 1.95 × 10⁵ s, and 2.4 × 10⁵ s, respectively. The mean traveling time for each calculation was 3.29 × 10³ s, 4.31 × 10³ s, 5.21 × 10³ s, 1.27 × 10⁴ s, and 2.19 × 10⁴ s, respectively. Through dividing the total traveling times without stress, the normalized total travel times with stress ratio 1, 2, 3, and 5 are 1.3, 1.58, 9.89, and 12.18.

The total and average travel times for all nuclides leaking out from the fractured rock masses increased with the increasing of stress ratio. The reason for this phenomenon was that fractures were compressed to be narrow and it was difficult for nuclides to leak.

(2) With Matrix Diffusion. In the case with matrix diffusion, the porosity of the matrix was 0.316%. For simplicity, advection, dispersion, and matrix diffusion were considered in this case, while sorption in the matrix was not allowed ( m³/kg). The effective diffusion coefficient of the nuclide in the matrix is 1.0 × 10⁻¹¹ m²/s. The breakthrough curves with different stress applied considering the effect of matrix diffusion were shown in Figure 15. From Figure 15, it was shown that the total traveling times with different stress boundaries were 3.03 × 10⁷ s, 7.98 × 10⁷ s, 1.086 × 10⁸ s, 8.02 × 10⁸ s, and 8.42 × 10⁸ s, respectively. The mean traveling time for each calculation was 4.88 × 10⁶ s, 1.77 × 10⁷ s, 2.89 × 10⁷ s, 5.2 × 10⁷ s, and 7.7 × 10⁷ s, respectively. Through dividing the total traveling times without stress, the normalized total travel times with stress ratio 1, 2, 3, and 5 are 2.63, 3.58, 26.5, and 27.8. The total and average travel times for all nuclides leaking out from the fractured rock masses increased with the increasing of stress ratio.

Comparing the breakthrough curves in Figure 14 with those in Figure 15, it could be seen that the total travel time for all particles outflowing from fractured rock masses without the matrix diffusion was much shorter than that with matrix diffusion. Thus, the effect of matrix diffusion on the nuclide particle transport was significant.

7. Conclusion

Although the nuclide transport with underground water in fractured rock masses had been studied for decades by many researchers, the coupled hydromechanical effect and the effect of matrix diffusion and stress were still not discussed systematically. Aiming at solving these problems, the time-domain random-walk (TDRW) method was firstly implemented into the discrete element method (DEM) in this study. During each time step of hydromechanical calculation, hydraulic aperture of each fracture and fluid pressure on fracture surface were calculated according to (15) and (14) to consider the coupled hydromechanical effect. To present the effect of matrix diffusion on nuclide transport, (24) and (27) were adopted to get the breakthrough curves. The implementation was then verified against the analytical solutions for nuclide transport in a single fracture and a single fracture network. After that, the implementation was applied to simulate the nuclide transport in a complex fracture network investigated in the DECOVALEX 2011 project. Finally, the effect of matrix diffusion and stress on nuclide transport was discussed in detail. Throughout this study, the following conclusions were drawn:(1)The implementation of the TDRW method into UDEC could well simulate the nuclide transport with underground water in fractured rock masses.(2)Through the hydromechanical calculation for DECOVALEX 2011, the closure of fractures increased with the increasing of applied stress and it caused decreasing of fluid flow in fractured rock masses. More time was needed for fluid or particle transport through a narrower path. This finally resulted in the fact that the total travel time of the nuclide particles with the external stress applied on the fracture network was longer than that without the applied stress.(3)According to the numerical results for the nuclide transport in the complex fracture network investigated in the DECOVALEX 2011 project, matrix diffusion has a significant effect on the nuclide transport. The travel time of the nuclide particle with the matrix diffusion modelled was much longer than that without the matrix diffusion.

Conflicts of Interest

The authors declare that the received funds do not lead to any conflicts of interest regarding the publication of this manuscript.

Acknowledgments

This work was supported by the Natural Science Foundation of China (nos. 51679172, 41130742, and 51574180), National Basic Research Program of China (973 Program) (no. 2014CB046904), and China Scholarship Council.

References

C.-F. Tsang, I. Neretnieks, and Y. Tsang, “Hydrologic issues associated with nuclear waste repositories,” Water Resources Research, vol. 51, no. 9, pp. 6923–6972, 2015.
View at: Publisher Site | Google Scholar
H. Funabashi, “Why the fukushima nuclear disaster is a man-made calamity,” International Journal of Japanese Sociology, vol. 21, no. 1, pp. 65–75, 2012.
View at: Publisher Site | Google Scholar
F. Cappa and J. Rutqvist, “Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO₂,” International Journal of Greenhouse Gas Control, vol. 5, no. 2, pp. 336–346, 2011.
View at: Publisher Site | Google Scholar
J. Rutqvist and O. Stephansson, “The role of hydromechanical coupling in fractured rock engineering,” Hydrogeology Journal, vol. 11, no. 1, pp. 7–40, 2003.
View at: Publisher Site | Google Scholar
G. Lei, P.-C. Dong, Z.-S. Wu et al., “A fractal model for the stress-dependent permeability and relative permeability in tight sandstones,” Journal of Canadian Petroleum Technology, vol. 54, no. 1, pp. 36–48, 2015.
View at: Publisher Site | Google Scholar
Y. Hao, P. Fu, and C. R. Carrigan, “Application of a dual-continuum model for simulation of fluid flow and heat transfer in fractured geothermal reservoirs,” in Proceedings of the 38th Workshop on Geothermal Reservoir Engineering, vol. SGP-TR-198, pp. 462–469, Stanford University, Stanford, Calif, USA, 2013.
View at: Google Scholar
M. Ahmadpour, M. Siavashi, and M. H. Doranehgard, “Numerical simulation of two-phase flow in fractured porous media using streamline simulation and IMPES methods and comparing results with a commercial software,” Journal of Central South University, vol. 23, no. 10, pp. 2630–2637, 2016.
View at: Publisher Site | Google Scholar
J. Douglas Jr., J. L. Hensley, and T. Arbogast, “A dual-porosity model for waterflooding in naturally fractured reservoirs,” Computer Methods in Applied Mechanics and Engineering, vol. 87, no. 2-3, pp. 157–174, 1991.
View at: Publisher Site | Google Scholar
A. R. Khoei, N. Hosseini, and T. Mohammadnejad, “Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model,” Advances in Water Resources, 2015.
View at: Publisher Site | Google Scholar
S. H. Lee, M. F. Lough, and C. L. Jensen, “Hierarchical modeling of flow in naturally fractured formations with multiple length scales,” Water Resources Research, vol. 37, no. 3, pp. 443–455, 2001.
View at: Publisher Site | Google Scholar
M. Oda, “An equivalent model for coupled stress and fluid flow analysis in jointed rock masse,” Water Resources Research, vol. 22, no. 13, pp. 1845–1856, 1986.
View at: Publisher Site | Google Scholar
J. D. Hyman, S. Karra, N. Makedonska, C. W. Gable, S. L. Painter, and H. S. Viswanathan, “DfnWorks: a discrete fracture network framework for modeling subsurface flow and transport,” Computers and Geosciences, vol. 84, pp. 10–19, 2015.
View at: Publisher Site | Google Scholar
J. Maillot, P. Davy, R. Le Goc, C. Darcel, and J. R. de Dreuzy, “Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models,” Water Resources Research, vol. 52, no. 11, pp. 8526–8545, 2016.
View at: Publisher Site | Google Scholar
S. Follin, L. Hartley, I. Rhén et al., “A methodology to constrain the parameters of a hydrogeological discrete fracture network model for sparsely fractured crystalline rock, exemplified by data from the proposed high-level nuclear waste repository site at Forsmark, Sweden,” Hydrogeology Journal, vol. 22, no. 2, pp. 313–331, 2014.
View at: Publisher Site | Google Scholar
M. Hu, J. Rutqvist, and Y. Wang, “A practical model for fluid flow in discrete-fracture porous media by using the numerical manifold method,” Advances in Water Resources, vol. 97, pp. 38–51, 2016.
View at: Publisher Site | Google Scholar
U. Svensson, M. Ferry, and H.-O. Kuylenstierna, “DarcyTools Version 3.4—Concepts, methods and equations,” Tech. Rep. SKB R-07-38, SvenskK€arnbr€anslehantering AB, Stockholm, Sweden, 2010.
View at: Google Scholar
S. P. Neuman, “Trends, prospects and challenges in quantifying flow and transport through fractured rocks,” Hydrogeology Journal, vol. 13, no. 1, pp. 124–147, 2005.
View at: Publisher Site | Google Scholar
Y.-Y. Jiao, H.-Q. Zhang, H.-M. Tang, X.-L. Zhang, A. C. Adoko, and H.-N. Tian, “Simulating the process of reservoir-impoundment-induced landslide using the extended DDA method,” Engineering Geology, vol. 182, pp. 37–48, 2014.
View at: Publisher Site | Google Scholar
L. Q. Choo, Z. Zhao, H. Chen, and Q. Tian, “Hydraulic fracturing modeling using the discontinuous deformation analysis (DDA) method,” Computers and Geotechnics, vol. 76, pp. 12–22, 2016.
View at: Publisher Site | Google Scholar
W. E. Morgan and M. M. Aral, “An implicitly coupled hydro-geomechanical model for hydraulic fracture simulation with the discontinuous deformation analysis,” International Journal of Rock Mechanics and Mining Sciences, vol. 73, pp. 82–94, 2015.
View at: Publisher Site | Google Scholar
F. Zheng, Y. Jiao, M. Gardner, and N. Sitar, “A fast direct search algorithm for contact detection of convex polygonal or polyhedral particles,” Computers and Geotechnics, vol. 87, pp. 76–85, 2017.
View at: Publisher Site | Google Scholar
F. Zheng, Y.-Y. Jiao, X.-L. Zhang, F. Tan, L. Wang, and Q. Zhao, “Object-oriented contact detection approach for three-dimensional discontinuous deformation analysis based on entrance block theory,” International Journal of Geomechanics, vol. 17, no. 5, 2017.
View at: Publisher Site | Google Scholar
H. Zheng, F. Liu, and C. Li, “Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method,” Applied Mathematical Modelling, vol. 39, no. 2, pp. 794–808, 2015.
View at: Publisher Site | Google Scholar | MathSciNet
Y. Wang, M. Hu, Q. Zhou, and J. Rutqvist, “Energy-work-based numerical manifold seepage analysis with an efficient scheme to locate the phreatic surface,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 15, pp. 1633–1650, 2014.
View at: Publisher Site | Google Scholar
S. Salimzadeh and N. Khalili, “Fully coupled XFEM model for flow and deformation in fractured porous media with explicit fracture flow,” International Journal of Geomechanics, vol. 16, no. 4, Article ID 04015091, 2016.
View at: Publisher Site | Google Scholar
T. Mohammadnejad and J. E. Andrade, “Numerical modeling of hydraulic fracture propagation, closure and reopening using XFEM with application to in-situ stress estimation,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 40, no. 15, pp. 2033–2060, 2016.
View at: Publisher Site | Google Scholar
T. Mohammadnejad and A. R. Khoei, “Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 37, no. 10, pp. 1247–1279, 2013.
View at: Publisher Site | Google Scholar
Z. Zhao, L. Jing, I. Neretnieks, and L. Moreno, “Numerical modeling of stress effects on solute transport in fractured rocks,” Computers and Geotechnics, vol. 38, no. 2, pp. 113–126, 2011.
View at: Publisher Site | Google Scholar
L. Quansheng, W. Yuexiu, and L. Bin, “Discrete element analysis of effect of stress on equivalent permeability of fractured rockmass,” Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 1, pp. 176–183, 2011 (Chinese).
View at: Google Scholar
J. Bodin, F. Delay, and G. de Marsily, “Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms,” Hydrogeology Journal, vol. 11, no. 4, pp. 418–433, 2003.
View at: Publisher Site | Google Scholar
J. Bodin, F. Delay, and G. de Marsily, “Solute transport in a single fracture with negligible matrix permeability: 2. Mathematical formalism,” Hydrogeology Journal, vol. 11, no. 4, pp. 434–454, 2003.
View at: Publisher Site | Google Scholar
F. Delay, P. Ackerer, and C. Danquigny, “Simulating solute transport in porous or fractured formations using random walk particle tracking: a review,” Vadose Zone Journal, vol. 4, no. 2, pp. 360–379, 2005.
View at: Publisher Site | Google Scholar
E. Stalgorova and T. Babadagli, “Field-scale modeling of tracer injection in naturally fractured reservoirs using the random-walk particle-tracking simulation,” SPE Journal, vol. 17, no. 2, pp. 580–592, 2012.
View at: Publisher Site | Google Scholar
M. Bechtold, J. Vanderborght, O. Ippisch, and H. Vereecken, “Efficient random walk particle tracking algorithm for advective-dispersive transport in media with discontinuous dispersion coefficients and water contents,” Water Resources Research, vol. 47, no. 10, Article ID W10526, 2011.
View at: Publisher Site | Google Scholar
B. A. Robinson, Z. V. Dash, and G. Srinivasan, “A particle tracking transport method for the simulation of resident and flux-averaged concentration of solute plumes in groundwater models,” Computational Geosciences, vol. 14, no. 4, pp. 779–792, 2010.
View at: Publisher Site | Google Scholar
G. Srinivasan, E. Keating, J. D. Moulton, Z. V. Dash, and B. A. Robinson, “Convolution-based particle tracking method for transient flow,” Computational Geosciences, vol. 16, no. 3, pp. 551–563, 2012.
View at: Publisher Site | Google Scholar
P. K. Kang, T. Le Borgne, M. Dentz, O. Bour, and R. Juanes, “Impact of velocity correlation and distribution on transport in fractured media: field evidence and theoretical model,” Water Resources Research, vol. 51, no. 2, pp. 940–959, 2015.
View at: Publisher Site | Google Scholar
S. Geiger, A. Cortis, and J. T. Birkholzer, “Upscaling solute transport in naturally fractured porous media with the continuous time random walk method,” Water Resources Research, vol. 46, no. 12, Article ID W12530, 2010.
View at: Publisher Site | Google Scholar
F. Delay and J. Bodin, “Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks,” Geophysical Research Letters, vol. 28, no. 21, pp. 4051–4054, 2001.
View at: Publisher Site | Google Scholar
V. Cvetkovic, A. Fiori, and G. Dagan, “Solute transport in aquifers of arbitrary variability: A time-domain random walk formulation,” Water Resources Research, vol. 50, no. 7, pp. 5759–5773, 2014.
View at: Publisher Site | Google Scholar
J. Bodin, “From analytical solutions of solute transport equations to multidimensional time-domain random walk (TDRW) algorithms,” Water Resources Research, vol. 51, no. 3, pp. 1860–1871, 2015.
View at: Publisher Site | Google Scholar
J. Bodin, G. Porel, F. Delay, F. Ubertosi, S. Bernard, and J.-R. de Dreuzy, “Simulation and analysis of solute transport in 2D fracture/pipe networks: The SOLFRAC program,” Journal of Contaminant Hydrology, vol. 89, no. 1-2, pp. 1–28, 2007.
View at: Publisher Site | Google Scholar
Itasca Consulting Group, UDEC User’s Guide, Itasca Consulting Group, Minneapolis, Minn, USA, 2004.
K. R. Bradbury and M. A. Muldoon, “Effects of fracture density and anisotropy on delineation of wellhead-projection areas in fractured-rock aquifers,” Applied Hydrogeology, vol. 2, no. 3, pp. 17–23, 1994.
View at: Publisher Site | Google Scholar
R. Parney and L. Smith, “Fluid velocity and path length in fractured media,” Geophysical Research Letters, vol. 22, no. 11, pp. 1437–1440, 1995.
View at: Publisher Site | Google Scholar
Y.-J. Park and K.-K. Lee, “Analytical solutions for solute transfer characteristics at continuous fracture junctions,” Water Resources Research, vol. 35, no. 5, pp. 1531–1537, 1999.
View at: Publisher Site | Google Scholar
J. P. Sauty, “An analysis of hydrodispersive transfer in aquifers,” Water Resources Research, vol. 16, no. 1, pp. 145–158, 1980.
View at: Publisher Site | Google Scholar
S. R. Yates, “An analytical solution for one‐dimensional transport in heterogeneous porous media,” Water Resources Research, vol. 26, no. 10, pp. 2331–2338, 1990.
View at: Publisher Site | Google Scholar
D. H. Tang, E. O. Frind, and E. A. Sudicky, “Contaminant transport in fractured porous media: analytical solution for a single fracture,” Water Resources Research, vol. 17, no. 3, pp. 555–564, 1981.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2017 Yuexiu Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1045

Downloads

1165

Citations