Radionuclide Transport in Fractured Rock: Numerical Assessment for High Level Waste Repository

Abstract

Deep and stable geological formations with low permeability have been considered for high level waste definitive repository. A common problem is the modeling of radionuclide migration in a fractured medium. Initially, we considered a system consisting of a rock matrix with a single planar fracture in water saturated porous rock. Transport in the fracture is assumed to obey an advection-diffusion equation, while molecular diffusion is considered the dominant mechanism of transport in porous matrix. The partial differential equations describing the movement of radionuclides were discretized by finite difference methods, namely, fully explicit, fully implicit, and Crank-Nicolson schemes. The convective term was discretized by the following numerical schemes: backward differences, centered differences, and forward differences. The model was validated using an analytical solution found in the literature. Finally, we carried out a simulation with relevant spent fuel nuclide data with a system consisting of a horizontal fracture and a vertical fracture for assessing the performance of a hypothetical repository inserted into the host rock. We have analysed the bentonite expanded performance at the beginning of fracture, the quantified radionuclide released from a borehole, and an estimated effective dose to an adult, obtained from ingestion of well water during one year.

1. Introduction

With the advent of nuclear technology, countries that employ it received many benefits from the energy production increase, but, as occurs to any activity related to industrialization, there is inevitably a waste production, but the volume of radioactive waste produced is very small, when compared with other industrial residues. The relatively long half-lives of some radionuclides present in this material, as well as their high activity and toxicity, pose a burden for the present and future generation. Part of this waste, the high level waste resulting from reprocessing and spent fuel disposal (HLW—high level waste), remains hazardous for thousands of years and therefore should be carefully stored, using barriers that delay the migration of radioisotopes to the biosphere.

Thus, the greatest challenge involves the management of high level radioactive waste, which must be isolated from the environment, and a disposal site for this waste, chosen to ensure safety for many thousands of years, OECD [1]. Although it is impossible to establish a safe depth, for which the biosphere could be effectively isolated from the buried radioactive wastes, it is recognized that a deep repository will have less risk of radionuclide releases to the biosphere due to human intrusion, animals, plants or other mechanisms, so the current trend is to have radioactive waste buried in deep and geologically stable host rock. In this context, it is very important to focus on the issue of groundwater movement, because the mechanism by which radionuclides in a fractured repository could return to surface is precisely groundwater movement. Thus, it appears that, besides the problems of possible hydrological area disposal of radioactive wastes, it is necessary to know how the lithological characteristics of the site can absorb/adsorb radionuclides transported and establish their mineralogical and chemical composition, so that the geological environment can act as a chemical barrier. Moreover, as most rocks allow for the passage of water, the chemical barrier that they can eventually present becomes more important than any physical barrier that may be provided, CoRWM [2].

Currently, the main design of final spent nuclear fuel repository is based on KBS-3 method, KBS-3 [3], which implies that the high level radionuclides are placed in a repository having a depth in the host rock of approximately 500 m, inside a copper/steel canister surrounded by compacted bentonite clay, as shown in Figure 1.

Figure 1 

The KBS-3 method for spent nuclear fuel.

The KBS-3 method for storage of spent nuclear fuel is based on a multibarriers system ensuring that, if any barrier should fail, another will mitigate radionuclide releases to the biosphere. Thus, the engineering barriers and geological environment surrounding the waste work together to provide a variety of safety functions, which control any potential releases of radioactivity from the repository and their movement through the host rock, Falck and Nilsson [4]. Focusing on the far field region, the host rock is fissured along its depth and groundwater flows into these fractures. Large and long zone, of the order of 1 km, generally vertical or almost vertical, are observed where the rock is more fissured than the rest of host rock and can have upward flow. These areas are called lineaments.

In evolution scenarios, the artificial barriers are assumed to fail in order to make predictions of radionuclide releases. The volume of rock that separates the repository from fracture zones and adjacent lineaments is of extreme importance to delay the radionuclide migration following fail in the engineering barrier system. Thus, the performance evaluation of the final spent nuclear fuel repository includes assessing the transport of radionuclides that accidentally escape through these barriers. The geosphere will serve as the outer barrier, which will slow the radionuclides migration to the biosphere. A thorough understanding of the migration mechanisms and the effect of spatial heterogeneity of rock properties on radionuclide transport in fractures are therefore essential for the evaluation of performance, Jacob [5].

The formulation in this study employed the model described in Li and Chiou [6] and Tang et al. [7] but inserting into the model the decay chain. For the numerical evaluation, this work makes use of some possible methods applied to partial differential equations that describe the radionuclides transport in a discrete fracture. The governing equations were discretized by finite differences technique, using the following schemes: fully explicit, fully implicit, and Crank-Nicolson. For each of these schemes, the advective term of the equation was discretized by the following numerical schemes: backward differences, centered differences, and forward differences and all of these discretizations were implemented in a Fortran 90 code.

To validate the numerical model of the decay chain proposed by finite differences technique, we used the analytical solution of Bauer et al. [8], who developed a solution in the Laplace domain using recursive formulas for different retardation factors in the chain. These recursive formulas can be used to construct transport solution for more complex multispecies, but the computational effort to obtain a solution increases as the chain gets longer.

