An Acad Bras Cienc (2020) 92(1): e20190427 DOI 10.1590/0001-3765202020190427
Anais da Academia Brasileira de Ciências | Annals of the Brazilian Academy of Sciences
Printed ISSN 0001-3765 I Online ISSN 1678-2690
www.scielo.br/aabc | www.fb.com/aabcjournal

ENGINEERING SCIENCES

Integral transform analysis of radionuclidetransport in variably saturated media using aphysical non-equilibrium model: application tosolid waste leaching at a uranium mining
installation
RENATO M. COTTA, CAROLINA P. NAVEIRA-COTTA, MARTINUS T. VAN GENUCHTEN,JIAN SU & JOÃO N.N. QUARESMA
Abstract: The Generalized Integral Transform Technique (GITT) was employed to simulatetransient 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 granularaquifer) 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 andunsaturated soil layers. The computer simulations were verified in great detail against
results obtained using the built-in routine NDSolve of the Mathematica platform andthe HYDRUS-1D software system. The present work reports the main results for 238 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 thevarious 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 initialconcentrations 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 technologicalexpertise 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 ofore is often only between 0.1% and 0.2%, which means that large amounts of ore have to be mined to
An Acad Bras Cienc (2020) 92(1)RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

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 thesame 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 ifthey contain certain minerals such as pyrite (FeS2 ), which then leads to continuous leaching of theore. 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 sterileore 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 withinformation 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 simulationsof 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 inenvironmental 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 ofthe migration of radioactive contaminants from liquid waste ponds at the same uranium mining facilityof INB (Orlande et al. 2006, Naveira-Cotta et al. 2013). The solution procedure was selected because ofits 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-porositytype) media. The solutions were verified against results obtained using the HYDRUS-1D softwaresystem (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,
An Acad Bras Cienc (2020) 92(1) e20190427 2 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

and multiple solutes in variably saturated media. The code solves the Richards equation forsaturated–unsaturated water flow and advection–dispersion type equations for heat and solutetransport.
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 channelof 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 initialconditions 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 andgeochemical properties of the soils and waste as obtained from batch and column experiments
for uranium dissolution, including validations against theoretical geochemistry models andmeasurements 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.08x106 tons)
originating from the top cover of the mine (approximately five meters deep), the leachate ore (onthe order of 2.11x106 tons, with granulometry below 10 mm) originating from piles assembled on theleaching platform, and sterile rocks (estimated at 11.13 x 106 tons). Figure 1a offers a detailed viewof 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 soilon 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 leachateore 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 initialconditions for the dissolved radionuclides. This washing with neutral water was interrupted when
production control identified a U3 O8 content of less than 20 ppm and a pH larger than 3 in the
An Acad Bras Cienc (2020) 92(1) e20190427 3 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

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

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

effluent aqueous solution. Only then were the piles removed from the leaching platform to thesolid 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 insolution.
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 ofdesorption and dissolution of contaminants during this percolation will cause an increase inradionuclide concentrations up to the foot of the pile, followed by dilution of the effluent withwater 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 hydraulicgradients 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).
An Acad Bras Cienc (2020) 92(1) e20190427 4 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Figure 2 illustrates the main steps in the simulation of radionuclide migration from the leachedore 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, qM , which percolates through the leached ore
pile, reaching the bottom drain with time-varying radionuclide concentrations. This flow mixes withwater percolating through the sterile rock with negligible contamination, which also flows into thepre-basin, thus causing a dilution of the contaminants but increasing the overall vertical infiltration
rate (qV ) 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 ofqH1 to a higher value of qH2 . The effects of advection, dispersion, retardation, radioactive decay andphysical 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 for238
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 aradiological impact assessment. The hybrid numerical-analytical GITT approach offers an accurate
and robust combination of analytical eigenfunction expansions, which eliminates the spatial variables
An Acad Bras Cienc (2020) 92(1) e20190427 5 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

through integral transforms, and a numerical methodology with automatic error control adopted forsolution of the resulting transformed ordinary differential system in the time variable.
As illustration of the formal integral transformation process, consider a transientadvection-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:∂ Tk (x, t)+ u(x, t, T` ).∇Tk (x, t) + Lk Tk (x, t) = Pk (x,t,T` ),∂tx ∈ V, t > 0, k,` = 1, 2, ..., mwk (x )

(1a)with initial and boundary conditions given respectively by
Tk (x, 0) = fk (x ), x ∈ V∂
αk (x ) + βk (x )Kk (x )T (x, t) = ϕk (x, t, T` ), x ∈ S, t > 0
∂n k

