Using Hydrus-HP1 to Estimate 226Ra Transport in Soils Following
the Use of Phosphogypsum in Agriculture
Camila R. Bezerra Coelho1, Márcia S. Batalha1, Maria Claudia Barbosa1,
Diederik Jacques2, and Martinus Th. van Genuchten3
1Department of Civil Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, RJ, Brazil
camila.rosabz@gmail.com, marciabatalha.nunes@gmail.com, mclaudia@coc.ufrj.br
2
Performance Assessments, Institute of Environment, Health and Safety, Belgian Nuclear Research Centre,
Mol, Belgium, djacques@sckcen.be3
Department of Mechanical Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro,
RJ, Brazil, rvangenuchten@hotmail.com

Abstract
Agricultural soils generally require the use of fertilizers and soil conditioners for optimal production.
Phosphate fertilizers produced from igneous phosphate rock often contain small amounts of naturalradionuclides (most notably 238U and 232Th), while the byproduct phosphogypsum (dihydrated calcium
sulfate) is enriched in 226Ra and 210Pb. Our study focused on the long-term (200 year) subsurface fate of226
Ra when phosphogypsum (PG) is applied annually for 27 years as an amendment to a typical Brazilian
Cerrado soil profile, after which PG applications ceased but agriculture continued. Expected water flow
and solute transport processes were estimated first using the standard HYDRUS-1D software package
assuming linear equilibrium transport, followed by simulations using the HP1 module of HYDRUS-1D.
Analyses with HP1 considered a range of geochemical and physical parameters affecting 226Ra transport
in the vadose zone, especially the effects of varying calcium concentrations resulting from the PG
applications. Several aspects related to the modeling approach, such using daily versus monthly averaged
weather data and assuming long-term steady-state flow conditions, are discussed also.

1.

Introduction

Natural radioactivity is common in most or all soils and rocks. This natural radioactivity may
change when fertilizers or amendments are used in agricultural operations to improve crop
production. For example, phosphate fertilizers often contain small amounts of radionuclides fromthe uranium (238U) and thorium (232Th) decay chains. During the phosphate production process,the naturally occurring U and Th decay series in the host rock are disrupted and some of theradionuclides tend to accumulate in both the fertilizer itself (approximately 86% of 238U and 70%
of 230Th) and the by-product phosphogypsum (about 80% of 226Ra and its daughter products)(Mazzilli et al., 2000; Papastefanou et al., 2006, Saueia and Mazzilli, 2006). For every ton ofphosphoric acid being produced (the basic constituent for the production of phosphate
fertilizers), about 4.8 tons of phosphogypsum (PG) is being generated. Because of its 226Ra
content, PG in many countries is considered a waste and stockpiled in very large quantities.
In Brazil, a major agricultural production area is the Cerrado region. Cerrado soils typically arehighly weathered, slightly acidic, leached, have a relatively low cation exchange capacity, andare poor in nutrients such as P, Ca and Mg. Natural gypsum has long been used as an amendment
or soil conditioner for such situations. Since PG is rich in P, Ca and Mg, many studies have1

