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

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 theMathematica 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 238U 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.


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 2 ), 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, 1994, 1998, Cotta & Mikhailov 1997, 2006, Cotta et al. 2016, 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 , 2003, Leal & Ruperti 2000, 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, 2012, Perez Guerrero & Skaggs 2010, 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). 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(Simunek et al. , 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 .
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
The deposit consisted mostly of three separate materials: the sterile soil (estimated at 1.08x10 6 tons) originating from the top cover of the mine (approximately five meters deep), the leachate ore (on the order of 2.11x10 6 tons, with granulometry below 10 mm) originating from piles assembled on the leaching platform, and sterile rocks (estimated at 11.13 x 10 6 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.

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 3 O 8 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 238 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. 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.

SOLUTION METHODOLOGY AND PROBLEM FORMULATION Formal Solution Procedure
The Generalized Integral Transform Technique (GITT) (Cotta 1990, 1994, 1998, Cotta & Mikhailov 1997, 2006, Cotta et al. 2016, 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: x ∈ V, t > 0, k, = 1, 2, ..., m with initial and boundary conditions given respectively by where the spatial operator involves both dispersion and linear dissipation: and n 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 u(x, t, T ) ≡ 0, P ≡ P(x, t), and φ ≡ φ(x, 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: with boundary conditions where the eigenvalues, μ ki , and related eigenfunctions, ψ ki (x) 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: where the symmetric kernelsψ ki (x) are given bỹ The integral transform of Eq. (1a) is then accomplished through the operator vψki (x)___dv, which after employing boundary conditions, Eqs. (1c) and (2b), yields the following transformed ordinary differential system: The initial conditions, Eq. (1b), are also integral transformed through the Equations (4) form an infinite ODE system of coupled nonlinear equations for the transformed potentials,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 (x, 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): where C(h)=dθ/dh presents the specific soil water content capacity as a function of the pressure head, θ 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: Boundary conditions (5c,d) are written in a general form so as to recover prescribed head conditions (β = 0), prescribed flux conditions (γ = 0), or free drainage conditions (γ = 0 and 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: where S e is effective saturation, θ r and θ 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: 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: subject to the filtered initial and boundary conditions: 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: 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, φ * (z, t).
Following the formalism in the GITT approach, the following fairly simple eigenvalue problem is chosen: where the eigenvalues μ i and respective eigenfunctions ψ i (z) are known in analytical explicit form. Equations (9a-c) allow definition of the following integral transform-inverse pair: where the normalized eigenfunctionsψ i (z) and the norms are given by: The filtered partial differential equation, Eq. (8a), is then operated on with L 0ψ i (z)___dz, thus eliminating the space variable and yielding the following transformed ODE system in the time variable, t: whereḡ An Acad Bras Cienc (2020) 92(1) e20190427 9 | 27 and G(z, t) The integral that defines the non-transformable nonlinear coupling term,ḡ 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) 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: The semi-analytical integration procedure now only requires analytical evaluation of two terms, as follows:ḡ Initial condition (8c) is similarly integral transformed through the operator L 0ψ i (z)___dz yielding: 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,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 238 U radioactive chain (Pontedeiro et al. 2010). The elements considered most relevant in the chain, in light of their individual half-lives, were 238 U, 234 U, 230 Th, 226 Ra and 210 Pb. The dual-porosity transport problem formulation is given as: with initial and boundary conditions: Again, filters are employed to make the boundary conditions homogeneous: which causes the filtered equations to be written as: in which the source terms are given by: The initial and boundary conditions are rewritten as: The integral transform method then starts with a proposed eigenvalue problem in the simple form: where the eigenvalues μ i and corresponding eigenfunctions ψ 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: Immobile: where the normalized eigenfunctionsψ i (z) and the normalization integrals are given by: Applying the operator L 0ψ i (z) ___dz on Eq. (18a,b), the transformed ODE system is obtained for the transformed potentials,C * m,i (t) andC * im,i (t): with transformed initial conditions with transformed initial conditions 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: where α=1.611x10 6 , θ s =0.287, θ r =0.075, β=3.96, K s =0.00944 cm/s, A=1.175x10 6 and γ=4.74, with initial and boundary conditions given by: 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.

Transport Problem
Verification of the transport problem solution was performed first for the physical non-equilibrium model of the 238 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 238 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. 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. Next, the hybrid GITT solution was verified against the well-established finite element Hydrus 1D code (Simunek et al. 2005(Simunek et al. , 2016, with special emphasis on radionuclide concentrations of the radioactive chain at the end of the solid waste pile. Figure 5 compares 238 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.

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.
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.  (Simunek et al. 2005) for parameters (5), (8), (9), (11). Experiments and direct measurements for the remaining parameters. Same data as case 3, but with initial concentration of Ra equal to 10 Bq/l 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 238 U and 234 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).
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.
Figures 8 and 9 present concentration evolutions for each element of the 238 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).
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 ( 238 U and 234 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.

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.