Integral transform analysis of radionuclide transport in variably saturated media using a physical non-equilibrium model: application to solid waste leaching at a uranium mining installation

Abstract

Abstract:

The Generalized Integral Transform Technique (GITT) was employed to simulate transient one-dimensional flow in variably saturated porous media, as well as radioactive waste transport within different layers (a solid waste pile, nearby soil, and a granular aquifer) towards the edge of a uranium mining installation under institutional control. Computational codes, written using the Mathematica software system, were implemented and tested for solution of the coupled advection-dispersion equations for an arbitrary number of daughter products within a radioactive chain migrating in saturated and unsaturated soil layers. The computer simulations were verified in great detail against results obtained using the built-in routine NDSolve of the Mathematica platform and the HYDRUS-1D software system. The present work reports the main results for ²³⁸U chain radionuclide transport using data extracted from a safety assessment of solid waste repositories at a uranium mining and milling installation in Caetité, state of Bahia, Brazil, operated by INB (Indústrias Nucleares do Brasil). Concentration evolutions of the various radionuclides obtained with the simulations were analyzed for five different cases to explore variations in the infiltration and recharge rates, the effect of assuming physical equilibrium or non-equilibrium transport conditions, and of different initial concentrations of some of the radionuclides.

Key words
uranium mining; radioactive waste leaching; unsaturated porous media; physical non-equilibrium transport; hybrid methods; Generalized Integral Transform Technique (GITT)

INTRODUCTION

Uranium mining and milling is an essential activity within the nuclear fuel cycle. Very few countries at present are self-sufficient in terms of uranium production and having the overall technological expertise for uranium enrichment. In 2017, Kazakhstan had the largest share of uranium production from mines (39% of the world supply), followed by Canada (22%) and Australia (10%) (WNA 2019). Conventional mines generally contain milling facilities where the ore is crushed, placed in leaching platforms, and then leached with sulfuric acid to dissolve the uranium oxides. The uranium content of ore is often only between 0.1% and 0.2%, which means that large amounts of ore have to be mined to produce a justifiable amount of uranium. Uranium mill tailings and leachates are normally dumped as a sludge in special ponds or as solid waste piles. The amount of waste produced is then nearly the same as that of the original ore before milling (about 99.9%). Except for the parcel of uranium extracted through leaching, the waste contains most of the original constituents of the ore. Since long-lived decay products such as thorium-230 and radium-226 are not removed, the waste may contain around 85% of the initial radioactivity of the ore, including small amounts of uranium that could not be extracted (Diehl 2011). The solid waste also still contains some of the leaching agents, despite washing operations at the leaching platform. Even if not, sulfuric acid can be formed within the waste piles if they contain certain minerals such as pyrite (FeS₂), which then leads to continuous leaching of the ore. Regulatory bodies for these reasons are always concerned about seepage from solid waste piles in view of the risk of contamination of adjacent soil layers and beyond (IAEA 2003).

The present work is part of a long-term study of the radiological impact of leachate and sterile ore piles at the Uranium Concentration Unit (URA) of INB (Indústrias Nucleares do Brasil) near Caetité in the state of Bahia, Brazil (Orlande et al. 2008). The URA/Caetité is an open pit mine, which at the time of this study was producing approximately 300 tons of uranium per year. Radioactive waste migration in soils near the leached and sterile ore piles and their environmental impact were analyzed based on input data obtained from documentation available at INB, from additional experiments and characterizations of the soils and tailings involved, and complemented when necessary with information from the literature.

In this paper we initially synthesize the site characterization and the adopted physical model. We then describe the methodology used for estimating the source term, as well as for the simulations of water flow and radionuclide transport in the leached ore pile and nearby soil layers. The hybrid numerical-analytical methodology known as the Generalized Integral Transform Technique (GITT) (Cotta 1990, 1993, 1994, 1998, Cotta & Mikhailov 1997, 2006, Cotta et al. 2016, 2017, 2018) was adopted for solution of the partial differential equations governing water flow and radionuclide transport. This hybrid method, which originates from classical integral transform methods (Mikhailov & Ozisik 1984), has been used extensively within thermal sciences and engineering, but more recently also in environmental engineering applications (Cotta et al. 1998, 2003, Leal & Ruperti 2000, 2001, Liu et al. 2000, Rocha & Cruz 2001, Heilbron Filho et al. 2002, Barros et al. 2006, Skaggs et al. 2007, Barros & Cotta 2007, Almeida et al. 2008, Perez Guerrero et al. 2009, 2010, 2012, Perez Guerrero & Skaggs 2010, Naveira-Cotta et al. 2013, Rubol et al. 2016). The GITT methodology was previously applied to studies of the migration of radioactive contaminants from liquid waste ponds at the same uranium mining facility of INB (Orlande et al. 2006, Naveira-Cotta et al. 2013). The solution procedure was selected because of its merits in analyzing long-term contaminant transport problems, which for certain applications may run for tens of thousands of years (such as the present one), and because it is one of the recommended methods for radioactive waste safety assessment by the International Atomic Energy Agency (IAEA 2004, Heilbron Filho et al. 2004).

We also present the verification of the integral transformation solutions as implemented here for fluid flow and contaminant transport in variably saturated physical non-equilibrium (dual-porosity type) media. The solutions were verified against results obtained using the HYDRUS-1D software system (Simunek et al. 2005, 2016), and other studies. HYDRUS-1D is a widely known and freely distributed finite-element model for simulating the one-dimensional movement of water, heat, and multiple solutes in variably saturated media. The code solves the Richards equation for saturated–unsaturated water flow and advection–dispersion type equations for heat and solute transport.

The developed GITT radionuclide migration code is employed subsequently to analyze a proposed conservative physical model that considers percolation from the leachate ore stack, dilution in a pre-basin (a natural reservoir formed by the confluence of a second bottom drain and the foot channel of the pile), vertical infiltration into soil below the pre-basin, leachate mixing with the underlying granular aquifer, and then lateral migration along the aquifer to the nearest stream (known as “Riacho das Vacas”, or in English, the “Cows´ Creek”). Using as much as possible local experimental data, results were obtained for five different cases to explore variations in the infiltration and recharge rates, by switching from physical equilibrium to non-equilibrium models, and for different initial conditions for the radionuclide concentrations in the leached ore pile. The simulation tool also permitted calculations of doses using well-established safety analysis scenarios (Pontedeiro et al. 2007, Naveira-Cotta et al. 2013).

Interested readers are referred to (Orlande et al. 2008) for details about the physical and geochemical properties of the soils and waste as obtained from batch and column experiments for uranium dissolution, including validations against theoretical geochemistry models and measurements obtained at the URA itself. An extensive experimental campaign of measurements was carried out of the solid waste pile, the pre-basin soil, and the granular aquifer, including analyses using inverse procedures for identification of unsaturated soil hydraulic and contaminant transport properties (Orlande et al. 2008, Moreira et al. 2016).

SITE DESCRIPTION AND PHYSICAL MODEL

Site Description

The deposit consisted mostly of three separate materials: the sterile soil (estimated at 1.08x10⁶ tons) originating from the top cover of the mine (approximately five meters deep), the leachate ore (on the order of 2.11x10⁶ tons, with granulometry below 10 mm) originating from piles assembled on the leaching platform, and sterile rocks (estimated at 11.13 x 10⁶ tons). Figure 1a offers a detailed view of the executed deposit in 2008, where the different modules have been identified: the sterile soil (blue pointer), the leachate ore (yellow pointer) and the sterile rocks (red pointer) modules. Figure 1b provides a view of the formation of each of the three modules while in execution, with the sterile soil on the left, the leachate ore in the center, and the sterile rock on the right side of the picture.

Figure 1a
General view of the solid waste deposit and location of the executed portion of the sterile soil (blue), leachate ore (yellow) and sterile rock (red) in 2008 (Google Earth).

Figure 1b
View of the constructed solid waste deposit in 2008 (sterile soil on the left, leachate ore in the center, and sterile rock on the right).

Physical Model for Radionuclides Migration

The physical model for radionuclides migration, starting with drainage and dissolution in the leachate ore piles to their transport to the pre-basin and into the granular aquifer, was constructed using several conservative assumptions as described below:

• The third and last washing of ore piles on the leaching platform established the initial conditions for the dissolved radionuclides. This washing with neutral water was interrupted when production control identified a U₃O₈ content of less than 20 ppm and a pH larger than 3 in the effluent aqueous solution. Only then were the piles removed from the leaching platform to the solid waste deposit. This initial condition is considered to be conservative since subsequent leaching, until the pile is removed, still allows some reduction of the uranium contents in solution.

• Drainage of the pile occurred mostly via bottom drains, whose original contours acted as natural paths of the infiltrated water in conjunction with the characteristic drainage of the rocky barren used for its construction. Secondary flows may settle towards the foot of the slope, going directly to the foot-collecting channel of the pile that flows into a pre-basin. However, infiltration into the clay soil of the original sloping terrain on which the leached ore was deposited, as well as along the drains, was conservatively neglected.