(1b)
(1c)

where the spatial operator involves both dispersion and linear dissipation:
Lk ≡ –∇.Kk (x )∇ + dk (x )(1d)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 areincorporated into the source terms. In any case, the nonlinear advection and source terms makethe 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 ofproblems, 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 suchthat 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 toan 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 situationthat admits exact solutions using classical integral transformation techniques, naturally leads toeigenvalue 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 ofproblem (1). Thus, the set of recommended auxiliary problems is given by:
Lk ψki (x ) = μ2ki wk (x )ψki (x ), x ∈ V

(2a)



∂
αk (x ) + βk (x )Kk (x )
ψ (x ) = 0 , x ∈ S∂ n ki(2b)with boundary conditions

An Acad Bras Cienc (2020) 92(1) e20190427 6 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

where the eigenvalues, μki , and related eigenfunctions, ψki (x ) are assumed to be known from analyticalexpressions 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:ZT̄k,i (t) = wk (x )ψ̃ki (x )Tk (x, t)dv, transforms(3a)v

Tk (x, t) =∞Xψ̃ki (x )T̄k,i (t), inverses

(3b)

i=1

where the symmetric kernels ψ̃ki (x ) are given byψ (x )
ψ̃ki (x ) = ki
1/2
Nki
Z
Nki = wk (x )ψ2ki (x )dv

(3c)
(3d)

v

R
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:dT̄k,i (t)
dt

+

∞
X

