Numerical modeling of radionuclide migration through a borehole disposal site

The migration of radionuclides from a borehole repository located about 20 km from the Akwapim fault line which lies in an area of high seismicity was analyzed for some selected radionuclides. In the event of a seismic activity, fractures and faults could be rejuvenated or initiated resulting in container failure leading to the release of radionuclides. A numerical model was solved using a two-dimensional finite element code (Comsol Multiphysics) by taking into account the effect of heterogeneities. Results showed that, the fractured medium created preferential pathways indicating that, fault zones generated potential paths for released radionuclides from a radioactive waste repository. The results obtained showed that variations in hydraulic conductivity as a result of the heterogeneity considered within the domain significantly affected the direction of flow.


Introduction
The site for the Borehole Disposal Facility is within the Accra Plains and is situated on the boundary between the Dahomeyan System and the Togo Series (Figure 1), both of which are of Precambrian age. The Dahomeyan is the major bedrock formation underlying the site and consists of gneiss and schist while the Togo Series are predominantly composed of quartzite and phyllite. The Akwapim hills and the Togo range mark the line of a major active fault zone that runs North-East and South-West to the coastal boundary fault (Muff & Efa 2006).
The whole of the site is covered by loose unconsolidated and weathered material that may reflect the presence of troughs formed by down faulted blocks which indicates the existence of seismic activity in the geologic past and it probably results from movements along the Akwapim fault line (Junner & Bates 1995).
The only major river near the Borehole Disposal Facility (BDF) site is River Onyasia, located about 1.3 km from the proposed facility and drains southwards through Achimota village to Accra. The Onyasia River has a depth of 0.6 m, a width of 6.8 m and a measured velocity of 0.8 m/s (2.5 × 10 7 m/y). The broad valley of the Onyasia River flanks the site on its eastern margin with swampy conditions generally found north-east of the site. Surface run-off in this area is very low, however, after heavy storms there is flow of water over the horizon below the top-soil (Darko et al. 1995).
Seismic surveys conducted in the area mapped out two weak lines suspected to have been caused as a result of faults or fractures. However, these have not been used in numerical modelling to follow and predict the migration of radionuclides at the site in case a seismic accident occurs (Essel et al. 2011).

Model definition
The hydrogeologic setting for groundwater flow at steady state shall be described by the conceptual model shown in Figure 2. Groundwater moves from a higher topographic surface to a lower topographic surface. The water table is a free surface across which there is vertical recharge, denoted by R. The base of the aquifer is considered impermeable.

Scenario development
In order to predict radionuclide release, the engineering barriers are assumed to fail in the event of a seismic activity. The system is thus, simplified into a two-dimensional conceptual model as shown in Figure 2.
The radionuclides are initially confined in the canister until a seismic activity occurs, leading to a crack of the barrier system such that the radionuclide inventory is released into the groundwater which is the major transport medium. It is assumed in the calculations performed that radionuclides start being released from the canister 30 years after closure of the repository.

Numerical illustrations
A two-dimensional numerical model was developed using Comsol Multiphysics (ver.3.4) similar to the proposed model in Figure 2. The lithology of the system was characterized by a porosity of 0.35, a transverse dispersivity of 0.005 m, a longitudinal dispersivity of 0.5 m and a hydraulic conductivity ranging from 10 −15 to 10 −5 m/s. The time-dependent solute transport equation as described by equation (1): Here, the dispersion tensor D L defines solute spreading by mechanical mixing and molecular diffusion. Equations for the tensor entries are: Where D ii , D jj are the principal components of the dispersion tensor based on the Darcy's velocity, D ij and D ji are the cross terms of the dispersion tensor, the  subscript "L" denotes longitudinal dispersivity, "T" the transverse dispersivity. v is the magnitude of the Darcy's velocity vector, D m represents the effective molecular diffusion in a saturated porous media and i, j are the spatial coordinates.

Governing equations
Fluid flow: assumptions ▪ On the scale simulated, the fracture system behaves as an equivalent porous medium ▪ The groundwater flow is assumed to be homogenous and subject to recharge ▪ Additionally, groundwater flows under steady state conditions. This means that, the velocity of flow is considered not to change with time since groundwater flow is naturally a slow process.