Until the validation of numerical modeling, this work considered that the matrix was infinite. Although Neretnieks [9] has demonstrated that diffusion in the matrix is an important mechanism that will retard the transport of radionuclides in the rock, considering it infinite may overestimate its retention capacity. Thus, a conservative assumption used for performance evaluation was to consider a zone change 1 cm thick bordering the fracture and a nonpermeable boundary at a distance of 10 cm. This analysis intentionally tends to be conservative, in such a way that errors in the process increase safety analysis margin.

To evaluate safety and system deposition performance, we present a model considering that the initial fracture region is filled with bentonite which expanded the engineering barrier system, in order to analyze how the bentonite can hinder the mobility of radionuclides and which of them are more susceptible. Finally, we present a migration analysis and quantification of radionuclides released from a borehole and an estimate of the effective dose to an adult by ingestion of well water during one year.

2. Conceptual Model

Initially, let us consider a single thin fracture located in a water-saturated porous rock. A contaminant source containing radionuclides is released instantaneously at the origin of the fracture, at time 𝑡=0. In order to allow for a one-dimensional analysis of transport along the fracture, coupled to one-dimensional molecular diffusion of the radionuclide from the fracture into the porous matrix, we will make several assumptions pertaining to the geometry and hydraulic properties of the system. The considerations made for geometry and hydraulic properties of the problem are summarized as follows.(1)The width of the fracture is much smaller than its length.(2)The fracture water is well mixed so that the concentration in the fracture is uniform in any plane that is perpendicular to the flow direction.(3)Transport along the fracture is much faster than transport within the rock matrix.(4)Transport in the porous rock matrix is controlled by molecular diffusion because the intact rock has low hydraulic conductivity.(5)The water flow in fracture is laminar.

In view of Assumptions 1 and 2, the mass transfer along the fracture can be described by a one-dimensional rectangular coordinate system. In view of Assumptions 3 and 4, the mass flux within the rock matrix can be taken in the direction perpendicular to the fracture plane. Assumption 5 allows us to use Darcy’s law, Chambre et al. [10]. Based on this law, the flow is described by a mean linear velocity which carries the dissolved substance (contaminants) by advection. From a geochemical point of view, an instant local chemical equilibrium is assumed between water and the solute.

The physical system is idealized as shown in Figure 2. It was considered that the rock contains a fracture plane with a width of 2𝑏. The surrounding rock is simplified as a fully saturated equivalent porous medium with an effective porosity of 𝜃. The water velocity is assumed constant with time and uniform in space within the fracture and zero within rock matrix. The processes considered are as follows:(1)advective transport along the fracture;(2)longitudinal hydrodynamic dispersion in the fracture, in the direction of the fracture axis;(3)molecular diffusion in the rock matrix;(4)adsorption onto the fracture surface;(5)adsorption in the rock matrix;(6)radioactive chain decay.

Figure 2 

Geometry of the model with a single, planar, infinite fracture.

Hydrodynamic dispersion means the spread of particles resulting from both mechanical dispersion and molecular diffusion. Mechanical dispersion is caused by velocity variations and molecular diffusion is caused by the concentration gradient. In the case of uniform flow and an isotropic medium, the hydrodynamic dispersion coefficient 𝐷_𝑓 can be expressed as 𝐷_𝑓=𝛼_𝐿𝑣+𝐷^∗, (1) where 𝛼_𝐿 is dispersivity in the directions of the fracture axis (L), 𝑣 is water velocity in the fracture (L/T), and 𝐷^∗ is coefficient of molecular diffusion in water (L²/T).

Inside rock, only the matrix diffusion is considered because the permeability is low. Also, longitudinal diffusion in the rock matrix is so slow that it is assumed not to contribute to the migration of solutes in a porous medium. Since the fracture is fully filled with groundwater, no sorption is considered in a fracture except for the one onto the interface of a fracture and a surrounding medium. Inside a rock, rather strong sorption occurs.

3. Governing Equations

The transport of each member of a decay chain can be described by two coupled, one-dimensional equations on the basis of the above assumptions, one for the fracture and one for the porous matrix. The coupling is provided by the continuity of concentration along the fracture-matrix interface.

Since the goal is to describe the transport of a decay chain, a second coupling is necessary, which relates the mass gain of the daughter product to the mass loss of its parent by means of the process of radioactive disintegration. The distribution of the chain of solute species is thus described by two systems of equations, one system representing advective transport along the fracture and the other system representing diffusive transport in the porous matrix.