akij (t, T` )T̄k,j (t) = ḡki (t, T` ), i = 1, 2, ..., t > 0, k, ` = 1, 2, ..., m

(4a)

j=1

R
The initial conditions, Eq. (1b), are also integral transformed through the operator v wk (x )ψ̃ki (x )___dv
to yield:ZT̄k,i (0) = f̄ki ≡ wk (x )ψ̃ki (x )fk (x )dv
(4b)
v

where
#
∂ Tk (x, t)∂ ψ̃ki (x )ḡki (t, T` ) = ψ̃ki (x )Pk (x, t, T` )dv+ Kk (x ) ψ̃ki (x )– Tk (x, t)
ds
∂n
∂n
v
S

(4c)

akij (t, T` ) = δij μ2ki + a∗kij (t, T` )

(4d)

Z

Z

"

with
δij =



 0, for i 6= j
 1, for i = jZ∗akij (t, T` ) = ψ̃ki (x )[u(x, t, T` ).∇ψ̃ki (x )]dv
v

(4e)
(4f)

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 therequired 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 Tk (x, t) in explicit
form.
An Acad Bras Cienc (2020) 92(1) e20190427 7 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORTFLOW 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):∂ h(z, t)
∂
∂h
∂h
C (h)
=
K (h)
– D(h) ,0 ≤ z ≤ L,
t>0
(5a)
∂t
∂z
∂z
∂z
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:h(z, 0) = f (z)γ(z) h(z, t) + β(z) K (h) (

(5b)

∂h
– 1) = φ(z, t),∂z

at z = 0 and z = L(5c,d)

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φ(z, t) = –β 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:

θs – θr1
θ(h) = θr + with m = 1 – ,n mn
1 + (αh)hin 2
K (h) = Ks Sle 1 – (1 – S1e/m )Se =(n > 1)(6a)(6b)

θ(h) – θr= (1 + (αh)n )mθs – θr(6c)where Se is effective saturation, θr and θs are the residual and saturated water contents, respectively,
Ks is the saturated hydraulic conductivity, m and n are quasi-empirical shape parameters, and l is apore 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:∂ h∗ (z, t)∂∂ h∗∂ h∗ (z, t)
C (h)
=
K (h)
– D(h)+ P(z, t) ,0 < z < L,
t>0
(8a)
∂t
∂z
∂z
∂z
where


∂ F (z, t)∂∂F∂ F (z, t)P(z, t) = –C (h)+K (h)– D(h)
,
∂t
∂z
∂z
∂z

0 < z < L,

An Acad Bras Cienc (2020) 92(1) e20190427 8 | 27

t>0

(8b)

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

subject to the filtered initial and boundary conditions:h∗ (z, 0) = f (z) – F (z, 0)γ(z) h∗ (z, t) + β(z) K ∗ (z)

(8c)

∂ h∗

= 0, at z = 0 and z = L(8d,e)∂zwhere K ∗ is a representative parameter of the boundary flux, for instance the linearized hydraulicconductivity at the two boundary positions, and where the nonlinear source terms to be eliminatedby the filter are given by:∂h+ β(z)K (h)(8f)∂zWhen 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, theproposed 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 ischosen:φ∗ (z, t) = φ(z, t) + β(z)[K ∗ (z) – K (h)]d2 ψi (z)+ μ2i ψi (z) = 0 ,
dz2

0<z<L

(9a)

dψi
= 0, z = 0 and z = L(9b,c)
dz
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:γ(z)ψi (z) + β(z) K ∗ (z)Transform:h̄∗i (t) =ZLψ̃i (z)h∗ (z, t) dz

(10a)

0

Inverse:∗h (z, t) =

∞
X
i=0

ψ̃i (z) h̄∗i (t)(10b)

where the normalized eigenfunctions ψ̃i (z) and the norms are given by:ψ̃i (z) =ψi (z)

M1i /2; Mi =ZL[ψi (z)]2 dz(10c,d)0

The filtered partial differential equation, Eq. (8a), is then operated on withRL0

ψ̃i (z)___dz, thus

eliminating the space variable and yielding the following transformed ODE system in the time variable,t:dh̄∗i (t)dt

= ḡi (t),t > 0, i = 1, 2, ...

(11a)

where
ḡi (t) =ZL

ψ̃i (z) G(z, t) dz

0

An Acad Bras Cienc (2020) 92(1) e20190427 9 | 27

(11b)

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

andG(z, t) =1 ∂∂ h∗ (z, t)
D(h) ∂ h∗ (z, t)P(z, t)[K (h)]–+
C (h) ∂ z
∂z
C (h)
∂z
C (h)

(11c)

The integral that defines the non-transformable nonlinear coupling term, ḡi (t), can involve
significant computational effort for large truncation orders, due to the oscillatory nature of theeigenfunctions. 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 intransforming nonlinear source terms. This is done by dividing the space variable domain in Nk
subintervals, and adopting a first order polynomial approximation for the non-oscillatory part of theintegrand, G(z, t), at each subinterval, such as:
G(z, t) ∼= ak z + bk ,zk–1 ≤ z ≤ zk

G(zk , t) = Gk = ak zk + bk and G(zk–1 , t) = Gk–1 = ak zk–1 + bk
ak =(12a-e)Gk – Gk–1
and bk = Gk – ak zk
zk – zk–1

Thus,

zk

ḡi (t) =

Nk Z
X

zk

Nk Z
X

ψ̃i (z) G(z, t) dz ∼=k=1zk–1

ψ̃i (z) (ak z + bk ) dz

(13)

k=1zk–1

The semi-analytical integration procedure now only requires analytical evaluation of two terms,
as follows:




zk
zk
Nk
Z
Z
X 



ḡi (t) ∼ak z ψ̃i (z) dz + bk ψ̃i (z) dz
(14)
=
k=1

zk–1

zk–1

Initial condition (8c) is similarly integral transformed through the operator

RL0h̄∗i (0) =ZL0ψ̃i (z) f (z)dz –ZL

ψ̃i (z) F (z, 0)dz

ψ̃i (z)___dz yielding:

(15)

0

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 beensolved numerically, the inverse formula, Eq. (10b), is again recalled to reconstruct the filtered potentialh∗ (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 andimmobile 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
An Acad Bras Cienc (2020) 92(1) e20190427 10 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

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:∂∂∂ Cm (z, t)[(θm + f ρKd )Cm (z, t)] =[θm Dm– qCm (z, t)] – β[Cm (z, t) – Cim (z, t)] + φm
∂t
∂z
∂z
∂
[(θ + (1 – f )ρKd )Cim (z, t)] = β[Cm (z, t) – Cim (z, t)] + φim∂ t imwith initial and boundary conditions:Cm (z, 0) = Cm,back (z); Cim (z, 0) = Cim,back (z)∂ Cm (z, t)= qCm,z0 (t)∂zz=0∂ Cm (z, t)∂ Cm (z, t)
qCm (L, t) + θm Dm= qCm,zL (t) or
=0
∂z
∂z
z =L
z=L
qCm (0, t) – θm Dm

(16a)
(16b)

(16c,d)(16e)(16f,g)Again, filters are employed to make the boundary conditions homogeneous:∗Cm (z, t) = Cm(z, t) + Fm (z, t)

(17a)

∗
Cim (z, t) = Cim(z, t) + Fim (z, t)(17b)which causes the filtered equations to be written as:∂∂∂ C ∗ (z, t)∗∗∗∗[(θm + f ρKd )Cm(z, t)] =[θm Dm m– qCm(z, t)] – β[Cm(z, t) – Cim(z, t)] + Pm (z, t)
∂t
∂z
∂z
∂
∗
∗
∗
[(θ + (1 – f )ρKd )Cim(z, t)] = β[Cm(z, t) – Cim(z, t)] + Pim (z, t)∂ t imin which the source terms are given by:

(18a)
(18b)

∂
∂
∂ Fm (z, t)[(θm + f ρKd )Fm (z, t)] +[θm Dm– qFm (z, t)] – β[Fm (z, t) – Fim (z, t)] + φm
∂t
∂z
∂z
(18c)
∂
Pim (z, t) = – [(θim + (1 – f )ρKd )Fim (z, t)] + β[Fm (z, t) – Fim (z, t)] + φim(18d)∂t
The initial and boundary conditions are rewritten as:
Pm (z, t) = –∗∗Cm(z, 0) = Cm,b (z) – Fm (z, 0) ; Cim(z, 0) = Cim,b (z) – Fim (z, 0)(18e,f)∗qCm(0, t) – θm Dm(18g)∗ (z, t)
∂ Cm

=0
∂z
z=0
∗ (z, t)∂ C ∗ (z, t)
∂ Cm
∗
qCm
(L, t) + θm Dm m
= 0 or
=0
∂z
∂z
z=L
z=L

(18h,i)

The integral transform method then starts with a proposed eigenvalue problem in the simpleform:
d2 ψi (z)+ μ2i ψi (z) = 0 ,
0<z<L
dz2
dψ
qψi (0) – θm Dm i=0dz z=0dψdψiqψi (L) + θm Dm i
= 0 or=0dz z=Ldz z=LAn Acad Bras Cienc (2020) 92(1) e20190427 11 | 27

(19a)
(19b)
(19c,d)RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

where the eigenvalues μi and corresponding eigenfunctions ψi (z) are known in explicit analyticalform, to allow the definition of the following integral transform-inverse pairs for the mobile andimmobile phases:
Mobile:
∗
C̄m,i
(t) =

Transform:ZL∗ψ̃i (z) Cm(z, t) dz

(20a)

0
∗
Cm
(z, t) =Inverse:∞X∗ψ̃i (z) C̄m,i
(t )

(20b)

∗
ψ̃i (z) Cim(z, t) dz

(21a)

i=0

Immobile:∗C̄im,i(t ) =

Transform:

ZL
0

∗
Cim
(z, t) =Inverse:

∞
X
i=0

∗
ψ̃i (z) C̄im,i
(t )

(21b)

where the normalized eigenfunctions ψ̃i (z) and the normalization integrals are given by:ψ̃i (z) =ψi (z)

N1i /2; Ni =ZL[ψi (z)]2 dz(22a,b)0

RL

Applying the operatorψ̃i (z) ___dz on Eq. (18a,b), the transformed ODE system is obtained for the0∗ (t) and C̄ ∗ (t):transformed potentials, C̄m,iim,iMobile:∞XAm,ij (t)∗ (t )
dC̄m,j

j=1

dt

=

∞
X
j=1

∗
∗
∗
Dm,ij (t)C̄m,j(t) + β[C̄im,i(t) – C̄m,i(t)] + p̄m,i (t)(23a)

with transformed initial conditions
∗
C̄m,i
(0) =

ZL

∗
ψ̃i (z) Cm,b(z)dz –0ZL

ψ̃i (z) Fm (z, 0)dz

(23b)

0

Immobile:∞
X

Aim,ij (t)∗ (t)dC̄im,j

j=1

dt

=

∞
X
j=1

0

∗
∗
∗
Aim,ij (t)C̄im,j(t) + β[C̄m,i(t) – C̄im,i(t)] + p̄im,i (t)(23c)

with transformed initial conditions∗C̄im,i
(0) =

ZL
0

∗
ψ̃i (z) Cim,b(z)dz –ZLψ̃i (z) Fim (z, 0)dz

0

An Acad Bras Cienc (2020) 92(1) e20190427 12 | 27

(23d)

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

where
Am,ij (t) =ZL0(θm + f ρKd )ψ̃i (z)ψ̃j (z) dz ;Am,ij (t) =0

Bm,ij (t) =ZL0dψ̃ (z) dψ̃j (z)(θm Dm ) i
dz ;
dz
dz

ZL

Cm,ij (t) =

q

0

0

Dm,ij (t) = –Aij (t) – Bij (t) + Cij (t) ;p̄m,i (t) =ZL

dAm,ij (t)dtdψ̃i (z)ψ̃j (z) dzdzψ̃i (z) Pm (z, t) dz(23e-m)0

Aim,ij (t) =ZL0(θim + (1 – f )ρKd )ψ̃i (z)ψ̃j (z) dz ;Aim,ij (t) =0p̄im,i (t) =ZL

dAim,ij (t)dtψ̃i (z) Pim (z, t) dz

0Routine NDSolve from the Mathematica system (Wolfram 2017) is again employed to numericallysolve Eqs. (23a-d), after which the inverses, Eqs. (20b) and (21b), are recalled to reconstruct the filtered∗ (z, t) and C ∗ (z, t), which are then summed to the filters, F (z, t) and F (z, t), to recoverpotentials, Cmmimim
the desired original potentials, Cm (z, t) and Cim (z, t).

VERIFICATION OF THE GITT SOLUTIONSSolutions of both the flow problem and of the radionuclides chain transport problem were
implemented using the mixed symbolic-numerical platform of the Mathematica system (Wolfram2017). Before proceeding to the physical analysis relevant to the radiological impact assessment, thesesolutions 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:

θ(h) =α(θs – θr )α + |h|β+ θr ; K (h) = Ks

A
A + |h|γ(24a,b)where α=1.611x106 , θs =0.287, θr =0.075, β=3.96, K s =0.00944 cm/s, A=1.175x106 and γ=4.74, with initial andboundary conditions given by:h(z, 0) = –61.5 cm(25a)h(40 cm , t) = htop = –20.7 cm(25b)h(0, t) = hbottom = –61.5 cm(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
An Acad Bras Cienc (2020) 92(1) e20190427 13 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

system truncation order, N, and the number of subintervals used in the semi-analytical integration,
Nk . 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 ofsubintervals up to Nk = 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 influenceof the overall truncation order on convergence for increasing order in the summation. Clearly, bothparameters 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
Nk =120) with results obtained by (Celia et al. 1990). The two independent approaches showed excellent
adherence.

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

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 beneaththe 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 theend of the solid waste pile. A maximum truncation order of N = 20 produced full convergence to fivesignificant digits in the concentrations, although graphical convergence was achievable already for
N<10.
Figure 4 compares 238 U concentration, again in [Bq/l], versus time (in years) at the bottom of thesolid 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 theAn Acad Bras Cienc (2020) 92(1) e20190427 14 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Table I. Convergence analysis of the GITT eigenfunction
expansion using a varying number of subintervals forthe semi-analytical integration (results for -h(15, 360),
N = 15).
Order

Nk = 60

Nk = 120

Nk = 1801

49.66049.73949.762

240.44640.58240.617341.84341.99442.033442.43542.65142.707541.15941.37841.434

6

41.54541.79541.861740.70540.97341.044840.70540.97341.044940.37540.66540.742

1039.98240.26840.34411

40.04540.33640.4131239.67439.96940.04813

39.61839.91039.98814

39.55739.85139.9301539.44939.74039.818

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

N = 20N = 25

N = 30149.44149.30649.2362

40.10539.89439.786

3

41.46841.23141.111

441.85741.50341.323

540.57040.21440.033640.85140.43340.221

739.92239.46439.2348

39.92239.46439.234939.49238.98638.73410

39.10238.60738.361

11

39.13138.61938.3651238.72438.20437.94813

38.68138.17237.922

14

38.58838.06737.81215

38.43337.92437.6761638.43337.92437.6761738.35537.84537.59918

38.34237.83637.59419

38.31537.80837.56620

38.32737.79837.5592137.80637.5682237.78537.5482337.80237.56424

37.80237.5642537.80837.56126

37.57227

37.57228

37.57429

37.57730

37.581An Acad Bras Cienc (2020) 92(1) e20190427 15 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Table III. Convergence of the GITT solution forthe mobile concentration of 238 U [Bq/l] at theend of the solid waste pile at different times[years]: Case 1.
Order

C(L,500)

C(L,5000)

C(L,50.000)2214.713.66093.57974

205.043.6782

3.5797

6

204.713.6812

3.57978

204.643.68203.579710204.623.6822

3.579712204.613.6823

3.579714204.603.6824

3.579716

204.603.6824

3.579718

204.603.6824

3.579720

204.603.6824

3.5797

maximum steps in the discretization of the NDSolve numerical solution were restricted to L/100 andtf /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, 2016), with special emphasis on radionuclide concentrations of theradioactive chain at the end of the solid waste pile. Figure 5 compares 238 U concentrations at the endof 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 4. Comparison of the GITT solutions(solid lines) and the NDSolve numericalsolution (Wolfram 2017) (symbols) for 238 U
concentrations in [Bq/l] at the end of threelayers: the solid waste pile (black curves),
the soil below the pre-basin (blue curves),
and the granular aquifer (red curves).

An Acad Bras Cienc (2020) 92(1) e20190427 16 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Figure 5. Comparison of GITT solutions (solid black line) against Hydrus 1D simulation results (red dots) for the 238 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 orderto 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 physicalequilibrium during transport within the solid waste pile (i.e., without preferential paths duringtransport within the pile, which otherwise would have resulted in lower concentrations leaving thewaste pile as exemplified in (Pontedeiro et al. 2010)). Case 1 is also conservative by considering physicalnon-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, andthe 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 rockpiles (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 andfinal 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.
An Acad Bras Cienc (2020) 92(1) e20190427 17 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Table IV. Input data for the radionuclide transport simulations.
VERTICAL SOIL HORIZONTAL SOIL
(Transition soil) (Granular aquifer)

PARAMETERSSOLID WASTE PILE

FLOW PROBLEM

unsaturated

unsaturated

saturated

1

Layer length

L

[m]

10

1.40

1488

2

Final simulation time

tfinal[years]50.00050.00050.000Infiltration 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

3

4

Natural water content

θ0

-

0.055

0.088

0.074

5

Residual water content

θr

-

0.050

0.043

0.046

6

Saturated water content

θs

-

0.317

0.290

0.388

7

Saturated hydraulic conductivity

Ks[cm/h]

1.06

3.11

6.48

8

Shape parameterα[1/cm]

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

Physicalequilibrium

Physicalnon-equilibrium

Physicalnon-equilibrium

0.106

0.290

0.388

0.121

0.290

0.388

TRANSPORT PROBLEM
Water content
12

-

θ

(Cases 1 and 2)(Cases 3 and 4)
3

13

Apparent density

ρb[g/cm ]

1.73

1.61

1.66

14

Dispersivity

αx

[m]

1.0

0.14

148.8

15

Dispersion coefficientD[m2 /year]

αx .q/θ

αx .q/θ

αx .q/θ

3

16

Distribution coefficient (U)Kd[cm /g]

44

298

122

17

Partition coefficient
Mass transfer coefficient(Cases 1 and 3)

f

-

0.999

β

[1/h]

0

0.64
0.058
0.079

0.64
4.6x10-4
5.2x10-4

-

-

f.θ

f.θ

18
19

Mobile water content

θm
2

20

Dispersion coefficient - MobileDm[m /year]-D.θ/θmD.θ/θm21

Dissolution constantK[1/year]

47.8

0

0

22

Equilibrium concentration

C eq[µg/dm3 ]

286.6

0

0

3

16,430

0

0

3

0

0

0

23

Initial condition - Mobile

C m (z,0) [µg/dm ]24

Initial condition - Immobile

C im (z,0) [µg/dm ]25

Boundary condition (z=0)

C m (0,t) [µg/dm3 ]26

Boundary condition (z=L)

-

-

2

nd

0

C 0v (t)C 0h (t)

kind

nd

nd

2

kind

2

kind

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

An Acad Bras Cienc (2020) 92(1) e20190427 18 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

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

5Same data as case 3, but with initial concentration of Ra equal to 10 Bq/l

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

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 highuranium 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 wherereactions 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
An Acad Bras Cienc (2020) 92(1) e20190427 19 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

to the other radionuclides, in particular within the first 5000 years. At later times, with the attainmentof non-equilibrium conditions, the uraniums remain with a radioactivity of one order of magnitudelarger 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, withthe 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 verticalsoil, and the granular aquifer). Results here are plotted for a much shorter total time to more clearlyvisualize 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 andPb, 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 chainat 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 thesame 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 fourcases become more pronounced. In particular, notice that the two cases assuming physical equilibriumin the two soil layers (cases 2 and 4), especially the granular aquifer, produce a significant delay inthe 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 tocases 1 and 3, as previously pointed out, the larger infiltration rate considered for case 3 produced amuch faster transport rate of the contaminants, albeit with lower maximum concentration values dueto the more pronounced dilution within the pre-basin. While for case 1 the maximum concentrationof the two uraniums occur at approximately 4000 years at values of about 34.6 Bq/l, for case 3 themaximum concentration is reached after approximately 2900 years but with a lower concentration
An Acad Bras Cienc (2020) 92(1) e20190427 20 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Figure 7. Case 1 concentrations versus time at the endof the three layers: (a) of 238 U; (b) of 234 U; (c) of 230 Th;(d) of 226 Ra; (e) of 210 Pb.

An Acad Bras Cienc (2020) 92(1) e20190427 21 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Figure 8. Concentrations versus time at the end ofwaste pile for cases 1 to 4: (a) of 238 U; (b) of 234 U; (c) of
230
Th; (d) of 226 Ra; (e) of 210 Pb.

An Acad Bras Cienc (2020) 92(1) e20190427 22 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

Figure 9. Concentrations versus time at the end of granular aquifer for cases 1 to 4: (a) of 238 U; (b) of 234 U; (c) of
230
Th; (d) of 226 Ra; (e) of 210 Pb; (f) also shows the total concentrations at the end of the granular aquifer.

An Acad Bras Cienc (2020) 92(1) e20190427 23 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

value of 33.4 Bq/l. For case 2, the maximum concentration was 24.3 Bq/l after 5700 years, while forcase 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 Raconcentration 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 forthe 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 transportof the other radionuclides. Results indicate the increased importance of the two radionuclides (Raand Pb), with a relatively fast decay of Ra leading to even larger activities of Pb. Please note thatat 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 arestill higher. Again, for large times within the equilibrium situation, the activities of the two uraniums
predominated.

Figure 10. Concentration evolution at the end of thewaste 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.

An Acad Bras Cienc (2020) 92(1) e20190427 24 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORTCONCLUSIONSThis 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 aradionuclide 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 granularaquifer, 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 numericalsolutions, 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.
Acknowledgments
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

COTTA RM. 1993. Integral transforms in computational heat andfluid flow. Boca Raton, FL: CRC Press.
COTTA RM. 1994. Benchmark results in computational heat and

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 forthree-dimensional steady turbulent dispersion in rivers andchannels. 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 tononlinear diffusion problems. Numer Heat Transfer, Part B 17:
217-226.

fluid flow: The integral transform method.
Transfer, Invited Paper 37: 381-394.

Int J Heat Mass

COTTA RM. 1998. The integral transform method in thermal andfluids sciences and engineering. New York, NY: Begell House.
COTTA RM, KNUPP DC & NAVEIRA-COTTA CP. 2016. Analytical heatand 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.

An Acad Bras Cienc (2020) 92(1) e20190427 25 | 27

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

COTTA RM & MIKHAILOV MD. 1997. Heat conduction: Lumped

LEAL MA & RUPERTI JR NJ. 2001. A hybrid solution for simulation

analysis, integral transforms, symbolic computation. New York,
NY: Wiley Interscience.

of 2-D contaminant transport in heterogeneous subsurface
systems. Hybrid Meth Eng 3(2): 111-129.

COTTA RM & MIKHAILOV MD. 2005. Semi-analytical evaluationof integrals for the generalized integral transform technique.

LIU C, SZECSODY JE, ZACHARA JM & BALL WP. 2000. Use of thegeneralized integral transform method for solving equations of

In: Proc. of the 4th Workshop on Integral Transforms andBenchmark Problems – IV WIT, CNEN, Rio de Janeiro, Brasil.

solute transport in porous media.
483-492.

COTTA RM & MIKHAILOV MD. 2006. Hybrid methods and symbolic
computations. In: Minkowycz WJ, Sparrow EM & Murthy JY (Eds),

MIKHAILOV MD & OZISIK MN. 1984.
Unified analysis andsolutions of heat and mass diffusion. New York, NY; also, Dover

Handbook of Numerical Heat Transfer, New York, NY: John Wiley,
2nd ed.

Publications, New York, NY, 1993: John Wiley.

COTTA RM, MIKHAILOV MD & RUPERTI JR NJ. 1998. Analysis ofradioactive waste contamination in soils: Solution via symbolic
manipulation. In: Proc. Int. Conf. on the Radiological Accident ofGoiâ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.

Contaminanttransport 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.
HEILBRON FILHO PFL, PEREZ-GUERRERO JS, PONTEDEIRO EM,
RUPERTI JR NJ & COTTA RM. 2002.
Analytical model forenvironmental impact due to solid waste disposal fromelectricity 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. Riode Janeiro, Brazil: E-Papers.
IAEA. 2003. Extent of environmental contamination by naturally
occurring radioactive material (NORM) and technologicaloptions for mitigation, Technical Report Series No. 419, Int.
Atomic Energy Agency, Vienna, Austria. http://www.pub.iaea.
org/MTCD/publications/Pubdetails.asp?pubId=6789.
IAEA. 2004. Safety assessment methodologies for near surface
disposal facilities, Vol. 1, Review and enhancement of safetyassessment approaches and tools, Int. Atomic Energy Agency,
Vienna, Austria.
LEAL MA & RUPERTI JR NJ. 2000. A numerical study for thetwo-dimensional solute transport in groundwater via integral
transform method. Hybrid Meth Eng 2(1): 111-129.

Adv Water Resources 23:

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 inheterogeneous media. Int J Heat Mass Transfer 52: 5029-5039.
NAVEIRA-COTTA CP, PONTEDEIRO EM, COTTA RM, SU J &VAN GENUCHTEN MT. 2013. Environmental impact assessmentof 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 deradionuclídeos e análise de impacto radiológico no sítio daspilhas 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 inhazardous 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 analyticalsolution 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
the advection-diffusion transport equation using

An Acad Bras Cienc (2020) 92(1) e20190427 26 | 27&
ofa

RENATO M. COTTA et al.

INTEGRAL TRANSFORM ANALYSIS RADIONUCLIDE TRANSPORT

change-of-variable and integral transform technique.

Int J

https://orcid.org/0000-0003-0965-0811

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 infinite 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 indual-porosity systems. Hybr Meth Eng 3(2&3): 1-18.
RUBOL S, BATTIATO I & BARROS FPJ. 2016. Vertical dispersion invegetated shear flows. Water Resour Res 15(10): 8066-8080.
SIMUNEKJ,

VAN

GENUCHTEN

Renato M. Cotta1,2

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 fortransport 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 irregulardomains. 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.
mining production, updated March 2019.

World uranium
http://www.

world-nuclear.org/information-library/nuclear-fuel-cycle/
mining-of-uranium/world-uranium-mining-production.aspx.
WOLFRAM S. 2017. Wolfram Mathematica 11.
Wolfram Research Inc.

Carolina P. Naveira-Cotta1
https://orcid.org/0000-0001-8073-5602

Martinus T. van Genuchten3,4
https://orcid.org/0000-0003-1654-8858

Jian Su3https://orcid.org/0000-0001-7853-1088João N.N. Quaresma5
https://orcid.org/0000-0001-9365-7498
1 LabMEMS,

Mechanical Eng. Dept., POLI & COPPE, Federal

University of Rio de Janeiro, Rua Muniz de Aragão, no. 360, Bloco
4, CT-2, Cidade Universitária, 21941-594 Rio de Janeiro, RJ, Brazil
2 General

Directorate for Nuclear and Technological
Development, DGDNTM, Brazilian Navy, Ed. 23 do AMRJ, Rua daPonte, s/n, Centro, 20091-000 Rio de Janeiro - RJ, Brazil
3 Nuclear Eng. Dept., POLI & COPPE, Federal University of Rio deJaneiro, UFRJ, Rua Horácio Macedo, Bloco G, Cidade Universitária,21941-450 Rio de Janeiro, RJ, Brazil
4 Centerfor Environmental Studies, CEA, São Paulo State
University, UNESP, R. 16 B, 99-193, Bela Vista, 13506-752 Rio Claro,
SP, Brazil
5 Chemical Eng. Dept., Federal University of Pará, UFPA, Institutode Tecnologia, Campus Guamá, 66075-110 Belém, PA, Brazil

Correspondence to: Renato Machado Cotta
E-mail: cotta@mecanica.coppe.ufrj.br

Authors contributions
RMC contributed on the physical modeling, solution
methodology, algorithm implementation and resultsinterpretation, CPNC contributed on the experimental input data
analysis, solution methodology and algorithm implementation,
MTVG contributed on the physical modeling and resultsinterpretation, JS contributed on the experimental input data
analysis and algorithm implementation, and JNNQ contributedon the solution methodology and experimental input data
analysis.

Champaign, IL:

How to cite
COTTA RM, NAVEIRA-COTTA CP, VAN GENUCHTEN MT, SU J & QUARESMA
JNN. 2020. Integral transform analysis of radionuclide transport invariably saturated media using a physical non-equilibrium model:
application to solid waste leaching at a uranium mining installation. An
Acad Bras Cienc 92: e20190427. DOI 10.1590/0001-3765202020190427.

Manuscript received on April 10, 2019; accepted for publication
on September 11, 2019

An Acad Bras Cienc (2020) 92(1) e20190427 27 | 27