advocated the use of PG as both a fertilizer and soil conditioner to reduce subsoil acidity and Al
toxicity, as well as to promote sodic soil reclamation and reduce soil erosion (Oster, 1982; Oates
and Caldwell, 1988; Embrapa, 2005; among others).
Detailed analyses are needed to understand the possible accumulation of radionuclides when fertilizers and/or phosphogypsum are being applied routinely over long periods of time to the same
field (Jacques et al., 2008; Batalha et al. 2011b). For this purpose we used the HYDRUS-1D
flow and transport simulator (Šimůnek et al., 2008) to predict the long-term fate and transport of226
Ra present in PG into and through the vadose zone to groundwater, assuming linear equilibrium transport. Rather than considering one single application, which should have very little impact on the environment (e.g., Rothbaum et al., 1979; Zielinski et al., 1997), the application is assumed to be repeated once every year for 27 years to a typical soil profile in the Cerrado regionof Brazil. The transport of 226Ra in the vadose soil is subsequently followed for another 170
years.
A recent study by Jacques et al. (2008) shows that radionuclide transport in soils is affected by arange of physical and geochemical processes. Their study focused on the transport of uranium
following the long-term use of phosphate fertilizers, as opposed to radium in this study. Since
radium and calcium released from applied PG may compete for the same sorption sites, several
simulations were performed using the HP1 (Jacques and Šimůnek, 2010) multicomponent
transport module of HYDRUS-1D. The module was obtained by coupling HYDRUS-1D with
the PHREEQC biochemical code (Parkhurst and Appelo, 1999). We were especially concerned
about the possible effects of varying calcium concentrations on 226Ra transport. Several aspects
related to the modeling approach, such as the use of daily versus monthly weather data and theeffect of assuming different transport regimes in the vadose zone (notably transient versus longterm steady-state flow) are discussed also. Additional details of the simulations are presented
elsewhere (Batalha et al., 2011, 2012).

2. Simulation Methodology
The scenario of concern is the repeated application of phosphogypsum to the same Cerrado soil
profile, and resultant 226Ra accumulation and transport in the vadose zone, and possibly into underlying groundwater. We considered a 2-m deep physically and chemically uniform soil profile,
without any preferential water flow. Since contaminant transport is very much affected by therates of water flow into and through the soil profile as a function of local weather conditions,transient precipitation and evaporation conditions were evaluated. Phosphogypsum was assumed
to be applied once each year for 27 years to the soil surface in the form of a short pulse withknown composition. Free drainage conditions were imposed at the bottom of the profile. In thisstudy we did not consider any uptake of water and solutes by plant roots.
Water flow in the soil profile is described using the standard Richards equation (Šimůnek et al.,2008):∂θ (h) ∂  ∂h  =
K ( h)  + 1  

∂z
∂z 
 ∂z  2(1)

where h is the soil water pressure head [L], θ is volumetric water content [L3L-3], t is time [T], z
is the spatial coordinate (positive upward in this study) [L], and K is the hydraulic conductivity[LT-1].
The volumetric water content (θ) and the hydraulic conductivity (K) in Equation (1) are strongly
nonlinear functions of the pressure head, h. These functions are described here using the constitutive relationship of van Genuchten (1980) as follows:

(

)

h< 0
θ −θ r  1 + α h== Se ( h )        
h≥0θ s −θ  
r

1

(2a)

K=
( h ) K s Sel 1 − (1 − Se1/ m )m 

(2b)

n −m



2

where Se is effective saturation [L3L-3], θr and θs are the residual and saturated volumetric soil
water contents, respectively [-], Ks is the saturated hydraulic conductivity [LT-1], and α [L-1], n>1[-], m=1-1/n [-] and l [-] are mostly empirical shape factors. Following van Genuchten (1980),
the pore-connectivity parameter l was fixed at a value of 0.5.
Radionuclide transport was simulated initially with HYDRUS-1D using the standard linear equilibrium advection-dispersion equation:∂ (θ RC ) ∂ 
∂C

=   
− qC  − µθ RC
θ D
∂t
∂z 
∂z


(3)

where C is the solution concentration [ML-3], R the retardation factor [-], D is the dispersion coefficient [L2T-1], q is the Darcy-Buckingham volumetric flux [LT-1], and μ is the decay coefficient [T-1] for the particular radionuclide being considered (in our case 226Ra). Since the half-lifeof 226Ra (1600 years) is much greater than the time period of our study (200 years), radioactive
decay will not be considered here (i.e., μ=0). The parameters R and D in Eq. (3) are given by:

R= 1+

ρb K dθq=
D   
Ddτ + λ

θ

(4)

(5)

