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

Authors: Yuexiu Wu; Quansheng Liu; Andrew H. C. Chan; Hongyuan Liu
Journal: Geofluids
 (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.


---

## Body

## 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 Section2; 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.


## 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 Section2; 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.Figure 1
Deformable blocks.


### 2.1. Mechanical Behavior

For any grid pointp in a triangle element, the total force is obtained as a sum of four terms:(1)Fp=Fpz+Fpc+Fpo+Fpg,where Fpo are the external applied loads, Fpz are the internal stresses contributed by zones adjacent to grid point p, Fpc are the contact forces from all the adjacent blocks, and Fpg are gravity forces. The internal stress is calculated using the following equation:(2)Fpz=∫Cσijnjds,where σij is element stress tensor and nj is the unit vector outward normal to the integral path C, which follows the closed polygonal line around the grid point p. The gravity force is calculated using the following equation:(3)Fpg=gpmp,where mp is the lumped gravitational mass, defined as the sum of one-third of the masses of triangles connected to the grid point, and gp is the gravitational acceleration.The contact forceFpc exists only when grid points are on block boundaries (like grid point n 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: Fpnc (normal contact force) and Fpsc (shear contact force):(4)Fpnc≔Fpnc-knΔunl,where Δun is the normal displacement increment for each contact, kn is joint normal stiffness, and l is the contact area, and(5)Fpsc≔Fpsc-ksΔusl,where Δus is the shear displacement increment and ks is joint shear stiffness. If the shear contact force exceeds the shear strength, it is modified according to the following equation:(6)Fpsc=sign⁡Δusc+tan⁡∅,where c is joint cohesion and ∅ is joint frictional strength.The motion of the grid point is obtained according to Newton’s second law:(7)du˙pdt=Fpmp,where u˙p is the velocity of grid point p and t 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:(8)Δεij=12ui,j+uj,iΔt,where Δt is time step.(9)Δσij=λΔεvδij+2μ-Δεij,where λ and μ- are Lame’s constant and Δεv 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 pressurep1 to a domain with pressure p2),(10)q=-kcΔp,where q is the flow rate, kc is a point contact permeability factor, and Δp is the pressure difference calculated using the following equation:(11)Δp=p2-p1+ρwgy2-y1,where ρw is the fluid density; g is the acceleration of gravity; and y2 and  y1 are y-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:(12)q=-kja3Δpl,where kj is a joint permeability factor (whose theoretical value is one-twelfth dynamic viscosity of fluid), a is the contact hydraulic aperture, and l 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):(13)uf=qA,where A 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 intoFpo. The external force is calculated using the following equation:(14)FpP=pniL,where p is the fluid pressure in a contact domain, ni is the normal direct of contact, and L is the contact length. With the motion of each block, the contact area will change following (15), which causes flow rate variation: (15)a=a0+Δun,where a0 is the initial fracture aperture and Δun is the normal displacement of contact area.


## 2.1. Mechanical Behavior