By the theory of mass balance and considerations made for the problem, the system of equations describing transport along the fracture is 𝜕𝐶_𝑓,1(𝑧,𝑡)−𝐷𝜕𝑡𝑅^𝑓,1_𝑓,1𝜕²𝐶_𝑓,1(𝑧,𝑡)𝜕𝑧²+𝑣𝑅_𝑓,1𝜕𝐶_𝑓,1(𝑧,𝑡)𝜕𝑧+𝜆₁𝐶_𝑓,1()+𝑞𝑧,𝑡¹(𝑧,𝑡)𝑅_𝑓,1𝑏=0,𝜕𝐶_𝑓,𝑘(𝑧,𝑡)−𝐷𝜕𝑡𝑅^𝑓,𝑘_𝑓,𝑘𝜕²𝐶_𝑓,𝑘(𝑧,𝑡)𝜕𝑧²+𝑣𝑅_𝑓,𝑘𝜕𝐶_𝑓,𝑘(𝑧,𝑡)𝜕𝑧+𝜆_𝑘𝐶_𝑓,𝑘()−𝑅𝑧,𝑡^{𝑓,𝑘−1}𝜆_𝑘−1𝐶_{𝑓,𝑘−1}(𝑧,𝑡)𝑅_𝑓,𝑘+𝑞_𝑘(𝑧,𝑡)𝑅_𝑓,𝑘𝑏=0for𝑧>0,𝑡>0, (2) where 𝐶_𝑓,𝑘(𝑧,𝑡) is the concentration of radionuclide 𝑘 in fracture water (M/L³), 𝑧 is the spatial coordinate along the fracture (L), 𝑡 is the time (T), 2𝑏 is the fracture aperture (L), 𝑣 is the groundwater velocity in the fracture (L/T), and 𝐷_𝑓,𝑘 is the longitudinal dispersion coefficient of solute 𝑘 (L²/T).

The fourth term on the left-hand side of (2) represents the mass loss of member 𝑘 by decay, and the fifth term in the 𝑘th equation of (2) denotes the mass gain of member 𝑘 by radioactive disintegration of its parent.

The radioactive decay constant, 𝜆_𝑘 (1/T), is defined by 𝜆_𝑘=ln2𝑡^𝑘_1/2, (3) where 𝑡^𝑘_1/2 is the half-life of the 𝑘th radioactive contaminant.

The retardation coefficient 𝑅_𝑓,𝑘 represents the effect of solute adsorption at the face of the porous rock. It is assumed to be governed by a linear adsorption isotherm which is defined as 𝑅_𝑓,𝑘𝐾=1+^𝑓,𝑘𝑏, (4) where 𝐾_𝑓,𝑘 is the surface distribution coefficient of a solute per unit area of the fracture rock interface over the unit volume (L).

The diffusive transport of the radionuclide chain in the porous matrix is governed by the following system of coupled diffusion equations: 𝜕𝐶_𝑝,1𝑧,𝑦,𝑡−𝐷𝜕𝑡𝑅^𝑝,1_𝑝,1𝜕²𝐶_𝑝,1𝑧,𝑦,𝑡𝜕𝑦²+𝜆₁𝐶_𝑝,1𝑧,𝑦,𝑡=0,𝜕𝐶_𝑝,𝑘𝑧,𝑦,𝑡−𝐷𝜕𝑡𝑅^𝑝,𝑘_𝑝,𝑘𝜕²𝐶_𝑝,𝑘𝑧,𝑦,𝑡𝜕𝑦²+𝜆_𝑘𝐶_𝑝,𝑘−𝑅𝑧,𝑦,𝑡_{𝑝,𝑘−1}𝜆_𝑘−1𝐶_{𝑝,𝑘−1}𝑧,𝑦,𝑡𝑅_𝑝,𝑘𝑦=0,for≥𝑏,𝑧>0,𝑡>0, (5) where 𝐶_𝑝,𝑘(𝑧,𝑦,𝑡) is the concentration of radionuclide 𝑘 in the porous matrix (M/L³) and 𝑦 is the spatial coordinate perpendicular to the fracture axis (L).

The matrix retardation coefficient 𝑅_𝑝,𝑘 represents the adsorptive loss of solute within the porous matrix, and it is defined as 𝑅_𝑝,𝑘𝜌=1+^𝑝𝐾_𝑝,𝑘𝜃, (6) where 𝜌_𝑝 is the density of the rock matrix, (M/L³), 𝜃 is the porosity of the rock, and 𝐾_𝑝,𝑘 is the distribution coefficient in the rock matrix, (L³/M).

The effective molecular diffusive coefficient 𝐷_𝑝,𝑘 (L²/T) is less than the free-solution value, 𝐷^𝑘_∗, because of the tortuous diffusion pathway. For an anisotropic medium, it can be represented by 𝐷_𝑝,𝑘=𝜏𝐷^𝑘_∗, (7) where 𝜏 is the tortuosity of the rock matrix.

The loss term 𝑞_𝑘(𝑧,𝑡) in the system of fracture equations (2) is equal to the diffusive flux of nuclide 𝑘 crossing the fracture-matrix boundary (M/L²T), and it is expressed in the form of Fick’s first law as 𝑞_𝑘(𝑧,𝑡)=−𝜃𝐷_𝑝,𝑘𝜕𝐶_𝑝,𝑘𝜕𝑦_𝑦=𝑏for𝑧>0,𝑡>0. (8)

This last equation provides a coupling between these two governing equations (2) and (5). The porosity 𝜃 is included in (8) to account for the portion of rock matrix that is open for diffusion owing to the presence of the solid phase.

The initial conditions for (2) and (5) are 𝐶_𝑓,𝑘()𝐶𝑧,0=0,for𝑧>0,_𝑝,𝑘𝑧,𝑦,0=0,for𝑦≥𝑏,𝑧>0. (9)

The boundary conditions in the fracture are given by 𝐶_𝑓,𝑘(𝑧,𝑡)=𝐶^𝑘₀𝐶for𝑧=0,𝑡>0,_𝑓,𝑘(∞,𝑡)=0,for𝑡>0, (10) where 𝐶^𝑘₀ is the concentration of 𝑘th radionuclide that is deposited instantaneously at the coordinate origin of the fracture (M/L³).

