Artículo en PDF
How to cite
Complete issue
More information about this article
Journal's homepage in redalyc.org
Sistema de Información Científica
Red de Revistas Científicas de América Latina y el Caribe, España y Portugal
Rev. Int. Contam. Ambie. 27(3) 215221, 2011
NUMERICAL SIMULATION OF
226
Ra MIGRATION AND DECAY IN
A SATURATED POROUS MEDIUM
Eduardo DE LA CRUZ
1
, Roberto GONZALEZ
1,2
, Jaime KLAPP
1
,
Luis Carlos LONGORIA
1
, Estela MAYORAL
1
and Ricardo DUARTE
1
Instituto Nacional de Investigaciones Nucleares, Carretera MéxicoToluca s/n, La Marquesa, Ocoyoacac
52750, Estado de México, México. Correo electrónico: eduardo.delacruz@inin.gob.mx
2
Facultad de Ciencias, Universidad Autónoma del Estado de México, El Cerrillo Piedras Blancas 50200, Estado
de México, México.
(Recibido junio 2010, aceptado mayo 2011)
Key words: Simulation, soil pollution, transport, porous medium
ABSTRACT
Numerical simulation of transport and radioactive decay of
226
Ra as a constant product
of
238
U decay from a simulated storage facility containing uranium tailings is presented
in this work. The simulation domain consists of six subsoil layers through which one
of the uranium daughters
226
Ra is migrating as well as decaying. The methodology
consists in ±rst calculating the
226
Ra velocity ±eld through the computational domain
by solving the momentum and mass conservation equation using Darcy´s Law. The
transport equation is then solved to obtain the extent of the contaminant´s migration.
The solution is obtained with the ±niteelement method. Results showed that the esti
mated
226
Ra migration time through 200 meters of soil was about ±ve hundred years.
Palabras clave: simulación, contaminación de suelo, transporte, medio poroso
RESUMEN
En este trabajo se presenta una simulación numérica de la migración y del decaimiento
de
226
Ra como un producto constante del decaimiento del
238
U a través del subsuelo
de un almacén modelado que contiene jales de uranio. El dominio de simulación para
modelar el transporte de
226
Ra consiste en seis diferentes capas de tierra a través de
las cuales los productos del
226
Ra van migrando y decayendo. La metodología calcula
el campo de velocidades por medio de la ecuación de momento y de conservación de
masa usando la ley de Darcy. Posteriormente se utiliza este resultado para resolver la
ecuación de transporte. La solución es obtenida con el método de elementos ±nitos.
Los resultados muestran que el tiempo estimado de migración del
226
Ra a lo largo de
200 metros de tierra fue de alrededor de 500 años.
E. de la Cruz
et al.
216
INTRODUCTION
Soil contamination due to natural and anthropo
genic pollutants is a global problem.
Soil pollution produces physical, chemical and
biological imbalances affecting living beings.
Radioactive contaminants receive special atten
tion due to their radiotoxicity and decay time. The
decay of a radioactive element is a physical constant
which is generally deFned as the half life, which
is the time it takes for the mass fraction of a given
radionuclide to be reduced by a half.
Table I
shows the
226
Ra radioactive decay chain
and the half life of the various radionuclides involved.
This work is focused on the
226
Ra radionuclide
as a product of
228
U decay. The
226
Ra isotope is the
most important uranium daughter in soil migration
studies because it is the longest lived element and
therefore is the most common nuclide in uranium
tailings.
One of the greatest risks of soil contamination
is that solid pollutants reach underground water.
As an example, it is known that one Curie of
90
Sr
dissolved in water could affect about 1000 liters of
water making it unacceptable for public drinking
(Stanley 1966).
Numerical simulation is a very powerful tool
to understand the solute migration phenomena as
it is necessary to predict processes for hundreds of
years. This work is concerned with the modeling of
a contaminant migration and decay through a porous
saturated stratum. The hydrogeological conditions
considered for the modeling correspond to those
typically found in the north of México.
A differential equation system representing this
problem taking into consideration the boundaries
and initial conditions was established. The equations
of the system were solved using the Fnite element
method. The software used for this purpose was
COMSOL.
Results show that the evolution of a
226
Ra radio
active plume exhibits a bell shape when a source of
contaminant is released in a Fnite time. A similar
behavior could be found in the literature for a chemi
cal agent dispersed in a porous media (Goltz 1986).
METHODOLOGY
Two approaches for simulating the solute migra
tion are considered; the Frst one deFnes the solid
matrix of the underground soil layers and then the
laws of hydrodynamics are applied. In this case the
simulated area is very small but large enough to
consider pores as channels. In this case, it is difFcult
to represent both, the variation of hydrodynamic
properties of the media space, and the deFnition of
the same solid matrix of the porous media.
The second approach, also called Darcy´s law,
introduces a macroscopic concept that takes into ac
count hydrodynamic properties that depends on the
morphology and constitution of the medium. The
±ow of the ±uid is considered to be a continuous
process introducing an average value of the physical
properties that characterize the porous media. These
soil properties are considered into the Darcy´s equa
tion and other ±ow equations giving a very good
estimate of the phenomenon over small and large
dimensions.
The simulated diffusion of the contaminant dis
solved in water ±owing through the underground
is based on the velocity Feld as the ±uid migrates
through the stratum.
Transport and diffusion of the contaminant in the
subsoil is a three step process: 1) Molecular diffusion
and hydrodynamic scattering are considered together
as a hydrodynamic dispersion. Molecular diffusion
is determined by the contaminant concentration
gradients and considers longitudinal and transverse
±ows (Sudicky 1989), 2) Transport due to “drag” of
contaminant particles in underground water, and 3)
Radioactive decay.
TABLE I.
HAL²–LI²E O²
226
Ra AND ITS DAUGHTERS
(KAPL 1996)
Radionuclide
Halflife
timescale
Nuclear
decay product
226
Ra
1599 yrs
222
Rn
222
Rn
3.8235 d
218
Po
218
Po
3.10 min
214
Pb
218
At
218
At
1.5 s
214
Bi
218
Rn
218
Rn
35 ms
214
Po
214
Pb
26.8 min
214
Bi
214
Bi
19.9 min
214
Po
210
Tl
214
Po
0.1643 ms
210
Pb
210
Tl
1.30 min
210
Pb
210
Pb
22.3 yrs
210
Bi
210
Bi
5.013 d
210
Po
206
Tl
210
Po
138.376 d
206
Pb
206
Ti
4.199 min
206
Pb
206
Pb
Stable