For any grid pointp in a triangle element, the total force is obtained as a sum of four terms:(1)Fp=Fpz+Fpc+Fpo+Fpg,where Fpo are the external applied loads, Fpz are the internal stresses contributed by zones adjacent to grid point p, Fpc are the contact forces from all the adjacent blocks, and Fpg are gravity forces. The internal stress is calculated using the following equation:(2)Fpz=∫Cσijnjds,where σij is element stress tensor and nj is the unit vector outward normal to the integral path C, which follows the closed polygonal line around the grid point p. The gravity force is calculated using the following equation:(3)Fpg=gpmp,where mp is the lumped gravitational mass, defined as the sum of one-third of the masses of triangles connected to the grid point, and gp is the gravitational acceleration.The contact forceFpc exists only when grid points are on block boundaries (like grid point n 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: Fpnc (normal contact force) and Fpsc (shear contact force):(4)Fpnc≔Fpnc-knΔunl,where Δun is the normal displacement increment for each contact, kn is joint normal stiffness, and l is the contact area, and(5)Fpsc≔Fpsc-ksΔusl,where Δus is the shear displacement increment and ks is joint shear stiffness. If the shear contact force exceeds the shear strength, it is modified according to the following equation:(6)Fpsc=sign⁡Δusc+tan⁡∅,where c is joint cohesion and ∅ is joint frictional strength.The motion of the grid point is obtained according to Newton’s second law:(7)du˙pdt=Fpmp,where u˙p is the velocity of grid point p and t 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:(8)Δεij=12ui,j+uj,iΔt,where Δt is time step.(9)Δσij=λΔεvδij+2μ-Δεij,where λ and μ- are Lame’s constant and Δεv 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 pressurep1 to a domain with pressure p2),(10)q=-kcΔp,where q is the flow rate, kc is a point contact permeability factor, and Δp is the pressure difference calculated using the following equation:(11)Δp=p2-p1+ρwgy2-y1,where ρw is the fluid density; g is the acceleration of gravity; and y2 and  y1 are y-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:(12)q=-kja3Δpl,where kj is a joint permeability factor (whose theoretical value is one-twelfth dynamic viscosity of fluid), a is the contact hydraulic aperture, and l 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):(13)uf=qA,where A 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 intoFpo. The external force is calculated using the following equation:(14)FpP=pniL,where p is the fluid pressure in a contact domain, ni is the normal direct of contact, and L is the contact length. With the motion of each block, the contact area will change following (15), which causes flow rate variation: (15)a=a0+Δun,where a0 is the initial fracture aperture and Δun 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: (1) fracture, (2) matrix, and (3) 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]:(16)∂cf∂t+λcf=-ufRf∂cf∂x+1Rf∂∂xDf∂cf∂x+DfbRf∂cm∂yy=b,where cf is the concentration of nuclide in fracture, uf is the fluid velocity in fracture, b is half of the aperture, cm is the concentration of nuclide in rock matrix, λ is the decay constant of nuclide, the term λcf represents the decay effect, Df is the hydrodynamic dispersion coefficient in fracture, Rf 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 x and y 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.Figure 2
Fracture with an aperture of 2b.The retardation factorRf is defined using the following equation [42]:(17)Rf=1+Kfb,where Kf is the surface sorption coefficient of nuclide onto the fracture surface. The nuclide transport equation in rock matrix can be written as follows [42]:(18)∂cm∂t+λcm=Da∂2cm∂y2,where Da is the apparent-diffusion coefficient of the nuclide in the matrix, which is calculated as Da=De/(θm+ρmKm); De is effective diffusion coefficient; θm is matrix porosity; ρm is matrix density; and Km 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: (1) advection-dispersion-sorption (ADS) transport equation, (2) decay effect, and (3) matrix diffusion effect. The ADS equation can be rewritten as (19) [42], considering the change of variable (∂cf/∂t=∂cf/∂x∂x/∂t=∂cf/∂xuf/Rf):(19)∂cf∂x=-Rfuf2∂∂tuf+∂Df∂xcf+Rf2uf3∂2∂t2Dfcf.By applying an equivalent change of variable, (19) can be written with Fokker-Planck formalism:(20)∂φ∂x=-1v2∂vφ∂t+12v3∂2∅φ∂t2.According to the solution of random-walk approach, the mean and variance of nuclide transport time distribution for a displacement of dx are obtained using (21a) and (21b), respectively [42]:(21a)μtdx=Rfuf2uf+∂Df∂xdx,(21b)σt2dx=Rf22Dfuf3dx.For Peclet number Pe=ufdx/Df larger than 10, it is shown that the travel time distribution is lognormal for a displacement of dx [30]:(22a)ln⁡Δtf=μln⁡+ZNσln⁡,(22b)μln=ln⁡μt1+σt2/μt2,(22c)σln2=ln⁡1+σt2μt2,where Δtf is the travel time of the nuclide for a displacement dx and ZN is a random number that follows a standard normal distribution.For very low flow velocities in short fractures (i.e., the Peclet numberPe is less than 10), the assumption of a lognormal distribution is not accurate. Bodin et al. [30] suggested that (22b) should be modified:(23)μln′=βμln=1-133Peln⁡μt1+σt2/μt2.The total travel time in a fracture of length dx can be written as(24)Δtf=exp⁡μln⁡+ZNσln⁡⁡forPe>10.When Peclet number Pe is less than 10, the variable μln is replaced by μln′.