• The source term for ²³⁸U chain radionuclide contamination was due to percolation of rainwater, estimated from a hydrologic balance, through the leached ore pile. The combined effects of desorption and dissolution of contaminants during this percolation will cause an increase in radionuclide concentrations up to the foot of the pile, followed by dilution of the effluent with water that infiltrated through the sterile pellets, but also finding its way to the pre-basin.

• From the pre-basin, having an area corresponding to that of the natural basin and a surrounding sedimentation basin, vertical infiltration occurred in transitional soil until encountering a granular aquifer on top of a fissured rock. Transport subsequently followed natural hydraulic gradients towards a downgradient creek (“Riacho das Vacas”). This step in the model was similar to that adopted in an analysis of liquid pond tailings at the same site (Orlande et al. 2006, Naveira-Cotta et al. 2013).

Figure 2 illustrates the main steps in the simulation of radionuclide migration from the leached ore pile of the solid waste deposit towards the granular aquifer. First, a water balance at the top of the ore pile provided a value for the infiltration rate, q_M , which percolates through the leached ore pile, reaching the bottom drain with time-varying radionuclide concentrations. This flow mixes with water percolating through the sterile rock with negligible contamination, which also flows into the pre-basin, thus causing a dilution of the contaminants but increasing the overall vertical infiltration rate (q_V ) into the soil just below the pre-basin. Once this infiltrated water reaches the underlying aquifer, a mixing zone is created in which the horizontal flux increases from its background value of q_H1 to a higher value of q_H2 . The effects of advection, dispersion, retardation, radioactive decay and physical non-equilibrium then govern the transport of the contaminants in the transition soil below the pre-basin and the horizontal aquifer. The governing flow and transport equations are presented in the next section.

Figure 2
Physical model for ²³⁸U chain radionuclide migration at the URA/Caetité solid waste site.

SOLUTION METHODOLOGY AND PROBLEM FORMULATION

Formal Solution Procedure

The Generalized Integral Transform Technique (GITT) (Cotta 1990, 1993, 1994, 1998, Cotta & Mikhailov 1997, 2006, Cotta et al. 2016, 2017, 2018) was adopted for solving the time-dependent partial differential equations governing both water flow and radionuclide transport, and subsequently used in a radiological impact assessment. The hybrid numerical-analytical GITT approach offers an accurate and robust combination of analytical eigenfunction expansions, which eliminates the spatial variables through integral transforms, and a numerical methodology with automatic error control adopted for solution of the resulting transformed ordinary differential system in the time variable.

As illustration of the formal integral transformation process, consider a transient advection-dispersion problem for m coupled potentials (for instance, heads, velocity components, concentrations, or temperatures). These potentials are defined within region V with boundary surface S, with transport considering nonlinear advection and source terms, as follows:

$$\begin{array}{cc}\hfill & w_k(𝐱)\frac{\partial T_k(𝐱,t)}{\partial t}+𝐮(𝐱,t,T_{\ell }).\nabla T_k(𝐱,t)+L_k{T}_k(𝐱,t)=P_k(𝐱,t,T_{\ell }),\hfill \end{array}$$ (1a)

$$\begin{array}{cc}\hfill & 𝐱\in V,t>0,k,\ell =1,2,...,m\hfill \end{array}$$

with initial and boundary conditions given respectively by

$$\begin{array}{cc}\hfill & T_k(𝐱,0)=f_k(𝐱),𝐱\in V\hfill \end{array}$$ (1b)

$$\begin{array}{cc}\hfill & [{\alpha }_k(𝐱)+{\beta }_k(𝐱)K_k(𝐱)\frac{\partial }{\partial 𝐧}]T_k(𝐱,t)={\phi }_k(𝐱,t,T_{\ell }),𝐱\in S,t>0\hfill \end{array}$$ (1c)

where the spatial operator involves both dispersion and linear dissipation:

$$\begin{array}{cc}\hfill & \hfill L_k\equiv -\nabla .K_k(𝐱)\nabla +d_k(𝐱)\end{array}$$ (1d)

and $𝐧$ represents the unitary normal vector pointed towards leaving boundary surface $S$.

The formulation of Eq. (1a) is general since any additional nonlinearities in the transient, dispersion and dissipation operators, as well as in boundary conditions (1c), can be merged into the equation and boundary source terms. Therefore, the linear coefficients that appear in Eq. (1a) can be interpreted as characteristic ones, while the remaining terms in the original formulation are incorporated into the source terms. In any case, the nonlinear advection and source terms make the problem a priori non-transformable. However, the ideas of the Generalized Integral Transform Technique (GITT) can be used to develop a hybrid numerical-analytical solution for this class of problems, with the integral transformation procedure leading to a coupled ordinary differential system to be solved numerically. In case the advection terms are null and the source terms are linear such that $𝐮(𝐱,t,T_{\ell })\equiv 0,P\equiv P(𝐱,t)$, and $\phi \equiv \phi (𝐱,t)$, this example would become a linear class I diffusion problem according to the classification in (Mikhailov & Ozisik 1984). This problem would then lead to an exact analytical solution that can be readily obtained using classical integral transform methods.

Following the solution path previously established for linear diffusion and advection-dispersion problems, the formal solution of the proposed nonlinear problem requires consideration of eigenfunction expansions for the associated potentials. The above-mentioned linear situation that admits exact solutions using classical integral transformation techniques, naturally leads to eigenvalue problems to be preferred also for nonlinear analyses. These, in turn, arise from direct application of separation of variables to the homogeneous and purely diffusive linear version of problem (1). Thus, the set of recommended auxiliary problems is given by:

$$\begin{array}{cc}\hfill & L_k{\psi }_{ki}(𝐱)={\mu }_{ki}^2{w}_k(𝐱){\psi }_{ki}(𝐱),𝐱\in V\hfill \end{array}$$ (2a)

with boundary conditions

$$\begin{array}{cc}\hfill & \hfill [{\alpha }_k(𝐱)+{\beta }_k(𝐱)K_k(𝐱)\frac{\partial }{\partial 𝐧}]{\psi }_{ki}(𝐱)=0,𝐱\in S\end{array}$$ (2b)

where the eigenvalues, ${\mu }_{ki}$, and related eigenfunctions, ${\psi }_{ki}(𝐱)$ are assumed to be known from analytical expressions or hybrid approaches for eigenvalue problems, including the GITT itself (Cotta 1993, Sphaier & Cotta 2000, Naveira-Cotta et al. 2009). Eqs. (2a,b), through the corresponding orthogonality property, allow for the definition of the integral transform-inverse pair:

$$\begin{array}{cc}\hfill & {\bar{T}}_{k,i}(t)={\int }_v{w}_k(𝐱){\tilde{\psi }}_{ki}(𝐱)T_k(𝐱,t)dv,\mathrm{\text{transforms}}\hfill \end{array}$$ (3a)

$$\begin{array}{cc}\hfill & T_k(𝐱,t)=\sum _{i=1}^{\infty }{\tilde{\psi }}_{ki}(𝐱){\bar{T}}_{k,i}(t),\mathrm{\text{inverses}}\hfill \end{array}$$ (3b)

where the symmetric kernels ${\tilde{\psi }}_{ki}(𝐱)$are given by

$$\begin{array}{cc}\hfill & {\tilde{\psi }}_{ki}(𝐱)=\frac{{\psi }_{ki}(𝐱)}{N_{ki}^{1/2}}\hfill \end{array}$$ (3c)

$$\begin{array}{cc}\hfill & N_{ki}={\int }_v{w}_k(𝐱){\psi }_{ki}^2(𝐱)dv\hfill \end{array}$$ (3d)

The integral transform of Eq. (1a) is then accomplished through the operator ${\int }_v{\tilde{\psi }}_{ki}(𝐱)\_\_\_dv$, which after employing boundary conditions, Eqs. (1c) and (2b), yields the following transformed ordinary differential system:

$$\begin{array}{cc}\hfill & \frac{d{\bar{T}}_{k,i}(t)}{dt}+\sum _{j=1}^{\infty }{a}_{kij}(t,T_{\ell }){\bar{T}}_{k,j}(t)={\bar{g}}_{ki}(t,T_{\ell }),i=1,2,...,t>0,k,\ell =1,2,...,m\hfill \end{array}$$ (4a)

The initial conditions, Eq. (1b), are also integral transformed through the operator ${\int }_v{w}_k(𝐱){\tilde{\psi }}_{ki}(𝐱)\_\_\_dv$ to yield:

$$\begin{array}{cc}\hfill & {\bar{T}}_{k,i}(0)={\bar{f}}_{ki}\equiv {\int }_v{w}_k(𝐱){\tilde{\psi }}_{ki}(𝐱)f_k(𝐱)dv\hfill \end{array}$$ (4b)

where

$$\begin{array}{cc}\hfill & {\bar{g}}_{ki}(t,T_{\ell })={\int }_v{\tilde{\psi }}_{ki}(𝐱)P_k(𝐱,t,T_{\ell })dv+{\int }_S{K}_k(𝐱)[{\tilde{\psi }}_{ki}(𝐱)\frac{\partial T_k(𝐱,t)}{\partial 𝐧}-T_k(𝐱,t)\frac{\partial {\tilde{\psi }}_{ki}(𝐱)}{\partial 𝐧}]ds\hfill \end{array}$$ (4c)