The boundary conditions in the rock matrix are 𝐶_𝑝,𝑘(𝑧,𝑏,𝑡)=𝐶_𝑓,𝑘()𝐶𝑧,𝑡,for𝑧>0,𝑡>0,_𝑝,𝑘(𝑧,∞,𝑡)=0,for𝑧>0,𝑡>0. (11)(12)

Equation (11) expresses the coupling between the concentration in the fracture and the concentration in the porous matrix.

4. Numerical Solution via Finite Differences

The governing equations presented in the previous section were numerically approximated by finite difference techniques. For both (2) and (5) the terms with second order spatial derivatives were discretized with centered difference schemes and derivatives in time were discretized with the well-known theta method. Besides, the advective term present only in (1) was discretized with backward, centered, and also forward differences. With these considerations, (2) can be rewritten as 𝐶^𝑛+1_{𝑓,𝑘,𝑖}−𝐶^𝑛_{𝑓,𝑘,𝑖}=(𝐷Δ𝑡1−Θ)^𝑓,𝑘𝑅_𝑓,𝑘𝐶^𝑛_{𝑓,𝑘,𝑖−1}−2𝐶^𝑛_{𝑓,𝑘,𝑖}+𝐶^𝑛_{𝑓,𝑘,𝑖+1}Δ𝑧²−𝑣𝑅_𝑓,𝑘𝐶^𝑛_{𝑓,𝑘,𝑖+𝑟}−𝐶^𝑛_{𝑓,𝑘,𝑖−𝑠}(𝑟+𝑠)Δ𝑧−𝜆_𝑘𝐶^𝑛_{𝑓,𝑘,𝑖}+𝑅_{𝑓,𝑘−1}𝜆_𝑘−1𝐶^𝑛_{𝑓,𝑘−1,𝑖}𝑅_𝑓,𝑘−1𝑅_𝑓,𝑘𝑏𝑞^𝑛_𝑘,𝑖𝐷+Θ^𝑓,𝑘𝑅_𝑓,𝑘𝐶^𝑛+1_{𝑓,𝑘,𝑖−1}−2𝐶^𝑛+1_{𝑓,𝑘,𝑖}+𝐶^𝑛+1_{𝑓,𝑘,𝑖+1}Δ𝑧²−𝑣𝑅_𝑓,𝑘𝐶^𝑛+1_{𝑓,𝑘,𝑖+𝑟}−𝐶^𝑛+1_{𝑓,𝑘,𝑖−𝑠}(𝑟+𝑠)Δ𝑧−𝜆_𝑘𝐶^𝑛+1_{𝑓,𝑘,𝑖}+𝑅_{𝑓,𝑘−1}𝜆_𝑘−1𝐶^𝑛+1_{𝑓,𝑘−1,𝑖}𝑅_𝑓,𝑘−1𝑅_𝑓,𝑘𝑏𝑞^𝑛+1_𝑘,𝑖, (13) where for 𝑟=𝑠=1 the advective term indicates centered differences, for 𝑟=0 and 𝑠=1 backward differences (upwind), and 𝑟=1 and 𝑠=0 forward differences. For diffusion in the rock matrix, (5) is rewritten in finite difference form as 𝐶^𝑛+1_{𝑝,𝑖,𝑗}−𝐶^𝑛_{𝑝,𝑖,𝑗}=(𝐷Δ𝑡1−Θ)^𝑝𝑅_𝑝𝐶^𝑛_{𝑝,𝑖,𝑗−1}−2𝐶^𝑛_{𝑝,𝑖,𝑗}+𝐶^𝑛_{𝑝,𝑖,𝑗+1}Δ𝑦²−𝜆𝐶^𝑛_{𝑝,𝑖,𝑗}+𝑅_{𝑝,𝑘−1}𝜆_𝑘−1𝐶^𝑛_{𝑝,𝑘−1,𝑖,𝑗}𝑅_𝑝,𝑘𝐷+Θ^𝑝𝑅_𝑝𝐶^𝑛+1_{𝑝,𝑖,𝑗−1}−2𝐶^𝑛+1_{𝑝,𝑖,𝑗}+𝐶^𝑛+1_{𝑝,𝑖,𝑗+1}Δ𝑦²−𝜆𝐶^𝑛+1_{𝑝,𝑖,𝑗}+𝑅_{𝑝,𝑘−1}𝜆_𝑘−1𝐶^𝑛+1_{𝑝,𝑘−1,𝑖,𝑗}𝑅_𝑝,𝑘. (14)

In (13) and (14), subscripts 𝑖 and 𝑗 denote spatial nodes in 𝑧 and 𝑦 directions, respectively, and index 𝑛 denotes discretization in time. The parameter Θ is a real constant that varies between 0 and 1, Anderson and Tannehill [11]. For Θ=0 and Θ=1, the explicit and implicit Euler methods, respectively, are obtained. For Θ=1/2, the Crank-Nicolson semi-implicit method is obtained.

5. Numerical Solution Performance

5.1. Validation of Numerical Method against an Analytical Solution

For validation of numerical modeling of decay chain, proposed by finite difference, a code in Fortran 90 was implemented, and the analytical solution developed in Bauer et al. [8] was used for comparison. We used the set of data presented in Table 1.