#### 3.2.2. Matrix Diffusion

To account for the interaction between ADS and matrix diffusion in a bond of lengthdx, the fluid velocity uf and dispersion coefficient Df in (22a) and (22b) were replaced by an apparent fluid velocity uf∗ and an apparent dispersion coefficient Df∗ [42]:(25)uf∗=Rdiffuf,Df∗=fuf∗,likeDf∗=αuf∗,orDf∗=Df,where α is a dispersivity constant and Rdiff is a retardation factor accounting for matrix diffusion:(26)Rdiff=1+Ω2σtRfexp⁡-ξ2⁡πerfc⁡ξ-ξ,where ξ=Ωdx/uf2/σt, Ω=θm+ρmKmDe/2b, θm is the matrix porosity, ρm is the bulk density of matrix, Km is the volumetric-sorption coefficient of nuclide in matrix, and De is the effective diffusion coefficient of nuclide in matrix.The modified (22a) and (22b) with uf∗ and Df∗ in place of uf and Df, respectively, were used to calculate the mean and variance of the log transform μln (or μln′) and σln2 according to (22b) and (22c), respectively. The total travel time for the coupled ADS and matrix diffusion is calculated using (27) for Pe>10 [42]:(27)Δtfm=exp⁡μln+ZNσln+Ωt0Rferfc-1⁡U012,where t0 is the particle residence time in the fracture for pure advection delayed by sorption reactions onto the fracture walls, expressed as t0=Rfdx/uf and  U01 is a random number drawn from a uniform distribution between 0 and 1. For Pe<10, the variable μln in (27) is replaced by μln′.


#### 3.2.3. Decay Effect

If the linear decay model is used, the nuclide mass decreases according to the following equation [42]:(28)mn+1p=mnpexp⁡-λtn+1-tn,where mn+1p and mnp represent the masses of the nuclide at node n+1 and n, 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:(29)Pij=Qj∑Q-,where Pij is the probability of particle transition from an inlet bond i to an outlet bond j; Qj is the flow rate in the outlet bond j; and ∑Q- 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.


## 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: (1) fracture, (2) matrix, and (3) 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]:(16)∂cf∂t+λcf=-ufRf∂cf∂x+1Rf∂∂xDf∂cf∂x+DfbRf∂cm∂yy=b,where cf is the concentration of nuclide in fracture, uf is the fluid velocity in fracture, b is half of the aperture, cm is the concentration of nuclide in rock matrix, λ is the decay constant of nuclide, the term λcf represents the decay effect, Df is the hydrodynamic dispersion coefficient in fracture, Rf 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 x and y 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.Figure 2
Fracture with an aperture of 2b.The retardation factorRf is defined using the following equation [42]:(17)Rf=1+Kfb,where Kf is the surface sorption coefficient of nuclide onto the fracture surface. The nuclide transport equation in rock matrix can be written as follows [42]:(18)∂cm∂t+λcm=Da∂2cm∂y2,where Da is the apparent-diffusion coefficient of the nuclide in the matrix, which is calculated as Da=De/(θm+ρmKm); De is effective diffusion coefficient; θm is matrix porosity; ρm is matrix density; and Km 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: (1) advection-dispersion-sorption (ADS) transport equation, (2) decay effect, and (3) matrix diffusion effect. The ADS equation can be rewritten as (19) [42], considering the change of variable (∂cf/∂t=∂cf/∂x∂x/∂t=∂cf/∂xuf/Rf):(19)∂cf∂x=-Rfuf2∂∂tuf+∂Df∂xcf+Rf2uf3∂2∂t2Dfcf.By applying an equivalent change of variable, (19) can be written with Fokker-Planck formalism:(20)∂φ∂x=-1v2∂vφ∂t+12v3∂2∅φ∂t2.According to the solution of random-walk approach, the mean and variance of nuclide transport time distribution for a displacement of dx are obtained using (21a) and (21b), respectively [42]:(21a)μtdx=Rfuf2uf+∂Df∂xdx,(21b)σt2dx=Rf22Dfuf3dx.For Peclet number Pe=ufdx/Df larger than 10, it is shown that the travel time distribution is lognormal for a displacement of dx [30]:(22a)ln⁡Δtf=μln⁡+ZNσln⁡,(22b)μln=ln⁡μt1+σt2/μt2,(22c)σln2=ln⁡1+σt2μt2,where Δtf is the travel time of the nuclide for a displacement dx and ZN is a random number that follows a standard normal distribution.For very low flow velocities in short fractures (i.e., the Peclet numberPe is less than 10), the assumption of a lognormal distribution is not accurate. Bodin et al. [30] suggested that (22b) should be modified:(23)μln′=βμln=1-133Peln⁡μt1+σt2/μt2.The total travel time in a fracture of length dx can be written as(24)Δtf=exp⁡μln⁡+ZNσln⁡⁡forPe>10.When Peclet number Pe is less than 10, the variable μln is replaced by μln′.