$$\begin{array}{cc}\hfill & a_{kij}(t,T_{\ell })={\delta }_{ij}{\mu }_{ki}^2+a_{kij}^{*}(t,T_{\ell })\hfill \end{array}$$ (4d)

with

$$\begin{array}{cc}\hfill & {\delta }_{ij}=\{\begin{array}{c}0,\mathrm{\text{for}}i\ne j\hfill \\ 1,\mathrm{\text{for}}i=j\hfill \end{array}\hfill \end{array}$$ (4e)

$$\begin{array}{cc}\hfill & a_{kij}^{*}(t,T_{\ell })={\int }_v{\tilde{\psi }}_{ki}(𝐱)[𝐮(𝐱,t,T_{\ell }).\nabla {\tilde{\psi }}_{ki}(𝐱)]dv\hfill \end{array}$$ (4f)

Equations (4) form an infinite ODE system of coupled nonlinear equations for the transformed potentials, ${\bar{T}}_{k,i}$. For computational purposes, system (4) is truncated at a large enough order for the required precision. Both commercial and public domain routines are available to numerically handle system (4) under automatic error control, such as routine NDSolve of the Mathematica platform (Wolfram 2017), as employed here. Once the transformed potentials are computed along the time variable, inversion formula, Eq. (3b), can be used to reconstruct to desired potentials $T_k(𝐱,t)$ in explicit form.

Flow Problem in Variably Saturated Porous Media

One-dimensional vertical flow of water in a variably-saturated rigid soil under isothermal conditions is modeled here using the transient one-dimensional version of the Richards equation, in terms of the pressure head h, given as (Simunek et al. 2005):

$$\begin{array}{cc}\hfill & C(h)\frac{\partial h(z,t)}{\partial t}=\frac{\partial }{\partial z}[K(h)\frac{\partial h}{\partial z}]-D(h)\frac{\partial h}{\partial z},0\le z\le L,t>0\hfill \end{array}$$ (5a)

where C(h)=d/dh presents the specific soil water content capacity as a function of the pressure head, $\theta$ is the volumetric water content, K is the hydraulic conductivity, D(h)=dK/dh is a soil water diffusivity type function, t is time, and z denotes the vertical distance from the soil surface downward.

The initial and boundary conditions are given by:

$$\begin{array}{cc}\hfill & h(z,0)=f(z)\hfill \end{array}$$ (5b)

$$\begin{array}{cc}\hfill & \gamma (z)h(z,t)+\beta (z)K(h)(\frac{\partial h}{\partial z}-1)=\phi (z,t),\mathrm{\text{at}}z=0\mathrm{\text{and}}z=L\hfill \end{array}$$ (5c,d)

Boundary conditions (5c,d) are written in a general form so as to recover prescribed head conditions ($\beta =0$), prescribed flux conditions ($\gamma =0$), or free drainage conditions ($\gamma =0$ and $\phi (z,t)=-\beta K(h)$ at $z=L$).

In order to describe the constitutive relations between water content and the pressure head, as well as between the hydraulic conductivity and the pressure head, which define the soil water retention and hydraulic conductivity relationships, the well-known van Genuchten model has been employed (van Genuchten 1980, Simunek et al. 2005) in the form:

$$\begin{array}{cc}\hfill & \theta (h)={\theta }_r+\frac{{\theta }_s-{\theta }_r}{{[1+{(\alpha h)}^n]}^m}\mathrm{\text{with}}m=1-\frac{1}{n},(n>1)\hfill \end{array}$$ (6a)

$$\begin{array}{cc}\hfill & K(h)=K_s{S}_e^l{[1-{(1-S_e^{1/m})}^n]}^2\hfill \end{array}$$ (6b)

$$\begin{array}{cc}\hfill & S_e=\frac{\theta (h)-{\theta }_r}{{\theta }_s-{\theta }_r}=(1+(\alpha h{)}^n{)}^m\hfill \end{array}$$ (6c)

where $S_e$ is effective saturation, ${\theta }_r$ and ${\theta }_s$ are the residual and saturated water contents, respectively, $K_s$ is the saturated hydraulic conductivity, m and n are quasi-empirical shape parameters, and l is a pore connectivity parameter.

In order to improve convergence of the eigenfunction expansion, the following filtering is proposed:

$$h(z,t)=h^{*}(z,t)+F(z,t)$$ (7)

where the filter F(z,t) satisfies the original boundary conditions so as to achieve homogeneous boundary conditions for the filtered potential, $h^{*}(z,t)$. The filtered problem is then written as:

$$\begin{array}{cc}\hfill & C(h)\frac{\partial h^{*}(z,t)}{\partial t}=\frac{\partial }{\partial z}[K(h)\frac{\partial h^{*}}{\partial z}]-D(h)\frac{\partial h^{*}(z,t)}{\partial z}+P(z,t),0<z<L,t>0\hfill \end{array}$$ (8a)

where

$$\begin{array}{cc}\hfill & P(z,t)=-C(h)\frac{\partial F(z,t)}{\partial t}+\frac{\partial }{\partial z}[K(h)\frac{\partial F}{\partial z}]-D(h)\frac{\partial F(z,t)}{\partial z},0<z<L,t>0\hfill \end{array}$$ (8b)

subject to the filtered initial and boundary conditions:

$$\begin{array}{cc}\hfill & h^{*}(z,0)=f(z)-F(z,0)\hfill \end{array}$$ (8c)

$$\begin{array}{cc}\hfill & \gamma (z)h^{*}(z,t)+\beta (z)K^{*}(z)\frac{\partial h^{*}}{\partial z}=0,\mathrm{\text{at}}z=0\mathrm{\text{and}}z=L\hfill \end{array}$$ (8d,e)

where $K^{*}$ is a representative parameter of the boundary flux, for instance the linearized hydraulic conductivity at the two boundary positions, and where the nonlinear source terms to be eliminated by the filter are given by:

$$\begin{array}{ccc}\hfill & \hfill & {\phi }^{*}(z,t)=\phi (z,t)+\beta (z)[K^{*}(z)-K(h)]\frac{\partial h}{\partial z}+\beta (z)K(h)\hfill \end{array}$$ (8f)

When eliminating the source terms of Eq. (8f), for the more general situation, the filter becomes implicit and has to be solved together with the transformed potentials. In its simplest form, the proposed filter is given by a straight line in the $z$ variable, defined by the two boundary conditions and the respective source terms, ${\phi }^{*}(z,t)$.

Following the formalism in the GITT approach, the following fairly simple eigenvalue problem is chosen:

$$\begin{array}{cc}\hfill & \frac{d^2{\psi }_i(z)}{d{z}^2}+{\mu }_i^2{\psi }_i(z)=0,0<z<L\hfill \end{array}$$ (9a)

$$\begin{array}{cc}\hfill & \gamma (z){\psi }_i(z)+\beta (z)K^{*}(z)\frac{d{\psi }_i}{dz}=0,z=0\mathrm{\text{and}}z=L\hfill \end{array}$$ (9b,c)

where the eigenvalues ${\mu }_i$ and respective eigenfunctions ${\psi }_i(z)$ are known in analytical explicit form. Equations (9a-c) allow definition of the following integral transform-inverse pair:

$$\begin{array}{cc}\hfill & \mathrm{\text{Transform:}}{\bar{h}}_i^{*}(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)h^{*}(z,t)dz\hfill \end{array}$$ (10a)

$$\begin{array}{cc}\hfill & \mathrm{\text{Inverse:}}{h}^{*}(z,t)=\sum _{i=0}^{\infty }{\tilde{\psi }}_i(z){\bar{h}}_i^{*}(t)\hfill \end{array}$$ (10b)

where the normalized eigenfunctions $\tilde{\psi }$_i(z) and the norms are given by:

$$\begin{array}{cc}\hfill & {\tilde{\psi }}_i(z)=\frac{{\psi }_i(z)}{M_i^{1/2}};M_i^{}=\underset{0}{\overset{L}{\int }}{[{\psi }_i(z)]}^{2}dz\hfill \end{array}$$ (10c,d)

The filtered partial differential equation, Eq. (8a), is then operated on with $\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)\_\_\_dz$, thus eliminating the space variable and yielding the following transformed ODE system in the time variable, $t$:

$$\begin{array}{cc}\hfill & \frac{d{\bar{h}}_i^{*}(t)}{dt}={\bar{g}}_i(t),t>0,i=1,2,...\hfill \end{array}$$ (11a)

where

$$\begin{array}{cc}\hfill & {\bar{g}}_i(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)G(z,t)dz\hfill \end{array}$$ (11b)

and

$$\begin{array}{cc}\hfill & G(z,t)=\frac{1}{C(h)}\frac{\partial }{\partial z}[K(h)\frac{\partial h^{*}(z,t)}{\partial z}]-\frac{D(h)}{C(h)}\frac{\partial h^{*}(z,t)}{\partial z}+\frac{P(z,t)}{C(h)}\hfill \end{array}$$ (11c)