Table 1 

Physical parameters used in simulation.

Figure 3 illustrates the best numerical method for the problem, which was implicit Euler with forward discretization of the advective term. These results do not supersede the other discretizations for the advective term, since the forward discretization has stability conditions that must be obeyed, Silveira [12]. The spatial and temporal meshes used were Δ𝑧=0,15 m and Δ𝑡=1 day. In this figure, we can see concentration distributions for a hypothetical three-member decay chain. It may be noted that the different members’ concentration successively rises to a maximum and then vanishes, and that the analytical solution confirms the numerical modeling.

Figure 3 

Concentration profiles 𝑐1, 𝑐2, and 𝑐3.

5.2. Performance Assessment of a Repository

The analysis was performed for a dual porosity model in geosphere for radionuclide flow and transport, so we estimated the performance of bentonite in the far field region. Some modifications in problem geometry were made, as shown in Figure 4.

Figure 4 

Modify geometry of proposed model.

The diffusive transport of the radionuclide chain in the bentonite extruded region is governed by the following system of coupled diffusion equations: 𝜕𝐶_𝑏,1(𝑧,𝑡)−𝐷𝜕𝑡𝑅^𝑏,1_𝑏,1𝜕²𝐶_𝑏(𝑧,𝑡)𝜕𝑧²+𝜆₁𝐶_𝑏,1(𝑧,𝑡)=0,𝜕𝐶_𝑏,𝑘(𝑧,𝑡)−𝐷𝜕𝑡𝑅^𝑏,𝑘_𝑏,𝑘𝜕²𝐶_𝑏,𝑘(𝑧,𝑡)𝜕𝑧²+𝜆_𝑘𝐶_𝑏,𝑘()−𝑅𝑧,𝑡_{𝑏,𝑘−1}𝜆_𝑘−1𝐶_{𝑏,𝑘−1}(𝑧,𝑡)𝑅_𝑏,𝑘=0,for𝑧>0,𝑡>0, (15) where 𝐶_𝑏,𝑘(𝑧,𝑡) is the concentration of radionuclide 𝑘 in the bentonite (M/L³) and 𝑧 is spatial coordinate along the fracture (L). The retardation bentonite coefficient 𝑅_𝑏,𝑘 represents the adsorptive loss of solute within the bentonite, and it is defined as 𝑅_𝑏,𝑘=1+1−𝜃𝜃𝜌𝐾_𝑏,𝑘, (16) where 𝜌 is the density of the bentonite, (M/L³), 𝜃 is the porosity of the bentonite, and 𝐾_𝑏,𝑘 is the distribution coefficient in the bentonite, (L³/M).

A conservative assumption used in the calculation was a 1 cm thick altered zone next to the fracture and a nonpermeable boundary at a distance of 10 cm. Table 2 (in the appendix) gives the porosity and the pore diffusion coefficient values, 𝜃 and 𝐷𝑝, extracted from Nykyri et al. [13], for the two rock layers used; one layer (1 cm thick) represents the rock matrix close to the fracture surface (subscripts 1) and the second (9 cm thick) represents the rock matrix further away from the surface (subscript 2). The entire rock matrix thickness available for sorption would thus be 10 cm. The physical parameters for each element near the fracture and the bentonite are also given in Table 2.

Table 2 

Physical parameters for elements in this study: distribution coefficient in the bentonite (𝐾_𝑏,𝑘), bentonite diffusion coefficient (𝐷_𝑏), bentonite porosity (𝜃_𝑏), distribution coefficient in the rock matrix (𝐾_𝑝,𝑘), pore diffusion coefficient (𝐷_𝑝), and porosity of the rock matrix (𝜃).

The model parameters depending only on environmental characteristics are shown in Table 3 (see the appendix). The water velocity in the fracture, the longitudinal dispersion coefficient of the solute, the retardation coefficient, and diameter of the fracture are presented in Li and Chiou [6].

Table 3 

Constant parameters for all nuclides in model: longitudinal dispersion coefficient (𝐷_𝑓), retardation coefficient (𝑅_𝑓), density of the rock (𝜌), groundwater velocity (𝑣), lineament velocity (𝑣_lin), and fracture aperture (2𝑏).

The choice of the discretized mesh was based on stability conditions of the model and its validation against analytical solution, so implicit Euler method with forward discretization in advective term was used. For every nuclide, the spatial meshes Δ𝑧=0.5 m and Δ𝑦=0.001 m were adopted, where the maximum lengths used in the simulation were 50 m for the 𝑧 direction and 0.1 m for the 𝑦 direction. A time step Δ𝑡=1 year was used. For the bentonite area, an extrusion length of 1 cm was used. In this region, only diffusion occurs and the spatial mesh Δ𝑧_bent=0.001 m was used.