in which ρb is the bulk density of the soil [ML-3], Kd is the distribution coefficient for linear sorption [L3M-1], λ is the longitudinal dispersivity [L], Dd is the diffusion coefficient in free water[L2T-1], and τ is the porous medium tortuosity factor [-] given by θ7/3θs-2 (Millington and Quirk,
1961).

3

Geochemical transport involving phosphogypsum was represented by a system of Nm reaction
equations for the primary (or master) species, in canonical form given byNm∑ν A = A
ji

j

(6)

i

j =1

and correspondent mass action equations assuming equilibrium:Nm

K i = γ i Ci ∏(γ j C j )− v ji(7)

j =1

where Nm the number of master aqueous species, Aj and Ai are the chemical equations for the
master and secondary species, respectively, νji are stoichiometric coefficients of the reaction, Ki
is the equilibrium constant of a reaction, γi is the activity coefficient of the aqueous complex, andγj is the activity coefficient for the jth master species in solution.
Transport of each species is represented in the above manner by a set of Nm multicomponent
reactive transport equations of the form (Jacques et al., 2008):∂ (θ C j ) ∂ ∂C j
=− qC j  + R0, j
θ D
∂t
∂z 
∂z


(8)

where R0,j is the reaction term representing the sources/sinks for each aqueous component, Cj isthe total concentration of master species jNsC=
c j + ∑ν ji ci
j

(9)

i =1

in which Ns the number of secondary aqueous species, and ci is the concentration of secondary
species i.
For 226Ra sorption we assumed a relatively simple equilibrium approach that lumps all of thesorption processes onto the solid phase into a single linear isotherm in which the distribution coefficient is a nonlinear function of the calcium concentration, as described in Batalha (2011b).

K d = K d 0CCap(10)in which Kd is the reference distribution coefficient for 226Ra when the calcium activity is 1, CCa
is the calcium activity in soil solution, and p is an empirical parameter whose value was fixed at-0.4 as deduced from data in USEPA (2004).

43. Cerrado Soils and PG Application Data
Flow boundary conditions at the soil surface were defined by meteorological data (rainfall andtemperature) obtained from a weather station in the western part of Minas Gerais state in Brazil(17o 25′ 11′′ S, 46o 88′ 64′′ W) for the period 1982-2005. These values were repeated many times
for the long-term (200 year) simulations. Daily potential evapotranspiration rates were calculatedwith the Hargreaves equation (e.g. Jensen et al., 1997) from available weather data and the geographic location of the site. Data of the rainwater composition were taken from Maier et al.
(1992) and adapted for use in HP1 module: pH 6.3; 5.59 10-5 mol/L of Ca and 6.05 10-5mol/L of
Mg. The soil chosen for the study is a typical Brazilian Cerrado oxisol, whose physical andchemical data were obtained from Oliveira (2008). Chemical parameters are listed in Table 1.
The granulometric composition of the soil was 16% sand, 63% fine sand, 5% silt and 16% clay.
Values of the soil hydraulic parameters needed for the water flow simulations were estimatedfrom the soil textural data using the Rosetta neural network based pedotransfer functions of
Schaap et al. (1998) as implemented in HYDRUS-1D. The soil bulk density was taken to be 1.50
g/cm3. For solute transport we used a diffusion coefficient in free water (Dd) of 1.0 cm2/d, whilea dispersivity (λ) of 20 cm was used for the simulations involving transient flow and 30 cm forthe steady-state flow regime.
Table 1. Soil chemical parameters of the Cerrado oxisol.
_______________________________________________________________________________________________
P
Ca2+
K+
Mg2+
Al3+
_______________________________________
pH
mmol/L______________________________________________________________________________________________________---5.3
0.081
0.500
0.410
0.100
0.333
_______________________________________________________________________________________________