The integral that defines the non-transformable nonlinear coupling term, ${\bar{g}}_i(t)$, can involve significant computational effort for large truncation orders, due to the oscillatory nature of the eigenfunctions. However, a semi-analytical integration procedure has been described in (Cotta & Mikhailov 2005, Cotta et al. 2013) which implies a drastic reduction of computation effort in transforming nonlinear source terms. This is done by dividing the space variable domain in $N_k$ subintervals, and adopting a first order polynomial approximation for the non-oscillatory part of the integrand, $G(z,t)$, at each subinterval, such as:

$$\begin{array}{cc}\hfill & G(z,t)\cong a_k{z}_{}+b_k,z_{k-1}\le z\le z_k\hfill \end{array}$$

$$\begin{array}{cc}\hfill & G(z_k,t)=G_k=a_k{z}_k+b_k\mathrm{\text{and}}G(z_{k-1},t)=G_{k-1}=a_k{z}_{k-1}+b_k\hfill \end{array}$$ (12a-e)

$$\begin{array}{cc}\hfill & a_k=\frac{G_k-G_{k-1}}{z_k-z_{k-1}}\mathrm{\text{and}}{b}_k=G_k-a_k{z}_k\hfill \end{array}$$

Thus,

$${\bar{g}}_i(t)=\sum _{k=1}^{N_k}\underset{z_{k-1}}{\overset{z_k}{\int }}{\tilde{\psi }}_i(z)G(z,t)dz\cong \sum _{k=1}^{N_k}\underset{z_{k-1}}{\overset{z_k}{\int }}{\tilde{\psi }}_i(z)(a_{k}z+b_k)dz$$ (13)

The semi-analytical integration procedure now only requires analytical evaluation of two terms, as follows:

$${\bar{g}}_i(t)\cong \sum _{k=1}^{N_k}{a}_k\left(\underset{z_{k-1}}{\overset{z_k}{\int }}z{\tilde{\psi }}_i(z)dz\right)+b_k\left(\underset{z_{k-1}}{\overset{z_k}{\int }}{\tilde{\psi }}_i(z)dz\right)$$ (14)

Initial condition (8c) is similarly integral transformed through the operator $\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)\_\_\_dz$ yielding:

$${\bar{h}}_i^{*}(0)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)f(z)dz-\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)F(z,0)dz$$ (15)

Substituting the inverse formula (Eq. 10b) into Eqs. (11) and (15), the transformed nonlinear ODE system is formed and then solved for the transformed potentials, ${\bar{h}}_i^{*}(t)$. Once this system has been solved numerically, the inverse formula, Eq. (10b), is again recalled to reconstruct the filtered potential $h^{*}(z,t)$, which is then added to the filter, $F(z,t)$ as in Eq. (7), to yield the head, $h(z,t)$.

Radionuclide Transport with Physical Non-equilibrium

The transport of radioactive contaminants through the solid waste piles and soil layers is analyzed using a physical non-equilibrium (dual-porosity) type model which accounts for both mobile and immobile water contents. A one-dimensional transient advection-dispersion formulation is adopted for the ²³⁸U radioactive chain (Pontedeiro et al. 2010). The elements considered most relevant in the chain, in light of their individual half-lives, were ²³⁸U, ²³⁴U, ²³⁰Th, ²²⁶Ra and ²¹⁰Pb. The dual-porosity transport problem formulation is given as:

$$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}[({\theta }_m+f\rho K_d)C_m(z,t)]=\frac{\partial }{\partial z}[{\theta }_m{D}_m\frac{\partial C_m(z,t)}{\partial z}-q{C}_m(z,t)]-\beta [C_m(z,t)-C_{im}(z,t)]+{\phi }_m\hfill \end{array}$$ (16a)

$$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}[({\theta }_{im}+(1-f)\rho K_d)C_{im}(z,t)]=\beta [C_m(z,t)-C_{im}(z,t)]+{\phi }_{im}\hfill \end{array}$$ (16b)

with initial and boundary conditions:

$$\begin{array}{cc}\hfill & C_m(z,0)=C_{m,back}(z);C_{im}(z,0)=C_{im,back}(z)\hfill \end{array}$$ (16c,d)

$$\begin{array}{cc}\hfill & q{C}_m(0,t)-{\theta }_m{D}_m{\frac{\partial C_m(z,t)}{\partial z}|}_{z=0}=q{C}_{m,z0}(t)\hfill \end{array}$$ (16e)

$$\begin{array}{cc}\hfill & q{C}_m(L,t)+{\theta }_m{D}_m{\frac{\partial C_m(z,t)}{\partial z}|}_{z=L}=q{C}_{m,zL}(t)or{\frac{\partial C_m(z,t)}{\partial z}|}_{z=L}=0\hfill \end{array}$$ (16f,g)

Again, filters are employed to make the boundary conditions homogeneous:

$$\begin{array}{cc}\hfill & C_m(z,t)=C_m^{*}(z,t)+F_m(z,t)\hfill \end{array}$$ (17a)

$$\begin{array}{cc}\hfill & C_{im}(z,t)=C_{im}^{*}(z,t)+F_{im}(z,t)\hfill \end{array}$$ (17b)

which causes the filtered equations to be written as:

$$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}[({\theta }_m+f\rho K_d)C_m^{*}(z,t)]=\frac{\partial }{\partial z}[{\theta }_m{D}_m\frac{\partial C_m^{*}(z,t)}{\partial z}-q{C}_m^{*}(z,t)]-\beta [C_m^{*}(z,t)-C_{im}^{*}(z,t)]+P_m(z,t)\hfill \end{array}$$ (18a)

$$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}[({\theta }_{im}+(1-f)\rho K_d)C_{im}^{*}(z,t)]=\beta [C_m^{*}(z,t)-C_{im}^{*}(z,t)]+P_{im}(z,t)\hfill \end{array}$$ (18b)

in which the source terms are given by:

$$\begin{array}{cc}\hfill & \hfill P_m(z,t)=-\frac{\partial }{\partial t}[({\theta }_m+f\rho K_d)F_m(z,t)]+\frac{\partial }{\partial z}[{\theta }_m{D}_m\frac{\partial F_m(z,t)}{\partial z}-q{F}_m(z,t)]-\beta [F_m(z,t)-F_{im}(z,t)]+{\phi }_m\end{array}$$ (18c)

$$\begin{array}{cc}\hfill & P_{im}(z,t)=-\frac{\partial }{\partial t}[({\theta }_{im}+(1-f)\rho K_d)F_{im}(z,t)]+\beta [F_m(z,t)-F_{im}(z,t)]+{\phi }_{im}\hfill \end{array}$$ (18d)

The initial and boundary conditions are rewritten as:

$$\begin{array}{cc}\hfill & C_m^{*}(z,0)=C_{m,b}(z)-F_m(z,0);C_{im}^{*}(z,0)=C_{im,b}(z)-F_{im}(z,0)\hfill \end{array}$$ (18e,f)

$$\begin{array}{cc}\hfill & q{C}_m^{*}(0,t)-{\theta }_m{D}_m{\frac{\partial C_m^{*}(z,t)}{\partial z}|}_{z=0}=0\hfill \end{array}$$ (18g)

$$\begin{array}{cc}\hfill & q{C}_m^{*}(L,t)+{\theta }_m{D}_m{\frac{\partial C_m^{*}(z,t)}{\partial z}|}_{z=L}=0or{\frac{\partial C_m^{*}(z,t)}{\partial z}|}_{z=L}=0\hfill \end{array}$$ (18h,i)

The integral transform method then starts with a proposed eigenvalue problem in the simple form:

$$\begin{array}{cc}\hfill & \frac{d^2{\psi }_i(z)}{d{z}^2}+{\mu }_i^2{\psi }_i(z)=0,0<z<L\hfill \end{array}$$ (19a)

$$\begin{array}{cc}\hfill & q{\psi }_i(0)-{\theta }_m{D}_m{\frac{d{\psi }_i}{dz}|}_{z=0}=0\hfill \end{array}$$ (19b)

$$\begin{array}{cc}\hfill & q{\psi }_i(L)+{\theta }_m{D}_m{\frac{d{\psi }_i}{dz}|}_{z=L}=0or{\frac{d{\psi }_i}{dz}|}_{z=L}=0\hfill \end{array}$$ (19c,d)

where the eigenvalues ${\mu }_i$ and corresponding eigenfunctions ${\psi }_i(z)$ are known in explicit analytical form, to allow the definition of the following integral transform-inverse pairs for the mobile and immobile phases:

Mobile:

$$\begin{array}{ccc}\hfill & \mathrm{\text{Transform:}}\hfill & \hfill {\bar{C}}_{m,i}^{*}(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)C_m^{*}(z,t)dz\end{array}$$ (20a)

