-------
Table 4-1. Sensitivities and Relative Sensitivities of Output F with Respect to the Six Input
Parameters, for Arbitrary Values of the Inputs.
Parameter
Sensitivity, S
Relative Sensitivity, Sr
a
e
b + cd + eln(c2f)
dc + 2e
dc + b + eln(c2f)
ac
e
cd
cd + b + eln(c2f)
b + cd
e£n(c2f) + b + cd
a
f
e£n(c2f) + b + cd
which is outside the domain of d. From these three input parameters (a,b,d), we can see the effects of linearity,
homogeneity and sign of the slope, as well as the domains of definition of the input parameters. Also in Figure 4-3 is
theplotof F(2,l,c,3,4,5) versus c. Because of the particular domain of c and the presence of ln(c4) inF, there is a
vertical asymptote at c = 0 in the graph of F versus c. In addition the linear-logarithmic variation of c in F results in
two zeros, one at a and one at p. The vertical asymptote shows up as a vertical asymptote in Fc but a zero in Frc.
The zeros (a, p) in F have no particular effect on Fc but become vertical asymptotes in Frc. The y zero in Fc carries
over into a zero of Frc. These zero/asymptotic relationships will exist between F, Fa, and Fra as long as we do not
run into indeterminant forms, 0/0 or °°/°°. For these cases, the asymptote and zeros in Fa, and Fra may disappear.
Finally, in Figure 4-4, we note that for the domain of definition in F(2,l,-l,3,e,5) there is a zero at a which produces
a vertical asymptote inFre. The quantity F(2,l,-1,3,4,1) has no zeros because of the restricted domain of definition
for f. One can also note the different characteristics in Fre and Frf at the right side of the domains of definition of e
and f, respectively. This is because the quantity F(2,l,-l,3,e,5) approaches a horizontal asymptote (F=3.219) for
sufficiently large e, while F(2,l,-l,3,4,f) approaches +°° in a logarithmic manner as f becomes large. Thus, by
studying these simple forms, one can better interpret the sensitivity and relative sensitivity results for the five models
being analyzed in this report.
32
-------
Table 4-2. Sensitivity and Relative Sensitivity of Output F to Individual Input Parameters with Reference to
the Base Case Given in Equation (4-10).
Parameter
a
b
c
d
e
f
Output F
Formula
Ui(5)- - a
L 2\
b 3
+ £n(25)
2 2
4 3c 1
ln(c ) + + ln(25) +
2 2
d 1
1 h fn(25)
1 1
4
- - + #j(25)
e
2 te(f ) - 1
Characterization
Linear and
Homogeneous in a.
Linear and Non-
homogeneous in b.
Linear-Logarithmic
and Nonhomo-
geneous in c.
Linear and Nonhomo-
geneous in d.
Nonhomogeneous in e,
Inverse Variation with e.
Logarithmic and
Nonhomogenous in f.
Sensitivity
s
1
9
1
2
4 3
p 9
\
2
A
e
2
—
f
Relative Sensitivity
sr
1
b
b - 3 + £n(625)
8
c + -
3
8/3 1 4
3 3
d
d - 1 - ln(625)
2
ln(5) e - 2
1
tn(f] -
2
33
-------
1.2
1.0
.8
F .6
1 ra
.4
.2
0
.4
-1012345
a
1.2
1.0
.8
"a .6
.4
.2
0
-1012345
0
F -2 -
rrb
.1 -
0 1/2 1 3/2 2
0 1/2 1 3/2 2
F
-02345
F
3 -
2 -
1 -
Figure 4-2.
0 1/2 1 3/2 2
Sensitivities and relative sensitivities of F with respect to a and b, with reference
to the base case in Equation (4-10).
-------
1 -
Frd
-2-
-3.
12345
-4 -3-2-10 12
-1.0
12345
d
3 -
F 2-
1 -
12345
d
Figure 4-3. Sensitivities and relative sensitivities of F with respect to c and d, with reference
to the base case in Equation (4-10).
35
-------
10
6 -
4 -
2 -
're o
-2 -
ia
a= llln (5)
1234567
e
'e 2-
1234567
1 234567
234567
f
Figure 4-4. Sensitivities and relative sensitivities of F with respect to e and f, with reference
to the base case in Equation (4-10).
36
-------
Section 5
Parameter Sensitivity Analysis: Hypothetical Modeling Scenario
The characteristics of a representative physical site are required to realistically study the sensitivity of the five
models proposed in this report for simulating the transport and fate of radionuclides in the unsaturated zone. Such a
site would be used to set the standards for base parameter choices in the sensitivity tests and for parameter range
determination. The ideal candidate site should have well characterized soil properties, and should have sufficient
chemical transport and fate data and climatological information for running the pertinent simulations. Once a site is
chosen, a conceptual model of that site is required so that the site data, the five proposed modeling codes, and the
intended application of the models are all consistent with one another. In what follows, we consider the site selection
process, selection of the candidate site, characteristics of the Las Graces Trench Site in New Mexico, development of
a conceptual model, and the base parameter selection.
5.1 Site Selection Process
The sensitivity studies reported in this document are part of an overall evaluation of the general applicability of
each model (i.e., HYDRUS, CHAIN 2D, FECTUZ, MULTIMED-DP 1.0, and CHAIN) for simulating the fate and
transport of radionuclides. The radionuclides of interest are the following five elements/isotopes (U.S. EPA, 2000b):
Element/Isotope Kd Default Values in ml/g
Plutonium (238Pu) 5
Strontium (90Sr) 1
Technetium ("Tc) 0.007
Tritium (3H) 0
Uranium (238U) 0.4 ,
where Kd is the soil/water partition coefficient for the linear Freundlich isotherm. These five elements/isotopes are
taken to be the parent species in the decay chains and the major segments of the chains are shown in Figure 5-1 along
with the corresponding half-lives of the species.
The quantity, Kd, in all five models being studied, is a lumped parameter representing known and unknown
phenomena. It tends to lose some of its meaning in the modeling world, while retaining its full meaning in the
laboratory. In the laboratory Kd is determined under carefully controlled conditions, but the real world cannot be
completely controlled or measured. Thus, the conditions surrounding the sorption phenomenon must be estimated,
and these estimates will only represent localized conditions in the vicinity of the sampling. Consequently, the large
heterogeneities of most, if not all, subsurface systems make it difficult to identify a single Kd value for the system.
Further, unless the model has a geochemical component, Kd will also represent everything one does not know about
the site's geochemistry. The relationships between Kd values, the chemical species, and a site's geochemistry are
discussed in U.S. EPA (1999ab, 2000b). However, in the current report, the default Kd values, given above, are
taken as the base values for the sensitivity analyses combined with an interval domain about these values, to partially
represent the heterogeneity found in the vadose zone.
Based on the above radionuclides, the site selection process consisted of identifying three features: site listing,
data availability, and existing field and modeling studies. By reviewing existing literature and databases, a list of
candidate sites was developed, see Table 5-1. These candidate sites, in general, were either contaminated with
radionuclides, or were current or future radioactive waste disposal sites. For example, Superfund sites with soil
contaminated by one or more of the above five radionuclides were identified in the U.S. EPA's VISITT database (i.e.,
the Vendor Information System for Innovative Treatment Technologies, Version 6.0). Many of these radionuclide
37
-------
PLUTONIUM,238Pu
88 yr. half-life
~v
URANIUM, 234U
2.4x 105yr. half-life
STRONTIUM, 9°Sr
29 yr. half-life
>.
YTTRIUM, 90y
64 hr. half-life
>.
ZIRCONIUM, 9°Zr
Stable
TECHNETIUM,"Tc
2.1 x 105yr half-life
V
RUTHENIUM, "Ru
Stable
TRITIUM, 3H
12 yr. half-life
•v
HELIUM-3,3He
Stable
URANIUM, 238U
4.5x 109yr. half-life
^
THORIUM, 234Th
24 day half-life
99.8%
PROTOACTINIUM,234m Pa
1 min half-life
URANIUM,234 U
2.4x 105yr. half-life
0.2%
PROTOACTINIUM, 234Pa
7 hr half-life
Figure 5-1. The major segments of the decay chains for the five elements/isotopes considered
to be parents in these analyses (U.S. EPA, 2000b).
contaminated sites were a result of past nuclear production and related activities, and are currently under the
environmental restoration program operated by the U. S. Department of Energy. With reference to low-level
radioactive waste storage sites, the most attractive areas were those having low annual precipitation, high
evapotranspiration, and thick unsaturated soils/sediments. An example of such a site is the Mojave Desert Waste-
Burial Site near Beatty, NV. In addition, the Nevada Test Site is a radionuclide contaminated site from nuclear
testing as well as an on-site/off-site low-level waste disposal facility. The only non-nuclear contaminated site and
non-nuclear storage site considered in this evaluation was the Las Graces Trench Research Site in New Mexico. The
field studies at this site have been used to provide data to test deterministic and stochastic models for water flow and
solute transport, Wierenga, et al. (1991), Hills, et al. (1991, 1994), and Rockhold, et al. (1996). The Las Graces
Trench Site was included in the list of candidate sites because it fit the above physical characteristics for waste
disposal facilities and because many detailed site characterization studies and field experiments have been conducted
at this location.
The availability of soil data and climatological records are essential for the proper simulation of water flow and
radionuclide transport in unsaturated soils. From the twenty-some sites in the original candidate list, four sites that
had the most available data were selected for further evaluation. These were the Las Graces Trench Site, the Beatty
Waste Storage Site, Idaho Falls National Laboratory and the Hanford Washington Site. In addition, the Perdido
Alabama Site was placed on the list; this was the site that was analyzed in the earlier sensitivity study of Nofziger, et
al. (1994). Actually, the Perdido Site, which is a benzene contaminated site in an area of high annual precipitation,
was assigned soil properties obtained from another site having the same soil series. Thus, no site-specific soil
properties were available for the Perdido site. This led to the elimination of the Perdido site from further
consideration. The existing field/modeling studies for the four remaining sites were reviewed and summarized. The
studies considered were water balance, recharge and infiltration experiments; tracer and radionuclide migration
studies; laboratory and field measurements of bulk density, Ks, and soil water retention curve data; and model
construction, testing, and simulation runs.
5.2 Selection of the Candidate Site
For the final four sites (Las Graces, Beatty, Idaho Falls, Hanford), additional factors and options were
considered in the choice of the best candidate site. Two of the important factors were the availability of required soil
data for the five pre-selected transport models in this study and the attainability of transport parameters from the
38
-------
Table 5-1. Partial List of Radionuclide Contaminated and Disposal Sites in the U.S. (U.S. EPA's VISITT Database).
SITE NAME
Location,
City
State
Contaminants
Site Characteristics
1 Rocky Flats Environ- Golden
mental Technology Site
Colorado
2 Idaho National
Engineering &
Environmental
Laboratory
3 Ethyl Corporation
4 DOE FUSRAP
St. Louis Site
Lreatability Study
5 Nevada Lest Site
Idaho Falls Idaho
Baton Rouge Louisiana
St. Louis Missouri
Mercury
Nevada
Plutonium, Uranium
Radioactive Materials,
Plutonium, Lritium
Uranium
Uranium
Plutonium, Lritium
6 Mojave Desert
Waste Burial Site
7 NL Industries, Inc.
8 Los Alamos Natl.
Laboratory
9 Las Cruces Lrench
Site
10 Sandia Natl. Lab.
11 West Valley Nuclear
12 Mound Demonstration
Project
13 Portsmouth Gaseous
Diffusion Plant (DOE)
14 EPASILE
Demonstration
15 Femald Feed Materials
Production Center
16 Apollo Fuel
Conversion Plant
17 Savannah River Site
18 INEL Pit 9 Pilot Project
19 K-25 Site
20 Hanford Site
21 DOE Morgantown
Energy Lech. Center
Beatty
Nevada
Various Radionuclides
Superfund Site. On-site waste
disposal
Radionuclide Contaminated Site:
Radioactive Waste Management
Complex was used for disposal
site for low-level and transuranic
radioactive waste
Superfund Site
Superfund Site
Hazardous Waste Site: Radio-
nuclide contaminated site from
nuclear testing: Low-level waste
disposal facility for both onsite
and off-site generated defense
low-level waste
Low-level radioactive waste and
hazardous chemical waste disposal
Pedricktown
Los Alamos
Las Cruces
Albuquerque
West Valley
Miamisburg
Piketon
Alliance
Fernald
Apollo
Aiken
Clem son
Oak Ridge
Richland
New Jersey
New Mexico
New Mexico
New Mexico
New York
Ohio
Ohio
Ohio
Ohio
Pennsylvania
South Carolina
South Carolina
Lennessee
Washington
Strontium
Plutonium, Uranium
None
Uranium, Plutonium
Strontium
Plutonium
Uranium, Lechnetium
Uranium
Uranium, Lechnetium
Uranium
Uranium
Uranium
Lechnetium, Uranium
Uranium, Strontium
Superfund Site
Superfund Site
Experimental Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site
Superfund Site. C
Morgantown West Virginia Various Radionuclides
disposal.
Superfund Site
39
-------
existing site modeling studies. Selection of the final candidate site was also dependent upon decision-making
options. For example, when verification of modeling results is an important issue (i.e., comparison of modeling
results with field and laboratory experimental data) in addition to data availability, the Las Cruces Trench Site was
the best candidate. However, when reliable model input data is the only issue and model verification is less
important (i.e., comparison among the SSG screening models is more important than model field verification or
comparison of the screening models with more comprehensive modeling of the area), then the Beatty Waste-Burial
Site is a good alternative candidate.
Four of the five proposed SSG models, HYDRUS, CHAIN 2D, MULTIMED-DP 1.0, and FECTUZ, require van
Genuchten soil-water retention parameters (Table 2-1) These parameters plus dispersivities were available at the Las
Cruces Site. For the Beatty Site, the van Genuchten parameters were not available and would have to be derived
from the raw soil-water content/pressure head data (Andraski, 1996). Furthermore, a great amount of model testing at
the Las Cruces Site has been conducted for the transport and fate of tritium, bromide and chromium in the
unsaturated zone. Considering the location and climatological features of the site, and the abundance of existing
field data in conjunction with significant modeling studies, the Las Cruces Trench Site was chosen as the candidate
site for the model evaluation in the SSG for radionuclides.
5.3 Characteristics of the Las Cruces Trench Site in New Mexico
In the soil screening level (SSL) process, generic SSLs for radionuclides can be calculated based on a number of
default assumptions chosen to be protective of human health for most site conditions. However, these are expected
to be more conservative than calculated site-specific SSLs. When site-specific SSLs are of interest, a simulation
using an appropriate model and site-specific data is required. Thus, in the current study, a typical physical site is
required for the evaluation of the HYDRUS, CHAIN 2D, FECTUZ, MULTIMED-DP 1.0, and CHAIN Codes in the
SSL process.
The Las Cruces Trench Research Site in New Mexico was chosen as the typical physical site because:
1. The site characteristics met the criteria of radioactive waste disposal areas (low annual precipitation,
high annual evapotranspiration, and a deep water table);
2. The site had been subjected to extensive testing of its soil physical and chemical properties and soil-
moisture distribution and movement in the unsaturated zone;
3. Results of tracer tests (chloride, bromide and tritium) and metal movement (chromium) were available.
This site lies in the Chihuahuan Desert Province of southern New Mexico (Bailey, 1980). This province is
mostly desert and the Rio Grande and the Pecos River and a few of their larger tributaries are the only perennial
streams. The area has undulating plains with elevations near 1200 m from which somewhat isolated mountains rise
600 m to 1500 m. There are washes which are dry most of the year that fill with water following a rain. Basins that
have no outlets drain into shallow playa lakes that dry up during rainless periods. Extensive dunes of silica sand
cover parts of the province, and in places, there are dunes of gypsum sand, the most notable being the White Sands
National Monument near Alamogordo (a town 100 km northeast of Las Cruces).
The climate of the Chihuahuan Desert is distinctly arid and the spring and early summer are extremely dry.
During July the summer rains usually begin, and they continue through October (Bailey, 1980). In general, these
summer rains are local torrential storms. Average annual temperature in the province ranges from 10°C to 18°C.
Summers are hot and long, and winters are short but may include brief periods when temperatures fall below
freezing. The characteristic vegetation of the Chihuahuan Desert is a number of thorny shrubs. These shrubs
frequently grow in open stands, but sometimes form low closed thickets. Short grass often grows in association with
the shrubs. On deep soils, mesquite is usually the dominant plant; creosote bush covers great areas in its
characteristic open stand and is especially common on gravel fans. Royo (2000) says the creosote bush is a desert
plant "par excellence," a true xerophyte. It is a drought-tolerant shrub with small dark green leaves and has an
extensive double root system - both radial and deep - to accumulate water from both surface and ground water.
These plants can tolerate up to two years without precipitation. The leaves are coated with a varnish-like resin which
reduces water loss by evaporation. The original creosote bush can live to about 100 years old, but it can produce
clones of the parent as the bush ages. These clones are produced in a circular pattern of genetically-identical plants,
expanding outward at the rate of about one meter every 500 years. The "King Clone" family on BLM land near
Victorville, California, is estimated to be 11,700 years old.
40
-------
The actual experimental site is located on the New Mexico State University college ranch some 20 or more
kilometers due north of the city of Las Graces, New Mexico. The site is on a basin slope of Mount Summerford at
the north end of the Dona Ana Mountains (Wierenga, et al, 1991; Defense Mapping Agency, 1987). Geologically,
these mountains are a domal uplift complex composed of younger rhyolitic and the older andesitic volcanics which
were intruded by monzonite. The covered trench that provides horizontal access to the experimental plots and is
used to provide soil samples for these plots is an evacuated earthen box with dimensions 26.4 m long, 4.8m wide,
and 6.0 m deep.
The published soil hydraulic properties for this site, given by Wierenga, etal. (1991), are listed in Table 5-2.
Some pertinent site characteristics obtained by Gee et al (1994) are given on Table 5-3. Table 5-2 shows that the
estimates for the hydraulic parameters were obtained for a uniform soil model and for a nine-layered soil model. The
layers in the layered soil model correspond to the nine soil layers identified at the site. The saturated hydraulic
conductivity for each soil layer was estimated by taking the geometric mean of the 50 laboratory-measured saturated
hydraulic conductivities obtained from each soil layer. Likewise, the water retention data from all 50 samples from a
given layer were used to estimate a,p,9r, and 9S for a single water retention curve for that layer. For the uniform soil
model, the geometric mean of 450 laboratory measured saturated hydraulic conductivities (nine layers with 50 per
layer) was used to estimate a uniform soil saturated hydraulic conductivity value. Likewise, the water retention data
for all 450 sample locations were simultaneously used to estimate single values for each of the parameters a,p,9r, and
9S in a least squares sense (Wierenga, et al., 1991).
In Porro and Wierenga (1993), the solute transport dispersivity (cm) was determined for six layers ranging over
0 < z < 500 cm. The values varied from 2.20 cm to 7.80 cm with an arithmetic average of 4.53 cm for the combined
layer of 500 cm. Minor adaptations of some of these data have been made in the conceptual model developed for the
sensitivity analyses. The details are given in the following subsection. Finally, Figure 5-2 shows the daily
precipitation and potential evapotranspiration (PET) at the Las Graces Site. The PET is calculated from daily
climatic data using Penman's general equation for a well-watered grass reference crop (Jensen, et al., 1990):
XEto = T • (Rn - G) + 6.43 (1 - T) W(e0 - e) , (5-1)
where the terms in Equation (5-1) are defined as:
'k = latent heat of vaporization in mega-joules per kilogram,
Eto = evapotranspiration rate (E^ from a well-watered grass reference crop, in kilograms per
meter squared per day,
^Eto = latent heat flux density in mega-joules per meter squared per day,
F = dimensionless parameter dependent upon surface elevation and air temperature, Table 6-1 of
Jensen, etal., 1990,
Pvj, = net radiation at the surface in mega-joules per meter squared per day,
G = heat flux density to the ground in mega-joules per meter squared per day,
W = a + b u2 = wind function in meters per second,
a,b = positive constants,
u2 = wind speed at 2m above surface in meters per second,
e0 = saturated vapor pressure of air at some height z in kPa,
e = water vapor pressure in air at height z in kPa,
PET = potential evapotranspiration rate in mm/d, given by Eto -^water density in kilograms per meter
squared per millimeter.
For the eighteen year period (1983-2000) shown in Figure 5-2, the annual precipitation ranged over the values from
11.5 cm/yrto 30.8 cm/yr, with an annual average of 22.5 cm/yr. For this same period, 28% of the precipitation came
in the January to June period, and 72% came in the July to December period. This is consistent with the
climatological description given by Bailey (1980).
5.4. Development of a Conceptual Model
The conceptual model for SSLs using detailed site-specific data is developed in a manner which is theoretically
and operationally consistent with the simplified methodology described in Section 2 of the Soil Screening Guidance
41
-------
Table 5-2. Soil Hydraulic Properties at the Las Graces Trench Site for SSG Model Evaluation Study
(reprinted from Water Resources Research, 1991, by PJ. Wierenga, R.G. Hills, and D.B. Hudson
with the permission of the American Geophysical Union, Washington, DC).
Layers Depth
(cm)
all 0 - 600
Saturated
Water
Content
(cnf/cnf)
0.321
Residual
Water
Content
(cnf/cnf)
Uniform Soil
0.083
van Genuchten
Alpha
Coefficient, a,
(cm1)
Model
0.055
van Genuchten
Beta
Coefficient, ft,
(—)
1.509
Saturated
Hydraulic
Conductivity, Ks,
(cm/d)
270.1
Layered Soil Model
1 0-15
2 15-140
3 140 - 205
4 205 - 250
5 250 - 305
6 305 - 370
7 370 - 460
8 460 - 540
9 540 - 600
0.348
0.343
0.336
0.313
0.302
0.294
0.310
0.325
0.306
0.095
0.091
0.085
0.071
0.072
0.090
0.073
0.083
0.078
0.042
0.062
0.060
0.068
0.040
0.070
0.027
0.041
0.047
1.903
1.528
1.574
1.537
1.550
1.711
1.418
1.383
1.432
539
250
267
300
250
334
221
172
226
Table 5-3. Characteristics of the Las Graces Trench Site for SSG Model Evaluation Study
(reprinted from Soil Sci. Soc. Am. J., 1994, by G.W. Gee, PJ. Wierenca, B.J. Andraski, M.H.
Young, MJ. Payer, and M.L. Rockhold with the permission of Soil Sci. Soc. Am. J., Madison,
WI).
Annual
Precipi-
tation
(cm/yr)
23
Annual
Potential
(Pan)
Evapora-
tion (cm/yr)
239
Annual
Potential
Recharge
(cm/yr)
8.7
Average
Daily Max.
Air Temper-
ature (°C)
28
Average
Daily Min
Air Temper-
ature (aC)
13
Elevation
(m)
1357
Depth to
Water
Table
(m)
60
Geology
Alluvial
Typical
Soil
Type
Berino fine
loamy sand
Typical
Vegeta-
tion
Creosote
bush
42
-------
"£•100
>§, 80
o 60
4°
o
£
a.
20
0
,1 UkiLiil
Las Cruces, NM
lllkJ I
Li
1983 1986 1989 1992
Calendar Year
1995
1998
2001
Figure 5-2. Daily precipitation and potential evapotranspiration (PET) at Las Cruces Site, NM. PET is
calculated from daily climate data using Penman's equation (Jensen et al., 1990).
for Radionuclides: Technical Background Document (U.S. EPA, 2000b). In so doing, it is assumed that the Las
Cruces Trench Site in New Mexico has been used as a waste disposal/storage facility where radionuclides from tank
leaks or improper waste disposal were released to the soil surface for a period of time with a specified total amount
of release set for each radionuclide (e.g., 10 mg/cm2 of 238U was released). Thus, a finite radionuclide contaminated
source is assumed. The driving force to send this material on a downward migration to the water table is the
infiltrated rainfall which produces a net annual recharge rate of 87 mm/yr. (Table 5-3). The site-specific soil
hydraulic properties given in the "all-layer" row of Table 5-2 and the mean layer dispersivity of 4.53 cm, obtained
for tritium transport through Berino fine loamy sand, are also used in the current analyses.
Further assumptions are:
The unsaturated zone is homogeneous although HYDRUS, CHAIN 2D, FECTUZ and MULTIMED-
DP 1.0 are capable of simulating layered soils and the hydraulic properties for layered soils are
available at the site.
There is no significant vapor pressure for most radionuclides (except for radon, which is not considered
in this study) and the dimensionless Henry's Law constant is assumed to be zero.
Only vertical flow and transport are considered; horizontal flow and transport are ignored, even though
CHAIN 2D is capable of simulating two-dimensional flow and transport.
There is no chemical or biological degradation in the unsaturated zone.
There is radioactive decay (Figure 5-1) in the unsaturated zone; however, since the published data on
decay rates (or half-lives) for radionuclides are reasonably accurate and precise, no sensitivity analysis
will be made on decay rates.
Complexation, oxidation-reduction, dissolution and precipitation, and ion-exchange are not considered
because these processes are not implemented in the five models being evaluated.
There is no facilitated transport (e.g., colloidal transport, preferential flow in fractures or root channels,
fingering pathways) of the radionuclides in the unsaturated zone.
The aquifer lying below the unsaturated zone is unconsolidated and unconfmed. However, flow and
transport of the radionuclides in the saturated zone are not considered and only leachate contributions
to the ground water at the center of the disposal area are of interest.
Initial concentrations of radionuclides in the soil are zero.
-------
In order to specify a total amount of radionuclide and a reasonable level of a radionuclide concentration release
from the waste source, a literature survey of radionuclide contamination in soils was conducted. Based on the survey
information, the radionuclide concentration released from the hypothesized waste source and the duration of
radionuclide release were determined. Thus, the total amount of radionuclide release can be obtained from the
product of the recharge rate, source concentration, and the duration time of waste release. The depth of radioactive
contaminated soil at the termination point of waste release from the source can be determined from the product of the
pore velocity times the duration time of waste release.
The decay series of the radionuclides in the sensitivity analyses in this report are based on the chain segments
given in Figure 5-1. The chain segments for plutonium-238 and uranium-238 are not complete but are sufficient for
current purposes because of the long half-life of uranium-234. Further, due to the time discretization in the five
models and the relatively short half-lives of the intermediate species in the strontium-90 and uranium-238 chains, the
following parent-daughter chains are only of importance in this report:
5.5 Base Parameter Selection
Table 5-4 provides the base values of selected input parameters to be used in the sensitivity analyses given in
later sections. These values are those that are basically used in the analyses of HYDRUS and CHAIN 2D. Subsets
of these values will be used in the analyses of FECTUZ, MULTIMED-DP 1.0 and CHAIN. Justification and
rationale for the use of these specific base parameters are presented in the following paragraphs.
The area of the disposal facility is arbitrarily taken as 400 m2, whose length and width are both 20 m. The one-
dimensional vertical models are assumed to be located at the center of the square, thus eliminating edge effects in the
unsaturated zone simulation. The total duration of the release of the radionuclide mass was arbitrarily chosen to be
1000 days; however, this value was reasonably consistent with the survey information mentioned in the previous
subsection.
The potential evapotranspiration (PET) from the surface of well-watered short grass was calculated from Las
Graces climatic data using Penman's equation (Jensen, et al. 1990). The weather/climatic data required for this
equation are temperature, relative humidity, wind speed, and solar radiation, along with the estimated albedo
coefficient of 0.21 and the 1357 m elevation of the site. Using the daily climatic data for the years 1983 to 2000, the
estimated value of PET is 204 cm/yr and the average precipitation at the site is 22.5 cm/yr. These values are
consistent with those reported by Gee, et al. (1994), see Table 5-3; thus, we use the annual potential recharge of 8.7
cm/yr given by Gee, et al. In Table 5-2, the all-layer residual moisture content is given as 0.083 and the all-layer
saturated moisture content is given as 0.321. It is expected that the "uniform" annual recharge of 8.7 cm/yr will
produce a soil moisture somewhere between 9r and 9S. Thus, we arbitrarily chose an initial water content, 9, equal
to the geometric mean of 9r and 9S, (9r9s )1/2. This value to two significant figures is 9 = 0.16.
As stated in the previous subsection, the total mass of the individual radionuclides released from the hypothetical
Las Graces disposal/storage site was chosen to be consistent with radionuclide releases from real sites throughout the
country. The relationship between the mass released and the concentration of the radionuclide in the recharge
water from the waste source is given by:
2 Duration (d) 1
Mass (mg/cm ) = Recharge (cm/yr) • • • Concentration (mg/L),
365 (d/yr) 1000 (cm'/L)
= (8.7) (fl)^) Concentration (mg/L), (5-2)
= 0.024 (Concentration, mg/L),
44
-------
Table 5-4. Base Values of Input Parameters for Unsaturated Zone Radionuclide Models (from Wierenga, et
al, 1991; Gee, et al., 1994; U.S. EPA, 2000ab; and U.S. EPA VISITT Database).
Parameters
Values
Source-Specific Parameters
Area of disposal facility (m2)
Width of disposal facility (m)
Length of disposal facility (m)
Mass release of Radionuclide 238U (mg/cm2)
Concentration of 238U in recharge water from waste source (mg/L)
Mass Release of Radionuclide "Tc (mg/cm2)
Concentration of "Tc in Recharge Water from Waste Source (mg/L)
Mass Release of Radionuclide 90Sr (mg/cm2)
Concentration of 90Sr in Recharge Water from Waste Source (mg/L)
Mass Release of Radionuclide 238Pu (mg/cm2)
Concentration of 238Pu in Recharge Water from Waste Source (mg/L)
Mass Release of Radionuclide 3H (mg/cm2)
Concentration of 3H in Recharge Water from Waste Source (mg/L)
Duration of Waste Source Being Completely Released (days)
Potential Recharge Rate (mm/yr)
Initial Water Content (cm3/cm3)
Soil Properties in Unsaturated Zone
Saturated Hydraulic Conductivity, Ks, (cm/d)
Porosity (-)
Residual Water Content (cm3/cm3)
Saturated Water Content (cm3/cm3)
Bulk Density (g/cm3)
van Genuchten Alpha Coefficient, a, (cm-1)
van Genuchten Beta Coefficient, p, (-)
Depth to Water Table (m)
Solute Transport Parameters
Decay Coefficient for Parent Product "Tc (1/d)
Decay Coefficient for Daughter Product "Ru (1/d)
Distribution Coefficient for "Tc (ml/g)
Distribution Coefficient for "Ru (ml/g)
Dispersivity (cm)
Diffusion Coefficient in Free Water (cm2/d)
Apparent Molecular Dispersion Coefficient (cm2/d)
Dispersion Coefficient (cm2/d)
400
20
20
10
417
3 x 10-4
1.25 x 10-2
4.8 x 10-3
0.2
2.4 x 10-9
1.0 x 10-7
2.6 x 10-9
l.lxlO-7
1000
87
0.16
270.1
0.358
0.083
0.321
1.70
0.055
1.509
6
9 x 10-9
Stable
0.007
5.0
4.53
1.73
0.33
1.01
or
Concentration (mg/L) = 41.7 (Mass, mg/cm2).
Using Equation (5-3) and the mass releases of 238U, "Tc, 90Sr, 238Pu and 3H in Table 5-4 results in the
corresponding concentrations of these species given in the table.
(5-3)
In order to keep the initial concentration and the total amount of the radionuclide entering the soil fixed for
varying recharge rates, q, in the sensitivity analyses, the duration of the source emissions in Equation (5.2) was
45
-------
allowed to vary. That is, the following product in Equation (5-2) was held fixed at its base value as the time duration
of the source varied with the recharge rate, q:
(Duration) x (Recharge Rate) = 8700 cm-d/yr. (5-4)
For the range of q used in the sensitivity analyses, 5.11 cm/yr < q < 10.95 cm/yr, the source duration ranges over the
interval 795 d < duration < 1703 d. However, for the species considered in this report, the differences in source time
duration have a negligible effect on the output values of concern in this report. That is, since the total mass of a
radionuclide and its initial concentration are held fixed, the differences produced by the total release times (i.e., the
pulse width of release) have sufficient time to smooth out in the soil column before the major parts of the
breakthrough curves (BTCs) are seen at the bottom, the 6 m depth, of the soil columns.
Table 5-3 lists the depth to water table at the Las Graces Test Site as 60 m. However, the detailed soil moisture
data are given only for the first 6 m, as listed in Table 5-2. Thus, for the hypothetical modeling scenario of this
report, we chose the 6 m depth to be the top of the water table in our sensitivity analyses. It is felt that this
assumption will meet the project's objectives. As stated in the previous subsection, this 6 m layer is taken to be
homogeneous with the soil properties listed in the all-layer row of Table 5-2, namely: 9S = 0.321, 9r = 0.083, VG-a =
0.055 cm-1, VG-(3 = 1.509, andKs = 270.1 cm/d.
Wierenga, et al. (1991) found that the bulk densities of the nine soil layers in Table 5-2 range in values from
1.66 to 1.74 g/cm3, thus giving a geometric mean of the end points of 1.70 g/cm3. In Table 5-3, Gee et al. (1994)
stated that the typical soil type at the Las Graces Trench Site is a Berino fine loamy sand. With these two pieces of
information and two "rales of thumb" used by Eagleson (1970), we can roughly determine the porosity (n) and the
effective porosity (rig) at the site. The following "rales of thumb" were used by Eagleson to analyze the properties of
a Touchet silt loam (ps = 2.60 g/cm3), a Columbia sandy loam (ps = 2.67 g/cm3), and an unconsolidated sand (ps =
2.71 g/cm3), where ps is the density of the solid matrix and p is the bulk density:
n=l-p/ps, (5-5)
ne = n-9r. (5-6)
Thus, using a density, ps, of 2.65 g/cm3 for Berino fine loamy sand and a bulk density of 1.70 g/cm3 gives a porosity
of 0.358 and an effective porosity of 0.275 for the all-layer soil column in Table 5-2.
The radionuclide used in the sensitivity analyses reported in Section 6 and 7 is technetium ("Tc). This species
possesses a rather long half-life and is highly mobile in the soil column, as seen by its default value for Kd (0.007
ml/g, as given in U.S. EPA 2000b). The decay coefficient for "Tc in Table 5-4 is derived from the radionuclide half-
life given in Figure 5-1 The daughter product of "Tc (i.e., ruthenium, "Ru) is stable, and its decay coefficient is
zero. The default values of Kd for "Ru is 5 ml/g, which indicates that this radionuclide is not very mobile in the soil
column. Because of "Tc's high mobility and long half-life, it possesses many of the characteristics of a conservative
species as it moves through the soil column. Conversely, a species such as plutonium (238Pu), with a shorter half-life
and a low mobility, may decay before reaching a receptor if the soil column is sufficiently long and facilitated
transport does not exist.
The dispersion coefficient, D, in the unsaturated zone is given by Hills, et al. (1991) as (also see Equations 2-22
and 2-23):
D=TwDw+DLq/0, (5-7)
where 9 is taken as the initial water content of 0.16 and |q| is the infiltration rate of 8.7 cm/yr. Taking the
dispersivity of 4.53 cm given for tritium transport at the Las Graces Site by Porro and Wierenga (1993), and the
above values of 9 and |q|, gives a value of 0.68 cm2/d for DL |q| 0"1 . The diffusion coefficient in free water is
assumed to be 1.73 cm2/d and the tortuosity factor TW is taken as 0.19 (Tomasko, etal., 1989). Thus, the value of
TWDW in Equation (5-7) is given as 0.33 cm2/d. Consequently, the sum of the two terms in Equation (5-7) is given as
1.01 cm2/d, the value of the dispersion coefficient in Table 5-4.
The HYDRUS and CHAIN 2D Codes have the capacity for accounting for water uptake by plant roots,
while the other three models do not have this feature (see Table 2-1). Therefore, root water uptake was only
considered in Appendix F to show its impact on radionuclide movement using the HYDRUS Code. As we have
46
-------
previously stated, the creosote bush (Larrea tridentata) is the dominant plant at the Las Graces Site. Gile, et al.
(1998) indicated that the root depth system of this plant varies with the soil environment and the slope of the terrain;
the depth of roots can extend 5 m, or so. Jenkins, et al. (1988) reported that the creosote bush roots at the Las Graces
Site have a vertical distribution over a range of 0.5 mto 3 m. Thus, for demonstration purposes, a root distribution
of 0.5 m to 2.5 m was used in the current study (see Appendix F).
47
-------
Section 6
Parameter Sensitivity Analysis: Implementation and Results
The basic elements of the parameter sensitivity analysis for the five models considered in this report are
illustrated in Section 4 through the use of a simple model, y = F(a,b,c,d,e,f). In the current section, these procedures
will be implemented for the HYDRUS Code as applied to the hypothetical modeling scenario described in Section 5.
In what follows, we consider the general procedures for parameter sensitivity analysis as applied to the five models
under investigation, the input parameters for constant recharge rate and constant water content, the input parameters
for constant recharge rate and variable water content, the output variables of interest to SSL application, and the
sensitivity analysis results for the HYDRUS Code.
6.1 General Procedures for Parameter Sensitivity Analysis
The procedures for conducting the parameter sensitivity analysis reported in this section and in Section 7 are
analogous to those given in Section 4.2 (Equations 4-9 to 4-13) for the simple model, y = F(a,b,c,d,e,f). For the
current case, the input/output relationships are represented by expressions of the following form:
Qij=Fij(Pbl>Pb2>...,PpJ...>Pbll)> (6-1)
where Qy represents the j* output for the ith model, Pbk is the base value of the kth parameter, Pp is the variable p*
parameter, and Fy- (-,-, ...,•) is the functional relationship between the n-1 base parameters and the one variable
parameter determined by the 1th model for the j* output. This functional relationship represents empirical physical
laws (e.g., Darcy's Law), conservation laws (e.g., the continuity equations) and constitutive equations (e.g., van
Genuchten soil parameter model).
The specific procedures for conducting a parameter sensitivity analysis for the current SSL application are as
follows:
1. Select a radionuclide from the list: 3H, "Tc, 238U, 90Sr, 238Pu;
2. Select a model for analysis from the list: CHAIN (i = 1), MULTIMED-DP 1.0 (i = 2), FECTUZ
(i = 3), CHAIN 2D (i = 4), HYDRUS (i = 5);
3. Choose the set of input parameters for the specific model selected in Step 2;
4. Select the model outputs of interest, j = 1,2,. ..,k;
5. Select a particular variable input parameter Pp and select its domain of variation,
Ap
-------
6D
at e + pKd
-JAC,
(6-3)
where p is the bulk density of the soil (g/cm3), c is the concentration of the liquid phase (mg/L), Kd is the distribution
coefficient (ml/g), D is the dispersion coefficient (cm2/d), and (J, is the decay rate (1/d) which is assumed fixed for a
given radionuclide. That is, for a given radionuclide, the decay rate, (J,, is the same for the sorbed phase and the
dissolved phase, and |A is sufficiently well known that no sensitivity analysis will be run on it. As in Equation (5-7),
9D can be written as:
= 0TwDw+DL|q
(6-4)
where Tw is a dimensionless tortuosity factor taken as 0.19, Dw is the molecular diffusion coefficient in free water
(cm2/d) and DL is the longitudinal dispersivity of the radionuclide (cm). Using the base values given in Table 5-4,
9iw Dw is about one-half of the value of DL |q|. The reason that the diffusion component of D is so important is
because of the low annual average recharge rate represented by |q|.
The denominators in the second and third terms in Equation (6-3) represent the product of the soil moisture and
the retardation factor R (unitless):
OR = 6 + pKd .
For the default values of Kd listed in Section 5 and the base values of 9 and p given in Table 5-4, we get the
following ratios, pKd + 9R, for the five radionuclides:
(6-5)
Radionuclides
Tritium (3H)
Technetium ("Tc)
Uranium (238U)
Strontium (90Sr)
Plutonium (238Pu)
0.000
0.012
0.680
1.700
8.500
OR
0.160
0.172
0.840
1.860
8.660
TpKd •*• 6R1 100
0.00%
6.98%
80.95%
91.40%
98.15%
These ratios show the relative importance of the distribution coefficient Kd to the product 9R as Kd increases. For
example, p and Kd are relatively unimportant in the transport and fate of technetium (99Tc), but very important for
that of plutonium (238Pu).
Equations (6-3) to (6-5) indicate how the seven input parameters (Kd,q,9,p,D,DL,Dw) influence the transport and
fate of radionuclides in the unsaturated zone when q and 9 are assumed to be constant input parameters. For
comparison purposes, all five models were analyzed under these constant q and 9 conditions for sensitivity to
variations of these seven parameters. The boundary conditions for the liquid phase concentration, c, were the source
term at the top of the 6m column, as described in Section 5.5, and a zero gradient at the bottom of the column.
However, the models do not all possess dispersion coefficients of the type shown in Equation (6-4). In the CHAIN
Code, only a constant D is specified, while MULTDVIED-DP 1.0 and FECTUZ only specify DL. CHAIN 2D and
HYDRUS specify both DLand Dw. The sensitivity analyses conducted for these five models for the seven input
parameters under the constant q and 9 assumptions are summarized in Table 6-1.
6.3 Input Parameters for Constant Recharge Rate, but Variable Water Content
For a constant recharge rate, q, and a variable water content, 9, in a vertical soil column, the governing equations
are Darcy's Law,
q = -K
(6-6)
and the equation of continuity,
49
-------
Table 6-1. The Sensitivity Analyses Performed (•) for the Five Models Under the Assumption of Constant
Recharge Rate and Constant Water Content.
Model
Input
Parameter
Sensitivity Analyses Performed
MULTIMED-
CHAIN DP 1.0 FECTUZ CHAIN 2D
HYDRUS
q
e
p
D
DL
D,,,
(6-7)
where q is in units of (cm/d), K is the unsaturated hydraulic conductivity (cm/d), h is the soil water pressure head
(cm), the z-coordinate is positive upward with the origin at the bottom of the layer, and S is a sink term due to plant
water uptake (cm3/cm3d). Combining Equations (6-6) and (6-7) leads to the modified Richards Equation:
-S
The van Genuchten Model for the parameters K(h) and 9(h) is as follows (see Appendix A):
K(h) = KsKr(h)
(6-8)
(6-9)
0(h) =
e, -e
^r, h<0
es, h>o,
(6-10)
where the effective water content Se is given by
S. =
e - e '
(6-11)
and where Ks is the saturated hydraulic conductivity (cm/d), Kr is the unitless relative hydraulic conductivity, 9S is
the saturated water content, 9r is the residual water content, a is the inverse of the air-entry value or bubbling
pressure head (cm-1), (3 is the unitless pore-size distribution index, and m is defined by
m = 1 - 1/P .
(6-12)
50
-------
For h < 0, Equation (6-10) can be rearranged to give h as a function of the effective water content Se:
For the parameter sensitivity results in this section and Section 7, it is assumed that the sink term, S, due to plant
water uptake is identically zero. Thus, from the continuity equation and from the assumption of constant recharge
rate, q, the water content, 9, is found to be independent of time. When 9 is independent of time, the modified
Richards Equation reduces to Darcy 's Law in Equation (6-6), and the water content becomes a space-varying
quantity, 9(z), through Equation (6-10).
As indicated in Table 2-1, the flow modules of MULTIMED-DP 1.0 and FECTUZ are based on Darcy's Law,
while those for HYDRUS and CHAIN 2D are based on Richards Equation. As stated above, the parameter
sensitivity analysis for the van Genuchten parameters (Ks,9s,9r,a,B ) are based on a constant recharge rate, q, and a
variable water content, 9(z). The determination of 9(z) for MULTIMED-DP 1.0 and FECTUZ follows from the
solution of Darcy 's Law, which is a nonhomogeneous, nonlinear, first-order ordinary differential equations in the
water pressure head, h(z). To solve this system, the head has to be specified at the bottom of the soil column (z = 0)
and at the top of the column (z = L). These boundary conditions for h are determined from the specifications of 9 at
the two boundaries, namely:
(1) 9= 9S at z = 0, giving h=0 at z = 0;
(2) 9 = Base value in Table 5-4 at z = L, giving a value of h at z = L derived from Equations (6-11) and (6-13).
Once h(z) is known, 9(z) can be determined from Equation (6-10). Given 9(z), the influence on the liquid phase
concentration, c, is experienced through Equation (2-15), of which Equation (6-3) is a simplified version. As the five
input parameters (Ks,9s,9r,a,B) are varied, one at a time, for the sensitivity analysis, different water content profiles,
9(z), are generated leading to different responses for the breakthrough curves (BTCs) of the liquid phase
concentrations.
The determination of 9(z) for HYDRUS and the one-dimensional version of CHAIN 2D is realized from the
steady -state solution of the Richards Equation. Using the expressions for K and 9 in Equations (6-9) and (6-10),
respectively, Richards Equation is a homogenous, nonlinear, partial differential equation in h, requiring both
boundary conditions and an initial condition. The boundary conditions for h are obtained in the same way as those
for MULTIMED-DP 1.0 and FECTUZ. The initial condition is specified as that value of h corresponding to the base
value of 9 given in Table 5-4, as calculated from Equation (6-13). Thus, the water content profile, 9(z), is given as
the steady-state solution of this initial, boundary value system for h(z,t), or 9(z,t), as given by equation (6-10).
Again, as the van Genuchten parameters (Ks,9s,9r,a,B) are varied, one at a time, for the sensitivity analysis, different
profiles, 9(z), are generated and the resultant BTCs for c are varied. Table 6-2 summarizes the sensitivity analyses
performed for the four models (CHAIN is not included) under the assumption of constant recharge rates, q, and
variable water content, 9(z).
Table 6-3 lists the base values of the twelve input parameters (Kd,q,9,p,D,DL,Dw,Ks,9s,9r,a,B), along with the
range of variation for each of the twelve parameters. The ranges of variation of the twelve parameters were chosen
to be reasonably consistent with the variations shown over the nine-layer soil column at the Las Graces Trench Site
given in Table 5-2, the hydrologic records shown in Figures 5-2 and F-l, and the physical properties of the Berino
fine loamy sand at the site. These base values and ranges were used in the parameter sensitivity analyses for all five
models being analyzed in this report.
6.4 Output Parameters Evaluated
The time evolution of the radionuclide concentration in the leachate at the bottom of the unsaturated zone (i.e.,
at the entry point to the ground-water zone) is of interest in the Soil Screening Level (SSL) process. The curve of
concentration versus time at this entry point is called the breakthrough curve (ETC). For sensitivity analysis, three
characteristics of these curves are of importance:
1 . The peak concentration of the radionuclide, Cpeak,
2. The time to reach peak concentration, Tpeak, and
3 . The time when the concentration of radionuclide is high enough so that the resulting concentration
at a receptor well will exceed the MCL (i.e., the time to reach MCL, denoted by TMCL).
51
-------
Table 6-2. The Sensitivity Analyses Performed (•) for Four of the Five Models Under the Assumption of Constant
Recharge Rate and Variable Water Content, where "Base" Represents the Base Parameter Values
Given in Table 5-4, and the Water Content Profile, 9(z), Varies with the Changing van Genuchten
Parameters (Ks,9s,9r, a, B).
Sensitivity Analyses Performed
Model
Input MULTIMED-
Parameter CHAIN DP 1.0
Kd Base
q Base
9 9(z)
p Base
D
DL Base
DW
9°
a
B
Table 6-3. Input Parameters for Each Model and the
Value of the Parameter.
Parameter, Pp
p Modela Symbol
1 1,2,3,4,5 Kd
2 1,2,3,4,5 q
3 1,2,3,4,5 9
4 1,2,3,4,5 p
5 1 D
6 2,3,4,5 DL
7 4,5 Dw
8 2,3,4,5 Ks
9 2,3,4,5 9S
10 2,3,4,5 9r
11 2,3,4,5 a
12 2,3,4,5 P
a 1 = CHAIN, 2 = MULTIMED-DP 1.0, 3
b B2 = 0.032 for CHAIN;
c B3 = 0.26 for CHAIN;
d A, = 1 60 fnr CHAIN-
FECTUZ CHAIN
Base
Base
9(z)
Base
Base
Base
Base
9(z)
Base
Base
Base
I
•
•
2D HYDRUS
Base
Base
9(z)
Base
Base
Base
I
•
•
Range of Each Variable Parameter, Along with the Base
Units
ml/g
cm/d
cm3/cm3
g/cm3
cm2/d
cm
cm2/d
cm/d
cm3/cm3
cm3/cm3
I/cm
1/1
= FECTUZ, 4 =
AP R
p rbp "p
0.001
0.016
0.10
1.62d
0.40
3.73
0.93
175
0.29
0.06
0.051
1.43e
= CHAIN 2D, 5 =
0.007 0.019
0.024 0.030b
0.16 0.22C
1.70 1.78
1.00 2.20
4.53 5.33
1.73 2.53
270 365
0.32 0.35
0.08 0.10
0.055 0.059
1.51 1.59e
HYDRUS;
A12 = 1.45 and B12 = 1.57 for FECTUZ
52
-------
In the SSL process, the concentration at a receptor well is assumed to exceed the MCL whenever the
concentration at the ground water entry point exceeds the MCL times a dilution attenuation factor (DAF). A value of
20 for DAF is proposed in U.S. EPA(2000ab), and this is the value used in this report. Thus, the three output
quantities that will be considered in the sensitivity analyses of this section and Section 7 are Cpeak, Tpeak and TMCL.
6.5 Sensitivity Results for the HYDRUS Model
As indicated by Tables 6-1 and 6-2, parameter sensitivity analyses for the HYDRUS Model were carried out for
the eleven input parameters (Kd,q,9,p,DL,Dw,Ks,9s,9r,a,B). For the first six parameters (Kd,q,9,p,DL,Dw) the
analyses were conducted under the constraints of constant recharge rates, q, and constant moisture content, 9. For
the van Genuchten parameters (Ks,9s,9r,a,B), the analyses were executed under the constraints of constant recharge
rate, q, and variable water content, 9(z). The BTCs for "Tc concentrations at the bottom of the homogeneous, 6
meter soil column were recorded, as given in Figure 6-la to 6-lk. Three BTCs for each parameter were obtained,
one for the base value of the parameter in question, Pp, one for the upper limit of the range in Table 6-3, Bp, and one
for the lower limit, Ap. Data from these BTCs were used to calculate the sensitivities and relative sensitivities of the
peak concentraitons at the depth of 6 meters to the eleven input parameters, see Figures 6-2a to 6-2k; the sensitivities
and relative sensitivities of the times to reach peak concentration at the depth of 6 meters to the eleven input
parameters, see Figures 6-3a to 6-3k; and the sensitivities and relative sensitivities of the time to exceed MCL to the
eleven input parameters, Figures 6-4a to 6-4k. The relative sensitivity curves show the rates of increase (+) or
decrease (-) of the outputs, Cpeak, Tpeak, and TMCL, to increases in the eleven input parameters.
Table 6-4 summarizes the relative sensitivities of Cpeak, Tpeak, and TMCL with respect to the eleven basic input
parameters for HYDRUS, the relative sensitivities being measured at the base values of the input parameters. The
dominant input parameters for all outputs are those which have the greatest influence on the advection term in
Equation (6-3), followed by those which have the greatest influence on the diffusion term. The parameters, p and
Kd, have a very small influence on the transport and fate of "Tc since 9 » pKd in the coefficients of Equation (6-
3). As one would expect from the statements about Equation (6-4), the outputs for DL and Dw are nearly the same,
as are their sensitivities. Thus, the ratios of the relative sensitivities of the outputs to DL and Dw should be roughly
in the same proportions as the ratio of their base values, 4.53 + 1.73 = 2.62 , which is the case in Table 6-4 if one
allows for round-off error. For the first six input parameters, Cpeak is most influenced by 9 and the sum of DL and
Dw, while Tpeak, and TMCL are most influenced by q and 9. The direction of these influences (increasing/
decreasing) is seen by the algebraic signs in Table 6-4 and is easily justified by referring to Equation (6-3).
The last five input parameters listed in Table 6-4 are concerned with the van Genuchten Model for 9(h) and
K(h), as given in Equations (6-9) to (6-13). These five parameters influence the distribution of soil moisture in the
soil column, and thus influence both the advection and diffusion terms in Equation (6-3). Because of the low annual
recharge rate of 87 mm/yr, the least influential parameter is Ks, followed by the soil structure parameter a. The most
influential input parameter for all the outputs is the soil texture coefficient (3. The next most influential parameters
are the upper and lower moisture bounds of the van Genuchten Model, respectively 9S and 9r.
Table 6-4. Relative Sensitivities for Cpeak, Tpeak and TMCL with Respect to the Input
Parameters for HYDRUS, Measured at the Base Values of the Input Parameters.
Relative Sensitivity (a), Base Values of the Input Parameters
Parameters Cpeak Tpeak TMCL
Kd
q
9
p
DL
Dw
Ks
es
9r
a
P
-0.06
+0.12
-1.10
-0.05
-0.28
-0.10
+0.05
-0.66
-0.29
+0.13
+1.38
+0.06
-1.10
+0.80
+0.06
-0.02
-0.01
-0.04
+0.57
+0.22
-0.05
-1.00
+0.07
-0.98
+0.86
+0.07
-0.07
-0.03
-0.04
+0.64
+0.24
-0.05
-1.03
53
-------
(a) Sensitivity of "To Breakthrough to Kd
(b) Sensitivity of "To Breakthrough to q
(o) Sensitivity of "To Breakthrough to 8
|
Distribution Coefficient (ml/g) HYDRUS
0.001
0.007
0.019
Reoharae Rate (om/d)
0.016
0.024
0.030
2000 4000 6000 8000 10000
Time (days)
HYDRUS Water Content (onf/orf)
0.010
0.016
0.022
2000 4000 6000 8000 10000
Time (days)
2000 4000 6000 8000 10000
Time (days)
(d) Sensitivity of "Tc Breakthrough to p (e) Sensitivity of "To Breakthrough to DL (f) Sensitivity of "To Breakthrough to D.
0.015| 1 1 1 1—- 1 0.015| 1 1 1 1 1 0.015
0.010
f
HYDRUS Bulk Dimity (a/on?)
1.62
1.70
1.78
0.010
HYDRUS DliporjMty (am)
3.73
4.53
5.33
_ _
2000 4000 6000 8000 10000
Time date
_ _
2000 4000 6000 8000 10000
Time «ay»))
Diffusion Coefficient In Water (orf f")
HYDRUS 0.93
1.73
2.53
2000 4000 6000 8000 10000
Time (days)
(g) Sensitivity of "To Breakthrough to K,
(h) Sensitivity of "To Breakthrough to 6S
(0 Sensitivity of MTo Breakthrough to 8r
|
Saturated Conductivity (cm O HYDRUS
175
270
365
HYDRUS Saturated Water Content (orf/onr)
0.29
0.32
0.35
2000 4000 6000 8000 10000
Time (days)
HYDRUS Residual Water Content (orf/onr)
0.06
0.08
0.10
2000 4000 6000 8000 10000
Time (days)
2000 4000 6000 8000 10000
Time (days)
(I) Sensitivity of "To Breakthrough to a
(k) Sensitivity of "Tc Breakthrough to
van Genuohten Retention Parameter a (onT1)
HYORUS 0.051
0.055
0.059
2000 4000 6000 8000 10000
Time (days)
van Ganuchten Retantlon Parameter 0 (—)
HYDRUS 1.43
1.51
1.59
2000 4000 6000 8000 10000
Time (days)
Figure 6-1. Sensitivity of "Tc breakthrough (through the 6m layer) to the system parameters using the HYDRUS
Model: (a) distribution coefficient, (b) recharge rate, (c) water content, (d) bulk density, (e) dispersi-
vity, (f) diffusion coefficient in water, (g) saturated conductivity, (h) saturated water content, (i) residual
water content, (]) van Genuchten retention parameter a, (k) van Genuchten retention parameter (3.
54
-------
u
£ -0.10
'w
c -0.20
(/>
-jj -0.30
ct
-0.40
'^~~*~--^-___' ' (a)
^^^^^~~^~-~^
•— — .
-
-
i i i
u
£" -0.50
~
c -1.00
(fl
-i -i-so
ct
-2.00
(c)
,.
^^
N,
N.
\
\
!
n nn rt rtrt
— u.uu
-0.02
t" -0.04
'» -0.06
c
0}
1/5 -0.08
-0.10
-
__^__^— *— ' ° * " _
-
1 1 1
-0.02
£" -0.04
'« -0.06
c
03
M -0.08
-0.10
-
^^^
^^^^
~
I I
^ 0.020
\
|> 0.015
£ 0.010
o
.* 0.005
D
D_
0
HYDRUS "Tc
-
i n
D D
1 1 1
^_^ u.uzu
\
I1 0.015
£ 0.010
o
^ 0.005
a
cp
Q_
0
HYDRUS "Tc
-
.
^^^^
^^^.^
1 1
0 0.005 0.010 0.015 0.020 0 0.10 0.20 0.30
Distribution Coefficient (ml/g) Water Content (cm3 cm"3)
OF- r, r* * n
.DU
£. 0.40
£ 0.30
w
c
$ 0.20
£ 0.10
0
(b)
-
-
*^^\
i i i i
U. 1 U
-£- °
1 -0.10
c
& -0.20
^ -0.30
-0.40
(d)
-
-
1 1 1
0.15
£•
£ 0.10
'w
C
w 0.05
0
-
\.
1 1 1 1
U.UU 1 U
0.0006
i" 0.0002
• —
"w -0.0002
c
1/5 -0.0006
-0.0010
-
_
_
-
1 1 1
Oi-i i-ij-i ^ non
.UZU
\
|> 0.015
^ 0.010
o
^ 0.005
a
0}
a_
0
HYDRUS "Tc
-
i i i i
^^
-------
£" -0.10
^>
" -0.20
CD
(A
•3 -o.3o
tx.
-0.40
-0.0002
£" -0.0004
>
| -0.0006
to
1/1 -0.0008
-0.0010
I
I1 0.015
" 0.010
o
0
.* 0.005
D
CU
Q_
0
3
i" -0.10
^>
" -0.20
CU
(ft
-i -0.30
-0.40
-0.0002
i" -0.0004
_>
| -0.0006
o
1/1 -0.0008
-0.0010
I
|> 0.015
c 0.010
0
0
.* 0.005
D
03
0
0
(e)
" —~~~— "
-
HYDRUS "Tc
-
.0 3.5 4.0 4.5 5.0 5.5 6
Dispersivity (cm)
•~^-— (0
: - — - ;
HYDRUS "Tc
-
.5 1.0 1.5 2.0 2.5 3
£ 0.05
>
1 0
CD
CO
i, -o.os
ct
-0.10
0.0006
~T 0.0002
>
« -0.0002
CD
1/1 -0.0006
-0.0010
I
1" 0.015
c 0.010
o
o
^ 0.005
D
0>
D_
0
.0 C
i" -0.200
^>
" -0.400
o
(fl
•5 -0.600
-0.800
-0.005
>-
~ -0.010
'w
c
w -0.015
-0.020
I
g1 0.015
c 0.010
0
o
^ 0.005
D
0)
0
.0 0.
(g)
i i i i
-
HYDRUS "Tc
1 1 1 1
) 100 200 300 400 5C
Saturated Conductivity (cm d~1)
(h)
1 1 1
;___—;
HYDRUS "Tc
1 1 1
280 0.300 0.320 0.340 0.3
Diffusion Coefficient in Water (cm2 d~')
Saturated Water Content (cm3cm~3)
Figure 6-3. Sensitivity and relative sensitivity of peak concentrations at the depth of 6m to the system parameters
using the HYDRUS Model: (e) dispersivity, (f) diffusion coefficient in water, (g) saturated
conductivity, (h) saturated water content.
56
-------
0.000
-0.100
-0.200
-0.300
-0.400
-0.500
0)
-0.010
-0.020
-0.030
-0.040
0.020
0.015
0.010
0.005
0
0
-
-
I 1 1 1 I
HYDRUS "Tc
-
1 1 1 1 1
05 0.06 0.07 0.08 0.09 0.10 0.
Residual Water Content (cm3 cm"3
0.200
*• 0.150
1 0-100
1J 0.050
o:
0
0.020
0.015
~ 0.010
0
^ 0.020
_i
I1 0.015
c 0.010
o
o
.* 0.005
HYDRUS
l9Tc
2.00
^ 1.50
^>
1 1.00
CD
)
-i 0.50
cc
0
0.010
0.008
^ 0.006
t>
| 0.004
^ 0.002
0
^ 0.020
_i
g1 0.015
^ 0.010
o
o
^ 0.005
HYDRUS
'9Tc
1.40 1.45 1.50 1.55 1.60
van Genuchten Retention Parameter J3 (—)
0.050 0.052 0.054 0.056 0.058 0.060
van Genuchten Retention Parameter a (cm"1)
Figure 6-2. Sensitivity and relative sensitivity of peak concentrations at the depth of 6m to the system parameters
using the HYDRUS Model: (i) residual water content, (j) van Genuchten retention parameter a, (k)
van Genuchten retention parameter (3.
57
-------
£ 0.15
>
'w
s= 0.10
(D
(/)
i 0-05
ct
0
50000
40000
t" 30000
'| 20000
cu
1/1 10000
o" 8000
cj 6000
c
o
0
4000
o
o
°- 2000
o
CD 0
E C
\—
0
£ -0.5
>
'«
c -1.0
0>
c/o
^ "1-5
-2.0
0
-100000
i" -200000
'1 -300000
CD
^ -400000
W
o~20000
o 15000
o
° 10000
D
CD
°- 5000
o
CD 0
(a)
^^"^
^^
-
\ \ \
HYDRUS "Tc
;_-^— — — ;
_>. 0.80
*>
| 0-60
co 0.40
^ 0.20
0
50000
40000
£" 30000
'« 20000
CD
1/1 10000
cT 8000
-------
£. -0.02
1 -°-04
c
& -0.06
£ -0.08
-0.10
0
-20
t" -40
> —
« -60
OJ
1/1 -80
V -100
o" 8000
-a
6 6000
c
o
o
^ 4000
a
01
°- 2000
o
CD 0
£ 3
\—
0
s" -0.005
" -0.010
0)
i, -0.015
-0.020
0
-20
t" -40
'| -60
CD
1/1 -80
V -10°
i? 8000
d 6000
c
o
o
4000
D
°- 2000
o
. 0.80
1 0.60
c
M 0.40
£ 0.20
0
20000
15000
-fc-
p 10000
'w
(/> 5000
0
'S' 8000
a
. 6000
o
c
0
O
^ 4000
0
a_
o 2000
1 0
(g)
-
-
~s^
j^°^
HYDRUS "Tc
1 1 1 1
) 100 200 300 400 50
Saturated Conductivity (cm d"1)
(h)
-— —
-
1 1 1
HYDRUS "Tc
; __— — ;
_
i i i
.1 0.5 1.0 1.5 2.0 2.5 3.0
"- Diffusion Coefficient in Water (crrrd"1)
0.280 0.300 0.320 0.340 0.360
Saturated Water Content (cm3 cm"3)
Figure 6-3. Sensitivity and relative sensitivity of time to reach peak concentrations at the depth of 6m to the system
parameters using the HYDRUS Model: (e) dispersivity, (f) diffusion coefficient in water, (g) saturated
conductivity, (h) saturated water content.
59
-------
£ 0.30
'in
c 0.20
CD
-jj 0.10
0
50000
40000
t" 30000
'| 20000
cu
1/1 10000
^ °
o" 8000
-o
cj 6000
c
o
0
^ 4000
o
°- 2000
o
a> 0
E o
\—
0.00
>- -0.05
>
=5 -°'10
c
<3 -0.15
^ -0.20
-0.25
0
-2000
i" -4000
'" -6000
1/1 -8000
~ 10000
o" 8000
tj 6000
o
0
4000
CD
°- 2000
o
CD 0
(0
-
^-^^
_
-
-
1 I I I 1
HYDRUS "Tc
; „_- ;
_
05 0.06 0.07 0.08 0.09 0.10 0.
£" -0.50
'in
c -1.00
(D
(/)
•5 -1-50
-2.00
0
-2000
£" -4000
'| -6000
^ -8000
-^10000
D" 8000
-a
-------
u.^u
£" 0.15
£
'to
c 0.10
CD
(/)
-5 0.05
0
' ' ' ^ >
^^ ~
^^
^*^
/^
- /^
•^ i i
1 .UU
_>, 0.80
S 0.60
w
c
ff> 0.40
^ 0.20
0
°~—«— ' , .
^^^-^
-
-
! j
""nnnn i-j-n-n-n-i
DUUUU
40000
:5? 30000
'«> 20000
c
0)
1/1 10000
^ 0
-
-
-
1 1 1
JUUUU
40000
t" 30000
'" 20000
E=
CD
1/1 10000
^ 0
-
-
°~-«— »-^
-
!
« CO
o 8000
£ 6000
2
-c 4000
u
D
CD
^ 2000
o
CD 0
HYDRUS "Tc
-
o e
-
1 1 1
O ouuu
jj 6000
2
-s= 4000
u
a
CD
^ 2000
o
- 0.08
1 0-06
c
& 0.04
H 0.02
0
Hd)
- — — — -^-___^____^ "
. .,
-
-
01 nnn
-100000
t" -200000
^ -300000
c
CD
*" -400000
-500000
/
/^
/
-
iii
1 UVJU
800
^~ 600
'" 400
c
CD
" 200
^ 0
-
-
-
,
CO
Tonnnn ^* ««««
w zuuuu
>^
D
T3
^15000
(J
2 10000
^
o
D
£ 5000
o
- o
HYDRUS "Tc
_
_
"~^^— _
1 1 1 1
Q OUUU
•a
-i 6000
u
2
-c 4000
o
0
^ 2000
0
"*~
CD 0
HYDRUS "Tc
_
_
-
1 1 1
0 0.01 0.02 0.03 0.04 0.05 .E 1.60 1.65 1.70 1.75 1.80
.*- Po^hcii-.-iaPci + a/'j^m H~1 ^ ' Rnllynarioiti/^r-i / r* rrf5"!
Figure 6-4. Sensitivity and relative sensitivity of time to exceed MCL to the system parameters using the
HYDRUS Model: (a) distribution coefficient, (b) recharge rate, (c) water content, (d) bulk density.
61
-------
^ -0.10
^>
1 -0.20
| -300
o
<" -400
,-^ -500
w
-D
^J 6000
2
-n 4000
D
01
^ 2000
£
1 -0.10
QJ
(/)
^ -0.15
o:
-0.20
-50
>^
£ -100
*w
(?) -150
^ -200
«
~u
£ 6000
2
-n 4000
D
to
^ 2000
£
1 1 1 1
- _^—^ -
HYDRUS "Tc
iii
_>. -0.02
1 -0.04
(?! -0.06
^ -0.08
-0.10
-1.00
~ -2.00
en
(?) -3.00
^ -4.00
in
T3
-i 6000
o
2
-c 4000
D
0}
^ 2000
o
CD 0
.0 E (
f—
_>, 0.80
">
1 0.60
(0 0.40
£ 0.20
0
15000
~ 10000
(?) 5000
^ 0
t/i
~o
-i 6000
o
2
-c 4000
D
CU
^ 2000
o
o) 0
(g)
-!—---• '
-
HYDRUS "Tc
1 1 1 1
) 100 200 300 400 50
Saturated Conductivity (cm d"1)
(h)
1 1 1
-
HYDRUS "Tc
1 1 1
0.5 1.0 1.5 2.0 2.5 3.0
Diffusion Coefficient in Water (cm2 dH)
0.280 0.300 0.320 0.340 0.360
Saturated Water Content (cm3 cm"3)
Figure 6-4. Sensitivity and relative sensitivity of time to exceed MCL to the system parameters using the
HYDRUS Model: (e) dispersivity, (f) diffusion coefficient in water, (g) saturated conductivity, (h)
saturated water content.
62
-------
±r o.3o
1 0.20
CD
15 °-10
C£
0
20000
15000
-£•
~ 10000
c
c/) 5000
^ 0
g- 8000
"D
-i 6000
o
2
-c 4000
o
0
CD
K 2000
0
CD 0
E o.
i-
0
£" -0.05
1 -0.10
CD
(fl
-0.20
0
-2000
~ -4000
_>
« -6000
"> -8000
-10000
V 8000
D
^ 6000
O
2 4000
_c
u
J 2000
0
- 0
(i)
^^^
- ^^^
_
1 1 1 I I
^ .
3 *
1 1 1 1 1
HYDRUS "Tc
1 1 1 1 1
t" -0.50
1 -1.00
CD
CO
-5 -1-50
-2.00
0
-2000
±T -4000
~
» -6000
cu
1/5 -8000
-10000
"£ 8000
D
^ 6000
_i
o
S 4000
^
o
S 2000
Q^
0
- o
' (k)
-
- ___— - — -
_
- _— — — — -
_
-
HYDRUS 99Tc
— — '
i i
0500.060 0.070 0.080 0.090 0.100 0.110 "> 1.400 1.450 1.500 1.550 1.600
Residual Water Content (cm3cm~3) ~ van Genuchten Retention Parameter jS ( — )
- ^—-——®-
1 1 I I
_
- __^---^^ -
-
1 1 1 1
HYDRUS 99Tc
-
-
iii
<" 0.050 0.052 0.054 0.056 0.058 0.060
~ van Genuchten Retention Parameter a (cm~1)
Figure 6-4. Sensitivity and relative sensitivity of time to exceed MCL to the system parameters using the HYDRUS
Model: (i) residual water content, (j) van Genuchten retention parameter a, (k) van Genuchten retention
parameter (3.
63
-------
Section 7
Comparison of Sensitivity Results Between Models:
Illuminating Numerical Differences
This section is concerned with the comparison of the parameter sensitivity results derived for the five models
being studied in this report. The input parameters being considered form subsets of the twelve parameters
(Kd,q,9,p,D,DL,Dw,Ks,9s,9r,a,p) described in Sections 6.2 and 6.3 The output parameters (Cpeak,Tpeak,TMCL) are
those introduced in Section 6.4. In what follows, we consider the modeling codes and their differences, the
parameter sensitivity results, and other results illuminating numerical differences and errors between the models.
7.1 The Modeling Codes and Their Differences
Table 2-1 lists the flow and transport modeling components, and the type of boundary conditions available to
each of the five codes being tested in this report. From this table, the depth and breadth of application of each code
can be determined. However, for the comparison of parameter sensitivity given in this section, we consider a
constant flow (i.e., recharge rate) through a 6m, vertical homogeneous soil column, with a constant or spatially
varying water content. The boundary conditions for the transport of the radionuclide, "Tc, consist of a stepwise
constant source at the top of the soil column and a zero concentration gradient at the bottom of the column.
As we have previously stated, the CHAIN Code is an analytical transport model which assumes uniform flow
conditions in the homogeneous soil column; that is, the soil moisture and recharge rate are both constant. Even
though CHAIN possesses an analytical solution for the pollutant breakthrough curves (BTCs) at the 6m level of the
soil column, numerical methods are required to generate the graphs of the BTCs from the closed-form analytical
expressions. Thus, the CHAIN Code can exhibit numerical error, with reference to the exact analytical expressions,
due to truncation of the series of functions representing the analytical solution in the numerical approximation
method, and due to computer round-off error.
The MULTIMED-DP 1.0 and FECTUZ Codes can handle constant recharge rates with constant or spatially
varying water contents. With respect to the BTCs for pollutant transport, these codes are semi-analytical models
which approximate the inversion of complex analytical solutions in transform space using numerical inverse
transform modules. Thus, the numerical errors most prevalent in these codes are those due to the numerical solution
of Darcy's Law for the pressure head (leading to a variable water content in the soil column) and those due to the
numerical inversion (i.e., a numerical Laplace inversion) of the analytical solutions for the pollutant concentration in
the Laplace transform space. If these two codes use the same differencing scheme or one of comparable accuracy for
solving Darcy's Law and the same or comparable numerical inversion schemes for the pollutant BTCs, then one
would expect that the BTCs for a given boundary-initial value problem would be nearly the same for both codes.
The CHAIN 2D (the one-dimensional, vertical version) and the HYDRUS Codes handle constant recharge rates
with constant, spatially varying, or temporally and spatially varying water contents. With respect to pollutant
transport, these codes solve the advection-dispersion-decay partial differential equations using finite-difference or
finite-element methods. Thus, the numerical errors most prevalent in these two codes are those due to the numerical
schemes used to solve the Richards Equation for the temporally and spatially varying pressure head (leading to a
temporally and spatially varying water content), those due to deriving a spatially varying water content from the
numerically obtained steady-state solution of the Richards Equation, and those due to the numerical schemes used to
obtain the BTCs for pollutant transport from the advection-dispersion-decay partial differential equations. If CHAIN
2D and HYDRUS use the same numerical solution schemes or ones of comparable accuracy for the Richards
Equation and the advection-dispersion-decay equations, then one would expect (as with MULTIMED-DP 1.0 and
FECTUZ) that the BTCs for a given boundary-initial value problem would be nearly the same for both codes.
64
-------
From the above discussions, one would be tempted to place the derived BTCs into three categories for the
current parameter sensitivity modeling scenario (the scenario as described in Sections 5, and 6.1 to 6.3). These three
categories, based on the potential for numerical error, would be those BTCs derived by the CHAIN Code, those
BTCs derived by the MULTDVIED-DP 1.0 and FECTUZ Codes, and those BTCs derived by the CHAIN 2D and
HYDRUS Codes. As pointed out in Sections 6.2 and 6.3 (see Tables 6-1 and 6-2), there are two other discriminating
factors between the codes for the current modeling scenario; namely: (1) the specification of a constant or spatially
varying water content in the homogeneous, 6m soil column, and (2) the construction of the dispersion factor, 9 D, in
Equation (6-4). Table 7-1 is a composite of Tables 6-1 and 6-2, and points out the impact that these two
discriminating factors have on the current parameter sensitivity analyses.
For the first seven parameters (Kd,q,9,p,D,DL,Dw) in Table 7-1, the five models were run under the constraints
of a constant recharge rate, q, and a constant water content, 9, using the base input values and input ranges given in
Table 6-3. As previously stated and as indicated in Table 7-1, CHAIN specifies a constant value for D and does not
resolve 9D into components as given in Equation (6-4). MULTIMED-DP 1.0 and FECTUZ only specify constant
values for DL in Equation 6-4, while CHAIN 2D (the one-dimensional, vertical version) and HYDRUS specify
constant values for both DL and Dw. Since the recharge rate, q, in the current sensitivity modeling scenario is so
small, the quantity 9iwDw in Equation (6-4) is relatively important as compared to the dispersion term DL|q|. Thus,
for the base values listed in Table 6-3, one would expect less total dispersion effects (i.e., effects due to 9D) on the
BTCs derived by MULTIMED-DP 1.0 and FECTUZ than for those derived by CHAIN 2D and HYDRUS. The
fourth row of Table 7-2 gives the numerical justification of this statement. As this row of values indicates, the base
value of D in the CHAIN Code gives a value for 9D which is nearly the same as those for CHAIN 2D and
HYDRUS, while the value of 9D for MULTIMED-DP 1.0 and FECTUZ is about 2/3 of the values for the values for
Table 7-1. The Sensitivity Analyses Performed (•) for the Five Models Under the Assumption of Constant
Recharge Rate and Constant Water Content, and the Sensitivity Analyses Performed (O) for Four of
the Five Models Under the Assumption of Constant Recharge Rate and Variable Water Content.
Model Input
Parameter
Sensitivity Analyses Performed
z
i
OH
Q
Q
W
S)
_
<
u
q
e
p
D
DL
Ks
es
er
a
P
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
65
-------
Table 7-2. The Values of the Dispersion Factor, 9D, as Derived from Equation (6-4) and the Base Values of
D, DL, Dw, 9 and q and the Ranges of D, DL and Dw as given in Table 6-3.
D
cm2/d
0.40
1.00
1.00
1.00
1.00
1.00
2.20
DL
cm
4.53
3.73
4.53
4.53
5.33
4.53
4.53
Dw
cm2/d
1.73
1.73
0.93
1.73
1.73
2.53
1.73
Djql
cm2/d
0.109
0.090
0.109
0.109
0.128
0.109
0.109
exwDw
cm2/d
0.053
0.053
0.028
0.053
0.053
0.077
0.053
9D in cm2/d
S;
«)
0.064
0.160
0.160
0.160
0.160
0.160
0.352
1
1-
1:
g"
§!
0.109
0.090
0.109
0.109
0.128
0.109
0.109
M
£*4
[^
0.109
0.090
0.109
0.109
0.128
0.109
0.109
§
§
1
0.162
0.143
0.137
0.162
0.181
0.186
0.162
g
1
0.162
0.143
0.137
0.162
0.181
0.186
0.162
the other three models. Thus, one would expect the BTCs for the base values of the input parameters to be nearly the
same for the CHAIN, CHAIN 2D and HYDRUS Codes, and nearly the same for the MULTIMED-DP 1.0 and
FECTUZ Codes, where the peak concentrations for the latter two models should be higher than those for the former
three models. This is, in fact, the case as we will see in Section 7.2. The numerical values for 9D in the other six
rows of Table 7-2 can be used to interpret the various sensitivity results given in the next subsections.
For the last five input parameters in Table 7-1 (Ks,9s,9r,a,p), only four models (CHAIN does not allow variable
9) are evaluated for constant recharge rate, q, and variable water content, 9(z). Based on the above discussions, one
would expect the BTCs for CHAIN 2D and HYRUS to be nearly the same, and those for MULTIMED-DP 1.0 and
FECTUZ to be nearly the same, where the peak concentrations for the latter two models should be higher than those
for the former two models. As we will see below, this is the case.
7.2 The Parameter Sensitivity Results
As indicated in Table 7-1, parameter sensitivity analyses are performed on all five models for the first four
parameters (Kd,q,9,p). CHAIN is the only model for which sensitivity analyses are performed for D. Sensitivity
analyses are performed for the six parameters (DL,Ks,9s,9r,a,p) on four of the five models (all except CHAIN), while
only CHAIN 2D and HYDRUS consider the parameter Dw. The "Tc BTCs at the bottom of the 6m, homogeneous
soil column are given in Figures 7-1 to 7-12, where the breakdown is as follows:
Figures 7-1 to 7-4 compare the BTCs for each of the five models for the first four parameters (Kd,q,9,p),
respectively;
Figure 1-5 gives the BTCs for the CHAIN Code for the parameter D; this figure is given for the purpose of
completion;
Figure 7-6 compares the BTCs for each of the four models (CHAIN being excluded) for the dispersivity,
DL;
Figure 7-7 compares the BTCs for CHAIN 2D and HYDRUS for the diffusion coefficient in water, Dw;
Figures 7-8 to 7-12 compare the BTCs for each of the four models (CHAIN being excluded) for the five van
Genuchten parameters (Ks,9s,9r,a,p), respectively.
These twelve figures of "Tc BTCs fulfill the expectations set forth in Section 7.1. These results, or
expectations, are summarized in the following list:
66
-------
0.015
0.010
0.005
Distribution Coefficient
HYDRUS
CHAIN 2D
FECTUZ
o>
E
c
o
c
-------
0.015
0.010
0.005
Recharge Rate = 0.030 cm/d
HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP
CHAIN
0.015
0.010
"c 0.005
-------
0>
C
o
c
O)
0.015
0.010
0.005
OL
0
0.015
0.010
0.005
0.015
0.010
0.005
Water Content
HYDRUS
CHAIN 2D
FECTUZ
0.22 cnf/cnf
MULTIMED-DP
CHAIN
2000
4000
6000
8000
10000
Water Content = 0.16 cnf/cnf
HYDRUS MULTIMED-DP
CHAIN 2D CHAIN
FECTUZ
2000
4000
6000
8000
10000
Water Content
HYDRUS
CHAIN 2D
0.10 c n?/c rtf
MULTIMED-DP
CHAIN
FECTUZ
2000 4000 6000
Time (days)
8000
10000
Figure 7-3. Sensitivity of "Tc breakthrough (through the 6m layer) to the water content using the CHAIN,
MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS Models, where 9 = 0.22 cmVcm3, 0.16
cmVcm3, and 0.10 cmVcm3.
69
-------
U.UIO
0.010
0.005
0
C
On4 K.
.U1O
cr
n»
c Concentration (m<
o o
• •
o o
O -"•
o en o
s~ <
n rH «;
U.U1O
0.010
0.005
0
C
Bulk Density
HYDRUS
CHAIN 9D
_ FECTUZ
^*"v
/^
) 2000 4000
Bulk Density
HYDRUS
CHAIN 9D
_ FECTUZ
XN
<_ V
f\
) 2000 4000
Bulk Density
HYDRUS
CHAIN 9D
_ FECTUZ
XX
<•. v
/s
1 ^*^ 1
) 2000 4000
Time
= 1 .78 g/cn?
MULTIMED— DP
CHAIN
6000 8000 1 0000
= V.70 g/cn?
MULTIMFD DP
pUAIkl
6000 8000 1 0000
= V.62 g/cn?
MULTIMED— DP
CHAIN
l^feh^ l
6000 8000 1 0000
(days)
Figure 7-4. Sensitivity of "Tc breakthrough (through the 6m layer) to the bulk density using the CHAIN,
MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS Models, where p = 1.78 g/cm3, 1.70 g/
cm3, and 1.62 g/cm3.
70
-------
0.015
0.010
0.005
Dispersion Coefficient = 2.2 (cm2/d)
CHAIN
o>
E
c
o
c
o>
o
c
o
o
o
01
o>
0.015
0.010
0.005
2500
5000
7500
Dispersion Coefficient = 1.0 (cm2/d)
CHAIN
0.015
0.010
0.005
2500
5000
7500
Dispersion Coefficient = 0.4
CHAIN
2500 5000 7500
Time (days)
10000
10000
10000
Figure 7-5. Sensitivity of "Tc breakthrough (through the 6m layer) to the dispersion coefficient using the CHAIN
Model, where D = 2.20 cmVd, 1.00 cmVd, and 0.40 cmVd.
71
-------
0.015
0.010
0.005
o>
E
c
o
o>
0.015
0.010
0.005
0.015
0.010
0.005
HYDRUS
CHAIN 2D
FECTUZ
Dispersivity = 5.33 cm
MULTIMED-DP
J_
2000
4000
6000
8000
10000
HYDRUS
CHAIN 2D
FECTUZ
Dispersivity = 4.53 cm
MULTIMED-DP
2000
4000
6000
8000
10000
HYDRUS
CHAIN 2D
FECTUZ
Dispersivity = 3./3 cm
MULTIMED-DP
2000 4000 6000
Time (days)
8000
10000
Figure 7-6. Sensitivity of "Tc breakthrough (through the 6m layer) to the dispersivity using the MULTIMED-DP
1.0, FECTUZ, CHAIN 2D, and HYDRUS Models, where DL = 5.33 cm, 4.53 cm, and 3.73 cm.
72
-------
0.015
0.010
0.005
Diffusion Coefficient in Water
HYDRUS
o>
E
c
o
c
at
0.015
0.010
0.005
0.015
0.010
0.005
= 2.5I3 err? a"1
CHAIN 2D
2000
4000
6000
8000
10000
Diffusion Coefficient in Water
HYDRUS
CHAIN 2D
= 1 .3 en? a"1
2000
4000
6000
8000
10000
Diffusion Coefficient in Water
HYDRUS
CHAIN 2D
= 0.^3 err? of1
2000 4000 6000
Time (days)
8000
10000
Figure 7-7. Sensitivity of "Tc breakthrough (through the 6m layer) to the diffusion coefficient in water using the
CHAIN 2D and HYDRUS Models, where Dw = 2.53 cm2/d, 1.73 cm2/d, and 0.93 cm2/d.
73
-------
0.015
0.010
0.005
Saturated Conductivity
HYDRUS
CHAIN 2D
FECTUZ
o>
E
c
o
c
a>
o
c
o
o
o
o>
o>
0.015
0.010
0.005
0.015
0.010
0.005
T~Je5~
cm d"1
MULTIMED-DP
2000
4000
6000
8000
10000
Saturated Conductivity
HYDRUS
CHAIN 2D
FECTUZ
= 2*70
cm d"1
MULTIMED-DP
2000
4000
6000
8000 10000
Saturated Conductivity
HYDRUS
CHAIN 2D
FECTUZ
= i?5 cm d^
MULTIMED-DP
2000 4000 6000
Time (days)
8000
10000
Figure 7-8. Sensitivity of "Tc breakthrough (through the 6m layer) to the saturated conductivity using the
MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS Models, where Ks = 365 cm/d, 270 cm/d,
and 175 cm/d.
74
-------
c
o
c
0
o
c
o
o
u
o>
o>
U.UIO
0.010
0.005
0
C
Ort 4 K.
.U1O
0.010
0.005
C
On* c
.Ul O
0.010
0.005
0
C
Saturated Water Content
HYDRUS
CHAIN 2D
_ FECTUZ
f
. J
) 2000 4000
Saturated Water Content
HYDRUS
CHAIN 9D
_ FECTUZ
£ \
r
) 2000 4000
Saturated Water Content
HYDRLJ^
CHAIN 2D
_ FECTUZ
A
f\
, J,
) 2000 4000
Time ((
= 0.35 cnr/cnf
MULTIMED DP
-
\
\
6000 8000 1 0000
= 0.32 cnr/crrf
-
V. '
6000 8000 1 0000
= 0.29 cnr/cnf
uiji TIMFD DP
-
6000 8000 1 0000
days)
Figure 7-9. Sensitivity of "Tc breakthrough (through the 6m layer) to the saturated water content using the
MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS Models, where 9S = 0.35 cmVcm3, 0.32
cmVcm3, and 0.29 cmVcm3.
75
-------
U.UIO
0.010
0.005
0
C
0/"\4 C
.UlO
X
£ 0.010
c
o
D
"c 0.005
-------
0.015
0.010
0.005
o>
E
~c
o>
0.015
0.010
0.005
0.015
0.010
0.005
van Genuchten a = 0.059 ci
HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP
2000
4000
6000
8000
10000
van Genuchten a = 0.055 cm"1
HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP
2000
4000
6000
8000 10000
van Genuchten a = 0.051 ci
HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP
2000 4000 6000
Time (days)
8000
10000
Figure 7-11. Sensitivity of "Tc breakthrough (through the 6m layer) to the van Genuchten retention parameter a
using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS Models, where a = 0.059cm-1,
0.055 cm-1, and 0.051cm-1.
77
-------
U.UI3
0.010
0.005
0
C
0.015
a*
o>
^ 0.010
c
JO
2
"c 0.005
Q>
O
C
O
0
£ °
« C
w 0.015
0.010
0.005
0
van Genuchten 0 = 1.59
HYDRUS MULTIMED DP
CHAIN 2D
_ FECTUZ
/^
f V
) 2000 4000 6000 8000 10000
van Genuchten /? = 1.51
HYDRUS MLJLTIMED DP
CHAIN 2D
_ FECTUZ
/~\
£ V
) 2000 4000 6000 8000 10000
van Genuchten /? = 1.43
HYDRUS MULTIMED DP
CHAIN 2D
_ FECTUZ
£\
2000 4000 6000
Time (days)
8000
10000
Figure 7-12. Sensitivity of "Tc breakthrough (through the 6m layer) to the van Genuchten retention parameter (3
using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS Models, where (3 = 1.59, 1.51,
and 1.43.
78
-------
The BTCs for the CHAIN 2D and HYDRUS Codes are nearly the same for all base parameter values and
upper and lower bounds of the parameter ranges, except for 9 = 0.10 in Figure 7-3; but even for this
exception, the peak concentration for CHAIN 2D is about 95% of that for the HYDRUS Code;
• For the first four parameters (Kd,q,9,p), the CHAIN BTCs are fairly close to those for CHAIN 2D and
HYDRUS, except for q = 0.016 cm/d in Figure 7-2; but even for this exception, the peak concentration for
CHAIN is about 87% of those for the CHAIN 2D and HYDRUS Codes;
• The BTCs for the MULTIMED-DP 1.0 and FECTUZ Codes are almost identical for all values of all
parameters considered.
The peak concentrations for the MULTIMED-DP 1.0 and FECTUZ BTCs range from 5% to 35% higher
than those for CHAIN 2D and HYDRUS, depending on the specific input parameter and its particular
value;
• For the base values of (Kd,q,9,p,DL), the concentration peaks of the MULTIMED-DP 1.0 and FECTUZ
BTCs are about 14% higher than those for CHAIN 2D and HYDRUS;
For the base values of the van Genuchten parameters (Ks,9s,9r,a,p), the concentration peaks of the
MULTIMED-DP1.0 and FECTUZ BTCs are about 21% higher than those of CHAIN 2D and HYDRUS.
The sensitivities and relative sensitivities of the output parameters (Cpeak, Tpeak, TMCL) to the twelve input
parameters (Kd,q,9,p,D,DL,Dw,Ks,9s,9r,a,p), corresponding to the BTCs given in Figures 7-1 to 7-12, are graphically
presented in Figures 7-13 to 7-24. The breakdown of these figures is as follows:
Figures 7-13 to 7-16 compare the output parameters and their sensitivities and relative sensitivities for each
of the five models for the first four input parameters (Kd,q,9,p), respectively;
Figure 7-17 gives the output parameters and their sensitivities and relative sensitivities for the CHAIN Code
for the input parameter, D; this figure is given for the purposes of completion;
Figure 7-18 compares the output parameters and their sensitivities and relative sensitivities for each of the
four models (CHAIN being excluded) for the dispersivity, DL;
Figure 7-19 compares the output parameters and their sensitivities and relative sensitivities for CHAIN 2D
and HYDRUS for the diffusion coefficient in water, Dw;
Figures 7-20 to 7-24 compare the output parameters and their sensitivities and relative sensitivities for each
of the four models (CHAIN being excluded) for the five van Genuchten parameters (Ks,9s,9r,a,p),
respectively.
As shown in Section 4 for the simple model, y = F(a,b,c,d,e,f), the relative sensitivity is a good measure for
comparing output parameter sensitivity to changing input parameters. For this reason, Table 7-3 was constructed
from the information contained in Figures 7-13 to 7-24. This table is a summary of the relative sensitivities for each
of the three output parameters for the "Tc BTCs with respect to all the pertinent input parameters for each of the
five models being tested. The relative sensitivities listed in the table are referenced to the base values of the input
parameters given in Table 6-3. The results evident from Table 7-3 are as follows:
The relative sensitivities are nearly the same for MULTIMED-DP 1.0 and FECTUZ for all output and input
parameters;
The relative sensitivities are nearly the same for CHAIN 2D and HYDRUS for all output and input
parameters;
The relative sensitivities are nearly the same for all models for the output parameters, Tpeak and TMCL, and
all input parameters;
79
-------
oo
o
-0.00
-0.02
£• -i
I
-0.10
0.020
0.015
0.010
0.005
0.04-
0.06- —
0.08 •
•no
HYDROS --
CHAIN 20
FECTUZ
MULTIMED-DP
CHAIN
0.20
£ 0.15
J °'
-5 0.05
0
50000
40000
* 30000
i 20000
I
1 10000
(b)
o 8000
S
d 6000
" 4000
°- 2000
5
0.005 0.010 0.015
Distribution Coefficient (ml/g)
0.020E
0.20
£ 0.15
| 0.10
OT
•5 0.05
0
50000
40000
I" 30000
| 20000
*» 10000
g- 8000
gJ 6000
2
4000
2000
(c)
0.005 0.010 0.015 0.0201
Distribution Coefficient (ml/g) p
0 0.005 0.010 0.015 0.020
Distribution Coefficient (ml/g)
Figure 7-13. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the distribution coefficient using the CHAIN, MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and
HYDRUS Models.
-------
(a)
(b)
(o)
u.ou
£. 0.40
| 0.30
> 0.20
1 0.10
o
09A
0.15
£•
£ 0.10
V) 0.05
0
g* 0.015
*-*
g 0.010
o
o
-at 0.005
i
°" 0
(
1 1 1 1
.
.
:7^^,
1 1 - J 1
••To
"x
X.
X"x
*v
^\,
1 — 1 L L _L-1_LJ '
1 1 1 1
— HYDRUS - - MULTIMED-DP
CHAIN 2D — CHAIN
-- FECTUZ
.
„, , j,»i — — •
•
i i i i
> 0.01 0.02 0.03 0.04 0.
U
£ -0.5
£
i -1-0
I'1-5
-2.0
0,
-100000
^ -200000
~
« -300000
* -400000
^
r^ 4 nnnn
. 8000
o
c
3 6000
D 40OO
o 2000
• 0
05 E C
P
i i i i
N /
S%— ^
i i i i
-
/
/
/
S
'
I I I I
I I I I
\
-
•
i i i i
) 0.01 0.02 0.03 0.04 0.
u
£* -0.50
s
I-1-00
•= -1.50
K
-2.00
0,
-100000
£ -200000
~
1 -300000
OT -400000
^ 8000
O
2 6000
D 4OOQ
o 2000
• 0
05 E (
P
i i i i
/
x^
1 1 1 1
/'*
-
-
1 1 1 1
1 1 1 1
-
Nk
^Kl^^
^^^^^
-
i i i i
> 0.01 0.02 0.03 0.04 0.
Figure 7-14. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the recharge rate using the CHAIN, MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS
Models.
-------
(a)
(b)
(o)
oo
K)
-0.25
£
£ -0.50
g -0.75
w -1.00
£ -1.25
-1.50
_ft ftft,
-0.02
£ -0.04
£
g -0.06
W -0.08
-0.10
g» 0.015
g 0.010
O
-g 0.005
(
NT--
- ^\ -
\
. . ^rt^r. ....^rr^*
X*"
X'
1 1
— HYDRUS - - MULTIMED-DP
CHAIN 2D — CHAIN
-- FECIUZ
*'** 0.10 0.20 0.
Wotor Content i o w o nT^)
_g- 0.80
| 0.60
c
M 0.40
• 0.20
0
CftftAft
40000
^ 30000
| 20000
OT 10000
d 6000
§
" 4000
°- 2000
• 0
30 E (
^^tifttit.—
"^'"
'
•
.
.
-fr^^— « — .
•
• i
.**
^s^
•
> 0.10 0.20 0,
Wotor Contont ( o nt o nf^)
_g- 0.80
£ 0.60
c
J> 0.40
• 0.20
o
CftftftA
40000
^ 30000
| 20000
OT 10000
^ 0
J 6000
2
^ 4000
S
06 2000
• 0
30 E (
^"^^^sur1 • •
"'""'
•
•
.
.
-.sssss=s^-.- •
•
1 1
•
^>^^^
> 0.10 0.20 0.
30
Figure 7-15. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the water content using the CHAIN, MULTIMED-DP 1.0, FECTUZ, CHAIN 2D and HYDRUS
Models.
-------
(a)
(b)
(o)
oo
L.J
U.1U
f •
£ -0.10
c
<7> -0.20
1 -0.30
-0.40
0.0010
0.0006
£ 0.0002
| -0.0002
w -0.0006
-0.0010
^ 0.020
g> 0.015
g 0.010
S
Jt 0.005
D
°" 0
1
-
-
1 1 1
— HYDRUS - - MULTIMED-OP
CHAIN 20 — CHAIN
-- FECTUZ
-
-
60 1.65 1.70 1.75 1.
Bulk Density (g/onf)
U.2U
^ 0.15
c 0.10
•S 0.05
-------
(a)
(b)
(o)
u.u
t-0.2
-0.4
c
w -0.6
1 -0.8
-1.0
0,
-0.005
£•
| -0.010
l -0.015
-0.020
£ 0.015
g 0.010
| 0.005
Q.
0
(
. . / .
: ^^^ :
- ^-——~ .
CHAIN ^c
^^^
. ^— — „ .
> 0.5 1.0 1.5 2.0 2
nfanAr«!nn flnmlltelmnt {eim*/A\
— U.UUU
£. -0.020
£ -0.040
c
> -0.060
• -0.080
-0.100
0,
-50
g -100
OT -150
^ -200
§.
n annn
d 6000
0
" 4000
i
°- 2000
5
« 0
•5 E (
" ^^*^^^^«^ "
-
-
CHAIN »»Tc
-
-
) 0.5 1.0 1.5 2.0 2
nfanAntlnn Hnaff Infant ( nm* /A\
U
^ -0.10
| -0.20
V)
•= -0.30
at
-0.40
0,
-200
^ -400
» -600
* -800
*->-1000
J 6000
2
£ 4000
S
06 2000
5
• 0
* 1 (
-^__ _
^^^
•
CHAIN ^c
-
) 0.5 1.0 1.5 2.0 2
Hfaruinilnn Coafflnlant ^om*/H^
Figure 7-17. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the dispersion coefficient using the CHAIN Model.
-------
(a)
(b)
oo
(Jl
u
£ -0.10
c -0.20
•= -0.30
at
-0.40
-0.0002
^ -0.0004
g -0.0006
OT -0.0008
-0.0010
5;
g 0.015
•^
g 0.010
o
°" 0
3
^^^^_
.
: ^**^- :
--'""
— HYDRUS -- MULTIMED-DP
CHAIN 20 - - FECTUZ
•
.
.0 3.5 4.0 4.5 5.0 5.5 6
— u.uu
£. -0.02
i ~°-04
OT -0.06
• -0.08
-0.10
-20
t-40
c -60
--80
-100
g.
D 800O
o 6000
o
jj 4000
D
°- 2000
s.
• 0
. b ,
n'W^i
^iiiZ^^^
•
•
•^o
.
•
.
.0 3.5 4.0 4.5 5.0 5.5 6
ni*nAr*luIfu fftm\
(o)
-0.10-
-5 -0.30-
at
—0.4ol_
0 -
-100 •
^ -200 •
1 -300 '
- -400 •
-500 _
I
8000
J 6000
I 4000
06 2000
"To
3.0 3.5 4.0 4.5 5.0 5.5 6.0
Dlsperslvlty (om)
Figure 7-18. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the dispersivity using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS Models.
-------
(a)
(b)
(o)
oo
ON
u
£ -0.10
| -0.20
>
•= -0.30
Of
-0.40
0,
-0.0002
£* -0.0004
£
| -0.0006
* -0.0008
-0.0010
^ 0.020
g* 0.015
g 0.010
-g 0.005
__^^^
.
— — — — ' "
•
1 1 1 1
— HYDRUS
CHAIN 20
-
•
u
£ -0.01
i -°-°2
V)
•= -0.03
ae
-0.04
0,
-10
£
3~ -20
W -30
'$ -40
n annn
S
a 6000
o
" 4000
D
°- 2000
.8
• 0
"TTT"TT""~--^
~—
-
1 1 1 1
"To
-
-
— u.uu
>. —0 02
1 -0.04
c
M -0.06
£ -0.08
-0.10
0,
-50
t
£ -100
W -150
^ -200
S
d 6000
2
-g 4000
D
06 2000
5
• 0
•
1 1 1 1
"To
-
•
0.5 1.0 1.5 2.0 2.5 3.0
Diffusion Coefficient In Water (on? d"1)
0.5 1.0 1.5 2.0 2.5 3.0
Diffusion Coefficient In Water (on? eT1)
0.5 1.0 1.5 2.0 2.5 3.0
Diffusion Coefficient In Water (oirf a"4)
Figure 7-19. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the diffusion coefficient in water using the CHAIN 2D and HYDRUS Models.
-------
(a)
(b)
(o)
. Senslttvlfy
» P F
; o 5 £
•5 —u.iu
Of
-0.20
t\ nnm
0.0005
1 0
"«
M -0.0005
-0.0010
n noo
g* 0.015
g 0.010
o
o
Jt 0.005
C
; _ ;
•
•
— HYDRUS -- FECTUZ
CHAIN 20 -- MULTIMED-DP -
'
) 100 200 300 400 5C
Saturated Conductivity (cm d~*)
U.IU
1* 0.05
i °
•5 —0.05
ae
-0.10
2,
1
1 0
n nnnn
a 6000
g
" 4000
JC
°- 2000
• 0
10 E (
-
-
• ^^ •
^ ^ ^
"
) 100 200 300 400 5C
Saturated Conductivity (cm d"*)
U.IU
1* 0.05
i °
•; —0.05
-0.10
2,
1
1 0
^> -2
-J 6000
£ 4000
o
06 2000
• 0
K> E C
•
•
. *~~- .
'
) 100 200 300 400 5G
Saturated Conductivity (cm d"1)
Figure 7-20. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the saturated conductivity using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS
Models.
-------
(a)
(b)
(o)
oo
oo
u
1* -0.200
c -0.400
W
1 -0-600
-0.800
0,
-0.005
I -0.010
W -0.015
-0.020
n n9n
g* 0.015
g 0.010
o
o
-g 0.005
•
— ^
*To
. 22=====-"-^ .
— HYDRUS -- FECTUZ
CHAIN 20 -•-• MULTIMED-DP •
•
•
l.UU
£. 0.80
~ 0.60
C
V> 0.40
^ 0.20
0
9fWW>
15000
1 10000
OT 5000
n nnnn
o 6000
0
" 4000
°- 2000
&
• 0
•
•
•
.
•
l.UU
£. 0.80
~ 0.60
C
W 0.40
£ 0.20
0
onnAn
15000
1 10000
W 5000
>• nnnn
J 6000
2
-c 4000
06 2000
• 0
•
.
•
•
.
•
0.280 0.300 0.320 0.340 0.360
Saturated Water Content (crn'onf*)
0.280 0.300 0.320 0.340
Saturated Water Content (cnr
0.360
0.280 0.300 0.320 0.340 0.360
Saturated Water Content (cm'cnf*)
Figure 7-21. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the saturated water content using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS
Models.
-------
(o)
(b)
(o)
u
^ -0.100
"5
| -0.200
V>
•= -0.300
Of
-0.400
-0.010
s| -0.020
OT -0.050
-0.040
A O9A
g* 0.015
g 0.010
o
.* 0.005
i
"- o
~_^
' ^^Sfc.
••To
,*••'"""
^I££;^*^
'""- -
— HYDRUS -- FECTUZ
CHAIN 20 -- MULTIMED-DP -
-
-
U.4U
fo.30
g °-20
Ul
•x 0.10
o:
0
40000
30000
« 20000
OT 10000
V °
n mvin
•g 800°
d 6000
" 4000
°- 2000
• 0
.
....^^ffrff^r..,
•
•
•
"•'-'-•-.-,-
-
-
U.4U
fO.30
c 0.20
-5 0.10
at
0
15000
g 10000
i/i 5000
^-> o
>» QAnn
-J 6000
£ 4000
06 2000
• 0
_---
^,~£1£Z''^ff^"~"
-
'*i"'*'»
1 ~ "*•''-
' "..
-
-
0.05 0.06 0.07 0.08 0.09 0.10 0.11
Residual Water Content (cm3 citf)
0.050 0.060 0.070 0.080 0.090 0.100 0.110
Residual Water Content (cm3 citf5)
0.050 0.060 0.070 0.080 0.090 0.100 0.110
Residual Water Content (om'onT3)
Figure 7-22. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the residual water content using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and HYDRUS
Models.
-------
(a)
(o)
u.zuu
£ 0.150
|
c 0.100
-5 0.050
0£
0
0.020
0.015
I 0.010
°5
c
OT 0.005
0
^ 0.020
g 0.015
g 0.010
^:\
-
•
050 0.052 0.054 0.056 0.058 O.C
Genuohten Retention Parameter a ( a
u
^ -0.05
i -°-10
•s -0.15
et
-0.20
0
-2000
£ -4000
| -6000
OT -8000
^10000
|8000
J 6000
2
£ 4000
i
* 2000
5
« 0
160 E 0.
nf1) p van
• _—^-^ •
• ^^^-rrrr. '
o^""**
•
•
•
050 0.052 0.054 0.056 0.058 0.060
Genuohten Retention Parameter a (citT1)
Figure 7-23. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the van Genuchten retention parameter a using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and
HYDRUS Models.
-------
(o)
(b)
(o)
z.uu
1" 1.50
c 1.00
OT
•= 0.50
at
0
0.008
^ n nne
_>
g 0.004
* 0.002
0
fi n9n
^> o.ozu
g> 0.015
N^
g 0.010
o
jc 0.005
8
"• o
; ^^^;
•
"To
— HYDRUS -- FECTUZ
CHAIN 20 -- MULTIMED-DP -
•
u
^ -0.50
" -1.00
VI
•x -1.50
at
-2.00
-2000
^ —Af\nn
g -6000
5 -8000
*£r 10000
n ftnnn
S
d 6000
" 4000
g
°- 2000
j>
• 0
.,,:•/
• ______--—-
•
"To
•
u
^ -0^0
c -1.00
(ft
•s -1.50
at
-2.00
-2000
^ _4nnn
^
« -6000
5 -8000
^r 10000
>* fiAAA
Z
-J 6000
2
-c 4000
o
g
« 2000
_o
• 0
' „«=»———
•
"To
-_- . _ . _
•
1.400 1.450 1.500 1.550 1.600 E
van Genuohten Retention Parameter (I (—) *~
1.400 1.450 1.500 1.550 1.600 E
van Genuohten Retention Parameter /8 (—) *~
1.400 1.450 1.500 1.550 1.600
van Genuohten Retention Parameter /3 (—)
Figure 7-24. Sensitivity and relative sensitivity of (a) peak concentration at the depth of 6m, (b) time to reach peak concentration at the depth of 6m, and (c)
time to exceed MCL at the receptor well to the van Genuchten retention parameter (3 using the MULTIMED-DP 1.0, FECTUZ, CHAIN 2D, and
HYDRUS Models.
-------
Table 7-3. Summary of Relative Sensitivities for the Outputs Obtained from the "Tc Breakthrough Curves with Respect to All the Pertinent Input
Parameters for All Models, Referenced to the Base Values of the Input Parameters.
Model Input
Parameter
peak
Q
So
rt
5
FECTUZ
o
HYDRUS
peak
Q
So
rt
5
FECTUZ
U
HYDRUS
5
FECTUZ
<
ffi
0
a
8
Kd
q
e
p
D
DL
Dw
es
er
a
P
-0.06 -0.05 -0.05 -0.06 -0.06
+0.40 +0.00 +0.00 +0.12 +0.12
-1.10 -0.62 -0.62 -0.90 -1.10
-0.17 -0.05 -0.05 -0.05 -0.05
-0.45 ~~ ~~ ~~ ~~
~~ -0.36 -0.36 -0.29 -0.28
~~ ~~ ~~ -0.10 -0.10
~~ +0.04 +0.04 +0.05 +0.05
~~ -0.51 -0.51 -0.66 -0.66
~~ -0.20 -0.20 -0.29 -0.29
~~ +0.05 +0.06 +0.13 +0.13
~~ +0.86 +0.86 +1.30 +1.38
+0.06 +0.06 +0.06 +0.06 +0.06
-1.10 -1.10 -1.10 -1.10 -1.10
+0.80 +0.80 +0.81 +0.78 +0.80
+0.06 +0.06 +0.07 +0.06 +0.06
-0.02 ~~ ~~ ~~ ~~
~~ -0.02 -0.02 -0.02 -0.02
~~ ~~ ~~ -0.01 -0.01
~~ -0.05 -0.04 -0.04 -0.04
~~ +0.60 +0.60 +0.57 +0.57
~~ +0.22 +0.22 +0.22 +0.22
~~ -0.05 -0.08 -0.05 -0.05
~~ -1.00 -1.00 -1.00 -1.00
+0.07 +0.07 +0.07 +0.07 +0.07
-0.98 -1.00 -1.00 -0.98 -0.98
+0.85 +0.92 +0.92 +0.84 +0.86
+0.08 +0.07 +0.07 +0.07 +0.07
-0.10 ~~ ~~ ~~ ~~
~~ -0.09 -0.09 -0.07 -0.07
~~ ~~ ~~ -0.03 -0.03
~~ -0.05 -0.04 -0.04 -0.04
~~ +0.66 +0.66 +0.64 +0.64
~~ +0.26 +0.26 +0.24 +0.24
~~ -0.06 -0.09 -0.05 -0.05
~~ -1.13 -1.15 -1.03 -1.03
-------
For the output C k, there are some differences in the relative sensitivities between the CHAIN Model, the
MULTIMED-DP 1.0 and FECTUZ Models, and the CHAIN 2D and HYDRUS Models;
All output parameters for all models are rather insensitive to Kd and p, as one would expect for the "Tc
BTCs as discussed in Section 6-2 (i.e., 9 » p Kd);
Conversely, all output parameters for all models are rather sensitive to 9;
The output parameters, Tpeak and TMCL, for all models are also rather sensitive to q, while output Cpeak is
much less sensitive to q, the sensitivity being model-dependent;
For the van Genuchten parameters, the model output parameters are most sensitive to p\ followed by 9S,
while the output parameters are least sensitive to Ks and a.
7.3 Other Results Illuminating Numerical Differences/Errors Between the Models
For the sensitivity modeling scenario employed in this report, the major differences between the models arose
due to an error in coding, the use of different numerical inversion schemes (inversion from Laplace space to physical
space), and the treatment of the dispersion factor, 9D, in Equation (6-4).
7.3.7 Correction of the MULTIMED-DP 1.0 Code
When the van Genuchten Model for the water retention parameters (Ks,9s,9r,a,p) was considered in the
sensitivity analysis depicted in Table 7-1, an error in the originally distributed MULTIMED-DP 1.0 Code was
detected. This error arose from the solution of Darcy's Law for the pressure head, which gives the water content in
the soil column through the use of the van Genuchten Model. Figure 7-25 shows the water content profiles in the 6m
soil column derived from Darcy's Law for the MULTIMED-DP 1.0 and FECTUZ Codes, and from the steady-state
solution of Richards Equation for the CHAIN 2D and HYDRUS Codes. The water contents obtained through the
use of the originally distributed MULTIMED-DP 1.0 Code are inconsistent with respect to those obtained from the
other three codes. The reason for the inconsistency was found to be the incorrect use of the residual water content,
9r, in the van Genuchten Module. When this error was corrected, the new water content results became consistent
with those of the other three models, as shown in Figure 7-25. Figure 7-25 also shows that the steady-state solution
of the Richards Equation gives the same results as the solution of Darcy's Law, as we know it should from our
discussion in Section 6-3.
u
1 0
|2.0
£ 3.0
Q.
- ,. ,
0.20 0.30 0.40 0.
Water Content (cnr?/cir?)
Figure 7-25. Water content distributions predicted by the HYDRUS, CHAIN 2D, FECTUZ, and MULTIMED-
DP 1.0 models. Note that the water contents (•) obtained from the originally distributed
MULTIMED-DP 1.0 code are in error. The corrected code gives a consistent water content
distribution (A.) with the other three models.
93
-------
7.3.2 Comparison of the Stehfest andDeHoog Inversion Algorithms
All the major modules of the MULTIMED-DP 1.0 and FECTUZ Codes are the same (or nearly the same),
except the MULTIMED-DP 1.0 Code has four alternate inverse transform modules (analytical, convolution, Stehfest,
DeHoog), while the FECTUZ only uses the DeHoog Module. For the current uses of these codes, the analytical and
convolution modules are not applicable. Thus, one may ask: Are there any differences between the Stehfest and
DeHoog Algorithms? We found that the answer is YES to this question for the generation of "Tc BTCs using the
MULTIMED-DP 1.0 Code.
The Stehfest Algorithm (Stehfest, 1970) produced non-realistic oscillations in the "Tc BTCs for most of the
parameters studied in this report. The DeHoog Algorithm (DeHoog et al., 1982) is an improved inversion scheme
which was developed to eliminate errors in the use of the analytical and convolution algorithms and to reduce the
oscillatory conditions often produced by the Stehfest Algorithm. We compared the "Tc BTCs produced by the
Stehfest and DeHoog Algorithms in the MULTIMED-DP 1.0 Code for several different input parameters. The
following results were evident from this study:
For all input parameters tested, the Stehfest Algorithm produced oscillatory behavior and peak
concentrations lower than those produced by the DeHoog Algorithm;
For all input parameters tested, the DeHoog Algorithm produced "Tc BTCs almost identical to those
produced by the FECTUZ Code;
As an example of this comparison, Figure 7-26 shows the "Tc BTCs for 9 = 0.08, 0.16, 0.22, where (a)
illustrates the oscillatory behavior of the Stehfest Algorithm in MULTIMED-DP 1.0, (b) illustrates the
corrective behavior of the DeHoog Algorithm in MULTIMED-DP 1.0, and (c) indicates the near perfect
comparison between the FECTUZ and MULTIMED-DP 1.0 Codes, when the DeHoog Algorithm is used.
7.3.3. Increasing the Base Value of Dispersivity in the FECTUZ Code
As we have previously discussed, CHAIN only specifies the value for D in the dispersion factor, 9D, in Equation
(6-4), while MULTIMED-DP 1.0 and FECTUZ only specify the value for DL. CHAIN 2D and HYDRUS specify
both DL and Dw. Since the recharge rate, q, is so small in our sensitivity modeling scenario, there is a significant
difference between the base value of 9D for MULTIMED-DP 1.0 and FECTUZ, and that for CHAIN 2D and
HYDRUS (see the fourth row of Table 7-2). This effect was evident in the "Tc ETC curves shown in Figures 7-1 to
7-12, where the peak concentrations from MULTIMED-DP 1.0 and FECTUZ were higher than those of CHAIN 2D
and HYDRUS. To see if this difference was primarily due to 9D, we increased the base value of DL in FECTUZ to
6.53 cm. This gives a value of 9D for FECTUZ of 0.157 cm2/d, versus that of 0.160 cm2/d for CHAIN and 0.162
cm2/d for CHAIN 2D and HYDRUS. Figure 7-27 compares the "Tc BTCs for the base parameters for the
HYDRUS, CHAIN 2D and CHAIN Models versus the "Tc ETC for the FECTUZ Model with the base value of DL
= 4.53 cm replaced by the value DL = 6.53 cm. This comparison proves that the major differences between the
model BTCs in Figures 7-1 to 7-12 are due to the formulation of the dispersion factor, 9D, and not to other factors.
94
-------
(a)
(b)
(c)
U.UI3
CT
^^
o>
E
^ 0.010
c
JO
"5
L.
c
0>
o
§ 0.005
O
u
i-
m
Q
I 1 1 I 1
MULTIMED-DP 1.0 Water Content (crf/crf)
_, ,_ , A AQ
Stehfest u>uo
n 1 R
0.22
_
s~\
/ \
' / \
: y V'"^
• / A '\
: / :- 1 \ \
•: \ ~'~ 1 \ "^
/ / '•/ \ \
iL—J i A i v i N'^, i
0 2000N-/>OefO 6000 "§000 TOOOO 12000
Time (days)
Oni c
.UlO
cr
o>
E
^ 0.010
c
,0
|
"E
v-' 0.010
c
£
|
n^
0,22
-
; \ A
: / \ -^\
: / \ /
': / \/ \
7 \ ^-
^ / A \
-^ / '• '• \ \
•: / '• •'" V ^
i.-' _/ i X i > N i
0 2000 4000 6000 8000 10000 12000
Time (days)
Figure 7-26. Sensitivity of "Tc BTCs at 6m depth to water content, where (a) uses Stehfest Algorithm of
MULTIMED-DP 1.0, (b) uses DeHoog Algorithm of MULTIMED-DP 1.0, and (c) uses DeHooj
Algorithm of FECTUZ Code.
95
-------
o
0>
o
c
o
o
0.015
0.010
0.005
Uniform Water Content Distribution (Base Case)
HYDRUS
CHAIN 2D
FECTUZ (DL=6.53 cm)
CHAIN
2500 5000 7500
Time (days)
10000
Figure 7-27. Comparison of the "Tc BTCs for the HYDRUS, CHAIN 2D and CHAIN Models for the base values
of the input parameters with the ETC for the FECTUZ Model with the base value of DL = 4.53 cm
replaced by the value DL = 6.53 cm.
96
-------
Section 8
Summary and Conclusions
The purpose of this report was to evaluate and perform sensitivity analyses on select computer models for
simulating radionuclide fate and transport in the unsaturated zone. The work reported here was performed in support
of the U.S. Environmental Protection Agency's (EPA's) project on Soil Screening Guidance (SSG) for Radionuclides
(U.S. EPA, 2000ab). The U.S. EPA developed the SSG as a tool to help standardize and accelerate the evaluation
and cleanup of soil contaminated with radioactive materials at sites on the National Priority List (NPL) with future
residential land use. This guidance is intended for the appropriate environmental professionals to be able to calculate
risk-based, site-specific, soil screening levels (SSLs) for radionuclides in soil, thus allowing them to identify areas
needing further investigation and possible remediation at NPL sites. The models reviewed and evaluated in the
current report form a small subset of the models available to the public, and those which are considered here have not
met U.S. EPA approval for exclusive use in the SSL process. Other models may be applicable to the SSL effort,
depending on pollutant- and site-specific circumstances.
The model selection process for the radionuclide SSL effort is briefly described in Section 2.1, and resulted in
the following five models being chosen for evaluation and sensitivity analysis:
One-dimensional CHAIN Code (van Genuchten, 1985),
• One-dimensional MULTIMED-DP 1.0 Code (Liu et al., 1995; Sharp-Hansen et al., 1995; Salhotra et al.,
1995)
• One-dimensional FECTUZ Code (U. S. EPA, 1995ab),
• One-dimensional HYDRUS Code (Simunek et al., 1998),
Two-dimensional CHAIN 2D Code (Simunek and van Genuchten, 1994).
The CHAIN Code is an analytical model which requires the assumption of uniform flow conditions. The
MULTIMED-DP 1.0 and FECTUZ Codes are semi-analytic models which approximate complex analytical solutions
using numerical methods such as numerical inverse transform modules. Transient and steady-state conditions can be
accommodated in both the analytical and semi-analytical models, while layered media are usually only
accommodated by the semi-analytical codes. The HYDRUS and CHAIN 2D Codes solve the pertinent sets of partial
differential equations using finite-difference or finite-element methods. The resolution in space and time of these
numerical models depends on the physical characteristics of the site in question, the computational resources
available to the modeler, and the purposes of the simulation. Numerical codes are used when simulating time-
dependent scenarios under spatially varying soil conditions and unsteady flow fields. Table 2-1 lists, in a
comparative manner, the flow and transport model components of each of the five models. Further explanation of
the model components is given in Appendices A to D.
Computer model predictions are uncertain and erroneous because of uncertainties and errors that can occur at
various points in the modeling process and in the model's structure. Errors can occur in conceptualization or model
structure when processes and the assumptions represented by the model fail to represent the intended reality of the
physical problem at hand. Errors in experimental measurements of model input parameters and model output
parameters used for calibration and testing can also lead to uncertainties. For numerical computer codes, the
discretization in space and time introduces error; the computer, itself, introduces roundoff error; and truncation errors
are produced whenever analytical expressions are replaced by truncated series of elementary functions. Other errors
and uncertainties result from the imprecise specification of boundary and initial conditions, poor model
parameterization, and incorrect spatial scale transitions (see Section 1.3). These uncertainties and errors lead to a
three-tier sensitivity analysis:
97
-------
Sensitivity of simulated results to conceptual model selection (Section 3),
Model output parameter sensitivity to varying model input parameters (Sections 4 to 6),
Comparison of sensitivity results due to numerical differences and errors (Section 7).
The breadth and depth of the variability and heterogeneity of the phenomena that can occur in the fate and
transport of radionuclides in the unsaturated zone are briefly considered in Section 3 of this report. Phenomena
considered include the following: the selection of a physical problem's time/space domain, the choice of boundary
and initial conditions for the flow and transport variables, the potential impacts of hy steretic effects and density and
thermal gradients, the importance of preferential pathways and facilitated transport, the form of the transport
dispersion processes, the impacts of non-ideal sorption and reaction processes, the importance of decay processes
and the evolution of decay products, the influence of multicomponent and multiphase systems, and the effects of
highly heterogeneous media. Against this backdrop, the five selected models were compared (see Table 2-1). All of
the models are deterministic, and thus, are incapable of stochastic simulation, except for Monte Carlo simulations for
certain input parameter or boundary condition distributions. Four of the five models are one-dimensional and only
CHAIN 2D is two-dimensional. Thus, truly three dimensional features cannot be easily, or even correctly, simulated,
except for idealistic features such as those presented in Figure 3 -1. Density gradient problems cannot be addressed
by these models, nor can facilitated transport and preferential pathways. The scale dependency of hydraulic
conductivity and dispersivity due to soil heterogeneities is not considered by these models. Other processes and
features that are not addressed by these models are nonlinear reactions, nonlinear nonequilibrium sorption,
nonequilibrium site heterogeneity, and variable decay processes. To address precipitation/evapotranspiration/
infiltration processes on a storm event basis for the several-year periods required for some radionuclides will require
extensive computational resources for the models in this group which are capable of performing such simulations.
Appendices E to I present the sensitivity results of the various models to the following conceptualizations, given
respectively as: nonuniform moisture distribution versus uniform moisture distribution, daily variation of
precipitation/evaporation/infiltration versus annual recharge rate, layered soil column versus homogeneous soil
column, nonequilibrium sorption of pollutants versus equilibrium sorption of pollutants, and variations in fate and
transport of radionuclides with different distribution coefficients and half-lives.
The basic elements of parameter sensitivity analysis are defined and illustrated in Section 4. To carry out this
analysis for radionuclides in the unsaturated zone realistically requires the characteristics and properties from a
representative physical site. The site selection process and the choice of the candidate site are given in Sections 5.1
and 5.2, respectively. The candidate site, located in the semi-arid West, is the Las Graces Trench Research Site in
New Mexico. Some of the soil hydraulic properties at this site and other site characteristics are given in Tables 5-2
and 5-3 and Figure 5-2. The hypothetical modeling scenario for this site assumed a step-wise constant source of
radionuclides at the top of a 6m, vertical, homogeneous soil column with a zero concentration gradient at the bottom
of the column. The input parameters (the base values given in Table 5-4) and the output parameters used in the
sensitivity analyses are as follows:
Input Parameters
Kd = Distribution Coefficient Dw = Diffusion Coefficient in Water
q = Recharge Rate Ks = Saturated Hydraulic Conductivity
9 = Soil Water Content 9S = Saturated Water Content
p = Bulk Density 9r = Residual Water Content
D = Dispersion Coefficient a = van Genuchten Retention Parameter
DL = Dispersivity P = van Genuchten Retention Parameter
Output Parameters
Cpeak = P63^ Concentration of Radionuclide ETC at 6m.
Tpeak = Time to Reach Peak Concentration at 6m.
TMCL = Time When Radionuclide Concentration at 6m is High Enough to Exceed MCL at Receptor Well
Section 6 presents the parameter sensitivity analysis results for the "Tc breakthrough curves (BTCs) at the
bottom of the 6m, vertical, homogeneous soil column, where the BTCs were generated by the HYDRUS Code. The
BTCs for the eleven input parameters (Kd,q,9,p,DL,Dw,Ks,9s,9r,a,P) are given in Figure 6-1, where the BTCs for the
first six input parameters (Kd,q,9,p,DL,Dw) were obtained under the constraints of constant recharge rate, q, and
98
-------
constant moisture content, 9, while those for the last five parameters (Ks,9s,9r,a,p) were obtained under the
constraints of constant recharge rate, q, and variable soil water content, 9(z), see Sections 6.2 and 6.3, respectively.
The sensitivity and relative sensitivity results for these BTCs are given in the following figures:
Figures 6-2a to 6-2k; for Cpeak;
Figures 6-3a to 6-3k; for Tpeak;
Figures 6-4a to 6-4k; for TMCL.
Table 6-4 summarizes the relative sensitivities of (Cpeak, Tpeak, TMCL) to the eleven input parameters for the
HYDRUS Code, the sensitivities being measured at the base values of the input parameters. The dominant input
parameters for all three output parameters are those which have the greatest influence on the advection term in the
transport equations, followed by those which have the greatest influence on the dispersion term. The parameters, p
and Kd, have a very small influence on the transport and fate of "Tc since 9 » pKd in the coefficients of the
transport equations. For the first six input parameters (Kd,q,9,p,DL,Dw), Cpeak is most influenced by 9, while Tpeak
and TMCL are most influenced by both 9 and q. For the last five input parameters (from the van Genuchten Model,
Ks,9s,9r,a,p), the most influential input parameter for all three output parameters is p, followed by 9S, then followed
by9r, a,andKs.
Section 7 compares the sensitivity results between the five models, illuminating numerical differences and
errors. The numerical differences and errors separating the five models are summarized in Table 8-1. Each set of
sensitivity results are derived for the 6m, vertical, homogeneous soil column; and the first seven input parameters
(Kd,q,9,p,D,DL,Dw) are compared under the constraints of constant recharge rate, q, and constant water content, 9.
The last five input parameters (Ks,9s,9r,a,p) are compared under the constraints of constant recharge rate, q, and
variable water content, 9(z). Figures 7-1 to 7-12 compare the "Tc BTCs at the bottom of the 6m soil column for the
various models and input parameters, while Figures 7-13 to 7-24 compare the sensitivities and relative sensitivities of
(Cpeak' Tpeak, TMCL) to the twelve input parameters for the various models. The breakdown of these various figures
with reference to input parameters and models is given in Table 8-2. The results derived from these figures and
calculations are as follows:
The BTCs for the CHAIN 2D and HYDRUS Codes are nearly the same for all input parameter values
considered.
• For the input parameters (K d,q,9,p), the CHAIN BTCs are fairly close to those for CHAIN 2D and
HYDRUS.
• The BTCs for the MULTIMED-DP 1.0 and FECTUZ Codes are almost identical for all values of the input
parameters considered.
• For the base values of (K d,q,9,p,DL), the concentration peaks of the MULTIMED-DP 1.0 and FECTUZ
BTCs are about 14% higher than those for CHAIN 2D and HYDRUS, because of the differences in the
dispersion factor, 9D, as given in Table 8-1.
• For the base values of the van Genuchten parameters (K s,9s,9r,a,p), the concentration peaks of the
MULTIMED-DP 1.0 and FECTUZ BTCs are about 21% higher than those of CHAIN 2D and HYDRUS,
for the reason given above.
When the base value of D L is increased to compensate for the absence of 9iwDw in 9D (see Table 8-1),
FECTUZ and MULTIMED-DP 1.0 peak concentrations of the BTCs closely approach those of CHAIN 2D
and HYDRUS (see Figure 7-27).
• The relative sensitivities are nearly the same for MULTIMED-DP 1.0 and FECTUZ for all output
parameters and all input parameters.
The relative sensitivities are nearly the same for CHAIN 2D and HYDRUS for all output parameters and all
input parameters.
The relative sensitivities are nearly the same for outputs (T peak, TMCL) for all models and all input
parameters.
• For the output (C peak), there are some differences in the relative sensitivities for the three groups of models:
CHAIN, MULTIMED-DP 1.0 and FECTUZ, CHAIN 2D and HYDRUS.
All output parameters for all models are rather insensitive to K d and p, as one would expect for the "Tc
BTCs since 9 » pKd in the coefficients of the transport equations. This will not be the case for higher Kd
radionuclides.
Conversely, all output parameters for "Tc BTCs are rather sensitive to 9 for all models.
The output parameters (T k, TMCL) for all models are also rather sensitive to q, while Cpeak is much less
sensitive to q, the sensitivity being model-dependent.
99
-------
Table 8-1. The Possible Numerical Differences and Errors (•) Separating the Five Models Tested for Water Flow
and Radionuclide Transport through the 6m, Vertical, Homogeneous Soil Column.
Numerical Differences
and Errors
Models Tested in Report
MULTIMED-
CHAIN DP 1.0 FECTUZ CHAIN 2D HYDRUS
Computer Roundoff Errors
Analytical Solution Truncation
Errors for Transport Equations
Inverse Transform Scheme Errors
for Laplace Space Analytical
Solutions of Transport Equations
Finite Difference/Finite Element
Solution Errors for Transport
Equations
Constant Recharge Rate, q, and
Constant Water Content, 9
Constant Recharge Rate, q, and
Variable Water Content, 9(z)
Numerical Integration Errors of
Darcy'sLaw
Numerical Integration Errors of
Richards Equation
Numerical Errors in Deriving the
Steady State Solution of
Richards Equation
Specification of
9D = 9lwDw + DL|q|
in the Transport Equations
Possible Numerical Diffusion
in the Transport Equations
Specifies Specifies Specifies Specifies Specifies
only only only both both
D DL DL DLandDw DLandDw
For the van Genuchten parameters (K s,9s,9r,a,p), the model output parameters are most sensitive to p\
followed by 9S, while the output parameters are least sensitive to Ks and a.
Numerical diffusion for the current modeling scenario is rather insignificant (see Figure 7-27) since the D L-
adjusted FECTUZ BTCs (a model with no numerical diffusion) correspond closely to those for CHAIN 2D
and HYDRUS, models which may have numerical diffusion.
The strongest discriminating differences (see Table 8-1) between the five models are those due to the form
of 9D, and those due to the uniform/nonuniform water content. Other numerical differences/errors appear
to be of lower order effect for the current modeling scenario.
Of the five models tested in this report and for the given modeling scenario, there are basically three
groupings: CHAIN, MULTIMED-DP 1.0 and FECTUZ, CHAIN 2D and HYDRUS.
100
-------
Table 8-2. The Figures in Section 7 Comparing "Tc BTCs and Output Sensitivities with Respect to Model Input
Parameters and Modeling Codes.
Model
Input
Parameter
Kd
q
e
p
D
DL
Dw
Ks
es
er
a
3
CHAIN
7-1
7-2
7-3
7-4
7-5
Figures
MULTIMED-
DP 1.0
7-1
7-2
7-3
7-4
7-6
7-8
7-9
7-10
7-11
7-12
Comparing BTCs
FECTUZ
7-1
7-2
7-3
7-4
7-6
7-8
7-9
7-10
7-11
7-12
CHAIN 2D
7-1
7-2
7-3
7-4
7-6
7-7
7-8
7-9
7-10
7-11
7-12
HYDRUS
7-1
7-2
7-3
7-4
7-6
7-7
7-8
7-9
7-10
7-11
7-12
Figures Comparing Sensitivities
CHAIN
7-13
7-14
7-15
7-16
7-17
MULTIMED-
DP 1.0
7-13
7-14
7-15
7-16
7-18
7-20
7-21
7-22
7-23
7-24
FECTUZ
7-13
7-14
7-15
7-16
7-18
7-20
7-21
7-22
7-23
7-24
CHAIN 2D
7-13
7-14
7-15
7-16
7-18
7-19
7-20
7-21
7-22
7-23
7-24
HYDRUS
7-13
7-14
7-15
7-16
7-18
7-19
7-20
7-21
7-22
7-23
7-24
101
-------
Section 9
References
Adler, P.M. and J.-F. Thovert. 1993. "Fractal porous media." Transport in Porous Media. 13,41-78
Andraski, B. J. 1996. "Properties and variability of soil and trench fill at an arid waste-burial site." Soil Sci. Soc.
Am. J. 60, 54-66.
Andraski, BJ. 1997. "Soil-water movement under natural-site and waste-site conditions: A multiple year field study
in the Mojave Desert, Nevada." Water Resour. Res.. 33(8), 1901-1916.
Assouline, S.D., D. Tessier, and A. Bruand. 1998. "A conceptual model of the soil water retention curve." Water
Resour. Res. 34(2), 223-231.
Ayyub, B.M. 1998. Uncertainty Modeling and Analysis in Civil Engineering. CRC Press, Boca Raton, FL. 506pp.
Bai, M., D. Elsworth, H.I. Inyang and J.-C. Roegiers. 1997. "Modeling contaminant migration with linear sorption
in strongly heterogeneous media." Journal of Environ. Eng. 123 (11), 1116-1125.
Bailey, R.G. 1980. Description of the Ecoregions of the United States. United States Department of Agriculture.
Forest Service. Misc. Publ. 1391. Washington, DC. 77pp.
Bear, J. 1993. "Modeling flow and contaminant transport in fractured rocks." Flow and Contaminant Tranport in
Fractured Rock. Academic Press, Inc. New York, NY. 1-37.
Bear, J. and J.J. Nitao. 1995. "On equilibrium and primary variables in transport in porous media." Transport in
Porous Media. 18, 151-184.
Bear, J. and A. Verruijt. 1987. Modeling Groundwater Flow and Pollution. D. Reidel Publishing Company.
Boston, MA. 414pp.
Beckie, R., A.A. Aldama and E.F. Wood. 1994. "The universal structure of the groundwater flow equations."
Water Resour. Res. 30(5), 1407-1419.
Bedient, P.B.,H.S. Rifai and CJ. Newell. 1999. Ground Water Contaminant: Transport and Remediation. Second
Edition. Prentice Hall PTR. Upper Saddle River, NJ. 604pp.
Berglund, S. and V. Cvetkovic. 1996. "Contaminant displacement in aquifers: Coupled effects of flow
heterogeneity and nonlinear sorption." Water Resour. Res. 32(1)23-32.
Boek, E.S. and P. van der School. 1998. "Resolution effects in dissipative particle dynamics simulations." Inter. J.
of Modern Physics C. 9(8), 1307-1318.
Brewer, K.E. and S. W. Wheatcraft. 1994. "Including multi-scale information in the characterization of hydraulic
conductivity distributions." Wavelets in Geophysics. Academic Press, Inc. New York, NY. 213-248.
Brokate, M. and J. Sprekels. 1996. Hysteresis and Phase Transitions. Springer-Verlag, Inc. New York, NY. 357
pp.
102
-------
Brooks, R.H. and A.T. Corey. 1964. Hydraulic Properties of Porous Media. Hydrology Paper No. 3. Colorado
State University. Fort Collins, CO.
Brusseau, M.L. and P.S.C. Rao. 1989. "Sorption nonideality during organic contaminant transport in porous
media." CRC Grit. Rev. Environ. Control. 19, 33-99.
Carerra, J. 1993. " An overview of uncertainties in modelling groundwater solute transport." Journal of Contain.
Hydrology. 13,23-48.
Charbeneau, R.J. 2000. Groundwater Hydraulics and Pollutant Transport. Prentice Hall. Upper Saddle River, NJ.
593 pp.
Chen, W. and RJ. Wagenet. 1995. "Solute transport in porous media with sorption-site heterogeneity." Environ.
Sci. Technol. 29, 2725-2734.
Chiles, J-P. and G. de Marsily. 1993. "Stochastic models of fractured systems and their use in flow and transport
modeling." Flow and Contaminant Transport in Fractured Rock. Academic Press, Inc. 169-236.
Clemo, T. and L. Smith. 1997. "A hierarchical model for solute transport in fractured media." Water Resour.
Res33(8), 1763-1783.
Corapcioglu, M.Y. and H. Choi. 1996. "Modeling colloid transport in unsaturated porous media and validation with
laboratory column data." Water Resour. Res. 32(12), 3437-3449.
Corwin, D.L., K. League, and T.R. Ellsworth. Editors. 1999. Assessment of Non-Point Source Pollution in the
Vadose Zone. Geophysical Monograph 108. American Geophysical Union. Washington, D.C. 369 pp.
Coveney, P.V. andK.E. Novik. 1996. "Computer simulations of domain growth and phase separation of two-
dimensional binary immiscible fluids using dissipative particle dynamics." Physical Review E. 54(5), 5134-
5141.
Crawford, J.W., P. Baveye, P. Grindrod and C. Rappoldt. 1999. "Application of fractals to soil properties,
landscape patterns, and solute transport in porous media." Assessment of Non-Point Source Pollution in the
Vadose Zone. Geophysical Monograph 108. American Geophysical Union, Washington, DC. 151-164.
Dagan, G. 1989. Flow and Transport in Porous Formations. Springer-Verlag. New York, NY. 465pp.
Defence Mapping Agency. 1987. Official Road Map: U.S. Army White Sands Missile Range. New Mexico.
Defence Mapping Agency. Hydrographic/Topographic Center. White Sand Missile Range, MM.
De Hoog, F.R., J.H. Knight and A.N. Stokes. 1982. "An improved method for numerical inversion of Laplace
Transforms". SIAM J. SCI. STAT. COMPUT. 3(3), 357-366.
Domenico, P.A. and F.W. Schwartz. 1990. Physical and Chemical Hydro geology. John Wiley and Sons. New
York, NY. 824pp.
Drew, D.A. and S.L. Passman. 1999. Theory of Multicomponent Fluids. Springer-Verlag. New York, NY. 308
pp.
Eagleson, P.S. 1970. Dynamic Hydrology. McGraw-Hill Book Co. New York, NY. 462pp.
Eeckhout, E.V. 1997. Use of Data Fusion to Optimize Contaminant Transport Predictions. LA- UR-98-845. U.S.
DOE, Los Alamos National Laboratory. Los Alamos, MM. 16 pp.
Faybishenko, B., P.A. Witherspoon and S.M. Benson. Editors. 2000. Dynamaics of Fluids in Fractured Rock..
Geophysical Monograph 122. American Geophysical Union. Washington, DC. 400 pp.
Fetter, C.W. 1999. Contaminant Hdyrogeology. Second Edition. Prentice Hall. Upper Saddle River, NJ. 500pp.
103
-------
Flury, M. and H. Fltihler. 1995. "Modeling solute leaching soils by diffusion-limited aggregation: Basic concepts
and application to conservative solutes." Water Resour. Res. 31(10), 2443-2452.
Gao, Y. and M.M. Sharma. 1994. "A LGA model for fluid flow in heterogeneous porous media." Transport in
Porous Media 17, 1-17.
Gee, G.W.., C.T. Kincaid, RJ. Lenhard and C.S. Simmons. 1991. "Recent studies of flow and transport in the
vadose zone." U.S. National Report to International Union of Geodesy and Geophysics. 1987-1990. American
Geophysical Union, Washington, DC. 227-239.
Gee, G.W., P J. Wierenga, BJ. Andraski, M.H. Young, M. J. Payer, and M.L. Rockhold. 1994. "Variations in water
balance and recharge potential at three western desert sites." Soil Sci. Soc. Am. J. 58(1), 63-72.
Gelhar, L.W. 1993. Stochastic Subsurface Hydrology. Prentice-Hall. Englewood Cliffs. NJ.
Gile, L.H., R.P. Gibbens, and J.M. Lenz. 1998. "Soil-induced variability in root systems of creosotebush (Larrea
tridentata) and tarbush (Flourensia cernua)." J. Arid Environ. 39(1), 57-78.
Gouyet, J-F. 1996. Physics and Fractal Structures. Springer-Verlag, Inc. New York, NY. 234pp.
Gray, W.G., A. Leijnse, R.L. KolarandC.A. Blain. 1993. Mathematical Tools for Changing Spatial Scales in the
Analysis of Physical Systems. CRC Press, Boca Raton, FL. 232 pp.
Groot, R.D. and P.B. Warren. 1997. "Dissipative particle dynamics: Bridging the gap between atomistic and
mesoscopic simulation. J. Chem. Phys. 107(11), 4423-4435.
Guymon. G.L. 1994. Unsaturated Zone Hydrology. PTR Prentice Hall. Englewood Cliffs, New Jersey. 210pp.
Hamed, M.M. and P.B. Bedient. 1999. Reliably-Based Uncertainty Analysis of Groundwater Contaminant
Transport and Remediation. EPA/600/R-99/028. National Risk Management Research Laboratory, ORD, US.
EPA. Cincinnati, OH.
Helton, J.C. 1993. "Uncertainty and sensitivity analysis techniques for use in performance assessment for
radioactive waste disposal" Reliability Eng. and System Safety 42. 327-367.
Hills, R.G., PJ. Wierenga, D.B. Hudson, and M.R. Kirkland. 1991. "The second Las Graces trench experiment:
experimental results and two-dimensional flow predictions." Water Resour. Res. 27(10), 2707-2718.
Hills, R.G., K.A. Fisher, M.R. Kirkland, and PJ. Wierenga. 1994. "Application of flux-corrected transport to the
Las Graces trench site." Water Resour. Res. 30(8), 2377-2385.
Hornung, U. 1997. Homogenization and Porous Media. Springer-Verlag, Inc. New York, NY. 279pp.
Huyakorn, P.S., and IE. Buckley. 1987. Finite Element Code for Simulating One-Dimensional Flow and Solute
Transport in the Vadose Zone. Technical Report prepared for U.S. EPA, Athens, GA.
Jang, Y-S, N. Sitar and A. der Kiureghian. 1994. "Reliability analysis of contaminant transport in saturated porous
media." Water Resour. Res. 30(8), 2435-2448.
Jenkins, M.B., R. A. Virgina, and W.M. Jarrell. 1988. "Depth distribution and seasonal populations of mesquite-
nodulating rhizobia in warm desert ecosystem." Soil Sci. Soc. Am. J. 52(6), 1644-1650.
Jensen, M.E.,R.D. Burman, and R.G. Allen. 1990. "Evapotranspiration and irrigation water requirements." ASCE
Manuals and Reports on Engineering Practice. No. 70. American Society of Civil Engineers, NY. 332 pp.
Jury, W. A. and R.L. Valentine. 1986. "Transport mechanisms and loss pathways for chemicals in soil." In: S.C.
Hern and S.M. Melancon (eds.), Vadose Zone Modeling of Organic Pollutants. Lewis Publishers, Inc., Chelsea,
Michigan, pp. 37-60.
104
-------
Kaiser, C. 1997. "A directed percolation model for clogging in a porous medium with small inhomogeneities."
Transport in Porous Media 26. 133-146.
Knox, R.C., D.A. Sabatini, and L.W. Canter. 1993. Subsurface Transport and Fate Processes. Lewis Publishers.
Boca Raton, FL. 430pp.
Kretzschmar, R., K. Barmettler, D. Grolimund, Y-d. Yan, M. Borkovec and H. Stecher. 1997. "Experimental
determination of colloid deposition rates and collision efficiencies in natural porous media." Water Resour. Res.
33(5), 1129-1137.
Litaor, M.I., G. Earth, E.M. Zika, G. Litus, J. Moffitt and H. Daniels. 1998. "The behavior of radionuclides in the
soils of Rocky Flats, Colorado." J. Environ. Radioactivity. 38(1), 17-46.
Liu, S., W.B. Mills, andB. Baena. 1995. Multimedia Exposure Assessment Modeling Including Fate and
Tranformation Products (MULTIMED_DP 1.0): Implementation. Tests and User's Manual. Technical Report
prepared by U.S. EPA, Athens, GA.
Logan, J.D. 1996. "Solute transport in porous media with scale-dependent dispersion and periodic boundary
conditions." Journal of Hydrology. 184,261-276.
Logan, J.D. 2001. Transport Modeling in Hydrochemical Systems. Springer-Verlag. New York, NY. 223pp.
Luckner, L., and W.M. Schestakow. 1991. Migration Processes in the Soil and Groundwater Zone. Lewis
Publishers. Chelsea, MI. 485 pp.
Mackay, D. 2001. Mulimedia Environmental Models: The Fugacity Approach. Second Edition. Lewis Publishers.
Boca Raton, FL. 261 pp.
Mayer, S., T.R. Ellsworth, D.L. Corwin and K. League. 1999. "Identifying effective parameters for solute transport
models in heterogeneous environments." Assessment of Non-Point Source Pollution in the Vadose Zone.
Geophysical Monograph 108. American Geophysical Union. Washington, DC. 119-133.
McCarthy, J.F. 1995. "Comparison of fast algorithms for estimating large-scale permeabilities of heterogeneous
media." Transport in Porous Media. 19, 123-1137.
McCarthy, J.F. and J.M. Zachara. 1989. "Subsurface transport of contaminants: Mobile colloids in the subsurface
environment may alter the transport of contaminants." Environ. Sci. Technol. 23(5), 496-502.
McCord, J.T., C.A. Gotway and S.H. Conrad. 1997. "Impact of geologic heterogeneity on recharge estimation using
environmental tracers: Numerical modeling investigation." Water Resour. Res. 33(6), 1229-1240.
Mehl, S. andM.C.Hill. 2001. "A comparison of solute-transport solution techniques and their effect on sensitivity
analysis and inverse modeling results." Ground Water. 39(2), 300-307.
Meyer, P.O., M.L. Rockhold and G. W. Gee. 1997. Uncertainty Analysis of Infiltration and Subsurface Flow and
Transport for SDMP Sites. NUREG/CR-6565. PNNL-11705. United States Nuclear Regulatory Commission.
Division of Regulatory Applications. Washington, DC.
Miller, E.E. and R.D. Miller. 1956. "Physical theory for capillary flow phenomena." J. App. Phys. 27, 324-332.
Miralles-Wilhelm, F. and L.W. Gelhar. 1996. "Stochastic analysis of transport and decay of a solute in
heterogeneous aquifers." Water Resour. Res. 32(12), 3451-3459.
Miyazaki, T., S. Hasegawa and T. Kasubuchi. 1993. Water Flow in Soils. Marcel Dekker, Inc. New York, NY.
296 pp.
Morel-Seytoux, H.J., P.O. Meyer, M. Nachabe, J. Touma, M. Th. Van Genuchten and R.J. Lenhard. 1996.
"Parameter equivalence for Brooks-Corey and van Genuchten soil characteristics: preserving the effective
capillary drive." Water Resour. Res. 32(5). 1251-1258.
105
-------
Mulla, D. J. and T.M. Addiscott. 1999. "Validation approach for field-, basin-, and regional-scale water quality
models." Assessment of Non-Point Source Pollution in the Vadose Zone. Geophysical Monograph 108.
American Geophysical Union. Washington, DC. 63-78.
Murray, ID. 1989. Mathematical Biology. Springer-Verlag. New York, NY. 767pp.
Nachabe, M.H. 1996. "Macroscopic capillary length, sorptivity, and shaping factor in modeling the infiltration rate."
Soil Sci. Soc. Am. J. 60, 957-962.
Neuman, S.P. and V.D. Federico. 1998. "Correlation, flow, and transport in multiscale permeability fields." In
Scale Dependence and Scale Invariance in Hydrology. Cambridge University Press. New York, NY. 354-397.
Nirmalakhandan, N. 2002. Modeling Tools for Environmental Engineers and Scientists. CRC Press. Boca Raton,
FL. 312pp.
Nofziger, D.L., J-S. Chen, and C.R. Haan. 1994. Evaluation of Unsaturated/Vadose Zone Models for Superfund
Sites. EPA/600/R-93/184. U.S. EPA. R.S. Kerr Environmental Research Laboratory, Ada, OK. 188pp.
Oldenburg, C.M. and K. Pruess. 1996. "Mixing with first-order decay in variable-velocity porous media flow."
Transport in Porous Media. 22, 161-180.
Porro, I., and P.J. Wierenga. 1993. "Transient and steady state solute transport through a large unsaturated soil
column." Ground Water. 31(2), 193-200.
Preuss, K., B. Faybishenko and G.S. Bodcarsson. 1999. "Alternative concepts and approaches for modeling flow
and transport in thick unsaturated zones of fractured rocks." Journal of Contain. Hydrology. 38, 281-322.
Probstein, R.F. 1994. Physicochemical Hydrodynamics: An Introduction. Second Edition. John Wiley and Sons,
Inc. New York, NY. 400pp.
Ravi, V. and J. Williams. 1998. Estimation of Infiltration Rate in the Vadose Zone: Compilation of Simple
Mathematical Models. Volume I. EPA/600/R-97/128a. National Risk Management Research Laboratory. U.S.
EPA Cincinnati, OH.
Renard, Ph., and G. de Marsily. 1997. "Calculating equivalent permeability: a review." Advances in Water
Resources. 20 (5-6), 253-278.
Richards, L.A. 1931. "Capillary conduction of liquids through porous mediums." Physics. 1,318-331.
Rockhold, M.L. 1999. "Parameterizing flow and transport models for field-scale applications in heterogeneous,
unsaturated soils." Assessment of Non-Point Source Pollution in the Vadose Zone. Geophysical Monograph
108. American Geophysical Union. Washington, DC. 243-260.
Rockhold, M.L., R.E. Rossi, and R.G. Hills. 1996. "Application of similar media scaling and conditional simulation
for modeling water flow and tritium transport at the Las Graces trench site." Water Resour. Res. 32(3), 595-
609.
Rockhold, M.L., C.S. Simmons and M.F. Payer. 1997. "An analytical solution technique for one-dimensional,
steady vertical water flow in layered soils." Water Resour. Res. 33(4), 897-902.
Royo, A.R. 2000. "Desert Plant Survival." DesertUSA Newsletter. Http://www.desertusa.com/du_plantsurv.html.
5pp.
Sachdev, P.L. 2000. Self-Similarity and Beyond: Exact Solutions of Nonlinear Problems. Chapman and Hall/CRC.
Boca Raton, FL. 319pp.
Salhotra, A.M., P. Mineart, S. Sharp-Hansen, T. Allison, R. Johns, and W.B. Mills. 1995.MulimediaExposure
Assessment Model (MULIMED 2.0) for Evaluating the Land Disposal of Waste-Model Theory. Final Report
prepared for U.S. EPA. Athens, GA.
106
-------
Sharp-Hansen, S., C. Travers, P. Hummel, T. Allison, R. Johns, and W. Mills. 1995. A Subtitle D Landfill
Application Manual for the Multimedia Exposure Assessment Model (MULTIMED 2.OX Technical Report
prepared for U.S. EPA. Athens, GA.
Shurbaji, A.-R. and A.R. Campbell. 1997. "Study of evaporation and recharge in desert soil using environmental
tracers, New Mexico, USA." Environ. Geology. 29 (3/4). 147-151.
Simon, W., P. Reichert and C. Hinz. 1997. "Properties of exact and approximate traveling wave solutions for
transport with nonlinear and nonequilibrium sorption." Water Resour. Res. 33(5). 1139-1147.
Simunek, J. and M.Th. van Genuchten. 1994. The CHAIN_2D Code for Simulating Two- Dimensional
Movement of Water. Heat, and Muliple Solutes. Research Report No. 136. U.S. Salinity Laboratory. USD A,
ARS. Riverside, CA.
Simunek, J., K. Huang, and M.Rh. van Genuchten. 1998. The HYDRUS Code for Simulating the One-Dimensional
Movement of Water. Heat, and Multiple Solutes in Variable-Saturated Media. Research Report No. 144. U.S.
Salinity Laboratory. USD A, ARS. Riverside, CA.
Spanos, T.J.T. 2002. The Thermophysics of Porous Media. Chapman & Hall/CRC. Boca Raton, FL. 220 pp.
Sposito, G. (Ed.). 1998. Scale Dependence and Scale Invariance in Hydrology. Cambridge University Press. New
York, NY. 423pp.
Stehfest, H. (1970). "Numerical inversion of Laplace Transforms." Commun. ACM. 13(1). 47-49.
Stockman, H. W. 1997. "A lattice gas study of retardation and dispersion in fractures: Assessment of errors from
desorption kinetics and buoyancy." Water Resour. Res. 33(8), 1823-1831.
Thibodeaux, LJ. 1996. Environmental Chemodynamics: Movement of Chemicals in Air. Water, and Soil. 2nd
Edition. John Wiley & Sons, Inc. New York, NY. 593 pp.
Tomaski, D., M. Reeves, B.A. Kelley, and J.F. Pickens. 1989. "Parameter sensitivity and importance for
radionuclide transport in double-porosity systems." In: B.E. Buxton (ed.) Proceedings of the Conference on
Geostatistical. Sensitivity, and Uncertainty Methods for Ground-Water Flow and Radionuclide Transport
Modeling. San Fransisco, CA. September 15-17, 1987. Battelle Press. Columbus, OH. pp. 297-321.
Toride, N., F.J. Leij, and M.T. van Genuchten. 1993. "A comprehensive set of analytical solutions for
nonequilibrium solute transport with first-order decay and zero-order production." Water Resour. Res. 29(7),
2167-2182.
Turcotte, D.L. 1992. Fractals and Chaos in Geology and Geophysics. Cambridge University Press. New York,
NY. 221pp.
U.S. EPA. 1995a. EPA's Composite Model for Leachate Migration with Transformation Products. EPACMTP.
Background Document. Office of Solid Waste, U.S. EPA, Washington, DC.
U.S. EPA. 1995b. EPA's Composite Model for Leachate Migration with TransformationProducts. EPACMTP.
User's Guide. Office of Solid Waste, U.S. EPA, Washington, DC.
U.S. EPA. 1996. Summary Report for the Workshop on Monte Carlo Analysis. EPA/630/R-96/010. Risk
Assessment Forum, U.S. EPA, Washington, D.C.
U.S. EPA. 1999a. Understanding Variations in Partition Coefficients. Kd. Values: The Kd Model. Methods of
Measurement, and Application of Chemical Reaction Codes. Volume I. EPA 402-R-99-004A. Office of
Radiation and Indoor Air & Office of Solid Waste and Emergency Response. U.S. EPA. Washington, D.C.
U.S. EPA. 1999b. Understanding Variation in Partition Coefficient. Kd. Values: Review of Geochemistry and
Available Kd Values for Cadmium. Cesium. Chromium. Lead. Plutonium. Radon. Strontium. Thorium. Tritium.
ffl). and Uranium. Volume II. EPA402-R-99-004B. Office of Radiation and Indoor Air & Office of Solid
Waste and Emergency Response. U.S. EPA. Washington, D.C.
107
-------
U.S. EPA. 2000a. Soil Screening Guidance forRadionuclides: User's Guide. EPA/540-R-00-007. Office of
Radiation and Indoor Air and Office of Solid Waste and Emergency Response, U.S. EPA, Washington, D.C.
U.S. EPA. 2000b. Soil Screening Gudiance for Radionuclides: Technical Background Document. EPA/540-R-00-
006. Office of Radiation and Indoor Air and Office of Solid Waste and Emergency Response, U.S. EPA,
Washington, D.C.
Van de Yen, T.G.M. 1989. Colloidal Hydrodynamics. Academic Press. New York, NY. 582 pp.
VanderHeijde, P.K.M. 1994. Identification and Compilation of Unsaturated/Vadose Zone Models. EPA/600/R-
94/028. R.S. Kerr Environmental Research Laboratory, U.S. EPA, Ada, OK.
VanderHeijde, P.K.M. and O.A. Elnawawy. 1993. Compilation of Ground-Water Models. EPA/600/R-93/118.
R.S. Kerr Environmental Research Laboratory, U.S. EPA, Ada, OK.
VanderHeijde, P.K.M. and D. A. Kanzer. 1995. Ground-Water Model Testing: Systematic Evaluation and Testing
of Code Functionality. Performance, and Applicability to Practical Problems. IGWMC-GWMI 95-01.
International Ground Water Modeling Center,Colorado School of Mines, Golden, CO.
Van Genuchten, M. Th. 1980. " A closed-form equation for predicting the hydraulic conductivity of unsaturated
soils." Soil Sci. Soc. Am. J. 44, 892-898.
Van Genuchten, M. Th. 1985. "Convective-dispersive transport of solutes involved in sequential first-order decay
reactions." Computer & Geosciences. 11(2), 129-147.
Van Genuchten, M. Th., and RJ. Wagenet. 1989. "Two-site/two-region models for pesticide transport and
degradation: theoretical development and analytical solutions." Soil Sci. Soc. Am. J. 53, 1303-1310.
Visintin, A. 1994. Differential Models of Hysteresis. Springer-Verlag. New York, NY. 407pp.
Vogel, T., and M. Cislerova. 1988. "On the reliability of unsaturated hydraulic conductivity calculated from the
moisture retention curve." Transport in Porous Media. 3, 1-15.
Vogel, T., M. Cislerova, and J.W. Hopmans. 1991. "Porous media with linearly variable hydraulic properties.
Water Resour. Res. 27(10), 2735-2741.
Vollmayr, H., F. Kleint and G. Schuurmann. 1997. "Discrete modeling of water and pesticide movement in soil."
Water Resour. Res. 33(7), 1743-1747.
Wang, Z., J. Feyen, D.R. Nielson, and M.T. van Genuchten. 1997. "Two-phase flow infiltration equations
accounting for air entrapment effects." Water Resour. Res. 33(12), 2759-2767.
Wang, Z., J. Feyen and D.E. Elrick. 1998. "Prediction of fingering in porous media." Water Resour. Res. 34(9),
2183-2190.
Wang, Z., J. Feyen and C.J. Ritsema. 1998. "Susceptibility and predictability of conditions for preferential flow."
Water Resour. Res. 34(9), 2169-2182.
Wang, Z., J. Feyen, M.T. van Genuchten, and D.R. Nielsen. 1998. "Air entrapment effects on infiltration rate and
flow instability." Water Resour. Res. 34(2), 213-222.
Wang, Z., L. Wu, and QJ. Wu. 2000. "Water-entry value as an alternative indicator of soil water-repellency on
infiltration rate and flow instability." J. Hydrol. 231-232. 76-83.
Wang, Z., Q. J. Wu, L. Wu, C.J. Ritsema, L. W. Dekker and J. Feyen. 2000. "Effects of soil water-repellency on
infiltration rate and flow instability." J. Hydrol. 231-232, 265-276.
Warrick, A.W., G.J. Mullen and D.R. Nielsen. 1977. "Scaling field-measured soil hydraulic properties using a
similar media concept." Water Resour. Res. 13,355-362.
108
-------
Weber, W.J., Jr., P.M. McGinley, and L.E. Katz. 1991. "Sorption phenomena in subsurface systems: concepts,
models and effects on contaminant fate and transport." Water Res. 25, 499-528.
Wierenga, P.J., R.G. Hills, and D.B. Hudson. 1991. "The Las Graces trench site: characterization, experimental
results, and one-dimensional flow predictions." Water Resour. Res. 27(1), 2695-2705.
Williams, J.R., Y. Ouyang, J.-S. Chen, and V. Ravi. 1998. Estimation of Infiltration Rate in Vadose Zone:
Application of Selected Mathematical Models. Volume II. EPA/600/R-97/128b. National Risk Management
Research Laboratory. U.S. EPA. Cincinnati, OH.
Yao, T-M. and J.M.H. Hendrickx. 1996. "Stability of wetting fronts in dry homogeneous soils under low infiltration
rates." Soil Sci. Soc. Am. J. 60, 20-28.
Yortsos, Y.C. andK. Shing. 1999. Evaluation and Analysis of Microscale Flow and Transport During Remediation.
EPA/600/R-99/022. National Risk Management Research Laboratory, Office of Research and Development,
U.S. EPA, Cincinnati, OH. 160pp.
Zheng, C., and G.D. Bennett. 1995. Applied Contaminant Transport Modeling: Theory and Practice. VanNostrand
Reinhold. New York, NY. 440pp.
109
-------
Appendix A
Empirical Models of the Unsaturated Soil Hydraulic Properties
Which Are Used in the Various Models
One-dimensional soil-water movement in a partially saturated rigid porous medium is described in the HYDRUS
Code(Simunek etal, 1998) by a modified form of the Richards Equation:
(A-l)
where 9 is the soil-water volumetric content (L3L~3), h is the soil-water pressure head (L), a negative quantity for
unsaturated conditions, K is the unsaturated hydraulic conductivity function (UP1), t is time (T), z is the vertical
coordinate (L), positive upward with an origin at some finite depth in the soil, and S(h,z) is a sink term defined as the
volume of water removed from a unit volume of soil per unit time due to plant water uptake (L3L-3T-!). The
formulation of Equation (A-l) assumes that the air phase plays an insignificant role in the liquid flow process and
that water flow due to thermal gradients can be neglected. The general form of K is given by :
K = K(h,z) = Ks(z) K,.(h,z) , (A-2)
where Ks is the saturated hydraulic conductivity (LT-1) and K,. is the relative hydraulic conductivity (1).
The HYDRUS Code features three models for the unsaturated soil hydraulic properties, 9 and K. In
homogeneous media, these properties can be expressed as nonlinear functions of the soil-water pressure head, 9(h)
and K(h). The three models are analytical expressions denoted by BC, VG, and VC, respectively given by Brooks and
Corey (1964), van Genuchten(1980), and Vogel and Cislerova (1988).
The BC-Model is defined by the following equations:
S =•
ah" , h < -I/a
(A'3)
K = KsSfn + c + 2, (A-4)
where Se is the effective water content, 9S is the saturated water content, 9r is the residual water content, a is the
inverse of the air-entry value or bubbling pressure, n is a pore-size distribution index, and c is a pore-connectivity
parameter (c was assumed to be 2.0 in the original Brooks-Corey study). Thus, this model consists of six free
parameters (9r, 9S, Ks, n, c, a). The quantities (9r, 9S, Ks) can be measured in the laboratory or in the field, while the
parameters (n,c,a) are considered to be empirical coefficients affecting the shapes of the hydraulic functions.
The VG-Model is defined by the analytical expressions given below:
A-l
-------
0(h) =
-,h<0
h>0
(A-5)
K(h)=KsS:[l-(l-Sl/m)mJ,
where
m=l-l/n, n
(A-6)
(A-7)
In some of the codes, other than that of HYDRUS, the parameter m is denoted by y, and the parameter n is denoted by
P. The pore-connectivity parameter, c, has been estimated to be about 0.5 as an average for many soils (Mualem,
1976). As in the BC-Model, this model contains the six parameters (9r, 9S, Ks, n, c, a).
The third set of hydraulic equations implemented in HYDRUS are the defining equations of the VC-Model, which
are modifications of the VG-Equations to add flexibility in the description of the hydraulic properties near saturation.
These modified equations are as follows:
0(h) =
9m-9a
hh
(A-8)
K(h)=
KsKr(h) ,
Y ,(h-hk)(Ks
. S '
h hs _
(A-9)
where
K.
'F(9r)-F(9)f
.F(9r)-F(6k)J
(A-10)
F(0) =
1 -
(A-H)
(A-12)
This model allows for a nonzero minimum capillary height, h,, by replacing 9S in the VG-Model by an extrapolated
parameter, 9m, which is slightly more than 9S (Figure A-la). In general, this change from 9S to 9m has little to no
effect on the soil-water retention curve; but its effect on the shape and value of K may be considerable, especially for
fine-textured soils when n is small (i.e., 1.0 < n < 1.3). A further increase in the flexibility of the analytical expressions
of the VC-Model is experienced by replacing 9r in the equation for 9(h) by another extrapolated parameter 9a < 9r.
A-2
-------
o
.2
"c
o
U
0
I
(a)
1% 0
o
3
T3
O
U
.O
"5
D
T3
h. 0
Pressure Head, h Pressure Head, h
Figure A-l. Schematics of the soil-water retention curve (a) and the hydraulic conductivity
function (b) for the VC-Model (from Simuneketal, 1998).
This model retains the physical meaning of 9r and 9S as measurable quantities while adding two new fitting
parameters, 9a and 9m. Equation (A- 10) assumes that the predicted hydraulic conductivity function is matched to a
measured value of the function, Kk = K(9k), at some water content 9k satisfying
9k<9s, Kk 9m~ 9k ~ 9S= Kk-Ks
the VC-Model reduces to the VG-Model.
(A-14)
The question arises, "Which model should a user employ?" The claim is that the VC-Model is more versatile
than the VG-Model, but three more unknown parameters are required to use the VC-Model. Of the BC- and VG-
Models, which is the best to use? Morel-Seytoux et al. (1996) defined an equivalence which provides a simple way to
convert BC-parameters into VG-parameters, and vice versa, when one set of parameters is known and preservation of
the maximum value of a physical characteristic called the "effective capillary drive" is guaranteed. This quantity is
defined by:
r'-rwdhc, (A-15)
o
where k^, is the relative permeability (or conductivity) to water and hc is the capillary pressure head, a positive
quantity equal to minus h. This equivalence makes infiltration capacity calculations insensitive to the model used
(BC or VG) to represent the soil hydraulic properties. It is strictly a matter of convenience for the user which
expressions are used. As one would expect, this equivalence only applies for situations which lead to high water
contents in some part of the domain of interest. However, because of the equivalence criterion which preserves the
asymptotic behavior of capillary pressure at low water contents and thus preserves the hydraulic gradient, the
authors claim there are at least plausible grounds that the parameter equivalence will be reasonably satisfactory for
situations involving drainage and evapotranspiration. But, the application of this paper to these latter processes
should be further tested.
Further support for the equivalence relationship of Morel-Seytoux et al. (1996) is given by the work of Nachabe
(1996). The work of this author was based on the fact that infiltration tests in the field routinely involve
measurements of sorptivity and macroscopic capillary length, the latter being related to the effective capillary drive
defined in Equation A-15. In addition, these two factors are strongly related through a shape factor b, which is a
measure of the nonlinearity of the soil water hydraulic diffusivity D(9), (L2T-!):
A-3
-------
d0
The macroscopic capillary length Xs, which is equivalent to the wetting front suction head, is defined by
J D(0) d0 ,
0)-K(0n)
where 90 is the supply water content applied at the surface and 9n is the dry, antecedent water content in the soil,
ahead of the wetting front. The sorptivity, S ,is well approximated by
S2 = J(00+0-20n)D(0)d0, (A.18)
en
where S , (LT~1/2), is defined by the early stages of an infiltration test. In Nachabe's study, relationships were
developed between hs, S, and shape factor b, and the parameters of the BC- and VG-Models for K(9) and D(9).
Numerical simulations using a dimensionless form of Richards Equation showed that the predicted infiltration rate is
not very sensitive to small variations in the shape factor, b, provided ^s is the same in all the simulations. This result
is important because b is difficult to determine accurately in the field. Thus, this similarity of the dimensional
infiltration solutions under small variations inb shows that ^s is a scale-factor (see Appendix B) because:
(1) ^s renders predictions of infiltration rates fairly insensitive to the models used for K(9) and D(9); and
(2) ^s reduces the number of parameters needed to characterize infiltration.
Therefore, given the generalized infiltration solution from this scale analysis, an infiltration curve into a particular soil
at a particular location and time can be deduced by choosing the proper units in time and space (i.e., scaling) of this
generalized solution, thus minimizing the computations for obtaining infiltration curves over spatially varying
domains. Further, the small range of variation of the shape factor, b, and the insensitivity to this variation by \-
scaled infiltration rates based on the BC- and VG-Models used in the Richards Equation supports the use of linear
scaling in the HYDRUS Code, see Appendix B.
In spite of the results obtained by the two previous investigations, researchers are still looking for better soil-
water retention curve (WRC) models. Assouline et al. (1998) introduced a new conceptual model for WRC's. Before
introducing this new concept, the authors gave a brief, but comprehensive, survey of WRC construction from the
work of Brooks and Corey (1964) to recent studies based on fractal and multifractal soil structures. The proposed
concept and model for the WRC given by Assouline et al. (1998) are based on two assumptions:
(1) A soil structure resulting from a uniform random fragmentation process where the probability of
fragmentation of an aggregate of soil material is proportional to its size; and
(2) A power function describing a relationship between the volume of the aggregates and the corresponding
pore volumes.
This new model WRC was tested on 12 different soil types that presented a wide range of textures and
mechanical analysis. The model accurately reproduced the very different shapes of the respective WRC's, from the
step-function shape of a sand to the sigmoidal almost linear shape of a clay. Also, the new model accurately
reproduced measured data over the whole range of water contents, from saturation conditions to water contents at
the wilting point. Compared with the VG-Model and another popular WRC-Model, the new model of Assouline et al.
exhibited increased flexibility and reproduced the measured data in the high and low water content ranges better than
the two other models. Thus, the authors claim two major benefits for the new WRC-Model:
(1) For this new two fitting parameter model, there is a gain in accuracy and flexibility over competing two fitting
parameter models, such as the BC- and VG-Models.
A-4
-------
(2) The conceptual basis on which the WRC-Model is constructed forms a basis for further research toward a
more physically based relationship between soil structure and its corresponding hydraulic properties.
References
Assouline, S, D.Tessier and A. Bruand. 1998. "A conceptual model of the soil water retention curve." Water Resour.
Res.. 34 (2\ 223-231.
Brooks, R. H. andA.T. Corey. 1964. Hydraulic Properties of Porous Media. Hydrology Paper No. 3. Colorado State
University. Fort Collins, CO.
Morel-Seytoux, H.I, RD. Meyer, M. Nachabe, J. Touma, M. Th. van Genuchten and R. J. Lenhard. 1996. "Parameter
equivalence for Brooks-Corey and van Genuchten soil characteristics: Preserving the effective capillary drive."
Water Resour. Res.. 32 (5), 1251-1258.
Mualem, Y. 1976. " A new model for predicting the hydraulic conductivity of unsaturated porous media." Water
Resour. Res.. 12(3), 513-522.
Nachabe, M. H. 1996. "Macroscopic capillary length, sorptivity, and shape factor in modeling the infiltration rate."
Soil Sci. Soc. Am. J. 60,957-962.
Simunek, I, K.Huang and M.Th. van Genuchten. 1998. The HYDRUS Code for Simulating the One-Dimensional
Movement of Water. Heat. andMultiple Solutes in Variably-Saturated Media. Version6.0. Res. Rep. No. 144.
U.S. Salinity Laboratory, Ag. Res. Ser, USD A, Riverside, CA.
van Genuchten, M. Th. 1980. "A closed-form equation for predicting the hydraulic conductivity of unsaturated
soils." Soil Sci. Soc. Am. J.. 44,892-898.
Vogel, T andM. Cislerova. 1988. "On the reliability of unsaturated hydraulic conductivity calculated from the
moisture retention curve." Transport in Porous Media. 3,1-15.
A-5
-------
Appendix B
A Discussion on the Scaling of Field Soil-Water Behavior
The scaling features for soil-water behavior found in the HYDRUS and CHAIN 2D Codes are part of a bigger
picture in science which includes the concepts of symmetry, invariance, and conservation. Further, it is difficult to
meaningfully discuss scaling without involving some of these more fundamental concepts. For example, Skoglund
(1967) defines scaling as a technique for relating the variables of a "prototype" to those of a "model" at
corresponding points and times in the two systems. Further, he says that scaling often requires that the ratios of the
corresponding dimensional variables of the "prototype" and the "model" at corresponding points and times be
constant. These constants are called scaling factors and are commonly used in the design and interpretation of
physical and mathematical models. But Skoglund says that this often obscures the more fundamental process of
searching for symmetry in nature. In what follows, the roles between scaling, symmetry, invariance, and
conservation are explained.
B.I Symmetry in Nature
Rosen (1995) says a symmetric situation is one which possesses the possibility of change; and if such a change
occurs, some aspect of the situation remains unchanged. Thus, the situation is said to be symmetric under the
change, or transformation, with respect to that aspect. Consequently, there are two essential components of
symmetry in nature (Rosen, 1995):
1. Possibility of change - it must be possible to perform the change, although the change does not actually
have to be performed.
2. Immunity - some aspect of the natural situation remains unchanged if the change is actually performed.
If a change is possible but some aspect of the situation is not immune to it, asymmetry exists. If there is no
possibility of change, with reference to the physical reality of the situation, then the concepts of symmetry and
asymmetry are inapplicable.
In the geophysical sciences, as opposed to several situations in fundamental particle physics, pure mathematical
symmetry is not the rule; but approximate symmetry does exist. That is, in natural situations there are aspects which
are approximately immune to possible transformations. There is no approximation in the transformation or in its
possibility; it must indeed be possible to perform the transformation. The approximation is in the immunity; that is,
some aspect of the natural situation must change by only a "little" when some transformation is performed (Rosen,
1995). The concept of "little" is based on the variations of the aspect in question being much less in magnitude than
some characteristic length, time, mass, velocity, or what have you. Thus, the situation is said to be approximately
symmetric under the transformation with respect to that aspect. In spite of the lack of perfect symmetry in the
biological and geophysical sciences, the search for symmetry is a fruitful undertaking as it has been and is in the
more fundamental areas of physics (Rosen, 1995).
As a science matures its reproductive and predictive features lead to the formulation of fundamental laws which
are usually expressed in mathematical statements and equations given in some system of coordinates (i.e., a reference
frame). Hence, Rosen (1995) states that the searchfor symmetry in nature consists of transformations, immunity from
these changes, frames of reference, and non-immunity from the changes (or asymmetry). This search requires a
general formalism of symmetry. Rosen shows that this formalism is couched in the mathematical language of
transformations and symmetry, which is known as group theory. As geometry is the appropriate mathematics for
investigating space, and calculus is the appropriate mathematics for investigating motion, group theory is the
mathematics of symmetry. It is essential for the search of symmetry in all fields of science, including investigations in
soil science and subsurface geology.
B-l
-------
B.2 Similitude, Transformation Groups, Inspectional Analysis, Self-Similarity
The search for symmetry in physical systems invokes the use of concepts such as similitude, dimensional
analysis, similar physical systems, inspectional analysis, transformations, groups, and self similarity. Skoglund (1967)
says that similitude includes both the concepts of dimensional analysis and similar physical systems. Dimensional
analysis is the special mathematics of units. All scientific measurements are referred to four internationally
recognized fundamental or primary units of measurement: the kilogram, meter, second, and degree Kelvin. The
fundamental or primary concepts of physical things and systems are mass M, length L, time T, and temperature 0.
These units of measurement and the primary concepts, although somewhat arbitrarily chosen, are the minimum
number required for all macroscopic measurements of physical phenomena. From a microscopic point of view (i.e.,
individual molecules, atoms and subatomic particles), the fundamental units are the kilogram, meter and second.
Thus, both viewpoints, macroscopic and microscopic, can be treated as one with respect to similitude concerns.
Dimensions comprise a class of units for the same variable, where the primary dimensions are M, L, T and 0.
For example, the universal gas constant in units of cal/degree K has the physical dimensions of ML2/T2 0. The
fundamental principle of dimensional substitution (Skoglund, 1967) states that any "complete equation of physics"
is satisfied when the units or dimensions of the variables are substituted for the variables themselves. Thus, any
"complete equation" can be converted to a nondimensional form by dividing by a variable or group of variables with
the same units as each of the sides of the equation. This forms a corollary of the principle of dimensional
substitution, namely, that any "complete equation" can be converted into a nondimensional form. Therefore, any
"complete equation" is independent of units, or unit free. "Complete equations" are also said to be dimensionally
homogeneous. In essence, the above properties of a complete physical equation form a set of equivalent definitions
for such an entity.
Birkhoff (1955) states that it is often claimed that "only dimensional homogeneous equations can be regarded as
having fundamental physical significance." However, not all physically significant equations are true in all units! In
fact, there exist no known fundamental units with respect to which all known physical laws are unit free. Birkhoff
says no "Swift's Lilliput" can exist, for to do so would require changing some of the fundamental physical constants,
such as the speed of light and Planck's constant in quantum mechanics. Nevertheless, for the current applications
we are not limited in this regard. Birkhoff shows that for normal velocity distributions and temperature variations, the
Newtonian mechanics of fluids possess as independent, fundamental units the quantities of mass, length, time and
temperature. Thus, the above concepts of complete physical equations and dimensional substitution are valid for
our current applications.
Two physical systems are said to be similar to each other if their pertinent, corresponding variables (e.g.,
velocity, acceleration, mass distribution, force fields and temperature fields) are proportional at corresponding
locations and times. For example, we have previously introduced the terms prototype and model, ^prototype of a
given problem is the physical system itself. A model of the problem is any system that is similar to it.
Mathematically speaking, and this is the only way that similarity can be precisely stated, the first requirement for the
similarity of two systems is that their nondimensional governing equations be the same. Further, to guarantee that
two similar systems possess the same nondimensional solutions it is required that both systems be well-formulated
or well-posed (Flavin and Rionero, 1996) and that their nondimensional initial and boundary conditions be equivalent
at corresponding points and times in their respective geometric domains. Thus, the complete establishment of
similarity, its existence and its properties, dictates that these requirements be met. Anything less than this leads to a
similarity between systems that is less than whole.
The various levels of similarity of importance in continuum mechanics are described in Skoglund (1967). A
model is geometrically similar to a prototype if all the geometric features of the model can be related to those of the
prototype by a single scale factor, as illustrated in Figure B-1. Kinematic similarity is concerned with similar
"particles" (solid or fluid particles) possessing similar velocities and accelerations at corresponding points and times
in similar geometries. Dynamic similarity is similarly defined for systems involving masses, forces, velocities,
accelerations, geometries, and time. Finally, the thermodynamic similarity of a prototype continuum and its
continuum model is based on the requirements and results derived from the principles of the conservation of mass,
momentum and energy, along with pertinent initial and boundary conditions for the dependent variables.
Mathematically speaking (obviously, there is also a very strong experimental component of the search), the
search for symmetry in the fields of continuum mechanics, which includes subsurface processes, first consists of
writing down the governing equations and initial and boundary conditions on which the area of continuum
mechanics in question is based. Then this system is tested for its invariance (i.e., similarity) under groups of
B-2
-------
Figure B-l. Geometrically similar figures, where (a) is the reference figure with characteristic lengthL*, (b) is a
similar figure with characteristic length!^, L* = ocjLj, with scale factor ocj = 2, and (c) is a similar
figure with a2L2 = L*, a2 = 1A (reprinted from UnsaturatedZone Hydrology, 1994, by Gary L.
Guymon, with permission of Pearson Education, Inc., Upper Saddle River, NJ).
transformations of interest. This procedure is called inspectional analysis by Birkhoff (1955). Forexample, two fluid
motions denoted by O and O' are dynamically similar if they can be described by Newtonian coordinate systems (i.e.,
systems in which Newton's laws of motion are valid) which are related by the transformation of space, time, and mass
of the form:
x' = a x, y' = a y, z' = a z, t' = P t, m' = yin,
(B-l)
where (a, P, y) are constant scale factors, and x, y, z, t, m have their usual meanings of length, time and mass. If T,
T(O) = O', represents the transformation in Equation (B-l) with scale factors (a, P, y), then T'1, T'1 (O') = O,
represents its inverse transformation with scale factors (or1, P'1, y1). The successive application of T and T'1, or T'1
and T, leads back to the original fluid motion. Thus, the product T^T or TT'1 defines the identity transformation I
with scale factors (1,1,1). Suppose we have two transformations of the type given in Equation (B-l): T1 with scale
factors (ab pb YI) and T2 with scale factors (a2, P2, y2). The product of Tj and T2, apply Tl to O and T2 to Ob is
defined by: T2 (T10) = T2 Oj = O2, where the scale factors of T2Tj are given by (a2 a1; P2 P1; y2 YJ). Transformations
which satisfy the above properties form a group. Thus, as stated above, Newton's three laws of motion (for speeds
much less than the speed of light) are invariant under this group of transformations. And in general, if the
hypotheses of a given physical theory (e.g., governing equations and initial and boundary conditions) are invariant
under a given group G of transformations, then the conclusions (e.g., system's solutions) are invariant under G; and
these invariances are related to the system's symmetries.
Another important group of transformations in continuum mechanics is the ten parameter Galilei-Newton Group
defined by Birkhoff (1955):
1. The space translation subgroup,
x =x2 + c2,
(B-2)
B-3
-------
2. The time translation subgroup,
t' = t + c; (B-3)
3. The three parameter subgroup of rigid motions,
xi = aikxk , (B-4)
k=l
||aik || = most general 3x3 orthogonal matrix;
4. The three parameter subgroup of moving axes translated with constant velocity,
Xj=Xi-bit. (B-5)
Newton's three laws of motion, for speeds much less than the speed of light, are also invariant under the Galilei-
Newton Group as they are under the group in Equation (B-l). Under these transformations, the definitions of the
material constants (density, viscosity, etc.) are left unchanged for constant mass in the original and primed systems.
In Birkhoff 's classical work of 1955, he concluded that because there are no general uniqueness theorems for the
combined system of Navier-Stokes Equations, the continuity equation, the equation of state, and general initial and
boundary conditions, complete inspectional analysis of hydrodynamic modeling is not possible. Such analysis
awaits fundamentally new theorems about partial differential equations. Although much advancement has taken
place since 1955, this is still the case in many areas of hydrodynamics, especially in some areas of interest to
subsurface scientists, such as hysteresis (Visintin, 1994) and phase transitions (Brokate and Sprekels, 1996).
Self-similarity is a special condition of a single system (Skoglund, 1967). A system is serf-similar if there exist a
separable variable of the governing equation and the initial and boundary conditions of the system. This separable
variable is called the similarity variable, and it permits one to reduce the number of variables/parameters which need
to be considered in the system. Such variables are very valuable in the search for symmetric solutions of special
partial differential and integral equations which satisfy special initial and boundary conditions. Classical examples of
the application of self-similarity variables are the solutions of certain diffusion systems, the Prandll boundary -layer
equations of viscous fluids, and the coagulation equation describing interacting liquid and solid particles. Whether
or not a system is serf-similar is not obvious, and the discovery of similarity variables may be a tedious process.
However, the use of group-theoretic techniques and other mathematical tools reduces some of this tedium as
evidenced by the work of Rogers and Ames (1989), Hill (1992), and Sachdev (2000). Several examples of interest to
the subsurface scientist are detailed in these three books.
B.3 Scale Dependence and Scale Invariance in Hydrology
Sposito ( 1 998a) has edited a recent survey concerned with symmetry and asymmetry in the field of hydrology.
The symmetric features identified in the field are related to the concept of scale invariance, while the asymmetric
features are related to scale dependence. If a system's features are symmetric, then there are no intrinsic scales for
these features, be they length scales, time scales, velocity scales, or what have you. The symmetric features are
similar for "any positive numerical values" of the scale factors in question. Of course, the use of the word "any" is
only mathematical; based on realistic physical considerations, there are both upper and lower limits on the values of
the pertinent scale factors. If some feature is asymmetric or becomes asymmetric for some reason, we say that this
feature's symmetry is broken and the feature becomes scale dependent. An intrinsic scale exists for the feature.
Thus, as the pertinent scales vary, the feature of the system changes in a nonsimilar manner.
The review of Sposito (1998a) is concerned with a variety of topics such as: scale analyses for land-surface
hydrology, scaling of river networks, spatial variability and scale invariance in hy drologic regionalization, scaling of
field soil-water behavior, scaling invariance of Richards Equation and its application to watershed modeling, scale
issues of heterogeneity and stochastic modeling of scale-dependent macrodispersion in vadose-zone hydrology,
dilution of nonreactive solutes in heterogeneous porous media and scale effects in large-scale solute-transport
models, and scale effects in fluid flow through fractured geologic media and transport in multiscale permeability
fields. Many of these topics are of direct concern to investigators of soil-moisture distribution and solute transport
in the unsaturated, subsurface media. In what follows, we will briefly consider some of the topics in this review and
show how they are related to the scaling features found in the HYDRUS and CHAIN-2D Codes.
B4
-------
One of the first papers on the similarity of soil-water behavior in unsaturated porous media and a paper related to
the subject at hand is that of Miller and Miller (1956). Their work is based on the principle of geometric similitude in
the microscopic arrangement of pores and solid particles, as illustrated in Figure B-l. Such similar porous media will
have the same porosities and the same volumetric water content 9. Thus, these similarity assumptions resulted in the
following scaling relationships:
Soil-Water Property
Symbol and Dimensions
Scaling Relationships
Volumetric Water Content
Soil-Water Pressure Head
Hydraulic Conductivity
em
h[L]
K[LT-!]
9* = 9
h* = ah
K* = or2K
where the asterisk denotes the referenced quantities, or the scaled properties that are the same for all the Miller-Miller
similar porous media. The quantity a denotes the scale factor for the pores and solid particles combined (as shown
in Figure B-l). The Miller-Miller similitude has been successfully applied to laboratory studies of sandy soils but has
not been as successful in the analysis of field soils. The reason for this is that the saturated soil moisture content 9S
usually varies from place to place in a field such as shown in Figure B-2. In this figure, 9S will usually vary (and thus,
the 9-distribution) as one moves from vertical profile to vertical profile. That is to say, the soil profile p; has a
different internal geometry than that of soil profile PJ, i # j. Geometric similitude is not satisfied.
Warrick, et al. (1977) accounted for the fact that total porosity of a field soil is highly variable even within a given
soil mapping unit by replacing 9 in the Miller-Miller similitude by Seo = 9/9s, where the residual soil moisture content
9r in the effective water content Se is taken to be zero. In essence, these authors introduced another scale factor into
the Miller-Miller system, namely the saturated soil moisture content 9S which accounts for the different internal
geometries in field soils. The resulting scaling relationships for this new system are as follows:
Soil-Water Property
Symbol and Dimensions
Scalins Relationships
Relative Saturation
Soil-Water Pressure Head
Hydraulic Conductivity
h[L]
K[LT-i]
c * = c
•Jeo oeo
h* = ahh
K*=ak-2K
Soil Surface
TR"
Figure B-2. The depiction of a set of vertical soil profiles ppp2,... distributed over afield mapping unit, where z
represents the local variable within a soil profile and R. (i = 1,2,...) represent the horizontal vectors
in the xy-plane giving the global positioning of the vertical profiles, p,,p2,....
B-5
-------
where the asterisk denotes the referenced quantities which may be fieldwide averages, whereas the un-asterisked
quantities denote field-plot-specific profiles of the quantities in question. The scale factors ah and ak are, in general,
not equal as in the Miller-Miller similitude case.
Warrick, et al. (1977) avoided a search for a microscopic physical length, as required by geometric similarity, by
merely deriving values of a,, that minimized the sum of squares
N M
where N represents the number of macroscopic locations within a field soil and M and number of observations of h
as a function of Seo. Figures B-3abcd illustrate this procedure for observations taken over an agriculture field
consisting of panoche soil. Figure B-3a shows 840 measurements of (9,h) plotted as h vs Seo, where the soil samples
were taken at six soil depths at 20 sites within the field, giving an N = 120. Each of these 120 soil samples were
analyzed in the laboratory with seven values of h, giving various values of Seo for each h, thus M is equal to 7. The
result of minimizing S S in Equation (B -6) was the coalescing of the data in Figure B -3 a into the single curve h* (Seo)
shown in Figure B-3b, the h*-curve being a low order rational function in Seo. Similarly, the 2640 values of (K,9) in
Figure B-3c, plotted as Seo vs K, were coalesced and described by the regression curve indicated in Figure B-3d. This
regression relationship is a cubic polynomial in Seo for the logarithm to the base e of K*. Warrick, et al. found that the
scale factors for ah were not equal to those for scaling K,ak. Further, it should be noted that the values of h(9) scaled
in Figure B-3ab were those measured in the laboratory on soil cores removed from the field, and the values of K(9)
scaled in Figures B-3cd relied on the laboratory results for h(9) to obtain estimates of 9(t) based on field tensiometric
measurements.
1
D
l.U
0.8
0.6
0.4
IT>
1
*
1
i
i i
a
"
OBSERVATIONS
» • FROM 6 SOIL
. • * t DEpTHS
{ill
0 -200 -400
Pressure Head, h in cm
1.0
0.8
0.6
0.4
0.2
SCALED
OBSERVATIONS
I
0 -200 -400
Reference Head, h* in cm
0.2 0.4 0.6 0.8 1.0
Saturalion,
0.2 0.4 0.6 0.8 1.0
Saturation, Sso
Figure B-3. (a) Unsealed observations of S eo(h), (b) Scaled observations of Seo(h*), showing the reference
relationship as a solid curve, (c) unsealed observations of K(Seo), and (d) scaled observations of
K*(Seo) (reprinted from Water Resources Research, 1977, by A.W. Warrick, G.J. Mullen andD.R.
Nielsen with the permission of the American Geophysical Union, Washington, DC).
B-6
-------
Based on the early work of Miller and Miller (1956) and Warrick, et al. (1977), as well as the principles of
similitude and invariance, Vogel, et al. (1991) introduced and studied the system of transformations that form the
foundation of the scaling relationships in the HYDRUS and CHAIN 2D Codes. These linear transformations
account for the linear variability in the soil moisture parameters and are defined by the following:
h=ahh*, (B-7)
e(h)-er=ae[e*(h*)-e;], (&*>
K(h) = akK*(h*), (B-9)
where once again the asterisk denotes reference quantities and the three scale factors (ah,ae,ak) are assumed to be
independent of one another.
Vogel, et al. (1991) claim that these linear transformations will explain the variability due to soil structure within a
given soil textural class, but will not account for the nonlinear phenomena expressed by different soil textural classes.
For example, if the soil-water behavior is represented by one of the models considered in Appendix A, say the VG-
Model, then the a fitting parameter in the model relates to the linear scaling factors and is dominated by the soil
structure. The shape factors in the quantities 9(h) and K(h) of the VG-Model are expressed by the n fitting parameter
which is highly correlated with soil texture and represents phenomena not accounted for in Equations (B-7) to (B-9).
Thus, the scaling laws in Equations (B-7) to (B-9), with their scale factors (ah ,(Xe ,(Xk ), reference hydraulic
characteristics (h*, 9*, K*), and the parameters of the reference quantities derived from analytical expressions (e.g.,
the BC-, VG-, or VC-Models), identify a set of similar soil classes where each member of the set corresponds to a
given textural class. Hence, this linear transformation allows an investigator to resolve a set of soil measurements
(from one field mapping unit, or a set of field mapping units) into a set of similar soil classes where the similar class
distinction is given by the various soil textural classes found in the data. Within a given similarity class, the soils are
linearly related by the transformations in Equations (B-7) to (B-9), and this relationship is dominated by varying soil
structure within the given textural class. These variations in soil structure may be due to certain linearly changing,
heterogeneous layers in a vertical profile and/or due to global changes in soil structure as one changes horizontal
location within field mapping units (Figure B-2). Figure B-4 is a schematic which summarizes the similarity
classification process.
Other features of the Vogel, etal. (1991) similarity system include the resolution of the scaling factors
(ah ,ae ,ak ) into products of two subfactors. One subfactor accounts for the local variation within a given profile
(the z-variation of profiles p l, p2,... in Figure B -2) and the other subfactor accounts for the global variation between
profiles (the R1; R2,... locations inFigure B-2). However, the reference quantities (h*, 9*, K*) are constructed to be
independent of global variations, while being locally dependent on some type of mean-field profile. Further, the
authors consider the invariance of Richards Equation (Appendix A and the main text) for a set of soil profiles with
respect to the transformations in Equations (B-7) to (B-9), supplemented by linear transformations in the variable, z,
time, t, infiltration rate, v, cumulative infiltration, I, and cumulative outflow, Q. Richards Equation was shown to be
invariant under these transformations, thus giving invariant solutions to the system, provided certain restrictions
were met. These restrictions were placed on the local and global subfactors of the scaling factors, and on the initial
and boundary conditions of the system. The results of the derived invariances led to simple relationships between
the soil-moisture scaling factors (ah ,ae ,ak ) and the dynamic scaling factors (% az, o^, a:, aQ). Because of the
restrictions placed on the system's auxiliary conditions, these invariances were shown to hold only for the following
situations:
1. Constant head infiltration into a stratified soil profile.
2. Water redistribution within a soil-water system with zero flux boundaries.
3. Drainage flow with constant suction head boundary conditions, such as in laboratory outflow experiments.
Finally, Vogel, et al. (1991) outline the procedures for identifying the (cth ,(Xe ,(Xk ) scaling factors and reference
quantities (h*, 9*, K*) from measured soil hydraulic data, as well as the indirect identification procedures of these
quantities from measured dynamic characteristics of soil-water systems, such as those obtained from field infiltration
tests and laboratory outflow studies. The authors conclude that these linear transformation procedures and their
resultant dynamical similarities produce the following advantages in the analyses of soil-moisture processes:
B-7
-------
Nonlinear Variability in Soils
Linear Variability in Soils
Shape Parameters Correlate with Soil Texture
Scaling Parameters Dominated by
Soil Structure
Figure B-4. The application of linear scaling to a set of soil-moisture observations, resulting in a set of m similar
soil classes. The m* class of similar soils accounts for the soil structures sml, sm2,..., smq within
the soil textural class m. The number of similar soil classes corresponds to the number of soil
textural classes.
1. There exists savings in the experimental and computing efforts required in the analyses of heterogeneous
soils;
2. There exists simplifications in the analyses of certain phenomena, such as the
derivation of 9 and K distributions from outflow experiments and infiltration tests; and
3. There exists efficient and viable methods of relating laboratory scale results to field scale results.
Rockhold, et al. (1996) applied the basic principles found in the previous three studies to simulate the water flow
and the transport of tritium measured at an initially unsampled domain of the Las Graces Trench Site. The hydraulic
parameters in this simulation were obtained in the following manner. From water-retention data for core samples
collected at the site (but not from the domain in question), the scale factors ah in h* = ah h were determined from a
scale-mean BC-Model (Appendix A) and the measured moisture data. Parameters for the soil-water retention were
then used in the Burdine (1953) relative-permeability model to estimate K(9). Saturated values forthe hydraulic
conductivity Ks were measured in the field at nearly 600 locations. The Ks data were scaled according to
B-8
-------
K(0)=a^K*(0). The probability distribution for ah and ak was found to be lognormal. The geostatistical horizontal
variograms for the log-transformed scale factors (ah, ak) show similar spatial structures to lags of 4 to 6 meters, which
is an indication of a Miller-Miller similitude structure. Thus, Rockhold, et al. (1996) found that it was not necessary to
invoke three independent factors (ah,ae,ak)to condition the hydraulic properties of the field to run their simulations.
Rather they used a constant value for 9S, constant values for the slopes of Seo(h) and K(Seo), and the same
distribution for ah and ak to condition the hydraulic properties of the field in their simulations. Simulations of water
flow adequately agreed with those measured without any further "calibration" of the model. Hence, the authors
showed the power of scaling analyses in the simulation of unsaturated flow and transport processes.
In the Sposito (1998a) review, which was previously mentioned, there are six sections that are of particular
interest to the subject at hand. These are Sections 5 to 10. In Section 5, Nielsen, et al. (1998) gives a detailed survey
of the scaling of field soil-water behavior, from the work of Miller and Miller, to the use of fractal and multifractal
concepts for representing soil structures. In addition, they review the work which combines linear-variability scaling
with the inverse solution of the Richards Equation to estimate in situ hydraulic properties of the soil. Although many
of the techniques surveyed have been successfully applied in various problems, the authors conclude that there
exists no generalized theory at this time for the comprehensive scaling of the behavior of field soil-water regimes.
In Section 6, Sposito (1998b) considers the relationship between scaling invariance and Richards Equation (the
Richards Equation which is free of sources and sinks). Using group analysis (Lie Groups in particular), Sposito
shows that the invariance of the pore size distribution under scaling leads to the result that similar media will show
mathematical identical forms of Se (ah) but not for 9 (ah). This invariance is compatible with a broad range of particle
sizes and the nonuniform porosity of field soils. Further, the scale invariance is a more general concept than fractal
serf-similarity. The Lie Group analysis of the Richards Equation, using the transformation
(e',t',z')=(ne,5t,az), (B-10)
where (|J,,8,a) are independent scale factors, requires that the parameters D(9), K(9), and h(9) mustbe/www or
exponential functions if the Richards Equation is scale invariant. For the derivation of invariant solutions of the
transformed Richards Equation, compatible (invariant) initial and boundary conditions are also required. The scale
invariance of the Richards Equation implies the scale invariance of Se, implies that the BC-Function can be a model for
h(9), and also implies that h(9) can be represented by a generic fractal model Perrier, et al. (1996). In addition, Sposito
(1998b) considers the effects of broken symmetries, for example the effects of time being unsealed, and the effects of
length being unsealed.
In Section 7, Haverkamp, et al. (1998) consider the scaling of Richards Equation and its application to watershed
modeling. They make a strong plea for the use of dimensional analysis and inspectional analysis (i.e., dynamic
similarity) as opposed to just the scaling of a soil's static hydraulic characteristics. This latter procedure fails to scale
vadose-zone behavior in a general way. With this in mind, this section is concerned with cumulative-infiltration
curves I(t) analyzed for different surface boundary conditions and for subsurface water movement governed by the
Richards Equation. They found that for general field soils, there is no unique dynamic similarity for the behavior of
soil-water movement in the vadose zone (see Figure B -4). Only for two soils (Green and Ampt, 1911; and Gardner,
1958) does unique dynamic similarity exist. However, these two soils define the limits of the envelope for all possible
similarity classes that exist for general field soils. For the description of general infiltration, three scaling parameters
are required, while only two are required forthe two limiting cases. Further, because of these two limiting cases, in
practical applications such as watershed modeling, one can assume that there exists a unique dynamic similarity, thus
simplifying certain hydrologic analyses.
In Section 8, Yeh (1998) investigates the key scale issues of heterogeneity in vadose-zone hydrology. He
defines a field scale representative elementary volume (FSREV) and considers the effects on specifying the hydraulic
parameters when our domain of interest is greater than FSREV and when it is less. He considers the use of stochastic
methods to give better estimates of these parameters, but dogmatically states that only large amounts of data can
lessen the uncertainties in vadose-zone parameters and can make the stochastic results statistically meaningful.
In Section 9, Russo (1998) considers the problem of solute transport through partially saturated heterogeneous
porous formations, while inSection 10, Kapoor and Kitanidis (1998) are concerned with the dilution of nonreactive
solutes in heterogeneous porous media. The key results of these papers show that solute plumes may not spread as
much in the longitudinal direction as expected in homogeneous cases, that the travel distances required for the
longitudinal component of the macrodispersion tensor to approach its asymptotic value (i.e., Fickian behavior) can be
B-9
-------
exceedingly large, and that current standard approaches of estimating plume concentrations can severely
overestimate dilution in heterogeneous media.
References
Birkhoff, G. 1955. Hydrodynamics: A Study in Logic. Fact and Similitude. Dover Publications, Inc. New York, NY.
186pp.
Brokate, M. and J. Sprekels. 1996. Hysteresis and Phase Transitions. Springer-Verlag, Inc. New York, NY, 357 pp.
Burdine, N.T 1953. "Relative permeability calculations from size distribution data." Am. Inst. Min. Metal. Pet. Eng.
198,71-77.
Flavin, J.N. and S. Rionero. 1996. Qualitative Estimates for Partial Differential Equations: An Introduction. CRC
Press. BocaRaton,FL. 368pp.
Gardner, W.R. 1958. "Some steady-state solutions of the unsaturated moisture flow equation with application to
evaporation from a water table." SoilSci. 85. 228-232.
Green, W.H. and G. A. Ampt. 1911. "Studies in soil physics: I. The flow of air and water through soils." J. Agric. Sci.
4,1-24.
Guymon,G.L. 1994. Unsaturated Zone Hydrology. PTR Prentice Hall. Englewood Cliffs, NJ. 210pp.
Haverkamp, R., J.-Y Parlange, R. Cuenca, P. J. Ross and T S. Steenhuis. 1998. "Scaling of the Richards equation and
its application to watershed modeling." In Scale Dependence and Scale Invariance in Hydrology. Edited by G.
Sposito. Cambridge University Press. New York, NY. 190-223.
Hill, J.M. 1992. Differential Equations and Group Methods for Scientists and Engineers. CRC Press. BocaRaton, FL
201pp.
Kapoor, V and P. Kitanidis. 1998. "Dilution of nonreactive solutes in heterogeneous porous media." In Scale
Dependence and Scale Invariance in Hydrology. Edited by G. Sposito. Cambridge University Press. New York,
NY, 291-313.
Miller, E.E. and R.D. Miller. 1956. "Physical theory for capillary flow phenomena." J. Appl. Phys. 27. 324-332.
Nielson, D.R., J.W. Hopmans andK. Reichardt. 1998. " An emerging technology for scaling field soil-water behavior."
In Scale Dependence and Scale Invariance in Hydrology. Edited by G. Sposito. Cambridge University Press.
New York, NY. 136-166.
Perrier, E., M. Rieu, G. Sposito and G. de Marsily. 1996. "Models of the water retention curve for soils witha fractal
pore size distribution." Water Resour. Res. 32(10). 3025-3031.
Rockhold, M.L., R.E. Rossi and R.G. Hills. 1996. "Application of similar media scaling and conditional simulation for
modeling water flow and tritium transport at the Las Cruses Trench Site." Water Resour. Res. 32(3). 595-609.
Rogers, C. and W.F. Ames. 1989. Nonlinear Boundary Value Problems in Science and Engineering. Academic Press,
Inc. New York, NY. 416pp.
Rosen, J. 1995. Symmetry in Science: An Introduction to the General Theory. Springer-Verlag. New York, NY. 213
pp.
Russo, D. 1998. "Stochastic modeling of scale-dependent macrodispersion in the vadose zone." In Scale
Dependence and Scale Invariance in Hydrology. Edited by G. Sposito. Cambridge University Press. New York,
NY 266-290.
Sachdev, PL. 2000. Serf-Similarity and Beyond: Exact Solutions of Nonlinear Problems. Chapman &Hall/CRC. Boca
Raton, FL. 319pp.
B-10
-------
Skoglund, V.J. 1967. Similitude: Theory and Applications. International Textbook Company. Scranton, PA 320pp.
Sposito, G. 1998a. Editor of Scale Dependence and Scale Invariance in Hydrology. Cambridge University Press. New
York, NY 423 pp.
Sposito, G. 1998b. "Scaling invariance and the Richards equation." In Scale Dependence and Scale Invariance in
Hydrology. Edited by G. Sposito. Cambridge University Press. New York, NY. 167-189.
Visintin,A. 1994. Differential Models of Hysteresis. Springer-Verlag. New York, NY. 407pp.
Vogel, T.,M. Cislerova and J. W. Hopmans. 1991. "Porous media with linearly variable hydraulic properties." Water
Resour.Res. 27(10). 2735-2741.
Warwick, A.W, G.W. Mullen and D.R. Nielsen. 1977. "Scaling field-measured soil hydraulic properties using a similar
media concept." WaterResour. Res. 13. 355-362.
Yeh, T. -C.J. 1998. "Scale issues of heterogeneity in vadose-zone hydrology: In Scale Dependence and Scale
Invariance in Hydrology. Edited by G. Sposito. Cambridge University Press. New York, NY. 224-265.
B-ll
-------
Appendix C
An Explanation of the Hysteretic
Characteristics of Soil-Water Properties
The HYDRUS Code is the only code of the five considered in this report to introduce the phenomenological
concept of hysteresis in the formulation of the soil-water properties, 9(h) and K(h). This phenomenon is certainly not
unique to the field of soil physics, for it has a long history in many engineering and scientific systems such as
mechanical, electrical, chemical, biological, physical, and geophysical systems, and even in psychological systems.
The presence of hysteretic components in a system tends to make the system highly nonlinear, may produce
instabilities and discontinuities, and certainly introduces memory dependent processes. Thus, hysteresis will make
easy problems complex, and complex problems even more complicated. But hysteresis is a natural phenomenon that
does exist and must be dealt with in many areas of soil physics, as well as in other scientific areas.
There are two sets of inquiries which one is concerned with when dealing with such a phenomenon. One set
consists of the following questions: In what physical systems does it occur? Under what conditions is the
phenomenon important? And, how does one formulate the phenomenological component in mathematical terms and
incorporate it into the governing equations of the system in question? The other set of inquiries deals with the well-
posedness of the formulated systems containing the phenomenon. The questions in this set of inquiries include:
Does a solution exist for the formulated system? Is the solution obtained unique for the given governing equation
and initial and boundary conditions? And, is the solution continuous with respect to the system parameters, the
form of the phenomenological component, and the initial and boundary conditions?
With reference to the first set of inquiries for the phenomenon of hysteresis, one runs across such concepts as
hysteresis loops, hysteresis operators, hysteretic models, partial differential equations with memory, and
discontinuous hysteresis, as well as ideas from catastrophe theory. Because of the complexity introduced by
hysteresis operators and the memory components, it is necessary to analyze the resultant systems with respect to the
concepts of well-posedness, our second set of inquiries. The tools used to answer this set of inquiries are highly
mathematical and involve such fields as functional analysis, semigroups, variational inequalities, and differential
inclusions. In the following subsections we will only briefly cover some of the key points in the first set of inquiries
and will not discuss any of the mathematical tools used in the second set of inquiries. The subsections that follow
are designated as:
(1) The Origins and Application of Hysteretic Phenomena,
(2) Hysteresis Loops, Operators and Models,
(3) Hysteresis in Soil-Moisture Parameters.
Although this appendix does not cover the mathematical aspects of hysteretic phenomena, we list several
pertinent references for the mathematically-inclined reader:
(1) For the origins and applications of hysteresis, the formulation of hysteresis operators and models, and the
investigations of the well-posedness of the systems see Visintin (1994), and Brokate and Sprekels (1996).
(2) For the theory and application of functional analysis see Rudin (1991) and Edwards (1995).
(3) For the formulation and applications of catastrophe theory see Gilmore (1981), and Castrigiano and Hayes
(1993).
(4) For the application of variational inequalities in porous media and other mechanical systems see Chipot
(1984), andHlavaceketal. (1988).
(5) For the theory and applications of semigroups in linear and nonlinear systems see Goldstein (1985) and
Miyadera(1992).
(6) For the definition, use, and solution of differential inclusions see Aubin and Cellina (1984) and Visintin
(1994).
C-l
-------
C. 1 The Origins and Applications of Hysteretic Phenomena
Hysteresis is a phenomenon that occurs in rather different situations. It can be a byproduct of fundamental
physical mechanisms, such as phase transitions; or a consequence of a degradation or imperfection, like the play in a
mechanical system; or a deliberate construct of a system in order to monitor the system's behavior, as in the case of
heat control via thermostats (Brokate and Sprekels, 1996). In physics, we encounter hysteresis in plasticity, friction,
ferromagnetism, ferroelectricity, superconductivity, adsorption and desorption, and other processes. Recently, shape
memory effects have been observed and exploited in some new materials. Hysteresis also occurs in the engineering
disciplines, such as in porous media filtration, granular motion, semiconductors, spin glasses, and mechanical
damage and fatigue. In addition, the phenomenon appears in chemical, economical, psychological, and biological
processes (Visintin, 1994).
As indicated above, hysteresis effects are often caused by phase transitions which are accompanied by abrupt
changes of some of the involved physical quantities, as well as the absorption or release of energy in the form of
latent heat. The area of the hysteresis loop (Figure C-l) gives a measure of the amount of energy that has been
dissipated or absorbed during the phase transformation. For example, phase transitions possess hysteresis effects
when there is undercooling prior to nucleation and superheating prior to vaporization. Other phase transitions
possessing hysteresis effects are as follows (Visintin, 1994):
Phase Transition
ferromagnetic
ferroelectric
liquid-vapor
martinsite
phase separation
C.2 Hysteresis Loops, Operators and Models
Order Parameter
magnetization
polarization
reduced density
strain
concentration
External Field
magnetic field
electric field
pressure
stress
chemical potential differences
Figure C-1 shows a typical continuous, closed hysteresis loop for a system whose state is defined by the two
scalar variables (u, v). The variable u represents the input to the system (i.e., the independent variable) and v
represents the output (i.e., the dependent variable). Both u and v depend on time t, and evolve as time increases in a
CONTINUOUS HYSTERESIS LOOP
Figure C-l. A continuous hysteresis loop for a system whose state is given by the couple (u,v), where u
is the input and v is the output (reprinted from Differential Models of Hysteresis, 1994, by Augusto
Visintin, with permission of Springer-Verlag, GmbH & Co. KG, Heidelberg, Germany).
C-2
-------
manner dictated by the system in question. When u increases from Uj to u2 , then the state of the system (u, v)
evolves along the lower curve of the loop, ABC. If u increases beyond u2, then the state (u, v) follows the curve from
C to E. If u now decreases from its value at point E, then the state (u, v) moves along the curve from E to C. As u
decreases from u2 to ub the state (u, v) follows the upper curve of the loop CD A. If u decreases beyond ub then (u,
v) follows the curve from A to F. Now, if u begins to increase from the u at point F, the cycle begins to repeat itself.
If u is between Uj and u2 and the evolution of the state (u, v) is inverted, (u, v) will move into the interior of the
region bounded by ABCDA, (i.e., region R). In standard examples, the couple (u, v) can attain any interior point of R
by a suitable choice of the input function u. These last two properties follow from the memory aspect of hysteresis.
As early as (1905), Madelung attempted to formalize these memory attributes into a set of rules governing the
experimentally observed branchings and loopings of ferromagnetic hysteresis. The rules which were developed are
as follows (see Figure C-2):
(1) As the state of the system (u, v) moves down curve F from point C, it arrives at turning point A, where the
state (u, v) now moves up curve Fj. Rule 1 states that Tl is uniquely determined by the coordinates of A.
(2) Point B is now taken as any point on curve Tl where the motion of the state (u, v) is inverted, point B
becomes a new turning point. Rule 2 says that the evolution of the state is now along curve F2 and that
curve starts at B and terminates at A.
(3) If the input value u decreases beyond its value for point A, Rule 3 says the state (u, v) follows along curve
F from A to D.
These rules of Madelung can be applied at any point on the loop ABCDA in Figure C-l, and at any point in its
interior R if previous state evolutions have reached that interior point.
We assume that the evolution of v is uniquely determined by that of u, and such a result is made precise by
formulating the concept of a hysteresis operator, W: u —> v. However, we know from above, that whenever
MADELUNG'S RULES
v
0
Figure C-2. A defining sketch of Madelung's Rules for the memory attributes of ferromagnetic
hysteresis (reprinted from Hysteresis and Phase Transitions, 1996 by Martin Brokate, Jtirgen
Sprekels, with permission of Springer-Verlag, GmbH & Co. KG, Heidelberg, Germany).
C-3
-------
U!
-------
V '
b
(uc,vc;
-a
Figure C-3a. A relay with hysteresis, or a delayed relay, defined by the parameters (a, b, u c, vc) with repect to the
system defined by states (u,v) (reprinted from Differential Models of Hysteresis, 1994, by Augusto
Visintin, with permission of Springer-Verlag, GmbH & Co. KG, Heidelberg, Germany).
Figure C-3b. An approximation to a continuous hysteresis loop by a linear combination of a finite family of
delayed relays. The quantity Rf is the region inside the discontinuous loop formed by the finite
family of relays (reprinted from Different! al Models of Hysteresis, 1994, by Augusto Visintin, with
permission of Springer-Verlag, GmbH & Co. KG, Heidelberg, Germany).
C-5
-------
complemented with pertinent initial conditions. Optimal control problems for ordinary differential equations with
hysteresis require specific attention since, even though hysteresis operators are in some sense continuous, they are
usually not differentiable. Understanding such systems cannot be accomplished by merely setting up a "numerical
solution" scheme, programming it, and generating computer output. Mathematical analyses are required for these
complex problems.
Constitutive laws in continuum mechanics formulated in terms of hysteresis operators lead in a natural way to
the partial differential equations of the conservation of mass, momentum and interval energy being coupled with
hysteresis operators. Mathematically speaking, the most interesting and complicated situations are those where the
hysteresis operators appear in the principal part of the partial differential equations. For these systems, basic
existence and uniqueness results are strongly linked in non-obvious ways to the properties of the hysteresis
diagrams and of the memory structure. For such situations, Visintin (1994) considers two model problems: the heat
equation (comparable to the concentration equation) with hysteresis,
du , d r -, 32u 32u
— + — P(u) - — - — = f(x,z,t) , (Ola)
dt dt dx dz
and the wave equation with hysteresis,
a2u a
dt2 dx
(C-4b)
where P(«) is the Preisach Model of hysteresis or some other hysteresis model. Although not studied in detail,
Visintin (1994) introduced two systems of prime interest to the subsurface investigator. One of these systems is the
partial differential equation model describing the growth of bacteria in the presence of nutrients. This phenomenon
exhibits a pattern of growing rings in response to nutrient diffusion. The constitutive relation between the
concentrations of nutrients and the activity of the bacteria leads to a bistability region in the "nutrient space", which
in turn leads to the occurrence of hysteresis in the evolution of the concentration of the bacteria. The other system
considered by Visintin is of direct interest to the problem at hand. This problem is the study of unsaturated water
flow through porous media whose constitutive laws, 9(h) and K(h), account for compressibility, and the history of
drainage/wetting, thus leading to hysteretic effects. In a simplified setting, the "Richards Equation" is written as:
f- [e + P(6)]-|^-f [P(6)]=-S, (C-5)
dt dz dz
where P(9) is the hysteretic constitutive equation. Suitable initial and boundary conditions must be specified for
Equation (C-5), and S and the starting point of P(9) must be given.
C.3 Hysteresis in Soil-Moisture Parameters
In Volume 2 of the Annual Review of Fluid Mechanics, Philip (1970) says that the usual state of the unsaturated
zone is that wetting and drying are simultaneously occurring in different regions of the zone, and that this leads to
the occurrence and the importance of the phenomenon of hysteresis. This phenomenon is most pronounced for
media and moisture ranges where capillarity is dominant. Philip says that the study of this problem started in the
late 1800 's with the investigation of interfacial stability. However, it was not until the late 1920 's that these ideas were
systematically applied to subsurface phenomenon. During this period, it was recognized that many of the possible
equilibrium interface configurations in capillary-porous media are unstable and that the configuration changes are
often spontaneous and uncontrollable. Such changes involve irreversibility and produce hysteretic effects in the
soil-moisture parameters. In addition, Philip says that there is another source of the discontinuities seen in stable
interface configurations, that source is due to the noncontinuous lines of solid-liquid-gas contact which exist due to
the geometrical impossibility of contact lines on some parts of the solid surface. The significance of this latter effect
was illustrated by Philip by making interface configuration calculations in a tube with a radius varying periodically
with the axial coordinate. Following up on the interfacial stability discussions of Philip are discussions given by
Lucknerand Schestakow (1991) in their book on subsurface processes. These discussions are concerned with
C-6
-------
interfacial tension, contact angles between gaseous-liquid interfaces and surface of the solid phase (i.e., a measure of
wettability), and the differential between contact angles at the advancing front of the wetting liquid phase and that at
the receding front (i.e., there exists a hysteresis of contact angles).
In some respects the above comments of Philip can be thought of as the more theoretical aspects of hysteresis in
soil-moisture parameters, while the more phenomenological aspects of the problem are framed in the language of soil
water potential and pressure head. This latter approach is more in line with the hysteresis components considered in
the HYDRUS Code. Water potential, symbolized by \|/, is the free energy of water, or its capacity to do work (Miller
and Donahue, 1995; Mauseth, 1998). When \|/ > 0, the water is capable of doing positive work, while \|/ < 0 means that
work must be done on the water to remove it from its given location and move it to a reference pool (e.g., moving
adsorbed water on soil particles to a free pool of water). The free energy of water can be increased by heating it,
putting it under pressure, and by elevating it. Conversely, the energy of water is reduced by cooling it, reducing its
pressure, and lowering it with respect to gravity. Water's capacity to do work can also be decreased by its adherence
to a solid substance due to hydrogen bonding between the water molecules and the material.
Soil water potential is a combination of the following effects (Miller and Donahue, 1995):
1. The surface area of soil particles and small soil pores that adsorb water - matrix potential \|/m;
2. The effects of dissolved substances - solute or osmotic potential \|/s;
3. The atmospheric or gas pressure effects - pressure potential \|fp;
4. The position of the water with respect to a reference state, a state usually taken as zero - gravitational
potential \|f
In summary, the water potential for soils independent of a reference state is given by
Vw=Vm+Vs+Vp, (06)
and the total water potential is defined by
Vt=Vw+Vg- (C-7)
Water moves whenever there is a difference of water potential within a specific mass of water; this is true in
subsurface soils, within plants, and between soils and living cells. However, if the water potentials of two regions are
equal, the regions are in equilibrium and there is no net movement of water. Temperature differences are usually not a
factor in the subsurface because the solutions being compared are assumed to be at nearly the same temperatures.
Even the temperature differences between leaves and roots in a plant are usually not of significance for water
movement.
With reference to the components in Equations (C-6) and (C-7), the quantity \|/s is zero for pure water and is
always negative for solutions. The quantity \|/p can either be positive or negative and the matrix potential \|/m is a
measure of water's adhesion to nondissolved structures such as soil particles, cell walls and membranes. Since
adhesion can only decrease water's free energy, \|/m is always negative. In soils \|/m is very important because so
much soil water is tightly bound to soil particles; while in living cells, \|/m is usually much less important than \|/s and
\|/p. Miller and Donahue (1995) state that most productive soils do not have a depth of water standing on them (i.e.,
|\|/p is small) and most soils have few salts (i.e., |\|/s is small), thus
Vt=Vw=Vm- (08)
So for most situations, \|/m is the dominant component of \|/t and represents about 95% or more of \|/t. For living cells,
the cell potential \|/c is fairly well approximated by
Vc=ys+yP^o. (c-9)
The description of hysteresis given by Guymon (1994) starts with how water is distributed in a pore and how \|/m
varies through the cross section of that pore. An example of such a pore is given in Figure C-4. In this figure are
shown areas which are drained by the pull of gravity, areas where water is held by capillary forces, and areas where
C-7
-------
Figure C-4. A cross section of a soil pore and the solid soil particles that make up its walls, showing areas
drained by the pull of gravity, areas where water is held by capillary forces, and areas where water
is held by surface forces (e.g., van der Waals forces) (reprinted from Soils in Our Environment,
1995,7th Ed., by R. W. Miller and R.L. Donahue, with permission of Pearson Education, Inc., Upper
Saddle River, NJ).
water is held by surface forces. Point a in the figure represents a point where air begins to enter the pore; for fine
grained soils, the air entry pore water pressure head is about -0.1 atm. Point b is located in the area of water being
held under the influence of capillary forces and may be represented by a water pressure head on the order of -3 to -5
atm. The concept of field capacity is defined as the pore water pressure at which gravity drainage is negligible and
for agriculture soils it occurs at about -2 to -3 atm. Capillary forces dominate surface forces over the range of
pressures from the air entry pressure to -1 atm and beyond in the negative direction down to a pressure as low as
possibly -5 atm. The permanent wilting point is the point at which plants can no longer extract water from the soil;
this may occur at from -5 to -15 atm and is indicated by Point c in Figure C-4. Point d, within the area where the water
film is held by surface forces, may represent a location where the pressure head is as low as -60 atm. Within a few
water molecules of the surface of the soil particles, Point e in Figure C-4, the water pressure may be as low as -8000
atm. With reference to hysteresis, Guymon (1994) states that the soil-water hysteretic memory of soils is
predominantly a capillary effect (as was also stated by Philip, 1970). In the ranges of soil water pressures where
chemical and other surface forces predominate, hysteresis is not evident. In addition, soil memory also ceases when
and where the pressure head is zero or greater. Thus, Guymon says that hysteretic memory effects commonly occur
in the pore pressure range from -5 atm to 0. Of course, the lower limit is dependent upon the soil structure and
texture.
The main cause of hysteresis advanced by Guymon is what is called the "ink bottle" effect as illustrated in
Figure C-5. This figure (Figure C-5a) illustrates that during a drying cycle starting from saturated conditions (h,, = 0),
water is drained from the soil. However, because of the presence of narrow necks, the capillary forces at a given
negative pressure head, hj < h0 = 0, will tend to keep some of the lower large pores filled or partially filled with water
as shown in Figure C-5a. Prior to the initiation of a wetting cycle, it is assumed in Figure C-5b that the soil has
reached the residual moisture content with a soil pressure head, h2 « hj < ho = 0. When the wetting cycle begins,
the soil moisture increases and h2 —> hj. Since the larger pores are starting from the "empty position" with respect to
capillary forces, the water retained in them is due to the lower narrow necks as h2 —> \ and is not due to the upper
narrow necks. Therefore, at the given soil-water pressure head hb the water content during drying is higher than that
during wetting. A secondary cause of hysteresis identified by Guymon is that due to wettability, or the so-called
"rain drop" effect. This effect is described as follows: a wetted soil while drying retains more sorbed water than a dry
soil while being wetted. Thus, the tendency for a higher water content during drying than during wetting is
enhanced.
C-8
-------
THE INK BOTTLE EFFECT
Air, Vapors
and
Gases
(a) A Draining Pore
(b) A Wetting Pore
Figure C-5. The "ink bottle" effect demonstrating that draining/drying soils under the influence of capillary
forces retain more water at a given soil-water pressure than a wetting soil at the same water
pressure (reprinted from Unsaturated Zone Hydrology, 1994, by Gary L. Guymon, with permission
of Pearson Education, Inc., Upper Saddle River, NJ).
Guymon (1994), as well as Luckner and Schestakow (1991), claims that significant effects of hysteresis have been
observed in the soil retention curves, 9(h), while only minor hysteresis effects have been observed in the
unsaturated hydraulic conductivity function, K(h). These authors believe that for most practical problems the
hysteretic effects experienced by K(h) may be ignored, especially since there is such a large error associated with
determining K(h). In spite of these conclusions, the HYDRUS Code does consider hysteresis in both 9(h) and K(h).
A hypothetical soil-moisture hysteresis loop with a discontinuity as the pressure head approaches zero is
illustrated in Figure C-6. The discontinuity is displayed by the positive difference between the saturated moisture
content for the main drying curve and that for the main wetting curve:
e? - er > o.
(C-10)
The reason for this difference is entrapped air in some of the soil pores, a condition that may be removed with time as
the moisture redistributes. Mathematically speaking, the varying discontinuity violates the time invariance of the
hysteresis loop and presents a time-varying condition that is most difficult to verify and document in the field, and to
simulate by modeling.
Also shown in Figure C-6 are sets of primary and secondary scanning curves for wetting and drying cycles.
These scanning curves follow the rules of Madelung depicted in Figure C-2. Knowing the time history of the wetting
and drying cycles throughout a soil column allows an investigator to know at each time and space location where he/
she is at with respect to the moisture state (9,h) on the main drying and wetting curves or the moisture state (9,h) on
scanning curves within the main loop's interior. However, this is much easier to say than to precisely formulate (as
indicated by the previous subsection), and much easier to say than to verify and document in the field. As indicated
by the HYDRUS Code, one can set up numerical procedures which seem to parrot the soil-moisture memory effects in
a soil column; but one should be cautious of the simulated outputs with respect to their verification, or lack thereof,
and their field documentation.
The main drying curve shown in Figure C-6 may be experimentally determined, but is more often expressed in
terms of one of the analytical forms found in Appendix A. The parameters occurring in the analytical model are
evaluated from field or laboratory data. Given the main drying curve, the main wetting curve is obtained by changing
(according to experimental results) one or more of the parameters in the main drying curve. The primary and
secondary scanning curves are usually obtained by a scaling procedure applied to the main drying and wetting
C-9
-------
„
O ^
CD jjj
I (1)
-1- o:
Q) §
3 O
<" £?
e/j ®
0) z
Q_ |
CD
is
^^ °?
=0=1
CO
-ATM
Primary Scanning Curves
Secondary Scanning Curves
o e
Volumetric Moisture Content
Figure C-6. A hypothetical soil-moisture hysteresis loop which is discontinuous for pressure heads near zero,
showing the main drying and main wetting curves, and example primary and secondary scanning
curves (reprinted from Subsurface Transport and Fate Processes, 1993, by R.C. Knox, D.A.
Sabatini, andLW. Canter, with permission of CRC Press, Inc., Boca Raton, FL).
curves (Lucknerand Schestakow, 1991). In the case of the HYDRUS Code, the analytical VG- or VC-Model of
Appendix A is used to give the main drying curve with experimentally derived parameters. Allowance is made in the
code for 9d being equal to 0™, or 9d being greater than 0™. Usually, the following relationships for the main
curves are assumed:
6rd = 67 = 6r, nd = nw = n, ad < aw. (C-ll)
Thus, the main wetting curve is primarily differentiated from the main drying curve by the choice of ad < aw, and
possibly by the saturated moisture contents if 0™ < 9d. If 0™ = 9d, the loop of the main curves closes as the
pressure head approaches zero.
In the HYDRUS Code, the drying scanning curves are scaled from the main drying curve, while those for wetting
are scaled from the main wetting curve. The scale factor for the pressure head in this process is unity, ah = 1. As
previously mentioned, this code does introduce hysteretic effects into the relationship K versus h. As with soil
moisture, allowance is made for differences between dry and wet saturated hydraulic conductivities:
Ksw < Ksd.
The construction of the hysteresis loop for K(h) is similar to that of the moisture hysteresis loop.
(C-12)
C-10
-------
References
Aubin, J.-P. and A. Cellina. 1984. Differential Inclusions: Set-Valued Maps and Viability Theory. Springer-Verlag.
New York, NY. 342pp.
Brokate, M. and J. Sprekels. 1996. Hysteresis and Phase Transitions. Springer-Verlag. New York, NY. 357pp.
Castrigiano, D.P.L. and S.A. Hayes. 1993. Catastrophe Theory. Addison-Wesley Publ. Co., Reading, MA. 250pp.
Chipot, M. 1984. Variational Inequalities and Flow in Porous Media. Springer-Verlag. New York, NY. 118pp.
Edwards, R.E. 1995. Functional Analysis: Theory and Applications. Dover Publ., Inc., New York, NY. 783pp.
Gilmore, R. 1981. Catastrophe Theory for Scientists and Engineers. John Wiley & Sons. New York, NY. 666pp.
Goldstein, J. A. 1985. Semigroups of Linear Operators and Applications. Oxford University Press. New York, NY. 241
pp.
Guymon,G.L. 1994. Unsaturated Zone Hydrology. PTR Prentice Hall. Englewood Cliffs, NJ. 210pp.
Hlavacek, I., J. Haslinger, J. Necas and J. Lovisek. 1988. Solution of Variational Inequalities in Mechanics. Springer-
Verlag. New York, NY. 275pp.
Knox, R.C.,D. A. Sabatini andL.W. Canter. 1993. Subsurface Transport and Fate Processes. Lewis Publ., Boca
Raton, FL. 430pp.
Luckner, L. and WM. Schestakow. 1991. Migration Processes in the Soil and Groundwater Zone. Lewis Publ.
Chelsea, MI. 485 pp.
Madelung, E. 1905. " Uber Magnetisierung durch schnell verlaufende Strome and die Wirkungsweise des
Rutherford-Marconi-schenMagnetdetektors." Ann. Phys. 17. 861-890
Mauseth, J.D. 1998. Botany: An Introduction to Plant Biology. 2/e. Multimedia Enhanced Edition. Jones and
BartlettPubl. Boston, MA. 837pp.
Miller, R. Wand R.L.Donahue. 1995. Soils in Our Environment. 7th Edition. Prentice Hall, Englewood Cliffs, NJ. 649
pp.
Miyadera, I. 1992. Nonlinear Semigroups. Translations of Mathematical Monographs. Vol. 109. Amer. Math Soc.
Providence, RI. 230pp.
Philip, J.R. 1970. "Flow in porous media." In: Annual Review of Fluid Mechanics. Vol.2. Annual Reviews, Inc. Palo
Alto,CA. 177-204.
Rubin. W 1991. Functional Analysis. 2ndEdition. McGraw-Hill, Inc. New York, NY. 424pp.
Visintin,A. 1994. Differential Models of Hysteresis. Springer-Verlag. New York, NY. 407pp.
C-ll
-------
Appendix D
The First-Order Decay Chains Used in the Various Models
The one-dimensional version of the transport equations for solutes involved in sequential
first-order decay reactions consists of the following:
i
dt dz dz
6R2 ^ + q^ = eD2 + F2(c2)c2 +G2(Cl,s,s2), (D-2)
dt dz dz |_ oz
3s*
— J-=H(ci,s* )-(m +n[)s*+Yi, i=U, (D-3)
dt
where 9 is the soil moisture content, "1" is the index for the parent species, "2" is the index for the daughter species,
R is the retardation accounting for gaseous and liquid solute phases and the equilibrium solid phase, c is the
concentration of the liquid phase, q is the Darcian fluid flux density, D is the dispersion coefficient, s* represents the
nonequilibrium solid phase of the solute, /j, is the degradation rate of s*, \i' is the rate of transfer of material from
parent to daughter, and y is the production of s* . The function F(c) represents degradation of the solute from
gaseous, liquid and solid phases, transfer losses from parent passed on to daughter for the three phases, losses due
to plant uptake, and changes due to the temperature variations in the system. The function Gl ( s, ) accounts for
production of the solute in all three phases, losses due to plant uptake, and the amount of solid phase at the
nonequilibrium sorption sites. The function G2 (q ,s* ,s*2 ) accounts for the daughter material transferred from the
parent in the gaseous, liquid, equilibrium solid and nonequilibrium solid phases, accounts for the production of the
daughter species in all phases, and accounts for root uptake losses. Finally, the function H(ci5s*) could be thought
of as the first-order sorption rate of s* due to the difference between equilibrium concentrations and nonequilibrium
concentrations.
The system given in Equations (D-l) to (D-3) is basically that which is presented in the HYDRUS Model. The
CHAIN 2D version of these equations replaces the vertical coordinate (z) by the vertical slab coordinates (x, z), the
quantity q by components in the x and z directions, and the dispersion coefficient D by a 2 x 2 matrix of dispersion
coefficients. The FECTUZ Code considers the one-dimensional version of the above system, but without
nonequilibrium sorption, zero-order production teams, root uptake and accounting for temperature variations. The
MULTIMED-DP 1 .0 Code is similar to FECTUZ but it only accounts for linear equilibrium sorption; it can also change
rate constants to account for temperature. Finally, the CHAIN Code is similar to that of MULTIMED-DP 1.0 but no
account is made of temperature variations and the first-order decay of the parent becomes the zero-order production
of the daughter.
In the following paragraphs, we list the type of reaction chains that are considered in the various codes. Specific
species and products are not identified in this discussion, only the types of chains that can be formulated in the
system in question. Both HYDRUS and CHAIN 2D will allow the user to specify six solutes, either coupled in an
unidirectional chain; or totally uncoupled, where each species is independent of the other. Examples that have been
run have had the following forms:
D-l
-------
Radionuclides
>•
Nitrogen Compounds
Gas
. .: .. .. .
Products
Uninterrupted Chain —
One Reaction Path for Pesticides
Gas
'! 'I 'I '
Product Product Product
Products
Interrupted Chain-
Two Independent Reaction Paths for Pesticides
Gas
Gas
Product
Product Product
Product
The FECTUZ Code can handle up to seven different chemical species, either in straight, unidirectional chains, or
in branched chains. Systems for which solutions have been obtained include the following:
Straight Chains
N<7
Simple Splitting Chains
D = Number of Daughters < 6
r|ij = Splitting Factors
D-2
-------
Simple Converging Chains
N<7
Seven-Member Branching Chain
•H12+ 1113 = T124+ T125 =
= 1
The MULTIMED-DP 1.0 Code only considers parents, daughters and granddaughters; and within each
subgroup, the code only considers 3, 4, and 4 variables (species), respectively. Pathways that have been considered
in hydrolysis transformations include those given below:
Hydrolysis Pathways
2
<
The CHAIN Code, as reported in van Genuchten (1985), considers four species in a straight chain. The
radioactive decay chain example problem given in this paper consists of the following radionuclides:
238
234
230
226
Pu-
Th-
Ra
The author states that the most critical species is the last one in the chain, 226Ra. This species is a high biological
hazard and is highly mobile, the retardation factor R4 is rather low.
D-3
-------
References
van Genuchten, M.Th. 1985. "Convective-dispersive transport of solutes involved in sequential first-order decay
reactions." Computers and Geosciences. 11(2): 129-147.
D4
-------
Appendix E
The Impact of Using a Nonuniform Moisture Distribution
versus a Uniform Distribution
Constant values of q, 9, and h were assumed for all models in the sensitivity analyses for the input parameters
Kd, q, 9, p, D, DL and Dw. On the other hand, a nonuniform water content distribution was employed by MULTIMED-
DP 1.0, FECTUZ, CHAIN 2D and HYDRUS for the analyses involving Ks, 9r, 9S, a, and (3. One should assume that
the nonuniform distribution is the "standard" and the uniform distribution is the "approximation." However, when
the van Genuchten Model for the water retention parameters (Ks 9S, 9r, a, (3) was considered in the sensitivity
analyses reported in Sections 6 and 7, an error in the originally distributed MULTIMED-DP 1.0 Code was detected.
The error arose from the solution of the governing equations for the pressure head, h, and the moisture content, 9.
When corrected, all systems were defined for a constant recharge rate, q, and boundary conditions resulting in the
following conditions for 9: 9 = base value = 0.16, at the surface; and at 6 m, 9 = saturated value = 0.32. Figure E-l
shows the results of solving Darcy's Law for MULTIMED-DP 1.0 and FECTUZ, and solving the Richards Equation
for CHAIN 2D and HYDRUS.
The water contents obtained from using the originally distributed MULTIMED-DP 1.0 Code are inconsistent with
respect to those obtained from the other three codes. The reason for the inconsistency was found to be the incorrect
use of the residual water content in the van Genuchten Module. When this error was corrected, the new water
content results became consistent with the results of the FECTUZ Code, as one would expect, and also became
consistent with the steady state (remember q and the boundary conditions are fixed) solutions for HYDRUS and
CHAIN 2D. For the simulations conducted in this report (i.e., time durations up to 10,000 days, and times to peak
concentrations of from 3000 to 8000 days), the approach to steady state moisture distributions for both HYDRUS and
CHAIN 2D is very fast. For the purposes of this report, the variable moisture content obtained from the van
Genuchten Model and Richards Equation is a steady state distribution as shown in Figure E-l. Further, one would
expect that all four models should give the same distribution (as the figure indicates) if the MULTIMED-DP 1.0 Code
is properly given. The revised code for MULTIMED-DP 1.0 is the one that was used in all of the sensitivity analyses
and any other subsequent calculations using this code.
Before comparing the breakthrough curves (BTCs) of "Tc for a uniform moisture content (9 = 0.16) with those
for a nonuniform distribution (Figure E-l), one should note the degree of nonlinearity in the distribution shown in
Figure E-l. Because of the homogeneity of the 6m layer, one notes that the steady state distribution is nearly
uniform(9 ~ 0.16) for the first 4m of depth. At about 4m, the distribution begins to break away from this uniform
value. However, even at a depth of 5.5m, the distribution has only deviated about 33% of the total difference (0.32 -
0.16). Thus, about 67% of the total deviation of the distribution takes place in the last 0.5m of the layer.
Consequently, one would assume that changes from going from a uniform water content to a nonuniform water
content (of the type in Figure E-l) would not be too great.
This comparison between uniform and nonuniform distributions of water content is illustrated in Figure E-2
using the breakthrough curves (BTCs) of "Tc predicted by the four models: MULTIMED-DP 1.0, FECTUZ, CHAIN
2D and HYDRUS. The bottom graph in this figure gives the BTCs predicted by all five models, including the CHAIN
Model, for the base values of the parameters given in Section 6. In the top graph, all input parameters are kept at
their base values except for 9 which is replaced by the distribution given in Figure E-1. No CHAIN results are given
in the top graph since 9 can only be a constant in this model.
In comparing the BTCs for uniform water content versus the nonuniform values, one sees that the nonuniform 9
reduces the Cpeak and increases the Tpeak and TMCL for each model. This is reasonable if we remember that the
nonuniform 9-distribution increases the overall moisture in the soil column versus that for 9 = 0.16 throughout the
column. Increasing moisture reduces the advection term in the transport equation and increases the diffusion term.
E-l
-------
u
1 0
I 2'°
£ 3.0
0.
- ,
0.10 0.20 0.30 0.40 0.50
Water Content (cnf/cnf)
Figure E-l. Water content distributions predicted by the HYDRUS, CHAIN 2D, FECTUZ, and MULTIMED-DP
1.0 Models. Note that the water contents (•) obtained from the originally distributed MULTIMED-
DP 1.0 Code are in error. The corrected code gives a consistent water content distribution (*) with
the other three models.
This results in lower Cpeak values for the nonuniform BTCs, produces delays in the Tpeak values and TMCL values
with respect to those for the uniform case, and produces an increase in diffusive spread of the BTCs for the
nonuniform case over that of the uniform case, although this spread is only about 1.08 times that for the uniform
case.
To get a better quantitative effect of the impact of using the uniform approximation for the sensitivity analyses of
the first seven input parameters (Kd, q, 9, p, D, DL ,DW), we used the CHAIN 2D Model to derive the BTCs and
sensitivities for the parameters Kd and DL for both the uniform and nonuniform water contents. These comparative
BTCs are given in Figure E-3 (Kd comparison) and Figure E-4 (DL comparison). As in Figure E-2, the Cpeak values for
the nonuniform water content BTCs have been uniformly reduced over the ranges of the input parameters compared
to those for a constant 9, while the T k values and TMCL values have been increased uniformly over the input
ranges. Table E-l lists the output vames for the uniform and nonuniform water content cases, values which are
measured at the base values of the input parameters Kd and DL. In addition, the table gives the relative sensitivities
of these three outputs to the inputs Kd and DL, again referenced to the base values of Kd and DL. As one can see,
the relative sensitivities for the uniform and nonuniform cases are basically the same. The percent differences of the
output values, based on the nonuniform case as the standard, are as follows:
Percent Differences
Output
*~peak
T
1peak
TMCL
15. 5% Decrease
6.7% Increase
6.4% Increase
12.1% Decrease
5. 9% Increase
6.1% Increase
Based upon these results and the BTCs in Figures E-2 to E-4, we believe that the use of the uniform approximation for
9 for the analysis of the sensitivity of the three outputs to the input parameters Kd, q, 9, p, D, DL and Dw is both valid
and representative of the sensitivities that would be obtained from the nonuniform water content distribution.
Further, using this uniform approximation allows us to directly compare the results of the CHAIN Model with the
other four models, as illustrated in the bottom graphs shown in Figure E-2. However, it should be emphasized that
these results are for a homogeneous layer under steady state conditions.
E-2
-------
0.015
0.010
0.005
o>
E
c
o
c
o
0.015
o
o
o
r
0.010
0.005
Nonuniform Water Content Distribution
HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP 1.0
Uniform Water Content Distribution (Base Case)
A CHAIN HYDRUS
CHAIN 2D
FECTUZ
MULTIMED-DP 1.0
2500 5000 7500
Time (days)
10000
Figure E-2. Comparison of the breakthrough curves predicted by the CHAIN, HYDRUS, CHAIN 2D, FECTUZ,
and MULTIIMED-DP 1.0 Models for the base case given in Section 6. The top curves are for a
nonuniform water content and the bottom curves are for 9 = 0.16 throughout the soil column.
There are no CHAIN results in the top graph because 9 can only be constant in this model.
Table E-l Comparison of Results Derived from Figures E-3 and E-4 for the Distribution Coefficient Kd and the
Dispersivity DL, respectively. The Values of Cpeak, Tpeak and TMCL are Given for the Base Values of
Kd and DL, and the Relative Sensitivities of These Output Quantities to Kd and DL are Given, These
Values Also Being Taken at the Base Values of Kd and DL.
Output
c-
Tpeak
TMCL
Property
Base Value (mg/1)
Relative Sensitivity
Base Value (d)
RelativeSensitivity
Base Value (d)
Relative Sensitivity
i
Uniform
Water Content
0.00732
-0.06
4766
+0.06
3554
+0.07
din ml/g
Nonuniform
Water Content
0.00634
-0.05
5108
+0.06
3799
+0.07
DL in cm
Uniform
WaterContent
0.00725
-0.28
4737
-0.02
3552
-0.08
Nonuniform
WaterContent
0.00647
-0.25
5036
-0.02
3783
-0.07
E-3
-------
0.0101
"5; 0.008
o 0.006
§ 0.004
I
o 0.002
Distribution Coefficient (ml/g)
0.001 CHAIN 2D
0.007
0.019
(nonunifotm water content)
0 2500 5000 7500 10000
Time (days)
0.010
0.008
o 0.006
0.004
0.002
Distribution Coefficient (ml/g) CHAIN 2D
0.001
0.007
0.019
(Uniform water content)
2500 5000 7500
Time (days)
10000
Figure E-3. Sensitivity of "Tc breakthrough (through the 6 m layer) to the distribution coefficient using the
CHAIN 2D Model, for a nonuniform water content (top) and for a uniform water content (bottom).
0.015
0.010
o.oos
CHAIN 2D (nonuniform water content)
Dlsperalvlty (am)
3.73
4.53
5.33
4000 6000 8000 10000
Time (days)
0.015
0.010
o
o
t—
s
CHAIN 2D DtspersMty (cm)
(uniform water content) 3.73
5.33
2000 4000 6000 8000 10000
Time (days)
Figure E-4. Sensitivity of "Tc breakthrough (through the 6 m layer) to the dispersivity using the CHAIN 2D
Model, for a nonuniform water content (top) and for a uniform water content (bottom).
E4
-------
Appendix F
The Impacts of Using Daily Precipitation Rates and Daily
Evapotranspiration Rates versus an Annual Average Recharge Rate
In the sensitivity analyses of Sections 6 and 7, the impact of the "natural cycles" of precipitation/runoff/
infiltration/evapotranspiration was assumed to be a constant mean annual recharge rate, q, which is applied year-
after-year for a total simulation time of 25 to 40 years. It is the constant rate, q, which carries the radionuclides down
through the homogeneous 6 m layer to the hypothetical water table at its lower boundary. The result of this ideal
hydrologic mechanism is the classical "bell-shaped" breakthrough curve (ETC) shown throughout this report. These
"bell-shaped" BTCs are probably sufficient to demonstrate the sensitivities of the three output quantities to the
twelve input parameters given in Section 6. However, this "ideal mechanism" is at most a very crude approximation
of the real mechanisms which produce ground-water contamination from radionuclides passing through the
unsaturated zone.
The attributes of radionuclide movement through the vadose zone depend on the details of the time and space
variations of the "cycles" of precipitation/runoff/infiltration/ evapotranspiration. Closely linked to these cycles are
the mechanisms of flow hysteresis and flow through preferential pathways. In fact, Yao and Hendrickx (1996) state
that worldwide it is thought that a major mechanism for ground-water contamination is the passage of pollutants
through preferential flow paths, which is an event type phenomenon not currently well parameterized for large time-
step simulations. These preferential flow paths can be produced by macropores such as soil cracks, old root channels
or animal paths; by spatial variability of hydraulic conductivity; or by unstable wetting fronts. The occurrence of
unstable wetting fronts is especially important since this mechanism can produce preferential pathways in
homogeneous layers which do not contain macropores and significant changes in hydraulic conductivity, such as the
6 m layer used in the current report. In several of the mechanisms which produce wetting front instabilities,./?*™;
hysteresis (Appendix C) plays a major role, see Section 3.4.
An example of the impact of flow hysteresis on the results of single storm events is given in Eagleson (1970).
He reported the results of a numerical experiment on the redistribution of moisture in a soil profile after the cessation
of infiltration. In this simulation it was found that percolation with hysteresis was much slower than percolation
without hysteresis. This simulation occurred over a two hour period and the percolation was through sandy soil.
Such a result applied to a soil with vegetation would allow a longer time for evapotranspiration to act under
hysteretic flow than under nonhysteretic flow, thus producing a greater loss of subsurface moisture under flow
hysteresis.
Modeling codes required to simulate the above phenomena which lead to flow fingering and preferential
pathways for moisture and pollutant transport through the vadose zone to the ground water should be two-
dimensional in character and preferably three-dimensional, have the ability to handle event precipitation/runoff/
infiltration/evapotranspiration, have the ability to simulate air compressibility and flow hysteresis, and the ability to
vary soil properties in a significant manner. In addition, such simulations would require time steps on the order of an
hour or less, and space increments fine enough to simulate a two-dimensional pattern of flow fingers in the horizontal
field of interest. Such simulations are beyond the objectives of the current report and beyond the capabilities of the
five models under consideration, although the one-dimensional HYDRUS Code considers most of the pertinent
hydrologic processes required for such simulations. Thus, even though the above preferential pathway phenomena
are potentially very important for the transport and fate of radionuclides through the vadose zone, they will not and
cannot be considered further in this appendix. Further, the time and space increments required for an analyses of such
phenomena would be prohibitive for the current applications, where our simulations run from 10,000 days to 15,000
days.
F-l
-------
To compromise, the simulations considered in this appendix are based on daily precipitation amounts without
allowance for runoff, and on daily evapotranspiration rates based on two different root water uptake stress response
distributions (Equations 2-3 and 2-6). Considering the 18 years of daily precipitation amounts shown in Figure 5-2,
one notes that there is only one peak above 80 mm/d, two peaks above 40 mm/d, and about 30 peaks above 20 mm/d.
These three categories, respectively, convert to 0.33 cm/hr, 0.17 cm/hr, and 0.08 cm/hr. However, one should keep in
mind that these hourly rates are lower, and probably much lower, than storm event rates. Further the frequency of
occurrence of the higher rates is greater when one uses storm events as opposed to using the daily precipitation rates.
In the work of Yao and Hendrickx (1996), for New Mexico type soils and hydrologic conditions, it was found
that wetting front instabilities and fingering occurred for infiltration rates between 0.3 and 12 cm/hr, incomplete
wetting without distinct development of fingering occurred for infiltration rates between 0.12 and 0.3 cm/hr, and
stable wetting front migration occurred for infiltration rates below 0.12 cm/hr. Assuming these laboratory results are
convertible to the homogeneous 6 m layer used in the current simulations, one would expect wetting front instability
to occur on many days throughout the total record of 6570 days. Thus, the complexities of wetting front instability
for storm events may be inconsistent with the assessment of the impacts of daily precipitation/evapotranspiration on
the transport of radionuclides in the unsaturated zone. However, such complexities, as stated above, require finer
time and space scale simulations than those used in the current study.
The precipitation record used in the current study, as mentioned above, is the one shown in Figure 5-2. Figure F-
1 shows this same record in terms of annual precipitation amounts and monthly average amounts for the 18 years of
record at the Las Cruces Site. As one can see from Figure F-la, nine of the first ten years of record are above the
annual average of 22.5 cm; the last eight years of the record are equal to or less than the mean annual average. For
the monthly averages, each of the first five months has a monthly value less than 1 cm. The monthly averages begin
to rise sharply in June; and from July to October, the monthly averages are greater than 2 cm per month. November's
average drops below 2 cm, while December returns to just over 2 cm.
As given in Equations (2-3) and (2-6), the sink term S in the modified Richards Equation for the soil moisture
distribution in the soil profile is expressed as the product of a water stress response function a(h) times a potential
water uptake rate Sp(z, t). In turn, the quantity Sp(z, t) is a product of the normalized water uptake distribution
function b(z, t) and the potential transpiration rate Tp(t). The quantity b (z, t) is related to the root distribution
function and in the current study was taken as uniform over an interval from 50 cm below the surface to 250 cm
below the surface. This type of root distribution is probably reasonable for mesquite tree and creosote bush
vegetative cover over semi-dry regions, such as the Las Cruces Site (see Royo, 2000; and Miller and Donahue,
1995). The potential transpiration rate was calculated from the site's climate data using Penman's Equation (Jensen
et al, 1990). The water stress response function a(h) used in the current simulation was the Feddes Module of the
HYDRUS Code (Simunek, et al, 1998). This function is the simple trapezoid shown in Figure F-2. The quantity hj
is the value of the pressure head below which roots start to extract water from the soil (a value of -20 to -50 cm); h2
is the pressure head below which roots start to extract water at the maximum possible rate, a(h2) =1; h3 is the
pressure head below which roots can no longer extract water at the maximum rate; and finally h4 is the pressure head
below which root water uptake ceases, this being usually equal to the wilting point. Water uptake by roots can be
easily varied by changing one or more of the values of the triplet (h1; h2 h3). In the current simulations, two sets of
this triplet were used to give two different scenarios of root water uptake. Two modifications that were not
considered in these simulations were root growth and decline, and the characteristics that desert plants use to reduce
water uptake and evapotranspiration rates, such as leaf coatings and root shrinkage which produces root-soil air gaps
that are highly resistant to moisture transfer (Nobel, 1994). Such modifications are probably important in the first
five months of the year (Figure F-lb) and during the dry years (Figure F-la). Thus, the results that follow may be
giving higher actual evapotranspiration (ET) rates than occur in nature, and lower recharge rates.
Figure F-3 summarizes the results of the two root water uptake scenarios run on the daily precipitation and PET
data for the Las Cruces Site. These two sets of curves give the cumulative amounts of precipitation, actual ET and
net recharge (precipitation minus actual ET, no allowance for runoff) in centimeters. The HYDRUS Model was used
to calculate the actual ET values, and thus the recharge rates. The chief characteristics of these two sets of curves are
given in the following notes:
The records for both sets of curves consist of three pieces, 0 to a, a to b, and b to c, covering a total of
14,140 days.
F-2
-------
(A
**
O
0)
p
1
o
aBejeAv A|muo|/\|
I
S2
re
0)
p
|Bnuuv
60
2
g aa
a a
aE
CH
' '^
" ^
it
.&
'o
.8
H
60
'H -3
S 8
F-3
-------
0)
w
^ 1.0--
.
at
81
0) *•
" e
**
a
0.5--
Q
£
LU
Figure F-2.
Soil Water Pressure Head, h
The water stress response function for the Feddes Module of the HYDRUS Code (Simunek, et al..
1998).
1000
800
600
400
- 200
^
0
1000
800
600
400
200
Cumul. Net Recharge (Avg. 67 mm/yr)
Cumul. Actual ET
Cumul. Precipitation
b
Cumul. Net Recharge (Avg. 54 mm/yr)
Cumul. Actual ET
Cumul. Precipitation
b
3000
6000 9000
Time (days)
12000 15000
Figure F-3. Cumulative amounts in centimeters of precipitation, actual evapotranspiration (ET), and net
recharge (precipitation minus actual ET) during a HYDRUS Model simulation using daily variable
precipitation and potential ET rates at the surface. Cumulative net recharge and ET vary between
the two figures because of differences in the root water uptake scenarios, (hj, h2, h3). The
precipitation/PET segment from "a to b" is repeated from "b to c."
F-4
-------
The records for the interval 0 to a cover 1000 days and represent a constant recharge rate and precipitation
rate (no ET taking place) of 0.024 cm/hr. This period was used to let the radionuclides distribute
throughout the "upper layers" of the 6 mvadose zone under the action of the base values given in Section 6.
Once the source of radionuclides quit emitting material at 1000 days, the daily precipitation and PET
records in Figure 5-2 were applied to the system for the transport and fate of the radionuclides contained in
the "upper layers" of the 6 m zone. Thus, the second piece of the records in Figure F-3 runs from a = lOOOd
to b = 7570d.
When the simulation reached 7570 days, the daily precipitation and PET records in Figure F-2 were
repeated, giving the third piece of the records in Figure F-3, which runs from b = 7570d to c = 14,140d.
The average recharge rates, 67 mm/yr for the top set of curves in Figure F-3 and 54 mm/yr for the bottom
set, were obtained from the graphs by taking the cumulative recharge values at 12,000 days and converting
these values to mm/yr.
Considering the wettest part of the 18 year record, the first ten years, the mean annual recharge rate for the
top set of curves in Figure F-3 over the first 4650 days (10 years plus 1000 days) is about 84 mm/yr and that
for the bottom set is 59 mm/yr.
The high frequency "oscillations" in all the cumulative curves shown in Figure F-3 are due to daily and
weekly variations in the precipitation and PET records.
Ignoring these oscillations in the cumulative precipitation curves (which, of course, are the same in both sets
of curves), one sees a curve made up of five line segments:
1. 0 to lOOOd, 87 mm/yr recharge rate for initial spreading of pollutants,
2. lOOOd to 5015d, first eleven years of record, the wettest period,
3. 5015d to 7570d, last seven years of record, the driest period,
4. 7570d to ll,585d, repeat of the wettest period,
5. ll,585d to 14,140d, repeat of driest period.
Ignoring the daily and weekly "oscillations" in the bottom set of curves in Figure F-3, one sees that the root
water uptake scenario is such that the divergence between the cumulative precipitation and cumulative
actual ET in each of the four segments from lOOOd to 14,140d is roughly constant from segment to segment.
This result produces a net recharge cumulative curve which is roughly a straight line, which means that after
the first 1000 days the pollutants are being transported by a constant recharge rate of about 52 mm/yr.
Once again ignoring the daily and weekly "oscillations" in the top set of curves in Figure F-3, one sees that
the root water uptake scenario is such that the divergence between the cumulative precipitation and
cumulative actual ET in the segment from lOOOd to 5015d is greater than that for the bottom set of curves,
giving a greater recharge rate of about 83 mm/yr. In the next segment, 5015d to 7570d, the cumulative
precipitation curve is approximately parallel to the cumulative actual ET, giving a recharge rate which is
relatively small at 20 mm/yr. The last two segments, 7570d to ll,585d and ll,585d to 14,140d, are a repeat
of the previous two segments.
The impacts that the variable recharge rates in Figure F-3 have on the transport and fate of radionuclides through
the 6 m of the vadose zone are shown by the "Tc breakthrough curves (BTCs) given in Figure F-4. The solid BTCs
in this figure corresponds to the daily computed recharge rates, while the dashed BTCs correspond to calculations
based on constant recharge rates (67 mm/yr and 54 mm/yr) for the entire time duration of the simulations. For both
sets of BTCs, the moisture distribution throughout the soil column is nonuniform (i.e., the modified Richards
Equation was solved). Observations that one can make about these results are as follows:
The solid curve in the bottom set of curves is basically the classical "bell-shaped" curve as one might expect
by the nearly constant recharge rate of 52 mm/yr after the first 1000 days of a rate of 87 mm/yr. The dashed
curve is the ETC for a constant recharge rate of 54 mm/yr from zero to the end of the simulation.
The solid curve in the bottom set of curves varies from that of the dashed curve because of the first 1000
days of pollutant distribution under a recharge rate of 87 mm/yr and because of daily variations of
F-5
-------
o
f
0.012
0.010
0.008
0.006
0.004
0.002
0
0.012
0.010
0.008
0.006
0.004
0.002
Variable Precip/Actural ET HYDRUS
Uniform Net Recharge
Average Net Recharge: 67 mm/yr)
Variable Precip/Actual ET HYDRUS
Uniform Net Recharge
Average Net Recharge: 54 mm/yr)
5000 10000
Time (days)
15000
Figure F-4. Comparison of predicted "Tc breakthrough curves (through the 6 m layer) using the variable
precipitation/actual ET versus uniform recharge rate in the HYDRUS Model. Average recharge
rate is calculated as the mean net amount of precipitation and actual ET from 0 to 12,000 days.
The net recharge varies between the two sets of curves due to the root-uptake scenario, (hbh2 ,h3).
precipitation and PET. The Cpeak of the solid curve is 1.16 times that of the dashed curve, and the Tpeak of
the solid curve is 1.01 times that of the dashed curve. At the scale shown in Figure F-4, the daily and
weekly "oscillations" are not apparent, but their cumulative effects are present.
The relationship between the dashed curve in the bottom set, which reflects a nonuniform moisture
distribution, and the corresponding uniform moisture content result is such that the Cpeak decreases by 14%
and the Tpeak increases 5.5% over those for the uniform case, based on the standard being the nonuniform
case. These percentages are consistent with those found for the parameters Kd and DL in Appendix E.
The solid curve in the top set of curves in Figure F-4 does not possess the classic "bell-shape". This is
primarily due to the nonlinear cumulative recharge shown in the top set of curves in Figure F-3. The
average recharge rate is 87 mm/yr for the first 1000 days, 83 mm/yr for the next 4015 days, followed by 20
mm/yr for the next 2555 days, then 83 mm/yr for 4015 days, and finally 20 mm/yr for the last 2550 days.
The first part of the solid curve in the top set is somewhat like that of the base case for HYDRUS given in
Section 6, except the variable precipitation Cpeak is 1.08 times higher and its T k 0.88 times smaller. After
5015 days, the average recharge rate drops sharply to 20 mm/yr, thus giving a long "heavy tail" to the
distribution. At 7570 days the larger recharge rate kicks in again, producing a cutoff of the tail at about
12,500 days, or so.
The dashed curve in the top set of Figure F-4 is the nonuniform moisture distribution ETC for a constant
recharge rate of 67 mm/yr. Comparing this curve to results for the uniform moisture case, one sees that the
Cpeak decreases by 13.6% and the Tpeak increases by 2.9% over those for the uniform case, based on the
standard being the nonuniform case. These results are consistent with those of Appendix E.
F-6
-------
References
Eagleson, P.S. 1970. Dynamic Hydrology. McGraw-Hill Book Co. New York, NY. 462pp.
Jensen, M.E.,R.D. Burman, and R.G. Allen. 1990. "Evapotranspiration and irrigation water requirements." ASCE
Manuals and Reports on Engineering Practice. No. 70. American Society of Civil Engineers, NY. pp.332.
Miller, R.W. and R.L. Donahue. 1995. Soils in Our Environment. 7th Edition. Prentice Hall, Englewood Cliffs, NJ.
649 pp.
Nobel, P.S. 1994. "Root-soil responses to water pulses in dry environments." In: Exploitation of Environmental
Heterogeneity by Plants: Ecophysiological Processes Above- and Belowground. M.M. Caldwell and R.W.
Pearcy, Eds. Academic Press. New York, NY. pp. 285-304.
Royo, A..R. 2000. "Desert Plant Survival." DesertUSA Newsletter. http://www.desertusa.com/du_plantsurv.html.
5pp.
Simunek, I, K. Huang, andM.Th. van Genuchten. 1998. The HYDRUS Code for Simulating the One-Dimensional
Movement of Water. Heat, and Multiple Solutes in Variably-Saturated Media. Research Report No. 144. U.S.
Salinity Laboratory. USD A, ARS. Riverside, CA.
Yao, T-M. and J.M.H. Hendrickx. 1996. "Stability of wetting fronts in dry homogeneous soils under low infiltration
rates." Soil Sci. Soc. Am. J. 60, 20-28.
F-7
-------
Appendix G
The Impact of Considering a Layered Soil Column
versus a Homogeneous Soil Column
In this appendix we are concerned with the impact on the flow and transport, and the sensitivity analyses of
replacing the 6 m homogeneous layer (the "all" layer in Table 5-2) by the 9-layer soil column given in Table 5-2. It
is recognized that these impacts will be considerably less than the impacts due to layering reported in Miyazaki et al.
(1993) and Rockhold et al. (1997). However, the objective of the current impact analysis is to assess the differences
of using the 9-layered soil column for the Las Graces Trench Site versus the single all-layer soil column, and to
decide if the use of one layer is comparable to the nine layers at this site.
As stated, Table 5-2 gives the soil hydraulic properties (KS,9S, 9r, a, P) at the Las Graces Trench Site for the first
6m of the soil profile. The 6 m homogeneous layer used in the sensitivity analyses in Sections 6 and 7 is the layer
with the hydraulic properties listed in the"all layer" row of the table. However, as seen in this table, the 600 cm layer
is broken down into 9 sublayers, having different thicknesses and hydraulic properties. The biggest outlier of all the
9 layers is, as one would expect, the 15 cm thick surface layer which contains the most organic matter and is the most
highly disturbed, both by natural and anthropogenic processes. The pathway through this surface layer represents
only about 2.5 percent of the total pathway through the 600 cm layer; thus, its effect on the total passage of "Tc
through the 6m vadose zone to the hypothetical water table is very minor. For the remaining eight layers, the range of
values for the five soil properties about their base values (Section 6) and the corresponding maximum percentage
deviation (based on the base value) are given below:
172 cm/d < Ks = 270 < 334 cm/d, 36.3 %,
0.294 <9S = 0.32 < 0.343, 8.1%,
0.071<9r = 0.08 < 0.091, 13.8%,
0.027 cm-1 < a = 0.055 < 0.070 cm-1, 50.9 %,
1.383
-------
c
o
o
c
o
o
0.010
0.008
0.006
0.004
0.002
HYDRUS Distribution Coefficient (ml/g)
Layered ° Uniform 0.001
Layered • Uniform 0.007
Layered * Uniform 0.019
2500 5000 7500
Time (days)
10000
Figure G-l. Sensitivity of "Tc breakthrough (through the unsaturated zone with a water table at a depth of 6 m)
to the distribution coefficients in a layered soil and in a uniform soil using the HYDRUS Model.
Cpeak anc* decrease T k and TMCL. In the layer from 305 to 370 cm, 9S and (3 tend to increase C k and decrease
Tpeak ancl TMCL> while 9r has the opposite effect. In the final three layers (370 to 600 cm), (3 tends to decrease Cpeak
and increase T k and TMCL, while the effects due to 9S and 9r are fairly close to those of the uniform case. From
this collection of "facts," it is difficult to dogmatically say how the layered and uniform cases should compare. Thus,
relying on the computations, one arrives at the following results, where percentages are based on the layered results
as being the standard and decreases/increases are related to changes from the uniform case:
Kd(ml/g)
0.001
0.007
0.019
(-*
^peak
9. 9% decrease
9. 5% decrease
9.0% decrease
T
peak
3. 4% increase
2.7% increase
2.7% increase
TMCL
2.4% increase
2.2% increase
2.0% increase
As suggested above, the differences between the results for the uniform soil and the layered soil are not too
great. And because of this, it seems reasonable to expect the sensitivity studies given in this report for the 6m
homogeneous layer are also applicable for the layered soil described in Table 5-2. The differences in Cpeak are
assumed to be within 10% of one another and those for Tpeak and TMCL are assumed to be around 3 to 4%, or less, of
one another. However, as initially stated, these conclusions will not be valid if the layering is of the type considered
by Miyazaki et al. (1993) and Rockhold et al. (1997), layering with significant changes in hydraulic properties from
layer to layer in the soil column.
References
Miyasaki, T., S. Hasegawa and T. Kasubuchi. 1993. Water Flow in Soils. Marcel Dekker, Inc., New York, NY,
296 pp.
Rockhold, M.L., C.S. Simmons and M.F. Payer. 1997. "An analytical solution technique for one-dimensional,
steady vertical water flow in layered soils." Water Resour. Res. 33(4):897-902.
G-2
-------
Appendix H
A Detailed Analysis of Nonequilibrium Sorption of Pollutants,
Mainly for the Radionuclide 90Sr
Thibodeaux (1996) states that the validity of the local equilibrium assumption (LEA) in solid-liquid soil
reactions depends on the degree of interaction between the macroscopic transport processes of water flow and
hydrodynamic dispersion, and the microscopic processes of molecular diffusion and sorbed-solute distribution in
conjunction with soil aggregate size. When the rate of change of solute mass during microscopic sorption processes
is fast relative to bulk flow, the interaction is nearly instantaneous, comparatively speaking, and it conforms with the
LEA. Deviations from local equilibrium occur as the interactions of the solute with the porous media become
increasingly time dependent with respect to the time scales of the bulk flow. This divergence also occurs as soil
aggregates increase in size and complexity, and the pore-class heterogeneity increases. In the present case, the
transport of radionuclides through the unsaturated zone is taking place under the influence of a very low annual
recharge rate, which leads to a hydrodynamic dispersion of only two times that of molecular diffusion. Under these
conditions, sites that were in local equilibrium at higher recharge rates will still be in local equilibrium. Some new
sites where nonequilibrium conditions existed or no liquid-solid reactions occurred at higher recharge rates will now
become under the LEA, while other new sites will pass from no action sites to nonequilibrium sites, and some will be
nonequilibrium sites under both the higher and lower flows. Thus, the adsorption-desorption cycle of a given
pollutant through a given soil profile has a variety of characteristic time scales dependent upon the type of liquid-
solid reaction mechanisms at the various sites in the soil matrix. Therefore, it seems reasonable that for "all"
recharge flows one considers nonequilibrium sorption as well as equilibrium conditions.
H.1 A Simplified Version of the Transport Equations
Although there are other models (see Section 3.6), the system used in CHAIN 2D and HYDRUS is a two-site
(equilibrium-nonequilibrium) model developed by van Genuchten and his colleagues (van Genuchten and Wagenet,
1989; Toride et al., 1993). The one-dimensional version of the transport equations using this model (i.e., the
HYDRUS Code) is given in Section 2 and Appendix D. It is this version that is used to generate results for the
current appendix. However, for purposes of discussion and argument, the following simplified set of transport
equations for the liquid phase pollutant concentration, c, and the nonequilibrium solid phase concentration, s, is
considered:
ac 0rwD +D
wwL
D/P '
at 6 + pfKd az 6 + pfKd az2 6 + pfKd
3s
— = D/P - US , (H-2)
3t
where ^ is the radioactive decay rate for both s and c, and D/P is a decay/production, sorption transfer function
defined by
D/P = 0)[(l-f)Kd C-s]. (H-3)
H-l
-------
The quantity, f, is a parameter fixed by the soil structure. When f = 1, there are no nonequilibrium sorption sites in
the soil and when f = 0, there are no equilibrium sorption sites. When 0 < f < 1, there are both equilibrium and
nonequilibrium sites. The quantity, co, is the first-order rate constant for nonequilibrium sorption in inverse time
units. Thus, we have two new input parameters for the nonequilibrium sorption process, (co,f).
When D/P = 0 (i.e., the nonequilibrium sites act as equilibrium sites), these transport equations reduce to the
following:
3c , q 3c _ OTWDW + DL q 32C
- TT - HC , (H-4)
3t 0 + pfKd 3z 0 + pfKd 3z
s = (1 - f) Kdc . (H-5)
Equations (H-4) and (H-5) define the situation when there are no nonequilibrium sites, which occurs when f = 1; or if
0 < f < 1, the equations define the situation when co is sufficiently large so as to dominate the effect of the decay rate,
H.2 The Transport of "Tc with a Kj = 0.007 ml/g
As a first example of Equations (H-l) and (H-2), let us consider the transport of "Tc through the 6 m soil
column. The initial distributions of (c,s) in the soil profile are taken as zero, and a source of c is postulated at the
surface. For "Tc, the decay rate, (J,, has a very small effect on the breakthrough curves (BTCs) since lifetimes for
these curves in the given system are from 10,000 days to 12,000 days for the default Kd of 0.0007 ml/g, while the
half-life of "Tc is 2.1 x 105 years. Thus the important decay/production in these two equations is D/P and not [i. As
the process begins, t > 0, pollutant c begins to be spread throughout the water in the soil profile. This then produces
a source of s via the D/P term; thus, s begins to be spread throughout the soil matrix in the column. As long as s is
less than (1 - f)Kdc, D/P is a sink term for c and a source term for s. If s is greater than (1 - f)Kdc, then D/P is a
source term for c and a sink term for s. Since the only other term in the s-equation is a very small sink term (i.e., |J,s),
s will approximately be equal to (1 - f)Kdc for sufficiently large co. Thus, under sufficiently high sorption rates co, we
have the "equilibrium" condition:
S - (1 - f)KdC . (H-6)
Equation (H-6) shows that for low Kd-pollutants, the amount of material adsorbed at the nonequilibrium sites for any
one time is very small, while for high Kd-pollutants, the opposite is usually true. Consequently, even though the low-
q problems considered in this report may partially violate the LEA, the resultant effect of nonequilibrium sorption
can still be very small if the pollutant in question (i.e., "Tc) has a very small Kd.
Using the HYDRUS Code with the nonuniform distribution for 9 and the base values of Section 6 for the other
10 standard input parameters, the BTCs of "Tc for different values of (co, f) were derived. One set of BTCs were
obtained for f = 0.47 and co = 0.003, 0.032, and 0.320^. Another set of BTCs were obtained for co = 0.032^ and f
= 0.27, 0.47, 0.67. These five curves were compared with the equilibrium case, f = 1 and co = 0, given by the solid
curve in the top set of curves in Figure E-2. It should be noted that when one is comparing such sets of curves,
Equations (H-l) and (H-2) are not exactly the ones solved by the HYDRUS Code since 9 is not constant but is
variable with depth as shown in Figure E-l. However, for order-of-magnitude arguments, Equations (H-l) to (H-6)
are adequate.
Considering Equations (H-l) and (H-2), one sees that the parameters (co, f) affect three processes (terms):
advection, diffusion, the decay/production term D/P. When s reaches its equilibrium value of (l-f)Kdc, the D/P term
is turned off and the s-distribution will track the c-distribution if c and s do not again get out of balance. This
tracking is such that when s = (l-f)Kdc, the c-distribution is the distribution that occurs when there are no
nonequilibrium sites. That is, the existence of active nonequilibrium sites requires an s-distribution which is not
equal to the equilibrium sorption condition given by (l-f)Kdc.
H-2
-------
For "Tc, with its default Kd of 0.007 ml/g and its large half-life, the c- and s-distributions should not get out of
balance once they are in balance. Further, the transfer of material between solution and the soil matrix is relatively
small for the (co, f) given above:
[(l-f)Kdcl = 0.51%ofc(t).
(H-7)
These results are verified when the "nonequilibrium" "Tc BTCs for the five combinations of (co, f) mentioned above
are compared with one another and with the equilibrium case for (co, f) = (0, 1). All of these BTCs are within 1% or
2% of one another at any given time. Thus, for a soil profile of 6 m, the nonequilibrium formulation given in
Equation (H-l) and (H-2) reduces to an equilibrium formulation for a radionuclide with a Kd of 0.007 ml/g and a
half-life of 76,650,000 days, if co is sufficiently large.
H.3 The Transport of "Tc with a Kd = 1.0 ml/g
For illustrative purposes, we next considered "Tc with a Kd equal to unity and (co, f) equal to (0.032d~1, 0),
(0.032d~1, 0.47) and (0, 1). The results of these numerical experiments were basically the same as the Kd of 0.007
ml/g results (i.e., nonequilibrium conditions quickly approach equilibrium conditions), even though the maximum
amount of material transferred between solution and soil matrix is much greater than that given in Equation (H-7):
[(l-f)Kdc]_ = 100%ofc(t),
(H-8)
Table H-l summarizes both sets of results (i.e.,K = 0.007, 1.0 ml/g) for three selected times and concentration
pairs (c, s) for the "Tc BTCs. The equilibrium cases are given by the pair (co, f) = (0,1), while the "nonequilibrium"
cases are given by the pairs (co, f) = (0.032, 0.47) and (0.032, 0). Columns headed by c and s are rounded off
calculated values from the HYDRUS Code, while the column headed by (l-f)Kdc is a check to see if c and s are in
equilibrium sorption balance. To round-off error, for both Kd values of 0.007 and 1.0 ml/g, all times, and a half-life
of 76,650,000 days, the equilibrium and "nonequilibrium" BTCs are the same and the c- and s-distributions are in
Table H-l. Comparison of Nonequilibrium and Equilibrium Results for "Tc BTCs for Kd Values of 0.007 and
1.0 ml/g.
Radionuclide Kd
ml/g
(d-i)
en
(d-i)
f
Time
(d)
c
mg/L
mg/kg
s
mg/kg
Technetium
"Tc 0.007 9.043xlO-9 0
0.007 9.043 xlO-9 0.032
1 9.043 xlO-9 0
1 9.043 xlO-9 0.032
1 9.043 xlO-9 0.032
0.47
0.47
3750
5005
6250
3750
5005
6250
40,000
46,500
60,000
40,000
46,500
60,000
40,000
46,500
60,000
0.00105
0.00646
0.00227
0.00111
0.00644
0.00222
3.67 x lO-4
6.22 x lO-4
1.38x10-4
3.67x10-4
6.12x10-4
1.38x10-4
3.67x10-4
6.02 x lO-4
1.38x10-4
0
0
0
4.12xlO-6
2.39xlO-5
8.24 x 10-6
0
0
0
1.94 x lO-4
3.24 x lO-4
7.31xlO-5
3.67 x lO-4
6.02 x lO-4
1.38 x lO-4
0
0
0
4.12xlO-6
2.39x10-5
8.19x10-6
0
0
0
1.90 x lO-4
3.24 x lO-4
7.35 x lO-5
3.63 x lO-4
6.01 x lO-4
1.40 x lO-4
H-3
-------
equilibrium sorption balance through the formula (l-f)Kdc = s. In addition, one sees a slight drop in the Cpeak (at t =
46,500d) for the Kd = 1.0 ml/g case from 6.22 x 10'4 mg/L for the equilibrium case (co, f) = (0, 1), to 6.12 x 10'4 mg/
L for (co, f) = (0.032.0.47) to 6.02 x 10'4 mg/L for (co, f) = (0.032, 0).
H.4 The Transport of 90Sr with a Kj = 1.0 ml/g
The third case that we considered was the base case scenario for90Sr; this radionuclide has a default Kd of 1.0
ml/g and a half-life of 10,585 days. Thus, one would expect significant decay over the lifetime of the BTCs for this
radionuclide. Figure H-la shows the BTCs for the concentrations of 90Sr in solution for the parameter pairs (co, f) =
(0, 1), (0.032d-!, 0.47), and (0.032d-!, 0). The BTCs for these three parameter pairs are basically equal, which
means that the two "nonequilibrium" cases, (co, f), (0.032d~1, 0.47), and (0.032CT4, 0), are the same as the equilibrium
case (co, f) = (0, 1). Further, Figure H-lb shows that the D/P term is essentially zero for all three parameter pairs (this
is trivially so for co = 0 and f = 1):
D/P = eo[(l -f)Kd c-s] = 0, for all time t. (H-9)
The s-distributions for f = 0 and 0.47 track those of the c-distributions through the relationship s = (1 -f)Kdc =
(1 -f)c.
The results in Table H-l and those shown in Figure H-lab support the claim that nonequilibrium sorption has
very little effect for the low recharge rates considered (0.014 < q < 0.032 cm/d), for the sorption reaction rates which
range over 0.003 < co < 0.32d~1, and for the specific radionuclides considered in this report, with Kd values that range
from 0 to 5 ml/g. This claim is further supported by the magnitude of the coefficients in Equations (H-l) and (H-2).
For example, in the 90Sr case, the coefficients for (co,f) = (0.032d~1, 0.47) are as follows:
(a)
1.5x10-*
o>
1.0x10-*
o 5. 0x10"*
O
HYDRUS MSr
Kd=1.0 ml/g
o = 0.032 f (1=0)
Equilibrium Sorption ff=1)
' (f=0.47)
o = 0.032
20000 40000
Time (days)
60000
80000
(b)
1.0x10-*
8. 0x10-*
6.0x10-*
4.0x10"*
2.0x10"*
Kd=1.0 ml/g
HYDRUS MSr
a = 0.032 dH (f=0)
Equilibrium Sorption (f=1)
a = 0.032 d"1 (f=0. 47)
16000 32000 48000 64000 80000
Time (days)
Figure H-l. The c- and s distribution of Equations (H-l) and (H-2) for 90Srfor (co,f) = (0,1),
(0.032 d'1, 0.47), (0.032d-!, 0). Distributions were derived by the HYDRUS Code, (a) gives the
concentration in solution and (b) gives the concentration on the soil matrix.
H-4
-------
Nonhomogeneous
Equation Radioactive Decay Sorption Decay Sorption Term
Eq. (H-l) -6.55xlO-5c -3.01 x 10'2c +5.67 x !Q-2s
Eq. (H-2) -6.55xlO-5s -3.20xlQ-2s +1.70 x 10'2c
From these coefficients, we note that the radioactive decay term is negligible as compared to the transfer of 90Sr
between solution and soil matrix. Thus, even though the radioactive decay greatly reduces the Cpeak values of the
90 Sr BTCs over their lifetimes of 65,000 to 80,000 days, this decay has little influence on the much faster adsorption/
desorption processes. Further, as shown in Figures H-la and H-lb, the nonequilibrium sorption process is
sufficiently fast that the D/P term goes to zero and the nonequilibrium sites, in reality, become equilibrium sites. A
more detailed analysis is given in the following paragraphs.
H.5 The Transport of 90Sr with a K,, = 1.0 ml/g, f = 0, and Varying Sorption Rates
We have found that the nonequilibrium sorption sites for 99Tc, in essence, acted as equilibrium sites for the low
recharge rates considered in this report when the sorption rates co are taken to be equal to or greater than 0.003d'1. In
addition, the nonequilibrium sorption sites for 90Sr acted as equilibrium sites for a sorption rate co of 0.032d-L.
Therefore, the question that arises is the following: At what sorption rates co will these systems truly possess
nonequilibrium sites that are not in equilibrium balance with the liquid phase concentration? We answered this
question for 90Sr, whose Kd value is taken as 1.0 ml/g and whose decay rate (J, is given as 6.55 x lO^d'1. However,
the reader should be aware that the values of co for which 90Sr possesses "true" nonequilibrium sites in this Las
Graces soil may be physically unrealistic. The current exercise was conducted to both understand this module of the
HYDRUS Code and to indicate the range of values for co that are required to unbalance equilibrium sorption
conditions for the very low recharge rates considered in this report.
For the radionuclide, 90Sr, withaKd= 1.0 ml/g and anf =0 (i.e., no equilibrium sorption sites, only
nonequilibrium sites), Equations (H-l) and (H-2) become:
— + v — = D—- - (AC - — co(c-s), (H-10)
at az az 6
3s .
— + (co + |A)S = coc, (H-ll)
at
where the coefficients (v,D) are defined by
v = q/6, D = TWDW + DL q /0 . (H-12)
Equation (H-ll) is a is a first order, linear, nonhomogeneous equation for s where s(z,0) is taken to be zero, as stated
in Section H.2. Solving Equation (H-ll) results in the following:
t
s(z,t) = co J exp[-(co + n)(t - TI)]C(Z,TI) dr| .
Combining Equations (H-10) and (H-13) produces the following linear, integrodifferential equation for the liquid
phase concentration c:
— +v— =D— +
at az az2
+
t
)[exp[-(co + jj,)(t -r|)l c(z,r|) dr| - c
(H-14)
col .-.-...- - .. -
0
H-5
-------
The integral term in Equation (H-14) radically changes the system defined by Equation (H-4) by introducing a
"history" dependence in the evolution of the liquid phase concentration. This integral term shows that the value of
c(z,t) at any point z in the soil column not only depends on the current time, but also depends on the entire history of
the evolution of c(z,r|) from the initial time r| = 0 to the current time r| = t at the point z in question. As we have seen
in Section H.2, this integral term loses its historic dependence if co is sufficiently large. But what happens if co is
sufficiently small? To answer this question, we let co equal a sequence of values from 6.5 x lO^d'1 to 6.5 x KHd'1,
remembering that [i equals 6.55 x lO^d'1.
Before considering the results obtained from the HYDRUS calculations, two comments are worth emphasizing:
(1) Equations (H-10) to (H-14) represent a "stripped-down" version of the HYDRUS Code used to
simulate the transport of 90Sr with parameters (Kd,f) set equal to (1.0, 0). The HYDRUS Code uses a
variable soil moisture, 9(z), in the simulation, while the "stripped-down" version does not. Equations
(H-10) to (H-14) are only used for explanatory and argumentative purposes.
(2) The "stripped-down" version of the HYDRUS Code (although f need not be zero) has been recently
solved by Drake et al. (2002), though the use of Laplace Transformations, the theory of residues, and
the Bromwich Integral in the complex plane. The z-domain was defined by 0 < z < L, the source at z =
0 was a variable stepwise-continuous function, and the condition at z = L was a variable flux of
pollutants. The exact solution is in the form of a single definite integral involving elementary functions
of the independent variables (z,t) and the system parameters.
Figures H-2a and H-2b depict the c- and s-distributions at the bottom of the 6 m soil column for the sequence of
co values listed in the above paragraph. Figure H-2a gives the usual BTCs for the concentration of 90Sr in solution at
the hypothetical water table (i.e., the bottom of the 6 m column). These BTCs show the usual "bell-shape" for co >
6.5 x lO-M'1 (keep in mind that Figure H-2a is a semi-log plot), but are radically different from the "bell-shape" for
co = 6.5 x lO^d'1 and 6.5 x 10-5cH. For these lowest two values of co, the "history" term in Equation (H-14) has the
greatest effect. Figure H-2b gives the corresponding concentration curves (CCs) for the solid phase at the water table
level of 6 m. As with the BTCs, the CCs show a characteristic "bell-shape" for co > 6.5 x lO^d'1, and a non "bell-
shape" for the lowest two cos. In fact, as co increases in size one can see, by careful analysis of the curves in Figures
H-2a and H-2b, that the BTCs approach the CCs. That is, as co increases, the nonequilibrium sites "convert" to
equilibrium sites, the "history" term in Equation (H-14) becomes less effective, and Equations (H-l) and (H-2)
change over to Equations (H-4) and (H-5).
Table H-2 reinforces the above conclusion. In this table we have listed the peak values of the BTCs and CCs at
the 6 m water table, C k and S k, respectively, and their corresponding times to Tpeak. The values for the BTCs
and CCs are quite different for the lowest two co values, while from co equal to 6.5 x KHd'1 and higher, the values for
the BTCs begin to quickly approach those for the CCs. For co > 0.032d'1, the BTCs and CCs tend to give the same
results. This means that the nonequilibrium sites become equilibrium sites (or act like equilibrium sites) for the low
recharge rates used in this analysis whenever the sorption rate co > 0.032d~1 and the pollutant is 90Sr. For other
pollutants with different Kds and decay rates, \i, this limiting sorption rate will vary, but the general trend is valid.
The final reinforcement of the above conclusion was made by considering the time evolution of the c- and s-
distributions throughout the entire 6 m soil column. Figures H-3 to H-6 give the liquid phase concentrations (the "a"
figures) and the solid phase concentrations (the "b" figures) throughout the soil column (0 for the surface and -600
cm for the water table) for various fixed times and various values of the sorption rate co. Figure H-3 is for a rate co of
6.5 x lO^d'1. It is obvious that the c- and s-distributions are not in balance, and nonequilibrium sorption conditions
do exist. An interesting observation is that the earlier time distributions for c have the classic "bell-shape," while the
later time distributions tend to lose this shape (also see the dotted curve in Figure H-2a). The reason for the non-bell-
shapes is that the "history" term in Equation (H-14) has had more time to accumulate the effects from earlier
distributions and has thus changed the classic diffusion equation. Figure H-4 is for a rate co of 6.5 x lO'M'1. These
distributions also indicate the existence of nonequilibrium conditions for the earlier time curves; while for the later
time curves, the c-distributions are rapidly approaching those of the solid phase, which indicates the approach of
equilibrium sorption at all sites in the soil matrix. Figure H-5 is for a rate co of 6.5 x lO^d'1. As one can easily see,
the shapes of the c-distributions are nearly the same as those for the s-distributions; however, the peaks of the s-
distributions are slightly greater than those of the c-distributions at the earlier times. As time increases the peaks for
both sets of distributions approach one another and the sorption sites approach equilibrium. In addition, both sets of
curves exhibit the classic "bell-shape" which indicates the reduced roll of the "history" term in Equation (H-14). The
H-6
-------
(a)
10-"
10-*
10"*
I/I
(b)
10-*
1.2x10-*
.
-
. Kd=1.0 ml/g
HYDRUS "Sr
B = 6. 5x 1 0-« d-'
a = 6.5x10-* d-(
a = 6. 5x 1 0-* d-1
a =6.5x10-* d"1
B = 6. 5x 1 0-* d'1
B = 6. SxlO-' d"1
f =0
f =0
f =0
f =0 -
f =0
f =0
1.0x10-*
8.0x10-
o>
£
I/I
E
*c
.a
= 6.0x10-*
4.0x10-
2.0x10-
20000 40000
Time (days)
60000
80000
HYDRUS
•°Sr
Kd=1.0 ml/g
= 6. 5x10-* dH
= 6. 5x10-= dH
= 6.5x10-* d-*
= 6.5x10-* d-*
= 6.
f =0
f=0
f=0
f=0
f=0
f =0
20000 40000
Time (days)
60000
80000
Figure H-2. (a) Breakthrough curves for the liquid phase concentration at the 6 m level, (b) Concentration
curves for the nonequilibrium solid phase at the 6 m depth. Curves for co = 6.5 x lO^d'1 and 6.5 x
KHd'1 are basically the same for both (a) and (b).
Table H-2. Liquid Phase and Solid Phase Peak Concentrations at the Hypothetical Water Table for a Sequence
of Sorption Rates, along with the Corresponding Times to Arrive at Those Peaks.
Sorption
Rate, CD
(d-1)
6.5 x 10-6
*6.5 x 10-5
6.5 x 10-4
6.5 x 10-3
**3.2xlO-2
6.5 x lO'2
6.5 x lO'1
Solid Phase Concentration
Speak (mg/kg) Tpeak (d)
6.5 x 10-4
6.7 x 10-4
4.1xlO-4
5.6 x 10-4
6.2 x 10-4
6.2 x 10-4
6.2 x 10-4
6,300
9,500
37,100
43,100
43,500
43,500
43,500
Liquid Phase Concentration
Cpeak(mg/L)
6.3 x lO'2
6.2 x lO'3
4.6 x lO'4
5.7 x lO'4
6.2 xlO'4
6.2 xlO'4
6.2 xlO'4
T ,(d)
peak v '
4,900
4,700
35,200
43,200
43,500
43,500
43,500
Decay Rate p for 90Sr is 6.55 x
Data from Figure H-1.
H-7
-------
final sets of distributions are given in Figure H-6, which are for co = 6.5 x 10-2cH; the curves for co = 6.5 x KHcH
were found to be "graphically" equivalent to those for co = 6.5 x lO^d'1. As one can see, for all times greater than
or equal to 2000 days, the c- and s-distributions are the same, and there are no nonequilibrium sites within the soil
column.
References
Drake, R.L., J-S. Chen and D.G. Jewett (2002). "An exact solution for the assessment of nonequilibrium sorption of
radionuclides in the vadose zone." Waste Management 2002 Conference. February 24-28, 2002, Tucson, AZ.
Available from www.wmsym.org/wm02 (click on WM'02 Proceedings).
Thibodeaux, L.J. 1996. Environmental Chemodynamics: Movement of Chemicals in Air. Water, and Soil. 2nd
Edition. John Wiley & Sons, Inc. New York, NY. 593pp.
Toride, N., F.J. Leij, and M.T. van Genuchten. 1993. "A comprehensive set of analytical solutions for
nonequilibrium solute transport with first-order decay and zero-order production." Water Resour. Res. 29(7),
2167-2182.
van Genuchten, M. Th., and RJ. Wagenet. 1989. "Two-site/two-region models for pesticide transport and
degradation: theoretical development and analytical solutions." Soil Sci. Soc. Am. J. 53, 1303-1310.
(a)
(b)
X. 8.0x10-*
ot
&
g 6.0x10-"
1
§ 4.0x10"*
Q
C
o
o
m 2.0x10-"
1.6x10-"
•c 1.2x10-"
.Q
"5
S 8.0x10-"
o
z
o 4.0x10""
HYDRUS Kd=1.0 ml/g
"°Sr a = 6.5x10"* d~* (f^))
2000 days
3000 days
4000 days
5000 days
-100 -200 -300 -400 -500 -600
Depth (cm)
HYDRUS Kd=1.0 ml/g
'"Sr a = 6.5x10-" f (f=0)
2000 days
3000 days
4000 days
5000 days
10000 days
20000 days
•r'-.T^rsgg™
-100 -200 -300 -400 -500 -600
Depth (cm)
r
HYDRUS Kd=1.0 ml/g
•"Sr a = 6.5x10-* dM (f=0)
2000 days
4000 days
6000 days
8000 days
10000 days
20000 days
-100 -200 -300 -400 -500 -600
Depth (cm)
5. OxIO-2
4. OxIO-2
3. CxIO"2
2.0x10-"
1.0x10-"
HYDRUS
"Sr a = e.
Kd=1.0 ml/g
* d-1 (f=0)
2000 days
4000 days
6000 days
8000 days
10000 days
20000 days
0 -100 -200 -300 -400 -500 -600
Depth (cm)
Figure H-3. Liquid phase concentration curves (a)
and solid phase concentration curves
(b) for 90Sr, for various times and for
co = 6.5 x 10-5d'1, where zero depth
is the surface and -600 cm is the
hypothetical water table.
Figure H-4. Liquid phase concentration curves (a) and
solid phase concentration curves (b) for
90Sr, for various times and for co = 6.5 x
10-4d-!, where zero depth is the surface
and -600 cm is the hypothetical water
table.
H-8
-------
(b)
HYDRUS Kd=1.0 ml/g
ir a = 6.5x10-" d"1 (f=0)
2000 days
4000 days
6000 days
BOOO days
20000 days
40000 days
(a)
-100 -200 -300 -400 -500 -600
Depth (om)
8.0x10-*! 1 1 1 1
HYDRUS Kd=1.0 ml/g
•°Sr a = 6.5X10-1 f (f=0)
2000 days
4000 days
6000 days
c 11 \ 8000 days
£ I \ 20000 days
£ I 40000 days
= 4.0x10-4 '
0 -100 -200 -300 -400 -500 -600
Depth (cm)
Figure H-5. Liquid phase concentration curves (a)
and solid phase concentration curves
(b) for 90Sr, for various times and for co
= 6.5 x lO-M'1, where zero depth is
the surface and -600 cm is the
hypothetical water table.
(b)
HYDRUS Kd=1.0 ml/g
"Sr o = 6.5x10-* d-1 (f=0)
2000 days
4000 days
6000 days
8000 days
20000 days
40000 days
-100 -200 -300 -400 -500 -600
Depth (om)
v? 8. Ox10-*r
oi | HYDRUS
Kd=1.0 ml/g
= 6.5x10-* dH (f=0)
2000 days
4000 days
6000 days
BOOO days
20000 days
40000 days
-100 -200 -300 -400 -500 -600
Depth (om)
Figure H-6. Liquid phase concentration curves (a) and
solid phase concentration curves (b) for 90Sr,
for various times and for co = 6.5 x lO^d'1,
where zero depth is the surface and -600 cm
is the hypothetical water table.
H-9
-------
Appendix I
Results from the Transport and Fate of Other
Radionuclides Not Considered in the Main Text
The fate and transport of radionuclides in a finite soil column (e.g., the five parent radionuclides listed in Section
5.1 passing through the 6 m soil column considered in this report) are strongly influenced by the distribution
coefficient, Kd, the radioactive decay rate, u, and the recharge rate, q. One could imagine that there exists some
decay-mobility scale (DMS) that represents the basic characteristics of a radionuclide as it passes through the finite
column. For example, for the smaller values of the DMS, the radionuclide could be long-lived and highly mobile in
the soil column, and thus have characteristics somewhat similar to a conservative species, Region I of the curve in
Figure 1-1. In this figure, the ordinate is the concentration, C, at the bottom of the soil column, normalized by the
concentration, C0, of the source at the top of the column, the ordinate being denoted as [C/C0]B. The values of
[C/C0]B for Region I of the curve, which are less than one, represent diffusive dilution but very little dilution due to
radioactive decay. For the larger values of DMS, the radionuclides could be short-lived and highly immobile in the
soil column, thus producing small values of [C/C0]B as shown in Region III of the curve in Figure I-1. The
radionuclides with intermediate half-lives and values of Kd would be characterized by the intermediate values of
DMS and Region II of the curve. However the exact form of the DMS is unknown, but its qualitative properties are
as represented in Figure 1-1. One reason a DMS is difficult to quantify is the fact that Kd is a lumped parameter
which represents many unknown geochemical processes in natural soil columns. Thus, a Kd value for a given
radionuclide is not unique but is highly dependent on soil structure, texture, and geochemistry (U.S. EPA 1999a,b;
2000a,b). Consequently, a DMS value for a given radionuclide would not be expected to be unique either. But, the
qualitative trends stated above, and shown in Figure 1-1, are potentially correct, and have value in relating diffusive
dilution and radioactive decay to pollutant travel time.
Decay-Mobility Chart
(Arbitrary Scale)
DMS
Figure 1-1. The normalized concentration of a radionuclide at the bottom of a soil column (the source being at the
top of the column) versus a hypothetical decay-mobility scale (DMS), where I represents highly mobile,
long-lived species, III represents highly immobile, short-lived species, and II represents species with
intermediate mobilities and half-lives.
1-1
-------
The five parent radionuclides given in Section 5.1, along with their default Kd values and their half-lives listed in
Figure 5-1, can be located "on the curve" in Figure 1-1 for this report's 6 m, homogeneous soil column. These
proposed locations are as follows:
"Tc and 238U are represented by Region I;
Travel time through the 6 m column for 238U is about 5 times greater than that for "Tc, thus diffusive
dilution is much greater for 238U than for "Tc;
3H is represented by the left-hand side of Region II;
90 Sr is represented by the right-hand side of Region II; and
238pu is represented by Region III.
This appendix reports the results of two sets of simulations illustrating the validity of the placement of the five
parent radionuclides on the curve in Figure I-1. Using the CHAIN and FECTUS Codes, the differences between the
transport and fate of a parent and that of a daughter are displayed for "Tc and its daughter, ruthenium "Ru, for the
6 m, homogenous soil column. The other set of calculations used the CHAIN Code to compute the Cpeak - and Tpeak-
values for the breakthrough curves (BTC's) at the 6 m level, for each of the five radionuclides under the influence of
a range of recharge rates, thus verifying the proposed locations of the radionuclides "along the curve" in Figure I-1.
LI The CHAIN Governing Equations
The CHAIN Equations covering the transport and fate for both these comparative exercises are as follows:
at e + pKd
eo a2c
+ pK, 3z2
at
9D
pic; az2
- u c
,
+
e + PKd
- ^-Hr
6 + pK*d
(1-2)
where c is the concentration of the parent in mg/L and the asterisked (*) quantities represent those for the daughter.
Using the base values for (q, 9, p, D) given in Section 6, the default values of Kd for the radionuclides given in
Section 5, and the half-lives given in Figure 5-1, the coefficients for the advection term, the diffusion term, and the
sink term in Equation (I-1) can be evaluated. These coefficients for the five radionuclides are given in Table I-1.
The relationship between \a and the half-life, t1/2, is given by:
H= (Ioge2)-t
1/2
(1-3)
Table 1-1 shows that as Kd increases from 0 to 5 ml/g, the seepage velocity (advection term coefficient)
decreases by a factor of 54, as does the diffusion term coefficient. Even though the diffusive mechanism decreases
greatly, it may still be very important because the great decrease in seepage velocity allows the diffusion to act over a
longer time, as exemplified by the following:
Radionuclide
3H
"Tc
238U
90Sr
238pu
(600 cm) -=- (Advection
Term Coefficient)
4,000 d
4,300 d
21,000 d
46,500 d
216,600 d
1-2
-------
Table 1-1. Coefficients of the Advection, Diffusion, and Sink Terms in Equation (1-1).
Radionuclide
Tritium 3H
Technetium "Tc
Uranium 238U
Strontium 90Sr
Plutonium 238Pu
Initial
Concentration
In Recharge
Water Cg
(mg/L)
l.lxlO-7
1.25 X ID'2
417
2.0 x ID'1
1.0 x 10-7
Advection Term*
Default
Value
Kd
(ml/g)
0
0.007
0.4
1
5
Half-Life
ty
(d)
4380
76,650,000
1.6425xl012
10,585
32,120
First-Order *
Decay Rate
M
(a-1)
1.5825xlO-4
9.043xlO-9
4.2201xlO-13
6.5484xlO-5
2.1580xlO-5
Coefficient
q
e + PKd
(cm/d)
0.15000
0.13962
0.02857
0.01290
0.00277
Diffusion Term*
Coefficient
6D
e + PKd
(cm2/d)
1.00000
0.93077
0.19048
0.08602
0.01848
q, 0, p, and D are assigned the base values in Section 6;
Kd is assigned values given in Column Three;
|i is obtained from half-life t1/2 using Equation (1-3).
Comparing the times, (600 cm + advection term coefficient), with the half-lives listed in Table 1-1, one observes
that the first-order decay rate [i will probably be most important for the transport and fate of 238Pu through the 6 m
layer, followed by that of 90Sr, followed by that of 3H. Because of the very long half-lives for 99Tc and 238U, their
decay rates are not very important over the 6 m layer with respect to the BTCs of these parents. However, with
respect to the production of the daughter, 99Ru from 99Tc and 234U from 238U, these very small decay rates may be
important.
1.2 Breakthrough Curves for 99Tc and It's Daughter, 99Ru
For the decay chain, 99Tc -> 99Ru, 99Tc has a very low decay rate \i which has negligible effect on the ETC for
99Tc over the 6 m depth of the current simulation. However, this low decay rate does produce a small amount of the
stable radionuclide, 99Ru, while the 99Tc moves from its surface source (of a 1,000 days duration) down through the
6 m soil column. The lifetime of the moving source of 99Ru equals the lifetime of the 99Tc BTCs for the given
recharge rate q = 0.0024 cm/d. This lifetime is equal to or less than 10,000 days while the half-life of 99Tc is about
76,650,000 days. The travel time for 99Ru through the layer is much greater than that for 99Tc because the Kd for
99Ru is 5 ml/g versus that of 0.007 ml/g for 99Tc. However, the travel time for 99Ru through the 6 m layer should be
on the same order as that of 238Pu, namely about 200,000 days, since the Kd values of 99Ru and 238Pu are taken to be
equal. However, 238Pu has a significant decay rate (238Pu -> 234U), while 99Ru is stable; |J* in Equation (1-2) is zero
for 99Ru. Since FECTUZ has a different 9D-value than CHAIN (see Equation 1-2), these two codes were used to see
the effects of diffusion on the BTCs of 99Ru for the 6 m layer. The diffusion term coefficient for FECTUZ is about
2/3 that of CHAIN.
The relative magnitude of the concentration c* of 99Ru produced by the decay of 99Tc can be obtained from
Equation (1-2) if one ignores all terms except the source term and the time derivative of c*:
C* = 1.8 x 10'
J c(t)dt .
(1-4)
The integral term in Equation (1-4) can be approximated from the ETC for 99Tc under the influence of a 0.024
cm/d recharge rate, as given in Figure I-2b. We can roughly approximate the area under this curve as given by 0.001
mg/L times 10,000 d. Thus, the concentration of 99Ru at 10,000 d is equal to 1.8 x 10'9 mg/L by Equation (1-4). The
FECTUZ/CHAIN result given in Figure I-2a is very nearly equal to this value at t = 10,000 d. Thus, over the first
10,000 d of the evolution of 99Ru, the dominant term in Equation (1-2) is the source of 99Ru being produced by the
decaying 99Tc.
1-3
-------
(a)
o>
E
c
o
c
-------
parent at the top of the column. However, the evolutions of parent and daughter species can be quite complex in
two-dimensional and three-dimensional, heterogeneous soil columns, as pointed out in Section 3.6 (Oldenburg and
Pruess, 1996).
1.3. The Sensitivity of the Five Parent Radionuclides to Recharge Rate
Figures 1-3 (a to e) illustrate the sensitivity of the transport and fate of the five parent radionuclides through the 6
m, homogenous soil column to a variable recharge rate, q, using the CHAIN Code for the simulation. The five sets
of curves are given in order of increasing default Kd values. As indicated by Equation (1-1), these Kd values govern
the rates of pollutant passage through the 6 m layer. In fact, a good representation of characteristic times of passage
(a)
5. On 10-"
•.3.75K10-"
2.5tt10-°
R.ohorg. Rat* (om/d)
0.014
CHAIN
Tim* (day*)
(d)
2. OilO
"Sr
(b)
0.010
I
I
'; 0.004
i
0.002
"To
Rtoharg* Rot* (om/d)
0.014
0.024
0.032
CHAIN
Tim* (days)
(e)
(c)
100.000
?
I
| 40.000
Rcoharg* Rat* (om/d)
0.014
0.024
0.032
CHAIN
2.0»fO
Raoharge Roto (om/d)
0.014
0.024
0.032
CHAIN
Tlm« (days)
Rwhargc Roto (om/d)
= 8:811
0.032
CHAIN
Tim. (days)
Tim* (days)
Figure 1-3. Sensitivity of radionuclide transport through the unsaturated zone to recharge rate (q) using the CHAIN
Model: (a) Tritium 3H; (b) Technetium "Tc; (c) Uranium 238U; (d) Strontium 90Sr; and (e) Plutonium
238Pu.
1-5
-------
can be obtained from the quotient 600 cm + v, where v is the seepage velocity, or coefficient of the advection term in
Equation (I-1):
v =
q
PKd
(1-5)
These characteristic times are listed in Table 1-2 for the five radionuclides and the three recharge rates q. Also listed
in this table are the normalized peaks of the BTCs shown in Figures 1-3 (a to e). These peaks are normalized by the
respective source concentrations of the radionuclides given in Table 1-1. Corresponding to these normalized C k
values are the times to peak concentration, T k, given in the last column. One should note the rather close
correspondence between the "characteristic time" values and the values of T k given in Table 1-2. The final column
to note in this table is the "type of decay." The descriptors given in the table can be made more explicit from the log-
log plot of Cpeak ^ C0 versus Tpeak in Figure 1-4.
To further illuminate the results in Table 1-2 and Figure 1-4, we introduce a relative diffusion factor (DIP) and a
relative decay factor (DEC). The value of DIP is defined as the quotient of the species Tpeak for q = 0.024 cm/d
divided by the Tpeak for the highly mobile radionuclide, 3H, at a q = 0.024 cm/d. The value of DEC is defined as the
quotient of the species Tpeak for q = 0.0024 cm/d divided by the species decay half-life. Table 1-3 lists these factors
for the five parent radionuclides. As DIP increases, the amount of diffusive dilution increases, while an increase in
DEC produces an increase in the amount of radioactive decay. Using these two factors, the locations of the five
curve segments in Figure 1-4 can be easily explained:
The diffusive dilutions of 3H and "Tc are nearly equal since their DIP values are nearly equal, and the
time ranges of their curve segments are nearly the same.
The curve segment for 3H is below that for "Tc because 3H has significantly decayed over the time
interval, while the decay of "Tc has been insignificant, although not zero as seen in Figure I-2a.
The curve segment for 238U, as that of "Tc, has only experienced insignificant decay, but is lower than
the curve segments for "Tc and 3H, and further to the right, because its DIP is almost five times that of
"Tc and 3H. The 238U's normalized concentrations are lower only because of diffusive dilution.
Table 1-2. Cpeak Normalized by Source Concentration C0 and Tpeak for Each Radionuclide and
for Three Recharge Rates.
Radionuclide
Tritium 3H
Technetium "Tc
Uranium 238U
Strontium 90Sr
Plutonium 238Pu
*Type
of Decay
Significant
Very Little
Very Little
Very Significant
Very Significant
*Type of Decay with respect **
to the 6 m soil column.
Recharge * * Characteristic Time
Rate q (600cm) + v
(cm/d) (d)
0.014
0.024
0.032
0.014
0.024
0.032
0.014
0.024
0.032
0.014
0.024
0.032
0.014
0.024
0.032
v = seepage velocity
6860
4000
3000
7370
4300
3225
36,000
21,000
15,750
79,715
46,500
34,875
371,145
216,600
162,375
q
9 + PK,
* "Normalized
peak
C *+C
peak o
1.72x10-'
3.23 x lO'1
4.18x10-'
4.49 x 10-1
5.64 x 10-1
6.31x10-'
9.98 x 10-2
1.29x10-'
1.48x10-'
4.47 x 10-4
3.28 xlO'3
7.44 x 10-3
1.66xlO-5
2.20 x 10-4
6.48 x 10-4
T ,
peak
(d)
7300
4400
3300
8000
4750
3560
35,500
21,000
15,900
66,100
43,100
33,200
287,000
194,000
152,000
*** C = Concentration at Source,
in Table I- 1.
1-6
-------
10 10 10 10
Time to Peak Concentration (Tpeak In days)
Figure 1-4. Cpeak -^ C0 versus Tpeak in days for five radionuclides, showing the effects of various Kd values and
various decay rates, where the time intervals of the curve segments are determined by the range of
discharge rates for the 6 m soil column.
The curve segments for 90Sr and 238Pu are both much lower and to the right of those 3H, "Tc and 238U
because their DEC values are much larger and their DIP values are also much larger. This means that
their normalized concentrations are lower than the other three species because of both greater diffusion
and greater decay.
The curve segment for 238Pu is lower and further to the right of that of 90Sr since its DIP is four and
one-half times that of 90Sr and its DEC is about one and one-half times that of 90Sr. This means that
over comparable time segments (i.e., the Tpeak values for the range of q for the 6 m soil column) more
of the 238Pu has diffused and decayed than that of 90 Sr.
Table 1-3. The Relative Diffusion Factors and the Relative Decay Factors for the Five Parent Radionuclides.
For Recharge Rate, q = 0.024 cm/d
Radionuclide
Tritium 3H
Technetium "Tc
Uranium 238U
Strontium 90Sr
Plutonium 238Pu
*Type
of Decay
Significant
Very Little
Very Little
Very Significant
Very Significant
Relative Diffusion
Factor
1.000
1.080
4.773
9.795
44.091
Relative Decay
Factor
1.005
0.000
0.000
4.072
6.040
*Type of Decay with respect to the 6 m soil profile.
1-7
-------
References
Oldenburg, CM. and K. Pruess. 1996. "Mixing with first-order decay in variable-velocity porous media flow."
Transport in Porous Media. 22:161-180.
U.S. EPA. 1999a. Understanding Variations in Partition Coefficients. Kd. Values: Model. Methods of
Measurement, and Application of Chemical Reaction Codes. Volume I. EPA 402/R-99/004A. Office of
Radiation and Indoor Air & Office of Solid Waste and Emergency Response. U.S. EPA. Washington, DC.
U.S. EPA. 1999b. Understanding Variations in Partition Coefficients. Kd. Values: Review of Geochemistry and
Available Kd Values for Cadmium. Cesium. Chromium. Lead. Plutonium. Radon. Strontium. Thorium. Tritium.
(3HX and Uranium. Volume II. EPA402/R-99/004B. Office of Radiation and Indoor Air and Office of Solid
Waste and Emergency Response. U.S. EPA. Washington, DC.
U.S. EPA. 2000a. Soil Screening Guidance for Radionuclides: User's Guide. EPA 540/R-00/007. Office of
Radiation and Indoor Air and Office of Solid Waste and Emergency Response. U.S. EPA. Washington, DC.
U.S. EPA. 2000b. Soil Screening Guidance for Radionuclides: Technical Document. EPA 540/R-00/006. Office
of Radiation and Indoor Air and Office of Solid Waste and Emergency Response. U.S. EPA. Washington, DC.
1-8
-------
S-EPA
United States
Environmental Protection
Agency
Please make all necessary changes on the below label,
detach or copy, and return to the address in the upper
left-hand corner.
If you do not wish to receive these reports CHECK HFRFlT
detach, or copy this cover, and return to the address in the
upper left-hand corner.
PRESORTED STANDARD
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
National Risk Management
Research Laboratory
Cincinnati, OH 45268
Official Business
Penalty for Private Use
$300
EPA/600/R-02/082
July 2002
-------