The radionuclide inventory in spent nuclear fuel varies depending on, for example, the fuel type, the burnup, the cooling time, and so forth. In this study, we considered the inventory of relevant elements: C¹⁴, Cl³⁶, Cs¹³⁵, Cs¹³⁷, I¹²⁹, Mo⁹³, Nb⁹⁴, Ni⁵⁹, Ni⁶³, Pd¹⁰⁷, Sm¹⁵¹, Se⁷⁹, Sn¹²⁶, Sr⁹⁰, Tc⁹⁹, and Zr⁹³ and the following decay chains: 4𝑛 (Pu²⁴⁰→U²³⁶), 4𝑛+1 (Cm²⁴⁵→Am²⁴¹→Np²³⁷→U²³³→Th²²⁹), 4𝑛+2 (Cm²⁴⁶→Pu²⁴²→U²³⁸→U²³⁴→Th²³⁰→Ra²²⁶), and 4𝑛+3 (Am²⁴³→Pu²³⁹→U²³⁵→Pa²³¹), in a deposition canister with spent nuclear fuel from type TVEL VVER-440 (PWR) with 1.44 tU and enrichment 4.0 wt%. The assumed burnup is 40 MWd/kgU. It was considered that the spent fuel cooling remained for 30 years. These data are in Anttila [14], where the composition and radioactive characteristics of the spent fuel were calculated using the code ORIGEN-S of the SCALE-5 software package.

For analysis of vertical migration, the lineament is regarded as an area of 3.5 m width with 7 fractures, each one with diameter equal to a half of horizontal fracture diameter and the velocity is calculated as an average value for the function given by Neretnieks [9]. We used in this lineament the parameters of the porous rock nearest horizontal fracture. In this region, the maximum lengths of 25 cm for diffusion and 475 m to the fracture system are used in the simulation. For the concentration at beginning of lineament we considered the concentration of outlet horizontal fracture but considering that there was an imbalance between the outflow in fracture and the inflow in the lineament, that should be compensated by dilution, once the system is in a saturated medium. In this case, continuous mixing exists and, to calculate the final concentration, the volume must be replaced by the flow rate. Thus, we have 𝐶_lin()=𝐶1,𝑡^𝑓(50,𝑡)𝑄_𝑓𝑄_lin (17) because the calculation for the mixture is a point mass balance. The product between flow rate and concentration is the mass flow per unit time.

The choice of the discretized mesh was based on stability conditions, so implicit Euler method with backward discretization in the advective term (unconditionally stable) was used. For every nuclide, the spatial meshes Δ𝑦=0.5 m and Δ𝑧=0.01 m were adopted, where the maximum length used in the simulation was 475 m in 𝑦 direction and 0.25 m in the 𝑧 direction.

5.2.1. Results

Figure 5 shows the concentration profile for the first group of relevant nuclides in high level radioactive waste: C¹⁴, Cl³⁶, Cs¹³⁵, Cs¹³⁷, I¹²⁹, Mo⁹³, Nb⁹⁴, and Ni⁵⁹.

(a)

(b)

(c)

(d)

Figure 5 

Concentration profiles in extruded bentonite region for 102, 103, 104, and 105 years.

As we can see, these nuclides are little affected by sorption on bentonite. However, there is clearly a decrease in the concentration of nuclides of lower half-life, such as C¹⁴ (with a half-life equal to 5.7 · 10³ years), Mo⁹³ (with a half-life equal to 4.0 · 10³ years), and Cs¹³⁷ (with a half-life equal to 3.0 · 10¹ years). Cs¹³⁷ is exhausted for simulation times above 10000 years.

Figure 6 shows the concentration profile for the second group of relevant nuclides in high level radioactive waste: Ni⁶³, Pd¹⁰⁷, Sm¹⁵¹, Se⁷⁹, Sn¹²⁶, Sr⁹⁰, Tc⁹⁹, and Zr⁹³ in the extruded bentonite region.

(a)

(b)

(c)

(d)

Figure 6 

Concentration profiles in extruded bentonite region for 102, 103, 104, and 105 years.

Note that Sn¹²⁶ and Tc⁹⁹ are the nuclides most affected by sorption on bentonite, with their retardation in migration observed at all times of simulation.

The nuclides with long life Pd¹⁰⁷, Zr⁹³, and Se⁷⁹ were little affected by bentonite. However, there is a decrease in the nuclides concentration of short life nuclides, such as Ni⁶³ (with a half-life equal to 9.6 · 10¹ years), Sm¹⁵¹ (with a half-life equal to 9.0 · 10¹ years), and Sr⁹⁰ (with a half-life of 2.9 · 10¹ years). The amount of these nuclides becomes totally depleted in a time span of 10⁵ years.

Figure 7 shows the concentration profile for the relevant nuclides in high level radioactive waste of 4𝑛 chain: Pu²⁴⁰→U²³⁶, in bentonite extruded region.

(a)

(b)

(c)

(d)

Figure 7 

Concentration profiles in extruded bentonite region, for nuclides of 4𝑛 chain (Pu240→U236) for 102, 103, 104, and 105 years.

We can see a small retention for the U²³⁶ at the first two times span simulated, 10² and 10³ years; after that, 10⁴ and 10⁵ years, the bentonite is saturated in this nuclide.

However, the Pu²⁴⁰ presents a significant retardation in migration at all simulated times, showing the importance of the clay for this nuclide retention, noting that this element is the one with greater radiological impact.

Figure 8 shows the concentration profiles for the relevant nuclides in high level radioactive waste of 4𝑛 + 1 chain: Cm²⁴⁵→Am²⁴¹→Np²³⁷→U²³³→Th²²⁹, in the bentonite extruded region. This series presents the radioactive family with greater radiological impact (Cm²⁴⁵→Am²⁴¹→Np²³⁷).

(a)

(b)

(c)