$$\begin{array}{ccc}\hfill & \mathrm{\text{Inverse:}}\hfill & \hfill C_m^{*}(z,t)=\sum _{i=0}^{\infty }{\tilde{\psi }}_i(z){\bar{C}}_{m,i}^{*}(t)\end{array}$$ (20b)

Immobile:

$$\begin{array}{ccc}\hfill & \mathrm{\text{Transform:}}\hfill & \hfill {\bar{C}}_{im,i}^{*}(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)C_{im}^{*}(z,t)dz\end{array}$$ (21a)

$$\begin{array}{ccc}\hfill & \mathrm{\text{Inverse:}}\hfill & \hfill C_{im}^{*}(z,t)=\sum _{i=0}^{\infty }{\tilde{\psi }}_i(z){\bar{C}}_{im,i}^{*}(t)\end{array}$$ (21b)

where the normalized eigenfunctions ${\tilde{\psi }}_i(z)$ and the normalization integrals are given by:

$${\tilde{\psi }}_i(z)=\frac{{\psi }_i(z)}{N_i^{1/2}};N_i^{}=\underset{0}{\overset{L}{\int }}{[{\psi }_i(z)]}^{2}dz$$ (22a,b)

Applying the operator $\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)\_\_\_dz$ on Eq. (18a,b), the transformed ODE system is obtained for the transformed potentials, ${\bar{C}}_{m,i}^{*}(t)$ and ${\bar{C}}_{im,i}^{*}(t)$:

Mobile:

$$\begin{array}{cc}\hfill & \sum _{j=1}^{\infty }{A}_{m,ij}(t)\frac{d{\bar{C}}_{m,j}^{*}(t)}{dt}=\sum _{j=1}^{\infty }{D}_{m,ij}(t){\bar{C}}_{m,j}^{*}(t)+\beta [{\bar{C}}_{im,i}^{*}(t)-{\bar{C}}_{m,i}^{*}(t)]+{\bar{p}}_{m,i}(t)\hfill \end{array}$$ (23a)

with transformed initial conditions

$$\begin{array}{cc}\hfill & \hfill {\bar{C}}_{m,i}^{*}(0)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)C_{m,b}^{*}(z)dz-\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)F_m(z,0)dz\end{array}$$ (23b)

Immobile:

$$\begin{array}{cc}\hfill & \sum _{j=1}^{\infty }{A}_{im,ij}(t)\frac{d{\bar{C}}_{im,j}^{*}(t)}{dt}=\sum _{j=1}^{\infty }{A}_{im,ij}^{\prime }(t){\bar{C}}_{im,j}^{*}(t)+\beta [{\bar{C}}_{m,i}^{*}(t)-{\bar{C}}_{im,i}^{*}(t)]+{\bar{p}}_{im,i}(t)\hfill \end{array}$$ (23c)

with transformed initial conditions

$$\begin{array}{cc}\hfill & \hfill {\bar{C}}_{im,i}^{*}(0)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)C_{im,b}^{*}(z)dz-\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)F_{im}(z,0)dz\end{array}$$ (23d)

where

$$\begin{array}{ccc}\hfill & A_{m,ij}(t)=\underset{0}{\overset{L}{\int }}({\theta }_m+f\rho K_d){\tilde{\psi }}_i(z){\tilde{\psi }}_j(z)dz;\hfill & \hfill A_{m,ij}^{\prime }(t)=\frac{d{A}_{m,ij}(t)}{dt}\\ \hfill & B_{m,ij}(t)=\underset{0}{\overset{L}{\int }}({\theta }_m{D}_m)\frac{d{\tilde{\psi }}_i(z)}{dz}\frac{d{\tilde{\psi }}_j(z)}{dz}dz;\hfill & \hfill C_{m,ij}(t)=\underset{0}{\overset{L}{\int }}q\frac{d{\tilde{\psi }}_i(z)}{dz}{\tilde{\psi }}_j(z)dz\\ \hfill & D_{m,ij}(t)=-A_{ij}^{\prime }(t)-B_{ij}(t)+C_{ij}(t);\hfill & \hfill {\bar{p}}_{m,i}(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)P_m(z,t)dz\end{array}$$ (23e-m)

$$\begin{array}{ccc}\hfill & A_{im,ij}(t)=\underset{0}{\overset{L}{\int }}({\theta }_{im}+(1-f)\rho K_d){\tilde{\psi }}_i(z){\tilde{\psi }}_j(z)dz;\hfill & \hfill A_{im,ij}^{\prime }(t)=\frac{d{A}_{im,ij}(t)}{dt}\\ \hfill & {\bar{p}}_{im,i}(t)=\underset{0}{\overset{L}{\int }}{\tilde{\psi }}_i(z)P_{im}(z,t)dz\hfill \end{array}$$

Routine NDSolve from the Mathematica system (Wolfram 2017) is again employed to numerically solve Eqs. (23a-d), after which the inverses, Eqs. (20b) and (21b), are recalled to reconstruct the filtered potentials, $C_m^{*}(z,t)$ and $C_{im}^{*}(z,t)$, which are then summed to the filters, $F_m(z,t)$ and $F_{im}(z,t)$, to recover the desired original potentials, $C_m^{}(z,t)$ and $C_{im}^{}(z,t)$.

VERIFICATION OF THE GITT SOLUTIONS

Solutions of both the flow problem and of the radionuclides chain transport problem were implemented using the mixed symbolic-numerical platform of the Mathematica system (Wolfram 2017). Before proceeding to the physical analysis relevant to the radiological impact assessment, these solutions were verified against literature results and independent numerical simulations.

Flow Problem

For verification of the flow problem GITT solution, we selected the work of (Celia et al. 1990). The constitutive relations employed in (Celia et al. 1990) were implemented here in the form:

$$\theta (h)=\frac{\alpha ({\theta }_s-{\theta }_r)}{\alpha +{\left|h\right|}^{\beta }}+{\theta }_r;K(h)=K_s\frac{A}{A+{\left|h\right|}^{\gamma }}$$ (24a,b)

where $\alpha$=1.611x10⁶, $\theta$_s=0.287, $\theta$_r=0.075, $\beta$=3.96, K_s =0.00944 cm/s, A=1.175x10⁶ and $\gamma$=4.74, with initial and boundary conditions given by:

$$\begin{array}{cc}\hfill & h(z,0)=-61.5\mathrm{\text{cm}}\hfill \end{array}$$ (25a)

$$\begin{array}{cc}\hfill & h(40\mathrm{\text{cm}},t)=h_{top}=-20.7\mathrm{\text{cm}}\hfill \end{array}$$ (25b)

$$\begin{array}{cc}\hfill & h(0,t)=h_{bottom}=-61.5\mathrm{\text{cm}}\hfill \end{array}$$ (25c)

Before comparing our results with those reported in (Celia et al. 1990), a convergence analysis of the GITT solution was carried out. This was done by independently analyzing both the transformed system truncation order, $N$, and the number of subintervals used in the semi-analytical integration, $N_k$. Sample results from this analysis, shown in Tables I and II, indicate convergence of the pressure head to at least the third significant digit up to a truncation order of $N=30$ and a total number of subintervals up to $N_k=180$.

Table I examines the influence of the number of subintervals in the semi-analytical integration, on the local results for the pressure head at z=15 cm and t=360 s, while Table II details the influence of the overall truncation order on convergence for increasing order in the summation. Clearly, both parameters need to be examined simultaneously for error-controlled solution of the potentials.

Figure 3 graphically compares fully converged results to the third significant digit (N=30 and N_k =120) with results obtained by (Celia et al. 1990). The two independent approaches showed excellent adherence.

Figure 3
Comparison of GITT solution (black solid line) with numerical results by (Celia et al. 1990) (red dots) for vertical infiltration in a soil column (t = 360 s).

Table I
Convergence analysis of the GITT eigenfunction expansion for varying truncation order N (results for -h(15, 360), N_k = 120).

Order	$N_k$ = 60	$N_k$ = 120	$N_k$ = 180
1	49.660	49.739	49.762
2	40.446	40.582	40.617
3	41.843	41.994	42.033
4	42.435	42.651	42.707
5	41.159	41.378	41.434
6	41.545	41.795	41.861
7	40.705	40.973	41.044
8	40.705	40.973	41.044
9	40.375	40.665	40.742
10	39.982	40.268	40.344
11	40.045	40.336	40.413
12	39.674	39.969	40.048
13	39.618	39.910	39.988
14	39.557	39.851	39.930
15	39.449	39.740	39.818

Table II
Convergence analysis of the GITT eigenfunction expansion for varying truncation order N (results for -h(15, 360), N_k = 120).

Order	$N$ = 20	$N$ = 25	$N$ = 30
1	49.441	49.306	49.236
2	40.105	39.894	39.786
3	41.468	41.231	41.111
4	41.857	41.503	41.323
5	40.570	40.214	40.033
6	40.851	40.433	40.221
7	39.922	39.464	39.234
8	39.922	39.464	39.234
9	39.492	38.986	38.734
10	39.102	38.607	38.361
11	39.131	38.619	38.365
12	38.724	38.204	37.948
13	38.681	38.172	37.922
14	38.588	38.067	37.812
15	38.433	37.924	37.676
16	38.433	37.924	37.676
17	38.355	37.845	37.599
18	38.342	37.836	37.594
19	38.315	37.808	37.566
20	38.327	37.798	37.559
21		37.806	37.568
22		37.785	37.548
23		37.802	37.564
24		37.802	37.564
25		37.808	37.561
26			37.572
27			37.572
28			37.574
29			37.577
30			37.581