Simulations of the cumulative flux leaving the bottom of the 2m soil profile over a period of 27
years assuming transient water flow were used to define the recharge rate (0.152 cm/d) used forthe long-term simulations assuming steady-state flow. Equation (1) for this purpose must be
augmented with a set of initial and boundary conditions. For Equation (1) we assumed an initialpressure head of -100 cm, the presence of atmospheric boundary conditions at the soil surface(i.e., given precipitation and calculated potential evapotranspiration data), and a free drainage
condition (dh/dx=0) at the bottom of the profile.
For the phosphogypsum and 226Ra transport calculations we considered PG application rates thatare common in Brazil (1.000 kg/ha year). The annual application was pulsed in 1 cm of rainwaterin the wet season (November). Measurements of the amount of 226Ra in the PG were carried out
at the Poços de Caldas Laboratory of the Brazilian Nuclear Energy Commission (CNEN). Results indicated an activity of 252±26 Bq per kg of PG. Mass balance calculations produced values of 73.45 mmol/L of CaSO4 and 3.04 x 10-13 mol/L of 226Ra when applied in 1 cm of rainwater for the transient simulations, and values of 483.2 mmol/L of CaSO4 and 20 x 10-10 mmol/L of226
Ra when simulating the constant flow regime using an averaged flux of 0.152 cm/day. A distribution coefficient (Kd) of 2000 cm3/g was used for radium based on data provided by ISAM
(1999) as described in Batalha (2011a).

5

4. Results and Discussion
Figure 1 shows calculated net infiltration fluxes (actual evaporation minus precipitation) into thesoil profile for nine years using transient conditions of precipitation and potential evaporation.
Results indicates two main seasons (wet and dry) during the year. The simulations using daily
weather data from Minas Gerais state (Brazil) revealed no run-off.
2
0
-2
-4
-6
-8
-10
Figure 1. Net fluxes across the soil surface during the first 9 years of the simulations(negative fluxes are downward into the soil profile).

Cumulative drainage fluxes from the soil profile are showed in Figure 2. Results are for bothtransient flow using daily weather data (left) and monthly averaged precipitation and evaporation
data (right). The importance of using daily weather data can be evaluated by the slopes of thetwo plots. The results indicate an average recharge rate of 0.152 cm/day when using daily weather data, and a nearly 50% lower value of 0.079 cm/d when monthly averaged values for precipitation and potential evaporation are used. This means that time-averaging weather data (especially precipitation) can lead to significantly lower estimates of the recharge rate.
0

0

-200

-100

-400

-200

-600

-300

-800

-400

-1000

-500

-1200

-600

Figure 2. Calculated cumulative recharge rates from the soil profile during 27 years using daily(left) and monthly averaged (right) weather data.

6

Figure 3 compares calculated 226Ra concentration distributions in the soil profile assuming transient (left) and steady-state (right) flow conditions. Results for transient flow were obtained using observed daily weather data, while those for steady-state flow assumed a constant long-time
average recharge rate of 0.152 cm/d as deduced from the transient simulations (Fig. 2, left).
Since the transient regime produced a slightly more dispersed contaminant front, a higher longitudinal dispersivity (DL = 30 cm) for the steady flux scenario was used. The results in Figure 3
indicate that long-term predictions can be based equally well assuming steady-state water flow,rather than transient flow, provided the correct long-term deep drainage (recharge) rate is used
for the steady-state simulations.

0

0

-5

-5

-10

-10

-15

-15

-20

-20

-25

-25

-30

-30

-35

-35

-40

-40

-45

45
226

Figure 3. Calculated Ra concentrations after T0= 0, T1=9, T2=18, and T3=27 years assuming
transient weather data (left) and a constant long-term average recharge rate (right).