Fluid flow: domain equations and boundary conditions
Steady groundwater flow is expressed with a conservation equation formulated with Darcy's law and expressed as: Where K is hydraulic conductivity (m/s); v i is the Darcy velocity (m/s); x i is the spatial distance (m) and H is the hydraulic head (m); Hydraulic head is a function of pressure and gravitational potential and it is defined: Where, h p is pressure head (m); and y is elevation (m). The equations for groundwater flow and solute transport are linked by: Where R is the recharge rate (m/s). Figure 3 shows the boundary conditions for the groundwater flow problem and stated by equations (6-8).
A zero flux Neumann condition represents the symmetry boundary at x = 0 m and the impermeable boundary at y = 0 m as follows: Hydraulic head is specified at x = 14 m with a Dirichlet condition: The annual average precipitation rate is recorded as 800 mm/y equivalent to 2.537 × 10 −8 m/s. The water table receives about 10% of this value. From a topographic height of y = 500, an inflow boundary condition that satisfies Neumann's condition is assumed and is written in the form: Solute transport: assumptions Transport of radionuclides is based on the following assumptions: ▪ Radioactive decay is the only reaction considered in the model. It is assumed to occur throughout the model in the liquid phase; ▪ For the purposes of this work, no gaseous release is considered; ▪ Transport of radionuclides is assumed to occur in the saturated zone;

Solute transport: domain equations and boundary conditions
For the transport equation, the initial and boundary conditions for solute transport are shown in Figure 4. In this work, it is also assumed that the contaminant source is continuous over a period of thirty years. Dirichlet conditions are used at the water table, where C(x, h, t) = 0, except for the segment 4m ≤ x ≤ 8m in which: Where C 0 is the relative concentration through year 30(t 0 ). The Dirichlet condition at the left boundary is: The Neumann condition defining the zero flow boundaries is: Finally, the initial condition indicating that the subsurface on the site is free of contaminants at the beginning of the simulation is:

Results and discussions
Numerical simulations were considered for Cobalt-60 (short half-life), Cesium-137 (medium half-life) and Americium-241 (long half-life) migrating from the repository through the subsurface environment for a 2dimensional cross-section along the streamlines of groundwater flow. These considerations were made because of the high amount of activity contained in the waste packages and the variability of their half-life as shown in Table 1.

Computer simulations
The models describe the steady-state fluid flow and follows up with a transient solute transport simulation.
Two partial differential equations (PDE) were solved for and these were assigned in separate mathematics interfaces in Comsol Multiphysics (version 3.4). The first partial differential equation (PDE) is stationary and it finds a solution to the Darcy velocity while the second partial differential equation (PDE) is time-dependent and finds a solution to the solute transport equation.
To test the effect of heterogeneity on solute behavior, a two-dimensional numerical transport model was created to investigate solute transport under two hydrogeologic conditions: (1) homogeneous hydraulic conductivity in porous subsurface medium and (2) heterogeneous hydraulic conductivity in a fractured medium. Table 2 describes the physical parameters used in the saturated homogenous porous medium for the numerical model of flow and transport. Except for increasing the hydraulic conductivity within the fracture, all other parameters are the same as in the homogenous model.

Model simulation results
Throughout the models, the amount of contaminant is shown by the colour bar. The activity concentration degree is indicated by the various colours, with red indicating an intense concentration.

Evolution of 60 Co
Radionuclide movement (streamline plot) and concentration estimates (surface plots) of 60 Co simulated at various times are illustrated in Figures 5a to 6. 60 Co has a half-life of 5.27 years and an initial inventory of 2.22 × 10 8 Bq/m 3 . From the simulations carried out, it was evident that, for a homogenous porous medium in a low conductivity zone, 60 Co contaminant source diffuses with the groundwater which has a flow velocity of 3.834 × 10 −8 m/s. It was observed that, the radionuclide migrates in the same direction as the groundwater flow velocity while at the same time undergoing radioactive decay. However, when the hydraulic conductivity was increased ( Figure 6) to account for heterogeneity within the system; it was observed that 60 Co radionuclide plume does not rise towards the surface. Instead, the plume is diluted into deeper groundwater flow systems as it decayed away. This may be attributed to the short radioactive half-life as well as the high conductivity zone exhibited by the fractures. Figure 7 shows the evolution of 60 Co contaminant concentration as a function of time. Here, the diffusion process takes place much longer which is accounted for by the widening of the peak hence the radionuclide plume concentration appears broad in both models. For a homogenous porous medium, the concentration increase steadily until it reaches a maximum and then decreases exponentially. The high value of the initial activity concentration together with the diffusion process taking place, explains why the observed radionuclide plume extends over a very long range of time.