Transport Problem

Verification of the transport problem solution was performed first for the physical non-equilibrium model of the ²³⁸U radioactive chain. We used for this purpose the numerical routine NDSolve of the Mathematica system (Wolfram 2017), which implements a Method of Lines numerical algorithm for partial differential equations. Input data were extracted from Case 1 to be discussed in the next section, including the infiltration rate in the solid waste pile, vertical transport in the soil beneath the pre-basin, and subsequent transport in the near horizontal aquifer. First a convergence analysis was again performed as illustrated in Table III below, for the mobile concentrations of ²³⁸U at the end of the solid waste pile. A maximum truncation order of $N=20$ produced full convergence to five significant digits in the concentrations, although graphical convergence was achievable already for N<10.

Table III
Convergence of the GITT solution for the mobile concentration of ²³⁸U [Bq/l] at the end of the solid waste pile at different times [years]: Case 1.

Order	C(L,500)	C(L,5000)	C(L,50.000)
2	214.71	3.6609	3.5797
4	205.04	3.6782	3.5797
6	204.71	3.6812	3.5797
8	204.64	3.6820	3.5797
10	204.62	3.6822	3.5797
12	204.61	3.6823	3.5797
14	204.60	3.6824	3.5797
16	204.60	3.6824	3.5797
18	204.60	3.6824	3.5797
20	204.60	3.6824	3.5797

Figure 4 compares ²³⁸U concentration, again in [Bq/l], versus time (in years) at the bottom of the solid waste pile (black curves), at the bottom of the vertical soil (blue curves), and at the end of the granular aquifer (red curves), as computed using GITT (solid lines) and NDSolve (symbols) for Case 1. For the GITT solution we used a truncation order of N=20 in the eigenfunction expansion, while the maximum steps in the discretization of the NDSolve numerical solution were restricted to L/100 and t_f /1000. The excellent agreement between the two sets of independent results is evident.

Figure 4
Comparison of the GITT solutions (solid lines) and the NDSolve numerical solution (Wolfram 2017) (symbols) for ²³⁸U concentrations in [Bq/l] at the end of three layers: the solid waste pile (black curves), the soil below the pre-basin (blue curves), and the granular aquifer (red curves).

Next, the hybrid GITT solution was verified against the well-established finite element Hydrus 1D code (Simunek et al. 2005, 2016), with special emphasis on radionuclide concentrations of the radioactive chain at the end of the solid waste pile. Figure 5 compares ²³⁸U concentrations at the end of the waste pile (L=10 m), again with excellent agreement, for cases 1 and 3 to be detailed in what follows. The verification analyses for variably-saturated flow and radionuclide transport provided here offer enough confidence about the accuracy of the implemented hybrid numerical-analytical solutions to be employed in the complete radionuclide transport problem to be discussed next.

Figure 5
Comparison of GITT solutions (solid black line) against Hydrus 1D simulation results (red dots) for the ²³⁸U concentrations (in Bq/l) versus time (in years) at the bottom of the solid waste pile for Cases 1 (left) and 3 (right).

Analysis of Radionuclide Transport at URA/Caetité

Next, the physical aspects of radionuclide migration at the URA/Caetité site are explored by analyzing concentration evolutions for each radionuclide based on the data compiled in Table IV. Those data include parameters that were determined experimentally using laboratory and field measurements, as well as some that were extracted from the report (Orlande et al. 2008) as specified in the table. Table V summarizes four separate cases that were analyzed. A fifth case was also examined in order to evaluate uncertainty in the initial conditions of the daughter products.

Table IV
Input data for the radionuclide transport simulations.

PARAMETERS				SOLID WASTE PILE	VERTICAL SOIL (Transition soil)	HORIZONTAL SOIL (Granular aquifer)
FLOW PROBLEM				unsaturated	unsaturated	saturated
1	Layer length	L	[m]	10	1.40	1488
2	Final simulation time	tfinal	[years]	50.000	50.000	50.000
3	Infiltration and recharge rates	q	[m/year]
	(Cases 1 and 2)			0.471	35.2	50.4
	(Cases 3 and 4)			0.802	50.4	65.6
4	Natural water content	$\theta$ 0	-	0.055	0.088	0.074
5	Residual water content	$\theta$ r	-	0.050	0.043	0.046
6	Saturated water content	$\theta$ s	-	0.317	0.290	0.388
7	Saturated hydraulic conductivity	Ks	[cm/h]	1.06	3.11	6.48
8	Shape parameter	$\alpha$	[1/cm]	0.0316	0.0396	0.0359
9	Shape parameter	n	-	3.93	2.48	3.36
10	Shape parameter	m	-	1-1/n	1-1/n	1-1/n
11	Pore connectivity	l	-	0.5	0.5	0.5
TRANSPORT PROBLEM				Physical equilibrium	Physical non-equilibrium	Physical non-equilibrium
12	Water content	$\theta$	-
	(Cases 1 and 2)			0.106	0.290	0.388
	(Cases 3 and 4)			0.121	0.290	0.388
13	Apparent density	$\rho$ b	[g/cm3]	1.73	1.61	1.66
14	Dispersivity	$\alpha$ x	[m]	1.0	0.14	148.8
15	Dispersion coefficient	D	[m2/year]	$\alpha$x.q/$\theta$	$\alpha$x.q/$\theta$	$\alpha$x.q/$\theta$
16	Distribution coefficient (U)	Kd	[cm3/g]	44	298	122
17	Partition coefficient	f	-	0.999	0.64	0.64
18	Mass transfer coefficient (Cases 1 and 3)	$\beta$	[1/h]	0	0.058	4.6x10-4
					0.079	5.2x10-4
19	Mobile water content	$\theta$ m	-	-	f.$\theta$	f.$\theta$
20	Dispersion coefficient - Mobile	Dm	[m2/year]	-	D.$\theta$/$\theta$m	D.$\theta$/$\theta$m
21	Dissolution constant	K	[1/year]	47.8	0	0
22	Equilibrium concentration	Ceq	[µg/dm3]	286.6	0	0
23	Initial condition - Mobile	Cm (z,$0$)	[µg/dm3]	16,430	0	0
24	Initial condition - Immobile	Cim (z,$0$)	[µg/dm3]	0	0	0
25	Boundary condition (z=0)	Cm ($0$,t)	[µg/dm3]	0	C0v (t)	C0h (t)
26	Boundary condition (z=L)	-	-	2nd kind	2nd kind	2nd kind

• Source: Rosetta estimates (Simunek et al. 2005) for parameters (5), (8), (9), (11). Experiments and direct measurements for the remaining parameters.

Table V
Simulated cases in the radionuclide transport analysis.

Case	Infiltration [m/year]	Runoff Sterile rock	Waste Pile Model	Vertical Soil Model	Horizontal Soil Model
1	0.471	0	Equilibrium	Non-Equilibrium	Non-Equilibrium
2	0.471	0	Equilibrium	Equilibrium	Equilibrium
3	0.802	0.2	Equilibrium	Non-Equilibrium	Non-Equilibrium
4	0.802	0.2	Equilibrium	Equilibrium	Equilibrium
5	Same data as case 3, but with initial concentration of Ra equal to 10 Bq/l

Case 1 constitutes the base scenario by making use of a complete hydrologic balance to estimate the infiltration rate (Orlande et al. 2008) and adopting the conservative hypothesis of having physical equilibrium during transport within the solid waste pile (i.e., without preferential paths during transport within the pile, which otherwise would have resulted in lower concentrations leaving the waste pile as exemplified in (Pontedeiro et al. 2010)). Case 1 is also conservative by considering physical non-equilibrium (dual-porosity) transport in both the vertical soil layer and the granular aquifer. Case 2 considers the situation of physical equilibrium for all three layers (the waste pile, the vertical soil, and the aquifer) so as to illustrate the influence of the preferential flow paths on radionuclide migration.

Case 3 differs from case 1 in terms of the imposed infiltration rate into the waste and sterile rock piles (0.802 rather than 0.471 m/year). Case 1 accounts for a complete hydrologic balance of the region leading to a lower long-time infiltration rate (0.471 m/year), while case 3 is very conservative by using for the infiltration rate the total annual precipitation rate, discounted by 20% runoff at the sterile piles only. Case 4 is similar to case 2 by analyzing the effect of preferential paths on radionuclide transport but adopting the larger infiltration rate of case 3.

For the four cases above we assumed that the initial radionuclide concentration was zero, except for uranium, since the only information available concerned the aqueous solution from the third and final washing of the leaching platform. Therefore, an additional case 5 was used by employing an average value for the Ra concentration (10 Bq/l) as determined from the aqueous solution. For Th and Pb we maintained zero initial concentrations.