Figure 3 shows that after 27 years of PG application, essentially no 226Ra had reached a depth ofabout 30 or 40 cm. Subsequent transport after PG applications ceased, but while maintaining thesame water flow conditions, remained equally slow. This is shown in Figure 4 which presents
calculated 226Ra distributions in the soil profile after 27, 100 and 200 years. The plots show very
little movement into the lower part of the profile, with essentially no 226Ra reaching a depth of 1
m after 200 years. Cumulative amounts leached beyond 40 cm depth after 100 and 200 years
were estimated to be 3.58 x 10-13 and 1.85 x 10-12 mmol/cm2. Also notice that the maximum concentrations in the upper part of the profile decreased substantially with time.
Figure 5 shows the influence of the calcium concentration on 226Ra transport after 9, 18 and 27
years of phosphogypsum application. The plots compare simulations for a constant distribution
coefficient, Kd (left) and for a distribution coefficient that varies as a function of Ca2+ activity(right). Results indicate that Ca2+ resulting from the PG application has a major effect on boththe transport rate of 226Ra and the solution concentration. However, this trend is only present initially. Once the PG is no longer applied (after 27 years), calculated 226Ra curves slowly converged to the point that after 100 or 200 years they became very similar (results not further
shown here).

7

0
-20
-40
-60
-80
-100
-120
-140
-160
Figure 4. Calculated 226Ra concentrations after T0= 0, T1=27, T2=100, and T3=200 years (left).

0

0

-5

-5

-10

-10

-15

-15

-20

-20

-25

-25

-30

-30

-35

-35

-40

-40

45

45

Figure 5. Calculated 226Ra concentrations after T1= 9, T2=18, and T3=27 years, without (left) andwith (right) accounting for the effect of the Ca+ activity on the Kd of 226Ra.

5. Conclusions
In this paper we analyzed the fate and transport of 226Ra present in phosphogypsum when applied
for many years as an amendment to a typical Brazilian Cerrado soil used for crop production.
The main factors that influenced water flow and radionuclide transport were the deep drainage
flux determined by rainfall and evaporation and the use of daily precipitation and potential evaporation data. Results indicate that for the assumed long-term transport scenario (200 years andhaving relatively low 226Ra activities), steady-state flow conditions could be assumed without8

significantly affecting the predictions. Most of the 226Ra remained in the upper 40 cm of the soil
profile.
The effect of PG derived calcium was important by decreasing the effective distribution coefficient of 226Ra, and as such increasing the transport rate. All of the simulations in this study assumed equilibrium water flow and solute transport conditions. We did not consider preferential
flow or nonequilibrium transport, or any other non-ideal transport process (such as colloidfacilitated transport). These processes could all lead to much higher rates of transport in the soil
profile (e.g., Šimůnek and van Genuchten, 2008).
In conclusion, we believe that modeling efforts such as those in this study using HYDRUS-1Dand HP1 are essential for detailed environmental assessment of the long-term use ofphosphogypsum or other naturally occurring radioactive materials (NORMs). They are especially needed for scenarios that cannot be easily monitored and/or occur over relatively long periodsof time.
References
Batalha, M. S., Multicomponent transport of contaminants released in the environment due tophosphogypsum use, MS dissertation in Mechanical Engineering, COPPE/UFRJ, Rio de Janeiro, RJ,
Brazil., 86 pp., 2011a (in Portuguese).
Batalha, M. S., C. R. Bezerra, M. C. Barbosa, D. Jacques, E. M. Pontedeiro and M. Th. van Genuchten,
Multicomponent transport of 226Ra present in phosphogypsum applied to a Brazilian Cerrado soil. VI
Brazilian Symposium on Geosynthetics, REGEO2011, Belo Horizonte, MG, Brazil, Nov. 21-14, 2011b(in Portuguese).
Batalha, M. S., C. B. Bezerra, D. Jacques, M. C. Barbosa, E. M. Pontedeiro, and M. Th. van Genuchten,
Multicomponent transport predictions of 226Ra in soil following the use of phosphogypsum,
WasteEng 2012,14th Int. Conf. on Engineering for Waste and Biomass Valorization, Oporto, Portugal,
2012.
Embrapa, Uso de Gesso Agrícola nos Solos do Cerrado, Empresa Brasileira de Pesquisa Agropecuária,Planaltina, Circular Técnica 32, 18 p., 2005.
Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, Modeling coupled hydrological andchemical processes: Long-term uranium transport following phosphorous fertilization, Vadose Zone
J., 7,698-711, 2008.
Jacques, D. and J. Šimůnek, Notes on HP1, a software package for simulating variably-saturated water
flow, heat transport, solute transport, and biogeochemistry in porous media, HP1 Version 2.2., Open
Report, SCK•CEN-BLG-1068, Mol, Belgium, 129 p., 2010.
Jensen, D. T., G. H. Hargreaves, B. Temesgen, and R. G.Allen, Computation of ETo under nonideal conditions, J. Irrig. Drain. Eng., ASCE, 123(5), 394-400, 1997.
Maier, M. H., M. Takino, and L. L. Miyamaru, Alguns aspectos da composição da água da chuva, Acta
Limnol. Brasil., IV, 247-254, Brazil, 1992
Mazzili, B., V. Palmiro, C. Saueia, and M. B. Nisti, Radiochemical characterization of Brazilian
phosphogypsum, J. Environ. Radioactivity, 40, 113-122, 2000.
Millington, R. J., and J. M. Quirk, Permeability of porous solids, Trans. Faraday Soc., 57, 1200-1207,
1961.
Oates, K. M., and A. G. Caldwell, Use of by-product gypsum to alleviate soil acidity, Soil Sci. Soc. Am J.,
49, 915-918, 1988
Oliveira, K. A. P., Aplicação do fosfogesso na agricultura do Cerrado e suas implicações radiológicas,
Dissertação de Mestrado, CDTN/CNEN, Belo Horizonte, MG, Brazil, 2008.