(d)

Figure 8 

Concentration profiles in extruded bentonite region, for nuclides of 4𝑛+1 chain (Cm245→Am241→Np237→U233→Th229) for 102, 103, 104, and 105 years.

We can note the growth of Th²²⁹. This nuclide did not exist in the waste when disposed in the repository but was generated by the decay of U²³³. Th²²⁹ presents significant retardation in migration as well as Am²⁴¹ and Cm²⁴⁵. Both nuclides, Am²⁴¹ and Cm²⁴⁵, have similar input parameters but differ in half-life (Cm²⁴⁵ half-life equal to 8.5 · 10³ years and Am²⁴¹ has a half-life of 4.3 · 10² years), thereby providing clearance profiles over time, while achieving a balance at the end of 100,000 years.

For nuclide U²³³, the retention is not significant. Np²³⁷ exhibits similar behavior, showing complete saturation in the bentonite region at the end of simulation time.

Figure 9 shows the concentration profile for the relevant nuclides in high level radioactive waste for 4𝑛 + 2 chain: Cm²⁴⁶→Pu²⁴²→U²³⁸→U²³⁴→Th²³⁰→Ra²²⁶, in the extruded bentonite region.

(a)

(b)

(c)

(d)

Figure 9 

Concentration profiles in extruded bentonite region, for nuclides of 4𝑛+2 chain (Cm246→Pu242→U238→U234→Th230→Ra226) for 102, 103, 104, and 105 years.

Here we can note the growth of Ra²²⁶ nuclide that did not exist in the waste when disposed in the repository but was generated by decay of Th²³⁰, and although it is the nuclide with the smallest half-life (1600 years) in this chain, we note that it has no significant decrease in time, since the Ra²²⁶ is generated continuously by its precursor and presents a slight retention in the bentonite extruded region.

For the Cm²⁴⁶, with half-life of 4700 years, which is the first chain element, we can observe its decrease in time. The nuclides Cm²⁴⁶, Pu²⁴², and Th²³⁰ are greatly affected by bentonite, showing significant retention.

Figure 10 shows the concentration profile for the relevant nuclides in high level radioactive waste for 4𝑛+3 chain: Am²⁴³→Pu²³⁹→U²³⁵→Pa²³¹, in the extruded bentonite region.

(a)

(b)

(c)

(d)

Figure 10 

Concentration profiles in extruded bentonite region, for nuclides of 4𝑛+3 chain (Am243→Pu239→U235→Pa231) for 102, 103, 104, and 105 years.

We can observe a significant delay in the nuclides Am²⁴³ and Pu²³⁹ for the whole time of simulation. We can also note the mass gain of Pa²³¹ by decay of its precursor U²³⁵ and the slight retardation of this nuclide (U²³⁵). We can note their saturation in the region of extruded bentonite at the end of simulation time.

In this chain, the Am²⁴³ appears as the radionuclide of greatest importance not only because of its half-life (4.7 · 10³ years) but also for being the precursor of Pu²³⁹.

Figure 11 shows an estimate of the committed effective dose to an adult in Sv/a, obtained from ingestion of well water during one year containing relevant radionuclides of spent fuel that were released by the repository.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11 

Estimated values of the committed effective dose to an adult in Sv/a, obtained from ingestion of well water containing relevant radionuclides of spent fuel that were released by the repository.

In Figure 11(a) we present the first group of important nuclides: C¹⁴, Cl³⁶, Cs¹³⁵, Cs¹³⁷, I¹²⁹, Mo⁹³, Nb⁹⁴, and Ni⁵⁹; in Figure 11(b) the second group: Ni⁶³, Pd¹⁰⁷, Sm¹⁵¹, Se⁷⁹, Sn¹²⁶, Sr⁹⁰, Tc⁹⁹; and Zr⁹³, and in Figures 11(c), 11(d), 11(e), and 11(f) the relevant nuclides of 4𝑛 (Pu²⁴⁰→U²³⁶), 4𝑛 + 1 (Cm²⁴⁵→Am²⁴¹→Np²³⁷→U²³³→Th²²⁹), 4𝑛+2 (Cm²⁴⁶→Pu²⁴²→U²³⁸→U²³⁴→Th²³⁰→Ra²²⁶), and 4𝑛 + 3 (Am²⁴³→Pu²³⁹→U²³⁵→Pa²³¹) chains, respectively. All graphs are on the same scale and show the upper restriction dose limit of 10⁻⁴ Sv/a.

We can see in Figure 12 the release rate to the environment, in Bq/a, of relevant radionuclides of spent fuel. In Figure 12(a) we show the first group of important nuclides: C¹⁴, Cl³⁶, Cs¹³⁵, Cs¹³⁷, I¹²⁹, Mo⁹³, Nb⁹⁴, and Ni⁵⁹; in Figure 12(b) the second group: Ni⁶³, Pd¹⁰⁷, Sm¹⁵¹, Se⁷⁹, Sn¹²⁶, Sr⁹⁰, Tc⁹⁹, and Zr⁹³; and in Figures 12(c), 12(d), 12(e), and 12(f) the relevant nuclides of 4𝑛, 4𝑛+1, 4𝑛+2, and 4𝑛+3 chains, respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12 

Release rates to the environment in Bq/a of spent fuel relevant nuclides.