Figures 6 present the concentration evolutions according to case 1, for the five radionuclides in the considered radioactive decay chain at, respectively, the end of the waste pile, the end of the vertical soil layer, and the end of the granular aquifer. Results show initially still relatively high uranium concentrations (about 100 Bq/l for both ²³⁸U and ²³⁴U) at the bottom of the waste pile due to washing of the original leachate ore (Fig. 6, top left). This is followed by a much longer period where reactions and decay within the waste pile cause the uranium concentrations at the bottom of the pile to decrease to approximately 3.6 Bq/l for the two uraniums. Clearly, the uraniums dominate relative to the other radionuclides, in particular within the first 5000 years. At later times, with the attainment of non-equilibrium conditions, the uraniums remain with a radioactivity of one order of magnitude larger than that of the other radionuclides, in the following order of importance, Pb, Ra and Th (which all were assumed to have zero initial concentration).

Figure 6
Case 1 concentration evolutions of the ²³⁸U decay chain at: (a) the end of the waste pile; (b) the end of the vertical soil; (c) the end of the granular aquifer.

Results for the vertical soil below the pre-basin (Fig. 6, top right) show the important dilution due to the mixing at the pre-basin with water leached through the sterile rock pile, but with overall transient behavior similar to the waste pile, also in view of the relatively small length of the vertical soil layer (1.4 m). Results in Fig. 6 (bottom plot) show that the contamination plume did not immediately reach the end of the granular aquifer, but once there, the concentration increased remarkably due to the preferential paths as simulated with the physical non-equilibrium model. Reactions and decay subsequently lowered the uranium concentrations. With respect to the final values of the mobile water concentrations, the sum of the activities of all five radionuclides did produce a total concentration of 1.65 Bq/l, with the two uraniums totaling 1.47 Bq/l and the other three radionuclides contributing approximately 0.18 Bq/l (0.11 Bq/l for Pb, 0.06 Bq/l for Ra and 0.01 Bq/l for Th). The maximum concentration of all radionuclides at the edge of the aquifer occurred after about 4000 years, with the two uraniums contributing equally to a total of approximately 34.6 Bq/l.

Figures 7 illustrate, still for the base case (case 1), the comparative behavior of the concentration evolutions of the five radionuclides, at the end of each of the three layers (the waste pile, the vertical soil, and the granular aquifer). Results here are plotted for a much shorter total time to more clearly visualize the initial transients up to steady state. The data show that the plume already arrived at the end of the aquifer before most of the ore (notably uranium) was leached out of the waste pile, followed by near-equilibrium conditions throughout the system. The marked dilution that occurs at the pre-basin is also clearly noticeable in that concentrations reduced significantly at the bottom of the vertical soil profile just below the pre-basin. The behavior of the Th is then followed by Ra and Pb, which in this case were all assumed to have zero initial concentrations.

Figure 7
Case 1 concentrations versus time at the end of the three layers: (a) of ²³⁸U; (b) of ²³⁴U; (c) of ²³⁰Th; (d) of ²²⁶Ra; (e) of ²¹⁰Pb.

Figures 8 and 9 present concentration evolutions for each element of the ²³⁸U radioactive chain at the end of the waste pile and at the end of the granular aquifer, respectively, for cases 1 to 4. Results for the waste pile (Fig. 8) indicate no distinctions between cases 1 and 2 or 3 and 4, since the same model is used for this layer (i.e., physical equilibrium). Results for case 3 in Fig. 9 show a more pronounced advance of the contamination front since a much higher infiltration rate was used (i.e., the total average precipitation rate of 0.802 cm/year).

Figure 8
Concentrations versus time at the end of waste pile for cases 1 to 4: (a) of ²³⁸U; (b) of ²³⁴U; (c) of ²³⁰Th; (d) of ²²⁶Ra; (e) of ²¹⁰Pb.

Figure 9
Concentrations versus time at the end of granular aquifer for cases 1 to 4: (a) of ²³⁸U; (b) of ²³⁴U; (c) of ²³⁰Th; (d) of ²²⁶Ra; (e) of ²¹⁰Pb; (f) also shows the total concentrations at the end of the granular aquifer.

On the other hand, Fig. 9 refers to the end of the aquifer where the differences among the four cases become more pronounced. In particular, notice that the two cases assuming physical equilibrium in the two soil layers (cases 2 and 4), especially the granular aquifer, produce a significant delay in the arrival of the contamination plume, in comparison to the two cases that allow for preferential flow paths (cases 1 and 3). The maximum concentration values are equally affected. With respect to cases 1 and 3, as previously pointed out, the larger infiltration rate considered for case 3 produced a much faster transport rate of the contaminants, albeit with lower maximum concentration values due to the more pronounced dilution within the pre-basin. While for case 1 the maximum concentration of the two uraniums occur at approximately 4000 years at values of about 34.6 Bq/l, for case 3 the maximum concentration is reached after approximately 2900 years but with a lower concentration value of 33.4 Bq/l. For case 2, the maximum concentration was 24.3 Bq/l after 5700 years, while for case 4 the concentration was 22.9 Bq/l, which occurred after 4200 years.

Finally, case 5 was included as a variant of case 3, to account for an initial non-zero Ra concentration condition as measured on the aqueous solution of the third washing (approximately 10 Bq/l). This case allows for a possible sensitivity analysis by considering non-zero initial values for the concentrations of the other radionuclides. Figures 10 illustrate the effect on the concentration evolutions for Ra and Pb at the end of the waste pile, with comparisons again relative to the transport of the other radionuclides. Results indicate the increased importance of the two radionuclides (Ra and Pb), with a relatively fast decay of Ra leading to even larger activities of Pb. Please note that at a certain moment in time, the activity of Pb becomes even larger than that of the two uraniums (²³⁸U and ²³⁴U), although during the initial 1500 years, the summed activities of the two uraniums are still higher. Again, for large times within the equilibrium situation, the activities of the two uraniums predominated.

Figure 10
Concentration evolution at the end of the waste pile with the initial Ra concentration of 10 Bq/l (case 5): (a) for the radionuclides Ra and Pb; (b) for all radionuclides - linear scale; (c) for all radionuclides - logarithmic scale.

CONCLUSIONS

This work presents a hybrid numerical-analytical solution, based on integral transforms, of the partial differential equations governing transient fluid flow and physical non-equilibrium transport of a radionuclide decay chain in a variably saturated porous medium. The methodology is employed in the analysis of radionuclides migration at an actual uranium mining and milling installation in Brazil. The proposed physical model involves infiltration and dissolution along the leached ore pile, dilution in a pre-basin with water originating from a sterile rock pile, transport along a transitional vertical soil layer below the pre-basin, application of a mixing model at the interface with the underlying granular aquifer, and finally lateral transport along the granular aquifer. The methodology was implemented in the mixed symbolic-numerical Mathematica platform. The hybrid solution was thoroughly analyzed in terms of convergence behavior and verified against literature results and independent numerical solutions, in particular solutions via the Hydrus 1D software for physical non-equilibrium transport. A thorough analysis of the concentration evolutions for each element of the radionuclide chain was also performed by considering five cases that account for variations in the infiltration/recharge rates, the presence of preferential paths, and possible initial conditions for the radionuclides (notably Ra) within the leached ore pile. The proposed methodology appears particularly useful for evaluating radioactivity doses under different scenarios for regulatory purposes.

ACKNOWLEGMENTS

The authors are grateful for financial support provided by INB (Indústrias Nucleares do Brasil), under grant no. COPPETEC PEM-7957. This work was also partially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), all of them government sponsoring agencies in Brazil.

REFERENCES