Evolution of 137 Cs
137 Cesium has a half-life of 30 years and an initial activity concentration of 5.66 × 10 12 Bq/m 3 . Evolution of 137 Cs contaminant source are shown in Figure 8a and b at various times (t = 250 years to t = 1000 years). In this case, 137 Cs is migrating in a low conductivity medium causing the radionuclide to be deposited onto sediments which sinks steeply downwards through the groundwater flow system. The contaminant source diffuses with the flowing groundwater which has a measured flow velocity of 3.834 × 10 −8 m/s. The flow velocity readily  moves with 137 Cs radionuclide in water causing contamination of groundwater. Figure 9a and b show the influence on radionuclide migration when the hydraulic conductivity is increased. After 40 years (Figure 9a) a fraction of its initial inventory would have undergone radioactive decay. Therefore in the presence of fractures, the increase in the hydraulic conductivity would provide preferential pathways for the flow velocity leading to an increase in mobility for the migrating radionuclide. The flow in this case causes the radionuclide plume to rise towards the surface which may finally discharge into wells, rivers and lakes endangering the environment, population and biota. Figure 9b shows that, radioactive decay, diffusive flux and the high hydraulic conductivity continues to have an effect in the direction and migration of the radionuclide source. These parameters tend to reduce the plume concentration as well as cause the radionuclide source to move away from the surface towards deeper groundwater flow systems contaminating groundwater resources. Figure 10 shows the evolution of 137 Cs contaminant concentration as a function of time for the homogenous and heterogeneous model. For 137 Cs, the radionuclide plume concentration appears broad suggesting that, the diffusion process takes place much longer which is accounted for by the widening of the peak. Also, the concentration increase steadily until it reaches a maximum and then decreases exponentially. The variation in the curves is a result of the presence of fractures which creates preferential pathways accounting for the high mobility of the radionuclide in the heterogeneous model.   Yeboah et al. SpringerPlus 2014, 3:155 Page 7 of 11 http://www.springerplus.com/content/3/1/155

Evolution of 241 Am
The simulation was implemented for 241 Am contaminant source with an initial activity concentration of 3.5 × 10 7 Bq/m 3 , a half-life of 432 years and a flow velocity of 3.534 × 10 −8 m/s. Figure 11a and b shows evolution of 241 Am radionuclide activity concentration at various times in a homogenous porous medium having a low hydraulic conductivity of 3.17 × 10 −5 m/s. From the simulations, 241 Am was observed to sink steeply downward into the groundwater flow system after 200 years. At 250 years (Figure 11a), it continued to steeply decline into the groundwater flow system until about 1000 years ( Figure 11b). Given that the value of the flow rate is low due to the low value of the hydraulic conductivity, transport of the radionuclide was mainly dominated by diffusion.     (Figure 12a and b), the radionuclide source will not have undergone any appreciable decay process because of the long half-life (t 1/2 = 432 years). Hence, majority of the radionuclide would appear to move towards the surface due to the high conductivity layer. In addition, lateral gradients can cause rapid rise of water-table leading to the possibility that the movement of 241 Am attributable to diffusion is faster than the rate at which groundwater is moving. The reduction in the initial activity concentration of 241 Am from 3.5 × 10 7 Bq/m 3 to 3.256 × 10 4 Bq/m 3 can be attributed to mass loss by the diffusive flux. Figure 13 shows the evolution of 241 Am activity concentration as a function of time for a homogenous and heterogeneous model. The nature of the peak in both models suggests that, advection influences the distribution moderately with diffusion being the dominant transport mechanism. It is noticeable that, the peak concentration does not reach the initial activity value because of mass loss by radioactive decay. Similarly, the concentration reaches a maximum and then decreases exponentially. 241 Am activity concentration extends over a very long range of time due to the high value of the initial activity, the long half-life and the diffusion process taking place.

Conclusions
The migration of three radionuclides namely, 60 Co, 137 Cs and 241 Am have been simulated using a twodimensional finite element numerical model code (Comsol Multiphysics).
Neglecting heterogeneity, simulated results showed that, all three radionuclides ( 60 Co, 137 Cs, 241 Am) within the low conductivity medium sunk steeply downward into the groundwater flow system by diffusing into the flowing groundwater. This caused the flow velocity to move readily with the radionuclide source causing contamination of groundwater resources.
In the presence of fractures, preferential pathways were created which gave rise to a rapid increase of the water-table and this caused the flow velocity to sweep the radionuclides with medium ( 137 Cs) to long ( 241 Am) half-life toward the surface endangering human population, the environment and biota.
For 60 Co, the plume was not noticeably seen at the surface even in the presence of high hydraulic conductivity but was rather diluted into deeper groundwater flow systems as it decayed away. This was attributed to the short radioactive half-life.
The results obtained showed contamination to be more sensitive to variations in hydraulic conductivity as a result of the heterogeneity considered within the domain. However, impact on groundwater was still inevitable.

Recommendation
It is recommended that, proper structural geological mapping including the use of stereograms should be made to be able to determine the fractures before radioactive waste is disposed of in an area.