### 3.2.2. Matrix Diffusion

To account for the interaction between ADS and matrix diffusion in a bond of lengthdx, the fluid velocity uf and dispersion coefficient Df in (22a) and (22b) were replaced by an apparent fluid velocity uf∗ and an apparent dispersion coefficient Df∗ [42]:(25)uf∗=Rdiffuf,Df∗=fuf∗,likeDf∗=αuf∗,orDf∗=Df,where α is a dispersivity constant and Rdiff is a retardation factor accounting for matrix diffusion:(26)Rdiff=1+Ω2σtRfexp⁡-ξ2⁡πerfc⁡ξ-ξ,where ξ=Ωdx/uf2/σt, Ω=θm+ρmKmDe/2b, θm is the matrix porosity, ρm is the bulk density of matrix, Km is the volumetric-sorption coefficient of nuclide in matrix, and De is the effective diffusion coefficient of nuclide in matrix.The modified (22a) and (22b) with uf∗ and Df∗ in place of uf and Df, respectively, were used to calculate the mean and variance of the log transform μln (or μln′) and σln2 according to (22b) and (22c), respectively. The total travel time for the coupled ADS and matrix diffusion is calculated using (27) for Pe>10 [42]:(27)Δtfm=exp⁡μln+ZNσln+Ωt0Rferfc-1⁡U012,where t0 is the particle residence time in the fracture for pure advection delayed by sorption reactions onto the fracture walls, expressed as t0=Rfdx/uf and  U01 is a random number drawn from a uniform distribution between 0 and 1. For Pe<10, the variable μln in (27) is replaced by μln′.


### 3.2.3. Decay Effect

If the linear decay model is used, the nuclide mass decreases according to the following equation [42]:(28)mn+1p=mnpexp⁡-λtn+1-tn,where mn+1p and mnp represent the masses of the nuclide at node n+1 and n, 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:(29)Pij=Qj∑Q-,where Pij is the probability of particle transition from an inlet bond i to an outlet bond j; Qj is the flow rate in the outlet bond j; and ∑Q- 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.


## 3.2.1. Advection-Dispersion-Sorption Equation