• ALMEIDA GL, PIMENTEL LCG & COTTA RM. 2008. Integral transform solutions for atmospheric pollutant dispersion. Environ Model Assess 13(1): 53-65.
• BARROS FPJ & COTTA RM. 2007. Integral transforms for three-dimensional steady turbulent dispersion in rivers and channels. Appl Math Model 31: 2719-2732.
• BARROS FPJ, MILLS WB & COTTA RM. 2006. Integral transform solution of two-dimensional model for contaminants dispersion in rivers and channels with spatially variable coefficients. Environ Model Software 21: 699-709.
• CELIA MA, BOULOUTAS ET & ZARBA RL. 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour Res 26(7): 1483-1496.
• COTTA RM. 1990. Hybrid numerical-analytical approach to nonlinear diffusion problems. Numer Heat Transfer, Part B 17: 217-226.
• COTTA RM. 1993. Integral transforms in computational heat and fluid flow. Boca Raton, FL: CRC Press.
• COTTA RM. 1994. Benchmark results in computational heat and fluid flow: The integral transform method. Int J Heat Mass Transfer, Invited Paper 37: 381-394.
• COTTA RM. 1998. The integral transform method in thermal and fluids sciences and engineering. New York, NY: Begell House.
• COTTA RM, KNUPP DC & NAVEIRA-COTTA CP. 2016. Analytical heat and fluid flow in microchannels and microsystems. Springer, New York, NY: Mechanical Engineering Series.
• COTTA RM, KNUPP DC, NAVEIRA-COTTA CP, SPHAIER LA & QUARESMA JNN. 2013. Unified integral transforms algorithm for solving multidimensional nonlinear convection-diffusion problems. Num Heat Transfer, part A - Applications 63: 1-27.
• COTTA RM, KNUPP DC & QUARESMA JNN. 2017. Analytical methods in heat transfer. In: Kulacki FA et al. (Eds), Handbook of Thermal Science and Engineering, Chapter 1, Springer, International Publishing.
• COTTA RM & MIKHAILOV MD. 1997. Heat conduction: Lumped analysis, integral transforms, symbolic computation. New York, NY: Wiley Interscience.
• COTTA RM & MIKHAILOV MD. 2005. Semi-analytical evaluation of integrals for the generalized integral transform technique. In: Proc. of the 4th Workshop on Integral Transforms and Benchmark Problems – IV WIT, CNEN, Rio de Janeiro, Brasil.
• COTTA RM & MIKHAILOV MD. 2006. Hybrid methods and symbolic computations. In: Minkowycz WJ, Sparrow EM & Murthy JY (Eds), Handbook of Numerical Heat Transfer, New York, NY: John Wiley, 2nd ed.
• COTTA RM, MIKHAILOV MD & RUPERTI JR NJ. 1998. Analysis of radioactive waste contamination in soils: Solution via symbolic manipulation. In: Proc. Int. Conf. on the Radiological Accident of Goiânia; 10 Years Later, p. 298-308. Brazil, October 1997; also in IAEA Conf. Proceedings Series, Vienna, IAEA-GOCP.
• COTTA RM, NAVEIRA-COTTA CP, KNUPP DC, ZOTIN JLZ, PONTES PC & ALMEIDA AP. 2018. Recent advances in computational-analytical integral transforms for convection-diffusion problems. Heat Mass Transfer, Invited Paper 54: 2475-2496.
• COTTA RM, UNGS MJ & MIKHAILOV MD. 2003. Contaminant transport in finite fractured porous medium: Integral transforms and lumped-differential formulations. Ann Nucl Energy 30(3): 261-285.
• DIEHL P. 2011. Uranium mining and milling wastes: An introduction. http://wise-uranium.org/uwai.html
  » http://wise-uranium.org/uwai.html
• HEILBRON FILHO PFL, PEREZ-GUERRERO JS, PONTEDEIRO EM, RUPERTI JR NJ & COTTA RM. 2002. Analytical model for environmental impact due to solid waste disposal from electricity generation. Hybrid Meth Eng 4(1&2): 1-26.
• HEILBRON FILHO PFL, XAVIER AM, PONTEDEIRO EM & FERREIRA RS. 2004. Segurança nuclear e proteçāo do meio ambiente. Rio de Janeiro, Brazil: E-Papers.
• IAEA. 2003. Extent of environmental contamination by naturally occurring radioactive material (NORM) and technological options for mitigation, Technical Report Series No. 419, Int. Atomic Energy Agency, Vienna, Austria. http://www.pub.iaea. org/MTCD/publications/Pubdetails.asp?pubId=6789.
  » http://www.pub.iaea.
• IAEA. 2004. Safety assessment methodologies for near surface disposal facilities, Vol. 1, Review and enhancement of safety assessment approaches and tools, Int. Atomic Energy Agency, Vienna, Austria.
• LEAL MA & RUPERTI JR NJ. 2000. A numerical study for the two-dimensional solute transport in groundwater via integral transform method. Hybrid Meth Eng 2(1): 111-129.
• LEAL MA & RUPERTI JR NJ. 2001. A hybrid solution for simulation of 2-D contaminant transport in heterogeneous subsurface systems. Hybrid Meth Eng 3(2): 111-129.
• LIU C, SZECSODY JE, ZACHARA JM & BALL WP. 2000. Use of the generalized integral transform method for solving equations of solute transport in porous media. Adv Water Resources 23: 483-492.
• MIKHAILOV MD & OZISIK MN. 1984. Unified analysis and solutions of heat and mass diffusion. New York, NY; also, Dover Publications, New York, NY, 1993: John Wiley.
• MOREIRA PHS, VAN GENUCHTEN MT, ORLANDE HRB & COTTA RM. 2016. Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column. J Hydrol Hydromech 64: 1-15.
• NAVEIRA-COTTA CP, COTTA RM, ORLANDE HRB & FUDYM O. 2009. Eigenfunction expansions for transient diffusion in heterogeneous media. Int J Heat Mass Transfer 52: 5029-5039.
• NAVEIRA-COTTA CP, PONTEDEIRO EM, COTTA RM, SU J & VAN GENUCHTEN MT. 2013. Environmental impact assessment of liquid waste ponds in uranium milling installations. Waste Biomass Valoriz 4(1): 197-211.
• ORLANDE HRB, COTTA RM, JIAN S, COUTO P, NAVEIRA-COTTA CP, MOREIRA PHS & SABINO VS. 2006. Estudo da dispersão de rejeitos radioativos das células dos depósitos de rejeitos líquidos tratados (ponds) na unidade de concentrado de urânio na INB-Caetité. Technical report. COPPETEC Project PEM-7552 (COPPE – INB).
• ORLANDE HRB, COTTA RM, JIAN S, NAVEIRA-COTTA CP, MOREIRA PHS, FONSECA HM, SABINO VS, AYRES JVC, VAN GENUCHTEN MT & QUARESMA JNN. 2008. Caracterização, migração de radionuclídeos e análise de impacto radiológico no sítio das pilhas de minério lixiviado e estéril da URA, Caetité, Bahia. Technical report. COPPETEC Project PEM-7957 (COPPE - INB), September.
• PONTEDEIRO EM, HEILBRON FILHO PFL & COTTA RM. 2007. Assessment of the mineral industry NORM/TENORM disposal in hazardous landfills. J Hazardous Materials 139(3): 563-568.
• PONTEDEIRO EM, VAN GENUCHTEN MT, COTTA RM & SIMUNEK J. 2010. The effects of preferential flow and soil texture on risk assessments of a NORM waste disposal site. J Hazardous Materials 174: 648-655.
• PÉREZ GUERRERO JS, PIMENTEL LCG, OLIVEIRA JR JF, HEILBRON FILHO PFL & ULKE AG. 2012. A unified analytical solution of the steady-state atmospheric diffusion equation. Atmosph Environ 55: 201-212.
• PÉREZ GUERRERO JS, PIMENTEL LCG, SKAGGS TH & VAN GENUCHTEN MT. 2009. Analytical solution of the advection-diffusion transport equation using a change-of-variable and integral transform technique. Int J Heat Mass Transfer 52: 3297-3304.
• PÉREZ GUERRERO JS & SKAGGS TH. 2010. Analytical solution for one-dimensional advection-dispersion transport equation with distance-dependent coefficients. J Hydrology 390(1): 57-65.
• PÉREZ GUERRERO JS, SKAGGS TH & VAN GENUCHTEN MT. 2010. Analytical solution for multi-species contaminant transport in finite media with time-varying boundary conditions. Transp Porous Media 85: 171-188.
• ROCHA RPA & CRUZ MEC. 2001. Integral transform solutions for one-dimensional transient flow and contaminant transport in dual-porosity systems. Hybr Meth Eng 3(2&3): 1-18.
• RUBOL S, BATTIATO I & BARROS FPJ. 2016. Vertical dispersion in vegetated shear flows. Water Resour Res 15(10): 8066-8080.
• SIMUNEK J, VAN GENUCHTEN MT & SEJNA M. 2005. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat and multiple solutes in variably-saturated media, Version 3.0. Volume 1. Department of Environmental Sciences, University of California, Riverside.
• SIMUNEK J, VAN GENUCHTEN MT & SEJNA M. 2016. Recent developments and applications of the hydrus computer software packages. Vadose Zone J 6(15): 1-25.
• SKAGGS TH, JARVIS NJ, PONTEDEIRO EM, VAN GENUCHTEN MT & COTTA RM. 2007. Analytical advection-dispersion model for transport and plant uptake of contaminants in the root zone. Vadose Zone J 6(4): 890-898.
• SPHAIER LA & COTTA RM. 2000. Integral transform analysis of multidimensional eigenvalue problems within irregular domains. Num Heat Transfer Part B 38: 157-175.
• VAN GENUCHTEN MT. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 44(5): 892-898.
• WNA - WORLD NUCLEAR ASSOCIATION. 2019. World uranium mining production, updated March 2019. http://www.world-nuclear.org/information-library/nuclear-fuel-cycle/ mining-of-uranium/world-uranium-mining-production.aspx.
  » http://www.world-nuclear.org/information-library/nuclear-fuel-cycle/ mining-of-uranium/world-uranium-mining-production.aspx
• WOLFRAM S. 2017. Wolfram Mathematica 11. Champaign, IL: Wolfram Research Inc.

Publication Dates

• Publication in this collection
  06 Apr 2020
• Date of issue
  2020

History

• Received
  10 Apr 2019
• Accepted
  11 Sept 2019