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

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.


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 (HLWhigh 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.
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.

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 onedimensional 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.
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 where   is dispersivity in the directions of the fracture axis (L), V is water velocity in the fracture (L/T), and  * is coefficient of molecular diffusion in water (L 2 /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.

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 where  , (, ) is the concentration of radionuclide  in fracture water (M/L 3 ),  is the spatial coordinate along the fracture (L),  is the time (T), 2 is the fracture aperture (L), V is the groundwater velocity in the fracture (L/T), and  , is the longitudinal dispersion coefficient of solute  (L 2 /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 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 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: where  , (, , ) is the concentration of radionuclide  in the porous matrix (M/L 3 ) 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 where   is the density of the rock matrix, (M/L 3 ),  is the porosity of the rock, and  , is the distribution coefficient in the rock matrix, (L 3 /M).
The effective molecular diffusive coefficient  , (L 2 /T) is less than the free-solution value,   * , because of the tortuous diffusion pathway.For an anisotropic medium, it can be represented by 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 2 T), and it is expressed in the form of Fick's first law as 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 The boundary conditions in the fracture are given by where   0 is the concentration of th radionuclide that is deposited instantaneously at the coordinate origin of the fracture (M/L 3 ).
Equation (11) expresses the coupling between the concentration in the fracture and the concentration in the porous matrix.

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 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 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.

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.
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 threemember 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.

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.
The diffusive transport of the radionuclide chain in the bentonite extruded region is governed by the following system of coupled diffusion equations: where  , (, ) is the concentration of radionuclide  in the bentonite (M/L 3 ) 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 where  is the density of the bentonite, (M/L 3 ),  is the porosity of the bentonite, and  , is the distribution coefficient in the bentonite, (L 3 /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.
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].
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 ), 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 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.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 14
Figure 6 shows the concentration profile for the second group of relevant nuclides in high level radioactive waste: Ni 63 , Pd 107 , Sm 151 , Se 79 , Sn 126 , Sr 90 , Tc 99 , and Zr 93 in the extruded bentonite region.
Note that Sn 126 and Tc 99 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 107 , Zr 93 , and Se 79 were little affected by bentonite.However, there is a decrease in the nuclides concentration of short life nuclides, such as Ni 63 (with a half-life equal to 9.6 ⋅ 10 1 years), Sm 151 (with a half-life equal to 9.0 ⋅ 10 1 years), and Sr 90 (with a half-life of 2.9 ⋅ 10 1 years).The amount of these nuclides becomes totally depleted in a time span of 10 5 years.
Figure 7 shows the concentration profile for the relevant nuclides in high level radioactive waste of 4 chain: Pu 240 → U 236 , in bentonite extruded region.
We can see a small retention for the U 236 at the first two times span simulated, 10 2 and 10 3 years; after that, 10 4 and 10 5 years, the bentonite is saturated in this nuclide.
However, the Pu 240 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 245 → Am 241 → Np 237 → U 233 → Th 229 , in the bentonite extruded region.This series presents the radioactive family with greater radiological impact (Cm 245 → Am 241 → Np 237 ).We can note the growth of Th 229 .This nuclide did not exist in the waste when disposed in the repository but was generated by the decay of U 233 .Th 229 presents significant retardation in migration as well as Am 241 and Cm 245 .Both nuclides, Am 241 and Cm 245 , have similar input parameters but differ in half-life (Cm 245 half-life equal to 8.5 ⋅ 10 3 years and Am 241 has a half-life of 4.3 ⋅ 10 2 years), thereby providing clearance profiles over time, while achieving a balance at the end of 100,000 years.

Science and Technology of Nuclear Installations
For nuclide U 233 , the retention is not significant.Np 237 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 246 → Pu 242 → U 238 → U 234 → Th 230 → Ra 226 , in the extruded bentonite region.
Here we can note the growth of Ra 226 nuclide that did not exist in the waste when disposed in the repository but was generated by decay of Th 230 , 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 226 is generated continuously by its precursor and presents a slight retention in the bentonite extruded region.
For the Cm 246 , with half-life of 4700 years, which is the first chain element, we can observe its decrease in time.The nuclides Cm 246 , Pu 242 , and Th 230 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 We can observe a significant delay in the nuclides Am 243 and Pu 239 for the whole time of simulation.We can also note the mass gain of Pa 231 by decay of its precursor U 235 and the slight retardation of this nuclide (U 235 ).We can note their saturation in the region of extruded bentonite at the end of simulation time.
In this chain, the Am 243 appears as the radionuclide of greatest importance not only because of its half-life (4.7 ⋅ 10 3 years) but also for being the precursor of Pu 239 .
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.
In Figure 11( We can see in Figure 12 the release rate to the environment, in Bq/a, of relevant radionuclides of spent fuel.In Figure 12( 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 137 , Ni 63 , and Sr 90 . The Pa 231 (in Figure 12(f)) and Ra 226 (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 3 of water and that an individual drinks 600 liters of contaminated water annually.

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 14 (with a half-life equal to 5.7 ⋅ 10 3 years), Mo 93 (with a half-life equal to 4.0 ⋅ 10 3 years), Cs 137 (with a half-life equal to 3.0 ⋅ 10 1 years), Ni 63 (with a half-life equal to 9.6 ⋅ 10 1 years), Sm 151 (with a half-life equal to 9.0 ⋅ 10 1 years), and Sr 90 (with half-life equal to 2.9 ⋅ 10 1 years).These last four nuclides had their concentration exhausted at time of 10 5 years.With respect to retention in the bentonite region, we note that the Sn 126 and Tc 99 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 229 , a nuclide that was not in the waste when disposed in the borehole but that was generated by decay of U 233 as well as the growth of Ra 226 which was generated by decay of Th 230 , and albeit it was the nuclide of lower half-life in 4 + 1 chain (1600 years), we note that it does not decrease significantly

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

Figure 4 :
Figure 4: Modify geometry of proposed model.

Figure 5
shows the concentration profile for the first group of relevant nuclides in high level radioactive waste: C 14 , Cl 36 , Cs 135 , Cs 137 , I 129 , Mo 93 , Nb 94 , and Ni 59 .

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

Table 1 :
Physical parameters used in simulation.

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 ().