In TDRW method, the solution of (16) can be divided into three steps: (1) advection-dispersion-sorption (ADS) transport equation, (2) decay effect, and (3) matrix diffusion effect. The ADS equation can be rewritten as (19) [42], considering the change of variable (∂cf/∂t=∂cf/∂x∂x/∂t=∂cf/∂xuf/Rf):(19)∂cf∂x=-Rfuf2∂∂tuf+∂Df∂xcf+Rf2uf3∂2∂t2Dfcf.By applying an equivalent change of variable, (19) can be written with Fokker-Planck formalism:(20)∂φ∂x=-1v2∂vφ∂t+12v3∂2∅φ∂t2.According to the solution of random-walk approach, the mean and variance of nuclide transport time distribution for a displacement of dx are obtained using (21a) and (21b), respectively [42]:(21a)μtdx=Rfuf2uf+∂Df∂xdx,(21b)σt2dx=Rf22Dfuf3dx.For Peclet number Pe=ufdx/Df larger than 10, it is shown that the travel time distribution is lognormal for a displacement of dx [30]:(22a)ln⁡Δtf=μln⁡+ZNσln⁡,(22b)μln=ln⁡μt1+σt2/μt2,(22c)σln2=ln⁡1+σt2μt2,where Δtf is the travel time of the nuclide for a displacement dx and ZN is a random number that follows a standard normal distribution.For very low flow velocities in short fractures (i.e., the Peclet numberPe is less than 10), the assumption of a lognormal distribution is not accurate. Bodin et al. [30] suggested that (22b) should be modified:(23)μln′=βμln=1-133Peln⁡μt1+σt2/μt2.The total travel time in a fracture of length dx can be written as(24)Δtf=exp⁡μln⁡+ZNσln⁡⁡forPe>10.When Peclet number Pe is less than 10, the variable μln is replaced by μln′.


## 3.2.2. Matrix Diffusion

To account for the interaction between ADS and matrix diffusion in a bond of lengthdx, the fluid velocity uf and dispersion coefficient Df in (22a) and (22b) were replaced by an apparent fluid velocity uf∗ and an apparent dispersion coefficient Df∗ [42]:(25)uf∗=Rdiffuf,Df∗=fuf∗,likeDf∗=αuf∗,orDf∗=Df,where α is a dispersivity constant and Rdiff is a retardation factor accounting for matrix diffusion:(26)Rdiff=1+Ω2σtRfexp⁡-ξ2⁡πerfc⁡ξ-ξ,where ξ=Ωdx/uf2/σt, Ω=θm+ρmKmDe/2b, θm is the matrix porosity, ρm is the bulk density of matrix, Km is the volumetric-sorption coefficient of nuclide in matrix, and De is the effective diffusion coefficient of nuclide in matrix.The modified (22a) and (22b) with uf∗ and Df∗ in place of uf and Df, respectively, were used to calculate the mean and variance of the log transform μln (or μln′) and σln2 according to (22b) and (22c), respectively. The total travel time for the coupled ADS and matrix diffusion is calculated using (27) for Pe>10 [42]:(27)Δtfm=exp⁡μln+ZNσln+Ωt0Rferfc-1⁡U012,where t0 is the particle residence time in the fracture for pure advection delayed by sorption reactions onto the fracture walls, expressed as t0=Rfdx/uf and  U01 is a random number drawn from a uniform distribution between 0 and 1. For Pe<10, the variable μln in (27) is replaced by μln′.


## 3.2.3. Decay Effect

If the linear decay model is used, the nuclide mass decreases according to the following equation [42]:(28)mn+1p=mnpexp⁡-λtn+1-tn,where mn+1p and mnp represent the masses of the nuclide at node n+1 and n, 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:(29)Pij=Qj∑Q-,where Pij is the probability of particle transition from an inlet bond i to an outlet bond j; Qj is the flow rate in the outlet bond j; and ∑Q- 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 Figure3). 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.Figure 3
Flow path in fracture network.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 Figure4.Figure 4
Flow chart of TDRW method’s implementation.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 M0 and the total number of particles is N, mass of each particle is expressed as(30)mpi=M0N.If the collected particle number on the outlet boundary isNi in a period of time dt, the outlet concentration is written as(31)ci=mpiNi∑Qoutdt,where ∑Qout 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) intoN parts, in which the time interval is Δt.(3)
Collect particles (number isNi) 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:(32)cfc0=12erfc⁡x-uft2Dft-exp⁡ufxDferfc⁡x+uft2Dft,where c0 is the constant concentration injection; and erfc⁡(x) 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−4 m/s, hydrodynamic dispersion coefficient in the fracture was 2.0 × 10−5 m2/s, and the aperture of the fracture was 2.5 × 10−4 m. Thus, the Peclet number (Pe=ufL/Df) 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.Figure 5
Comparison between numerical and analytical solutions for nuclide particle transport in a single fracture.


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