NUMERICAL SIMULATION OF
226
Ra MIGRATION AND DECAY
217
The groundwater fow and solute transport equa
tions are coupled to Darcy´s di±±usion ±ormalism
which assumes that the di±±usion velocity is
u
= –
k
∇
H
,
where
k
is the permeability and
H
the hydraulic head.
The velocity
u
is the fuid fow velocity ±rom an in
²nitesimal sur±ace that represents both the solid and
the pore spaces. We ²rst calculate the Darcy velocity
²eld
u
=(
u
,
w
), where
u
and
w
are the velocities in
the x and z directions, respectively. Then, assuming
that the contaminant ±ollows the previously obtained
velocity ²eld, the contaminant transport equation
coupled to the nuclear isotope decay equation ±or the
226
Ra chain is solved.
According to Darcy, the fuid speci²c discharge
rate through a porous material is proportional to the
hydraulic head loss and is given by the equation
(Fetter 1994):
∇
•
(
K
∇
H
) +
R
= 0,
(1)
where R is the water fow rate [m/s].
Solute transport typically is a timedependent
equation, and is described by the equation
θ
s
(∂C/∂t) +
∇
∙(θs D
L
∇
C + uC) – λC = 0,
(2)
where θ
S
is the porosity,
D
L
the hydrodynamic disper
sion tensor [m³/s],
C
the dissolved concentration [kg/
mµ],
u
the Darcy´s speed [m/d],
t
the time [s], and λ
the radioactive decay rate [s
–1
].
The dispersion represents the solute spreading by
mechanical mixing and molecular di±±usion. The equa
tions ±or the tensor entries are (Zheng and Wang 1998)
θ
D
ii
=
α
L
(
u
i
2
/
u
) +
α
T
(
u
j
2
/
u
) +
D
m
,
(3)
θ
D
ii
= θ
D
ii
= (
α
L
–
α
T
)[(
u
i
u
j
)/
u
)]+
D
m
,
where θ
D
ii
are the principal components o± the dis
persion tensor based on the Darcy´s velocity, θ
D
ij
and θ
D
ji
are the cross terms o± the dispersion tensor,
the subscript “L” denotes longitudinal dispersivity,
“T” the transverse dispersivity,
u
= (
u
i
2
+
u
j
2
)
1/2
is
the magnitude o± the Darcy´s velocity vector, and
D
m
represents the e±±ective molecular di±±usion in
a saturated porous media.
The
226
Ra decay chain is shown in
Table I
. The
disintegration rate is proportional to the number o±
nuclei present in the sample:
dN/d
t
= – λN,
(4)
where N is the number o± nuclei o± the radioactive
species present at time
t
, and, λ
226
= ln2/1599 (yrs
–1
)
±or the
226
Ra →
222
Rn decay.
In this work is assumed that
226
Ra decays directly
to the stable isotope
206
Pb, because all hal±li±e tim
escales ±rom
222
Rn to
206
Ti are much shorter than ±or
226
Ra (see
Table I
).
THE INPUT MODEL
Figure 1
shows the con²guration o± the simulated
system having hydrogeological conditions similar to
those ±ound in the north o± México. The dimensions
o± the simulated domain are 80 m wide and a depth
o± 200 m.
The domain is constituted o± six horizontal layers
with di±±erent characteristics. The top layer is denoted
as K_CO5 which corresponds to the sur±ace soil
layer. This is the layer where the
226
Ra is entering
into the system and then ²lters to groundwater layers.
The infow velocity o± 334.4 mm/year represents the
rain±all per unit time in this site.
For the present simulation a saturated medium is
considered; the fuid is a combination o± water and
dissolved contaminant, assumed not reactive and
homogeneous.
The lithology o± the system was characterized by
di±±erent porosities ±or each layer within a range o±
0.20 to 0.35, permeability between 10
–9
to 10
–6
cm
2
,
a molecular di±±usion coe±²cient o± 1.34
x
10
–9
m
2
/s,
a transverse dispersivity o± 5 mm and a longitudinal
dispersivity equal to 0.5 m as shown in
Fig. 1
and
table II
.
Fig.1.
The six layers considered in the computational system.
The permeability coe±²cient o± each layer is denoted by
the constant k.
1500 m
1495 m
1485 m
1400 m
Boundary
Zero Fluz
Boundary
Zero Fluz
1310 m
1300 m
110 m
190 m
outflow boundary condition
K_CO1=10e–9 cm
2
K_CO2=10e–6 cm
2
K_CO3=10e–3 cm
2
K_CO6=10e–7 cm
2
K_CO4=10e–7 cm
2
K_CO5=10e–6 cm
2
Recharge 334 mm/year
E. de la Cruz
et al.
218
The left and right boundaries along the zdirection
are located at x=110 m and x=190 m, respectively,
where we impose the following Neumann zero Fow
condition:
n
•
(
k
∇
H
) = 0,
(5)
where
n
is a unit vector normal to the boundary. This
boundary condition implies that at the boundary the
velocity is parallel to the boundary and that the x
component of the velocity is zero, i.e. there is no Fuid
Fow across the left and right boundaries.
±rom the top boundary (z=1500 m) water Fows
into the system at a rate of 334 mm/year; an inFow
boundary condition that satis²es the Neumann condi
tion is assumed and is written in the form
n
•
(
k
∇
H
) = R.
(6)
The top boundary condition for the hydraulic head
corresponds to the Dirichlet condition
H
(
z
=1500,
t
)=
H
0
,
(7)
where
H
0
is the barometric hydraulic head.
In the bottom boundary (
z
=1300 m), an outFow
boundary condition was imposed, assuming the de
rivatives of the variables that characterize the Fuid
are zero at the boundary.
±or the transport equation, the initial and bound
ary conditions are given by Eqs. (8) to (11). In this
work it is also assumed that the contaminant source
is continuous over a period of ²fty years and zero
after this time. The solute transport condition for the
top boundary is
–
n
•
(–θ
sD
L
∇
C
+
uC
)=
C
0
,
(8)
where C
0
is the concentration of
226
Ra that enters the
system and has two values along of simulation: a)
C
0
= 1
x
10
–6
kg/m
3
for 0<
t
<50 years and
C
0
=0 for
t
>
50 years.
C
0
indicates that a small amount of Fuid
(water + radionuclide) is injected to the system with
a hydraulic recharge rate of 334.4 mm/year (product
rain) from the time
t
=0 to
t
= 50 years, i.e. the injec
tion of
226
Ra only occurs during the ²rst ²fty years,
and afterwards the source is removed.
Equation (8) is applied only in the border segment
140<
x
<160 meters (see
Fig. 1
); for the rest of the
border the Neumann condition C(x,h,t) = 0 is applied.
The Dirichlet condition on the right and left
boundaries is
C
(0,
y
,
t
) = 0.
(9)
The Neumann condition de²nes the zero Fow
boundaries:
–
n
• (–θ
sD
L
∇
C
) = 0.
(10)
±inally, the initial condition indicating that the
subsoil on the site is free of contaminants at the
beginning of the simulation is
C
(
x
,
y
, 0) = 0.
(11)
The initial concentration is zero elsewhere.
RESULTS AND DISCUSSION
Diffusion simulations of a contaminant dissolved
in water Fowing through the underground implies at
least two types of physical analysis. ±or this case the
hydro geological context is ²rst solved, then using
the velocity ²eld, the behavior of the Fuid through
the stratum is simulated.
The methodology consists in ²rst calculating
the
226
Ra velocity ²eld through the computational
domain by solving the momentum and mass conser
vation equation using Darcy´s Law. Then we solve
the contaminant transport equation by the procedure
shown in
fgure 2
, see below.
Figure 3
shows the pressure distribution and the
velocity ²eld in the system at two simulation times.
Pressures are given in [Pa] showing its variations by a
color table. Red corresponds to the maximum pressure.
±or the bottom panel of
fgure 3
the pressure reaches
6.651 [Pa]. The vectors represent the Fow velocity ²eld
due to the pressure ²eld. The magnitude is proportional
to the length of the vector; the direction and sense are
clearly indicated by the vector orientation.
It was also found that pressure and velocity
achieves a steady state after a lapse of about two
thousand years.
TABLE II.
TYPE O± SOILS AND GEOLOGIC CHARACTE
RISTICS
Strate
Composition
Porosity
CO1
Sand conglomerate
0.20
CO2
Limestone and sand
0.30
CO3
Sand gravel and limestone mixed
0.35
CO4
Sandstone
0.25
CO5
Limestone and sand
0.30
CO6
Sandstone
0.25
NUMERICAL SIMULATION OF
226
Ra MIGRATION AND DECAY
219
The top and bottom panels of
fgure 3
show a
similar pattern for the velocity and pressure distri
butions at two different times. It can be seen that
the ±uid lines follows a pattern due to the pressure
distribution and that the pressure depends on the
hydrogeological characteristics of each layer; for
example, the ±uid velocity is faster in the K_CO3
layer due to the hydraulic high permeability as
signed to this layer. This causes the pressure to
increase on the right side of the system, speci²cally
in the K_CO6 stratum.
Figure 4
shows the evolution of the
226
Ra plume
as a function of time. The contaminant (
226
Ra) enters
through a segment of the upper boundary of stra
tum K_CO5 for 140<
x
<160 meters, and
t
=0 to
t
=50
years. It can be seen that the
226
Ra is seeping through
underground layers.
In the same
fgure 4
we show the amount of con
taminant concentration in each point. Colors indicate
the
226
Ra concentration degree. Red indicates a more
intense concentration.
The concentration varies due to advection, diffu
sion, and additionally, suffers the effects of radioac
tive decay along its path.
Figure 4
(from top to bottom) shows that the ±ow
is going towards the left boundary, but it is reversed
as a result of the velocity boundary condition which
requires that the xcomponent of the velocity at the
boundary be zero.
The remaining evolution of the contaminant is
shown in the lower panels of
fgure 4
; the ±uid “fol
lows” jagged lines and widens in the xdirection.
Figure 5
shows the evolution of the contaminant
concentration versus time at x=165, z=1300, which
is located in the lower boundary of the system. It can
be seen that the concentration of the contaminant at
this point increases gradually with time until reach
ing a maximum at approximately 2.3
x
10
10
seconds
(729 years) and then decreases exponentially. The
enlargement of the plot (bottom panel of
Fig. 5
),
allows us to obtain a particular time estimated when
the contaminant concentration reaches a signi²cant
value in the lower boundary of the system. The esti
mated value is 1.6
x
10
10
seconds, which is equivalent
to 507 years.
For the ²rst ²fty years the
226
Ra concentration
keeps increasing because this time is much shorter that
the
226
Ra halflife timescale. After this time the
226
Ra
composition decays to zero. As mentioned before,
all halflife timescales from
222
Rn to
206
Ti are much
shorter than for
226
Ra (see
Table I
), and thus in the
timescale of
fgure 5,
the ²nal product is the stable
isotope
206
Pb.
Begin
Initial
conditions
(Darcy´s Law)
Initial conditions
(transport equation and
radiactive decay)
Compute
of
velocity
field
Results
pollutant
concentrations
END
NO
YES
Is correct?
Fig. 2.
Problem resolution diagram
Time=8.64e10
Surface: Pressure [Pa]
Arrow: Velocity field
Max: 6.651
1500
1400
1350
1450
1350
100 110
120 130 140 150 160
170 180 190
200
Min: –4.791e–16
6
5
4
3
2
1
0
1500
1400
1350
1450
1350
100 110
120 130 140 150 160
170 180 190
200
Min: –0.0277
Time=2.009059e10
Surface: Pressure [Pa]
Arrow: Velocity field
Max: 3.076
3
2.5
2
1.5
0.5
0
1
Fig. 3.
Pressure and velocity vector ²eld at times
t
= 637 yrs (left panel) and
t
= 2739 yrs (right panel). Time in
the ²gure is given in seconds
E. de la Cruz
et al.
220
Fig. 5.
Evolution of the contaminant concentration as a function of time [s] for the Fx point in space (x, z)=
(165, 1300)m. The right panel is an ampliFcation of the left panel
0.05
0.045
0.04
0.035
0.03
0.025
0.02
Concentration, c [kg/m
3
]
0.015
0.01
0.005
0
–0.005
0
0.5
1
1.5
2
2.5
3
Time
x10
10
3.5
4
4.5
6
5.5
1.1
1.2
1.3
Time
1.4
1.5
1.6
x10
–3
9
8
7
6
5
4
3
2
1
0
Concentration, c[kg/m
3
]
CONCLUSIONS
The migration and decay of a contaminant in
soil of the type found in the North of México has
been simulated, allowing a better understanding of
its underground behavior. The time required by the
contaminant to reach the bottom of the lower stratum
was estimated to be about Fve hundred years.
The particular ±ow obtained in the numerical
simulation is due to the soil characteristics of the
layers and the pressure distribution in the system.
Hydrogeological computer simulation studies will
enable to assess not only the long time behavior of the
contaminants, usually impossible to recreate by experi
ment, but also to evaluate the migration routes of radio
nuclides dissolved in water through the underground.
Fig. 4.
Evolution of the contaminant concentration C from the time
t
= 0 to
t
= 318 yrs. In the Fgure time
is given in seconds.
1500
1420
1400
1340
1320
1300
1480
1460
1440
1360
1380
1500
1420
1400
1340
1320
1300
1480
1460
1440
1360
1380
110
120 130 140 150 160
170 180 190
Min: –0.0844
110
120 130 140 150 160
170 180 190
Time=6.3072e7
Surface: Concentration, c [kg/m
3
]
Time=1.702944e9
Surface: Concentration, c [kg/m
3
]
Max: 1.161
Max: 6.433e–7
x10
–7
1
0.8
0.6
0.4
0.2
0
6
5
4
3
2
0
1
–1
Min: –1.4e–7
1500
1420
1400
1340
1320
1300
1480
1460
1440
1360
1380
110
120 130 140 150 160
170 180 190
Min: –2.60e–8
Time=4.068144e9
Surface: Concentration, c [kg/m
3
]
Max: 2.611e–7
x10
–7
2.5
2
1.5
1
0.5
0
Time=1.002845e10
Surface: Concentration, c [kg/m
3
]
1500
1420
1400
1340
1320
1300
1480
1460
1440
1360
1380
110
120 130 140 150 160
170 180 190
Max: 8.618e–8
x10
–8
8
7
6
5
4
3
2
0
1
Min: –5.659e–9
NUMERICAL SIMULATION OF
226
Ra MIGRATION AND DECAY
221
The model allows the estimation of the transport
time, the dissolved contaminants concentration and
their eventual arrival in aquifers with consequent
contamination. Furthermore, the
226
Ra concentration
curves follow a pattern similar to that reported in a
±eldscale experiment of the migration of contami
nants through a porous medium.
The numerical simulation presented in this work
has captured the basic physical mechanism of the
migration and decay of a contaminant as it moves
through an underground porous medium. It is fea
sible to apply similar studies to assess the suitability
of places for containment of contaminants as well
as to identify factors that determine the amount of
runoff water that could carry the contaminant, and
characterize the hydraulic operation of water tanks
or ²uids.
REFERENCES
Fetter C.W. (1994).
Applied Hydrogeology
. PrenticeHall.
Goltz M.N. and Roberts P.V. (1986). Interpreting organic
solute transport data from a ±eld experiment using physi
cal nonequilibrium models. J. Contam. Hydrol. 1, 7793.
KAPL (1996).
Chart of the Nuclides
, Fifteenth edition.
Knolls Atomic Power Laboratory and Lockheed Martin
Company.
Stanley N.D. (1966) Hydrogeology. Wiley, New York,
464 pp
Sudicky E. A. (1989). The Laplace transform Galerkin
technique: application to mass transport in groundwa
ter. Water Resour. Res. 25, 18331846.
Zheng C. and Wang P. (1998). MT3DMS: A modular three
dimensional multispecies transport model for simula
tion of advection, dispersion and chemical reactions
of contaminants in groundwater systems. university
of Alabama, 239 p.