9

Oster, J. D., Gypsum usage in irrigated agriculture: A review, Fertilizer Res., 3, 73-89, 1982.
Papastefanou, C., S. Stoulos, A. Ioannidou, and M. Manolopoulou, The application of phosphogypsum inagriculture and the radiological impact, J. Environ. Radioactivity, 89, 188-198, 2006.
Parkhurst, D. L., and C. A. J. Appelo, User's guide to PHREEQC (Version 2), A computer program forspeciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, Water Resources Investigations Report 99-4259, Denver, CO, 1999.
Rothbaum, H. P., D. A. McGaveston, T. Wall, A. E. Johnston, and G. E. G. Mattingly, Uranium accumulation in soils from long-continued applications of superphosphate, J. Soil Sci. 30, 147-153, 1979.
Saueia, C. H. R., and N. P. Mazzilli, Distribution of natural radionuclides in the production and use ofphosphate fertilizers in Brazil, J. Environ. Radioactivity, 89, 229-239, 2006.
Schaap, M. G., F. J. Leij, and M. Th. van Genuchten, Neural network analysis for hierarchical predictionof soil hydraulic properties, Soil Sci. Soc. Am. J., 62, 847-855, 1998.
Šimůnek, J., and M. Th. van Genuchten, Modeling nonequilibrium flow and transport processes using
HYDRUS, Vadose Zone J., 7(2), 782-797, 2008.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van Genuchten, The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media,
Version 4.0, HYDRUS Software Series 3, Department of Environmental Sciences, Univ. of California,
Riverside, CA, 2008.
U.S. Environmental Protection Agency, Understanding Variation in Partition Coefficient, Kd, Values.
Vol III, EPA 402-R-04-002C, EPA, Washington, DC, 2004
van Genuchten, M. Th., A closed-form equation for predicting the hydraulic conductivity of unsaturated
soils, Soil Sci. Soc. Am. J., 44, 892-898, 1980.
Zielinski, R.A., S. Asher-Bolinder, A. L. Meier, C. A. Johnson, and B. J. Szabo, Natural or fertilizerderived Uranium in irrigation drainage: A case study in southeastern Colorado, Appl. Geochem., 12, 921, 1997.

10