In these figures we can note that the largest peak release rate comes from inventory nuclides with higher activities and shorter half-life, for example, Cs¹³⁷, Ni⁶³, and Sr⁹⁰.

The Pa²³¹ (in Figure 12(f)) and Ra²²⁶ (in Figure 12(e)) are nuclides that do not have significant amounts initially in the waste but they arise from their precursors decay and their high rates of release represent an important contribution to dose.

Figure 13 shows the maximum dose rate to which a member of the public may be exposed by ingestion of water containing radionuclides considered in this study that migrated from repository to environment. The dose is obtained estimating the radionuclide flow entering the environment, that is, the output of lineament region multiplied by the conversion factor dose (Tables 4 and 5 in the appendix). Dose conversion factor considers that the annual releases from the repository to the biosphere are diluted in 100,000 m³ of water and that an individual drinks 600 liters of contaminated water annually.

Table 4 

Radionuclide chain decay data.

Table 5 

Radionuclide data.

Figure 13 

Maximum dose rate to an exposed member of public by ingestion of water with radionuclides considered in this study in Sv/a.

6. Conclusion

Performance assessment is the use of mathematical models to simulate the long-term behavior of engineering and geological barriers of a nuclear waste repository. Thus, analysis of radionuclide migration in fractured porous media is an important part of the safety assessment of a deep host rock for high level radioactive waste.

In this work, the numerical solution to a simple geometry, including decay chain, was developed to study the migration of radionuclides in the far field region of a hypothetical repository. To validate the numerical modeling of the decay chain proposed by finite differences, we used an analytical solution in the literature, since analytical solutions are useful for checking the results of numerical codes.

This study is a deterministic analysis, where the assumptions are purposely conservative, which means that they must ensure that the results, with a high degree of certainty, overestimate the radioactive releases. The analysis covers the releases of relevant spent fuel radionuclides and their arrival to the biosphere. The main quantitative results are expressed as indicative dose rates from consumption of well water and as activity releases from borehole to the geosphere and from the geosphere to the biosphere.

The radioactive decay tends to reduce the concentrations of nuclides of short half-life more efficiently. Thus, there is clearly a decrease in the concentration of nuclides of lower half-life, such as C¹⁴ (with a half-life equal to 5.7 · 10³ years), Mo⁹³ (with a half-life equal to 4.0 · 10³ years), Cs¹³⁷ (with a half-life equal to 3.0 · 10¹ years), Ni⁶³ (with a half-life equal to 9.6 · 10¹ years), Sm¹⁵¹ (with a half-life equal to 9.0 · 10¹ years), and Sr⁹⁰ (with half-life equal to 2.9 · 10¹ years). These last four nuclides had their concentration exhausted at time of 10⁵ years. With respect to retention in the bentonite region, we note that the Sn¹²⁶ and Tc⁹⁹ are the nuclides most affected by sorption on bentonite, with their retention observed at all times of simulation.

In the decay chains considered, we can see the growth of Th²²⁹, a nuclide that was not in the waste when disposed in the borehole but that was generated by decay of U²³³ as well as the growth of Ra²²⁶ which was generated by decay of Th²³⁰, and albeit it was the nuclide of lower half-life in 4𝑛 + 1 chain (1600 years), we note that it does not decrease significantly in time, because its generation compensates its loss by decay.

The nuclides Cm²⁴⁶, Pu²⁴², Th²³⁰, Am²⁴³, and Pu²³⁹ have a significant retention in bentonite region at all simulation times. This is also the case of Pu²⁴⁰, which presents significant migration retardation at all simulated times. This shows the importance of bentonite clay for this later nuclide retention, noting that it is one nuclide that has greater radiological impact, because the toxicity of plutonium is among the largest of any known element.

We can also note the mass gain of Pa²³¹ by decay of its precursor U²³⁵ and the slight delay of U²³⁵ and its saturation in bentonite extruded region at the end of simulation time. In this chain, the Am²⁴³ appears as the most important radionuclide, not only because of its half-life but also because it is the precursor of Pu²³⁹.

We can see that the highest peak release rate comes from nuclides with higher activity and short half-life, for example, Cs¹³⁷, Ni⁶³, and Sr⁹⁰. The Pa²³¹ and Ra²²⁶ are nuclides with minor concentrations in the waste, but they arise from the radioactive decay and have high rates of release and dose. Thus, radionuclides generated by decay are among those causing higher impact on the individual by drinking water containing radionuclides released from the repository. In this simulation, the Sr⁹⁰ is the nuclide with greater contribution to dose, but its value is below the permissible limit of 0.1 mSv/year. Thus, this work guides further studies, indicating that special attention should be given to this nuclide since it is a conservative study.

Finally, the present model constitutes an important guidance for the complete and detailed safety assessment of a high level radioactive waste repository. The results for the hypothetical case analyzed in this paper show that final spent nuclear fuel repository, based on KBS-3 model, can achieve an adequate level of protection, with a minor contribution to the public annual radiation dose.

Appendix

Data for Simulation

Tables 2, 3, 4, and 5 present the physical parameters for fracture, extruded bentonite region, porous matrix, and lineament.

Disclosure

The authors report that reference to SCALE-5 software package was only to inform the reader that data from a referenced paper, concerning composition and radioactive characteristics of spent fuel that used this package, were used in this work.