The simple fracture network was shown in Figure6. 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−6m. 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.Figure 6
A simple fracture network.Figure 7
Comparison between numerical and analytical solutions of particle transport in a simple fracture network.From the comparisons shown in Figures5 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.


## 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:(32)cfc0=12erfc⁡x-uft2Dft-exp⁡ufxDferfc⁡x+uft2Dft,where c0 is the constant concentration injection; and erfc⁡(x) 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−4 m/s, hydrodynamic dispersion coefficient in the fracture was 2.0 × 10−5 m2/s, and the aperture of the fracture was 2.5 × 10−4 m. Thus, the Peclet number (Pe=ufL/Df) 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.Figure 5
Comparison between numerical and analytical solutions for nuclide particle transport in a single fracture.


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

The simple fracture network was shown in Figure6. 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−6m. 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.Figure 6
A simple fracture network.Figure 7
Comparison between numerical and analytical solutions of particle transport in a simple fracture network.From the comparisons shown in Figures5 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 Figure8, 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 (σv=σh=0) 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 (σv/σh=1,2,3,5). 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.Figure 8
Calculation model studied in the DECOVALEX 2011 project.Figure 9
Hydromechanical boundary condition.The block material property wasK=54.23 GPa, and G=34.11 GPa.The normal behavior of the fractures in the calculation model was nonlinear, as shown in Figure10 [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 μm and a second moment of 1.0. The minimum and maximum apertures were 1 μm and 200 μm, respectively. The residual aperture was assumed to be 20% of the initial aperture.Table 1
Material parameters of the fracture in the calculation model.
Mechanical properties
Values
Fractures
Shear stiffness,Ks (GPa/m)
434
Friction angle (°)
24.9
Dilation angle (°)
5
Cohesion (MPa)
0Figure 10
Relationship between normal displacement and normal stress.The hydrodynamic dispersion coefficient in fractures was 1.0 × 10−5 m2/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 Figures11–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−8 m3/(s·m) and 4.06 × 10−9 m3/(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−8 m3/(s·m) and 3.13 × 10−8 m3/(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−7 m3/(s·m) and 8.02 × 10−8 m3/(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.Figure 11
Flow rate distribution of fractures on the bottom outlet boundary.Figure 12
Flow rate distribution of fractures on the right outlet boundary.Figure 13
Flow rate distribution of fractures on the left outlet boundary.


#### 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 × 104 s, 2.58 × 104 s, 3.12 × 104 s, 1.95 × 105 s, and 2.4 × 105 s, respectively. The mean traveling time for each calculation was 3.29 × 103 s, 4.31 × 103 s, 5.21 × 103 s, 1.27 × 104 s, and 2.19 × 104 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.Figure 14
Breakthrough curve with different stress ratios without matrix diffusion.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 θm was 0.316%. For simplicity, advection, dispersion, and matrix diffusion were considered in this case, while sorption in the matrix was not allowed (Km=0 m3/kg). The effective diffusion coefficient of the nuclide in the matrix is 1.0 × 10−11 m2/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 × 107 s, 7.98 × 107 s, 1.086 × 108 s, 8.02 × 108 s, and 8.42 × 108 s, respectively. The mean traveling time for each calculation was 4.88 × 106 s, 1.77 × 107 s, 2.89 × 107 s, 5.2 × 107 s, and 7.7 × 107 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.Figure 15
Breakthrough curve with different stress ratios with matrix diffusion.Comparing the breakthrough curves in Figure14 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.


## 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 Figure8, 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 (σv=σh=0) 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 (σv/σh=1,2,3,5). 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.Figure 8
Calculation model studied in the DECOVALEX 2011 project.Figure 9
Hydromechanical boundary condition.The block material property wasK=54.23 GPa, and G=34.11 GPa.The normal behavior of the fractures in the calculation model was nonlinear, as shown in Figure10 [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 μm and a second moment of 1.0. The minimum and maximum apertures were 1 μm and 200 μm, respectively. The residual aperture was assumed to be 20% of the initial aperture.Table 1
Material parameters of the fracture in the calculation model.
Mechanical properties
Values
Fractures
Shear stiffness,Ks (GPa/m)
434
Friction angle (°)
24.9
Dilation angle (°)
5
Cohesion (MPa)
0Figure 10
Relationship between normal displacement and normal stress.The hydrodynamic dispersion coefficient in fractures was 1.0 × 10−5 m2/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 Figures11–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−8 m3/(s·m) and 4.06 × 10−9 m3/(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−8 m3/(s·m) and 3.13 × 10−8 m3/(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−7 m3/(s·m) and 8.02 × 10−8 m3/(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.Figure 11
Flow rate distribution of fractures on the bottom outlet boundary.Figure 12
Flow rate distribution of fractures on the right outlet boundary.Figure 13
Flow rate distribution of fractures on the left outlet boundary.


### 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 × 104 s, 2.58 × 104 s, 3.12 × 104 s, 1.95 × 105 s, and 2.4 × 105 s, respectively. The mean traveling time for each calculation was 3.29 × 103 s, 4.31 × 103 s, 5.21 × 103 s, 1.27 × 104 s, and 2.19 × 104 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.Figure 14
Breakthrough curve with different stress ratios without matrix diffusion.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 θm was 0.316%. For simplicity, advection, dispersion, and matrix diffusion were considered in this case, while sorption in the matrix was not allowed (Km=0 m3/kg). The effective diffusion coefficient of the nuclide in the matrix is 1.0 × 10−11 m2/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 × 107 s, 7.98 × 107 s, 1.086 × 108 s, 8.02 × 108 s, and 8.42 × 108 s, respectively. The mean traveling time for each calculation was 4.88 × 106 s, 1.77 × 107 s, 2.89 × 107 s, 5.2 × 107 s, and 7.7 × 107 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.Figure 15
Breakthrough curve with different stress ratios with matrix diffusion.Comparing the breakthrough curves in Figure14 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.


## 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 Figures11–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−8 m3/(s·m) and 4.06 × 10−9 m3/(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−8 m3/(s·m) and 3.13 × 10−8 m3/(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−7 m3/(s·m) and 8.02 × 10−8 m3/(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.Figure 11
Flow rate distribution of fractures on the bottom outlet boundary.Figure 12
Flow rate distribution of fractures on the right outlet boundary.Figure 13
Flow rate distribution of fractures on the left outlet boundary.


## 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 × 104 s, 2.58 × 104 s, 3.12 × 104 s, 1.95 × 105 s, and 2.4 × 105 s, respectively. The mean traveling time for each calculation was 3.29 × 103 s, 4.31 × 103 s, 5.21 × 103 s, 1.27 × 104 s, and 2.19 × 104 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.Figure 14
Breakthrough curve with different stress ratios without matrix diffusion.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 θm was 0.316%. For simplicity, advection, dispersion, and matrix diffusion were considered in this case, while sorption in the matrix was not allowed (Km=0 m3/kg). The effective diffusion coefficient of the nuclide in the matrix is 1.0 × 10−11 m2/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 × 107 s, 7.98 × 107 s, 1.086 × 108 s, 8.02 × 108 s, and 8.42 × 108 s, respectively. The mean traveling time for each calculation was 4.88 × 106 s, 1.77 × 107 s, 2.89 × 107 s, 5.2 × 107 s, and 7.7 × 107 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.Figure 15
Breakthrough curve with different stress ratios with matrix diffusion.Comparing the breakthrough curves in Figure14 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.



---

*Source: 5940380-2017-08-15.xml*