&EPA
United States
Environmental Protection
Agency
Environmental Research
Laboratory
Athens GA 30613
EPA/600/3-89/048a
July 1989
Research and Development
Risk of
Unsaturated/Saturated
Transport and
Transformation of
Chemical
Concentrations (RUSTIC)
Volume I: Theory and
Code Verification
-------
EPA/600/3-89/048a
RISK OF UNSATURATED/SATURATED TRANSPORT
AND TRANSFORMATION OF CHEMICAL
CONCENTRATIONS (RUSTIC)
Volume I: Theory and Code Verification
by
J.D. Dean,1 P.S. Huyakorn,^ A.S, Donigian, Jr.,-*
K.A. Voos,1 R.W. Schanz,1 Y.J. Meeks1, and R.F. Carsel4
Woodward-Clyde Consultants1
Oakland, CA 94607
n
HydroGeologic''
Herndon, VA 22070
AQUA TERRA Consultants3
Mountain View, CA 94043
Contract No. 68-03-6304
Project Officer
Robert F. Carsel
Assessment Branch
Environmental Research Laboratory^
U.S. Environmental Protection Agency
Athens, GA 30613
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GA 30613
..., u, 6°604
-------
DISCLAIMER
The information in this document has been funded wholly or in part by the
United States Environmental Protection Agency under Contract No. 68-03-6304
to Woodward-Clyde Consultants. It has been subject to the Agency's peer and
administrative review, and it has been approved for publication as an EPA
document. Mention of trade names of commercial products does not constitute
endorsement or recommendation for use by the U.S. Environmental Protection
Agency.
ii
-------
FOREWORD
As environmental controls become more costly to implement and the penalties
of judgment errors become more severe, environmental quality management
requires more efficient analytical tools based on greater knowledge of the
environmental phenomena to be managed. As part of this Laboratory's research
on the occurrence, movement, transformation, impact, and control of
environmental contaminants, the Assessment Branch develops management or
engineering tools to help pollution control officials reach decisions on the
registration and restriction of pesticides used for agricultural purposes.
The pesticide regulatory process requires that the potential risk to human
health resulting from the introduction or continued use of these chemicals be
evaluated. Recently much of this attention has been focused on exposure
through leaching of pesticides to groundwater and subsequent ingestion of
contaminated water. To provide a tool for evaluating this exposure, the
RUSTIC model was developed. RUSTIC simulates the transport of field-applied
pesticides in the crop root zone, the unsaturated zone, and the saturated
zone to a drinking water well, taking into account the effects of
agricultural management practices. The model further provides estimates of
probable exposure concentrations by taking into account the variability in
the natural systems and the uncertainties in system properties and processes.
Rosemarie C. Russo, Ph.D.
Director
Environmental Research Laboratory
Athens, Georgia
iii
-------
ABSTRACT
This publication contains documentation for the RUSTIC model. RUSTIC links
three subordinate models in order to predict pesticide fate and transport
through the crop root zone, unsaturated zone, and saturated zone to drinking
water wells: PRZM, VADOFT, and SAFTMOD. PRZM is a one-dimensional finite-
difference model which accounts for pesticide fate and transport in the crop
root zone. This release of PRZM incorporates several features in addition to
those simulated in the original PRZM code: specifically, soil temperature
simulation, volatilization and vapor phase transport in soils, irrigation
simulation and a method of characteristics (MOC) algorithm to eliminate
numerical dispersion. PRZM is now capable of simulating fate and transport
of the parent compound and up to two daughter species. VADOFT is a one-
dimensional finite-element code which solves the Richard's equation for flows
in the unsaturated zone. The user may make use of constitutive relationships
between pressure, water content, and hydraulic conductivity to solve the flow
equations. VADOFT may also simulate the fate and transport of two parent and
two daughter products. SAFTMOD is a two-dimensional finite-element model
which simulates saturated solute flow and transport in either an X-Y or X-Z
configuration. The codes are linked together with the aid of a flexible
execution supervisor which allows the user to build models which are tailored
to site-specific situations. In order to perform exposure assessments, the
code is equipped with a Monte Carlo pre- and post-processor.
This report on the RUSTIC modeling system was submitted in fulfillment of two
separate work assignments under Contract No. 68-03-6304 by Woodward-Clyde
Consultants under the sponsorship of the U.S. Environmental Protection
Agency. RUSTIC was developed by a Project Team headed by Woodward-Clyde
Consultants as the prime contractor, with AQUA TERRA Consultants and
HydroGeologic, Inc. as subcontractors. This report describes the work
performed from January 1986 through September 1988. The final manuscript was
prepared for publication and printing by AQUA TERRA Consultants under EPA
Contract No. 68-03-3513.
IV
-------
TABLE OF CONTENTS
Page
Foreword iii
Abstract iv
Figures viii
Tables xiii
Acknowledgments xv
1.0 Introduction 1
1.1 Background and Obj ectives 1
1.2 Concept of Risk and Exposure Assessment 2
1. 3 Overview of RUSTIC 6
1.3.1 Overview of PRZM 9
1.3.2 Overview of the Vadose Zone Flow and Transport Model
(VADOFT) 11
1.3.3 Overview of the Saturated Zone Flow and Transport
Model (SAFTMOD) 12
1.3.4 Model Linkage 13
1.3.5 Monte Carlo Processor 14
1.3.6 Overview Summary 15
1.4 Overview of Volumes I and II 15
2.0 Pesticide Root Zone Model (PRZM Release II) 17
2 .1 Introduction and Background 17
2.1.1 Introduction 17
2.1.2 Background 18
2 . 2 Features and Limitations 18
2.2.1 Features 18
2.2.2 Limitations 20
2.3 Description of the Equations 22
2.3.1 Transport in Soil 22
2.3.2 Water Movement 30
2.3.3 Soil Erosion 37
2.3.4 Volatilization 38
2.3.5 Irrigation Equations 56
2.4 Numerical Solution Techniques 59
2.4.1 Chemical Transport Equations 60
2.4.2 Volatilization 62
2.4.3 Soil Temperature 65
2.4.4 Furrow Irrigation 66
2.5 Results of PRZM Testing Simulations 68
2.5.1 Transport Equation Solution Options 68
2.5.2 Testing Results of Volatilization Subroutines 69
v
-------
TABLE OF CONTENTS (continued)
Page
2.5.3 Testing Results of Soil Temperature
Simulation Subroutine 84
2.5.4 Testing of Daughter Products Simulation 85
3.0 Vadose Zone Flow and Transport Model (VADOFT) 95
3 .1 Introduction 95
3 . 2 Overview of VADOFT 95
3.2.1 Features 95
3.2.2 Limitations 97
3 . 3 Description of Flow Module 97
3.3.1 Flow Equation 97
3.3.2 Numerical Solution 103
3 .4 Description of the Transport Module 107
3.4.1 Transport Equation 107
3.4.2 Numerical Solution of the Transport Equation 108
3 . 5 Results of VADOFT Testing Simulations HO
3.5.1 Flow Module (Variably Saturated Flow Problems) no
3.5.2 Transport Module HI
3.5.3 Combined Nonlinear Flow and Transport Modules 115
4.0 Saturated Zone Flow and Transport Model (SAFTMOD) 135
4.1 Introduction 135
4. 2 Overview of SAFTMOD 135
4.2.1 Features 135
4.2.2 Limitations 138
4. 3 Description of Flow Module 138
4.3.1 Flow Equation 138
4.3.2 Numerical Solution of Flow Problems 140
4.4 Description of the Transport Module 144
4.4.1 Transport Equation 144
4.4.2 Numerical Solution of the Transport Problems 145
4.4.3 Treatment of a Decay Chain 146
4.4.4 Spatial and Temporal Discretizations 147
4.5 Results of SAFTMOD Testing Simulations 149
4.5.1 Flow Module 149
4.5.2 Transport Module 153
5.0 Linkage of VADOSE and Saturated Zone Models 167
5 .1 Temporal Model Linkage 169
5. 2 Spatial Linkages 169
5.2.1 PRZM and VADOFT 169
5.2.2 VADOFT and SAFTMOD 171
6.0 Uncertainty Preprocessor 177
6 .1 Introduction 177
6.2 Overview of the Preprocessor 177
6.2.1 Description of the Method 178
6.2.2 Uncertainty in the Input Variables 179
vi
-------
TABLE OF CONTENTS (concluded)
Page
6.3 Description of Available Parameter Distributions 180
6.3.1 Uniform Distribution 180
6.3.2 Normal Distribution 181
6.3.3 Log-Normal Distribution 181
6.3.4 Exponential Distribution 182
6.3.5 The Johnson System of Distributions 182
6.3.6 Triangular Distribution 183
6.3.7 Empirical Distribution 183
6.3.8 Uncertainty in Correlated Variables 185
6.3.9 Generation of Random Numbers 186
6.4 Analysis of Output and Estimation of Distribution
Quantiles 188
6.4.1 Estimating Distribution Quantiles 189
6.4.2 Confidence of u_ 191
7 .0 References 193
8 .0 Appendices 201
8 .1 Quality Assurance 201
8.1.1 Algorithms 201
8.1.2 Software 202
8.1.3 References 203
vii
-------
FIGURES
Page
1.1 Decision path for risk assessment 3
1.2 Time series plot of toxicant concentrations 5
1.3 Frequency distribution of toxicant concentrations 5
1.4 Cumulative frequency distribution of toxicant concentrations 5
1.5 Time series of toxicant concentrations with moving
average window of duration tc 7
1.6 Linked modeling system configuration 8
2.1 Pesticide root zone model 19
2.2 Schematic representation of a single chemical in a soil layer 23
2.3 Schematic of pesticide vapor and volatilization processes 39
2.4 Variability of infiltration depths within an
irrigation furrow 58
2.5 Schematic of the top two soil compartments and the overlying
surface compartment: (a) without plant canopy and (b) with
plant canopy 64
2.6 Comparison of simulation results at high Peclet number 70
2.7 Comparison of simulation results at low Peclet number 71
2.8 Comparison of volatilization flux predicted by PRZM
and Jury's analytical solution: test cases #1 and #2 74
2.9 Comparison of volatilization flux predicted by PRZM and
Jury's analytical solution: test cases #3 and #4 75
2.10 Sensitivity of cumulative volatilization flux to Kd
and decay rate 79
2.11 Effects of DELX on volatilization flux and pesticide decay 80
viii
-------
FIGURES (continued)
Page
2.12 Comparison of constant and two-step decay rates 82
2.13 Effects of two-step decay rates on volatilization flux
and pesticide decay 83
2.14 Comparison of soil temperature profiles predicted by
analytical and finite difference solutions 86
2.15 Comparison of soil temperature profiles predicted by
analytical and finite difference solutions 87
2.16 Schematic of a system of parent and daughter pesticides 89
2.17 Conversion of Ci to Cj to Co with no adsorption
and no decay 92
2.18 Conversion of C^ to €2 to Cj with decay but no
adsorption 93
2.19 Conversion of Ci to Cy to C? with decay and adsorption 93
2.20 Conversion of aldicarb to aldicarb sulfoxide to aldicarb sulfone.. 94
3.1 Cross-section view showing the solute migration path
in the vadose and saturated zones of an unconfined aquifer 96
3.2 Logarithmic plot of constitutive relations for clay, clay
loam, loam, and loamy sand: (a) saturation vs. capillary
head and (b) relative permeability vs. saturation 100
3.3 Logarithmic plot of constitutive relations for silt, silty
clay loam, silty clay and silty loam: (a) saturation vs.
capillary head and (b) relative permeability vs. saturation 101
3.4 Logarithmic plot of constitutive relations for sandy clay,
sandy clay loam, sandy loam, and sand: (a) saturation vs.
capillary head and (b) relative permeability vs. saturation 102
3.5 Standard plot of relative permeability vs. saturation for
clay, clay loam, loam, and loamy sand 103
3.6 Standard plot of relative permeability vs. saturation for
silt, silt clay loam, silty clay, and silty loam 104
3.7 Standard plot of relative permeability vs. saturation
for sandy clay, sandy clay loam, sandy loam, and sand 105
ix
-------
FIGURES (continued)
Page
3.8 Finite-element discretization of a soil column
showing node and element numbers 106
3.9 Simulated pressure head profiles for the problem
of transient upward flow in a soil column 113
3.10 Simulated profile of water saturation for the
problem of transient upward flow in a soil column 114
3.11 Simulated pressure head profiles for five cases of
the problem of steady infiltration in a soil column 116
3.12 Simulated profiles of water saturation for five
cases of the problem of steady infiltration in a
soil column 117
3.13 Simulated concentration profiles for the problem
of solute transport in a semi-infinite soil column 118
3.14 Simulated concentration profiles for two cases of
the problem of solute transport in a soil column of
finite length: (a) A = 0 d'1, and (b) A = 0.25 d'1 120
3.15 Simulated outflow breakthrough curve for case 1 of
the problem of solute transport in a layered soil
column 123
3.16 Simulated outflow breakthrough curve for case 2 of the
problem of solute transport in a layered soil column 124
3.17 One-dimensional solute transport during absorption
of water in a soil tube 126
3.18 Simulated profiles of water saturation during
absorption of water in a soil tube 128
3.19 Simulated concentration profiles for the problem
of one-dimensional solute transport during
absorption of water in a soil tube 129
3.20 Problem description for transient infiltration
and transport in the vadose zone 130
3.21 Infiltration rate vs. time relationship used in
numerical simulation 130
3 . 22 Simulated water saturation profiles 131
-------
FIGURES (continued)
Page
3 . 23 Simulated pressure head profiles 132
3.24 Simulated vertical Darcy velocity profiles 133
3.25 Simulated solute concentration profiles 134
4.1 Schematic description of pesticide migration in
unconfined aquifer system subject to pumping 136
4.2 Discretization of aquifer regions: (a) regular
region and (b) irregular region 148
4.3 Problem description of flow to a well in a two-aquifer
system with a confining aquitard 150
4.4 Finite element grids used in simulating well flow in
two aquifer system: (a) axisymmetric grid and
(b) areal grid 151
4.5 Simulated drawdown-time relationships for the pumped
aquifer 152
4.6 Problem description for flow to parallel drains in an
unconfined aquifer 153
4.7 Finite element grids used in the simulation of flow
through parallel drains: (a) areal grid, and
(b) cross-sectional grid 154
4.8 Simulated steady-state and transient water table profiles 155
4.9 Problem description of flow to a gallery 156
4.10 Finite element grid used in simulation of flow to a gallery 156
4.11 Simulated head profiles along the horizontal line
passing through the gallery 157
4.12 One dimensional dispersion in uniform flow:
(a) problem description and
(b) simulated concentration profiles 158
4.13 Areal dispersion in uniform flow: (a) problem
description, and (b) modeled region and boundary
conditions 159
xi
-------
FIGURES (continued)
Page
4.14 Simulated concentration profiles along the x-axis
for case 1 161
4.15 Simulated concentration profiles along the x-axis
for case 2 162
4.16 Areal dispersion of a solute slug in uniform groundwater
flow: (a) problem description, and (b) simulated
concentration profiles 163
4.17 Solute transport from an injection well fully
penetrating an aquifer overlain and underlain by
aquitards 164
4.18 Simulated concentration profiles in the aquifer:
(a) case 1 and (b) case 2 165
5.1 Linked modeling system configuration 168
5 . 2 Resolution of component model timesteps 170
5.3 Definition sketch for a VADOFT/SAFTMOD control
volume (x-y) configuration 172
5.4 Schematic of VADOFT/SAFTMOD x-z configuration 176
6.1 Triangular Probability Distribution 184
xii
-------
TABLES
Page
2-1 Summary of soil temperature model characteristics 48
2-2 Input parameters for the test cases - analytical solution 73
2-3 Trifluralin volatilization losses, amounts remaining in
soil, and estimated losses via other pathways for the
120-day field test 78
2-4 Input parameters for the test cases - Watkinsville site 78
2-5 Simulation results of using different compartment depth
(DELX) 81
2-6 Simulated soil temperature profile after one day for
different compartment thicknesses (time step = 1 day) 88
3-1 Soil properties and discretization data used in
simulating transient flow in a soil column Ill
3-2 Soil properties used in simulating steady-state
infiltration 112
3-3 Iterative procedure performance comparison 114
3-4 Values of physical parameters and discretization
data used in simulating one-dimensional transport
in a semi-infinite soil column 115
3-5 Concentration profile curves at t - 25 hr and
t = 50 hr showing comparison of the analytical
solution and results from VADOFT 119
3-6 Values of physical parameters and discretization
data used in simulating one-dimensional transport
in a finite soil column 121
3-7 Concentration profile curves showing comparison of the
analytical solution and VADOFT 122
3-8 Values of physical parameters used in the simulation
of transport in a layered soil column 124
xiii
-------
TABLES (concluded)
Page
3-9 Breakthrough curves (at z = 86.1 cm) computed using
the analytical solution and VADOFT: Case 1 125
3-10 Breakthrough curves (at z - 86.1 cm) computed using
the analytical solution and VADOFT: Case 2 126
3-11 Values of physical parameters and discretization
data used in simulating transport in a variably
saturated soil tube 127
3-12 Values of physical parameters and discretization data
used in simulating transient infiltration and contaminant
transport in the vadose zone 127
4-1 Values of physical parameters for problem 2 160
4-2 Values of physical parameters and discretization data
for the problem of areal dispersion of a solute slug in
uniform groundwater flow 164
4-3 Values of physical parameters for the problem of solute
transport from an inj ection well 166
xiv
-------
ACKNOWLEDGMENTS
A number of individuals contributed to this effort. Their roles are
acknowledged in the following paragraphs.
Mr. Voos of Woodward-Clyde Consultants (WCC) programmed the execution
supervisor and linked the models within the overall model architecture and
design of RUSTIC developed by Mr. Jack Kittle of Aqua Terra Consultants. The
model linkage was conceived by Mr. Dean and Dr. Atul Salhotra of WCC, Mr.
Kittle, and Dr. Huyakorn. Dr. Huyakorn and his staff wrote the time/space
bridging subroutines for the linkage.
Release II - PRZM was written by Dr. J. Lin and Mr. S. Raju of Aqua Terra
Consultants under the direction of Mr. Donigian. Mr. Schanz (WCC) and
Ms. Meeks (WCC) wrote the irrigation and MOC algorithms. Mr. Dean wrote the
daughter products algorithms which were implemented by Dr. Lin.
The VADOFT code was written and documented by Dr. Huyakorn, H. White, J.
Buckley, and T. Wadsworth of HydroGeologic. Mr. John Imhoff of Aqua Terra
performed many useful testing simulations with VADOFT. SAFTMOD was written
and documented by Dr. Huyakorn and J. Buckley.
The Monte Carlo pre- and post-processors were written by Dr. Salhotra, Mr.
Phil Mineart, and Mr. Schanz of Woodward-Clyde. Final model assembly of the
code and documentation and model testing was performed by Woodward-Clyde
Consultants. Drs. Salhotra and P. Mangarella (WCC) peer reviewed the final
documentation. The support of the editorial and graphics staff of WCC is
appreciated.
The final manuscript was prepared for publication, incorporating review
comments, by Ms. Susan Sharp-Hansen, Mr. S. Raju, and Ms. Dorothy Inahara of
Aqua Terra Consultants.
The authors would like to acknowledge the support of the U.S. Environmental
Protection Agency, and Mr. Lee Mulkey, Chief, Assessment Branch; Mr. Bob
Carsel, Project Officer; and Dr. Rudy Parrish for their suggestions, input,
and helpful comments.
xv
-------
SECTION 1
INTRODUCTION
This publication contains documentation for a linked model, known as RUSTIC,
for contaminant transport in the root, vadose and saturated zones. A brief
section on background and objectives for the model development effort
follows in this introduction (Section 1.1). Section 1.2 gives a synopsis of
risk and exposure assessment concepts. The reader who has sufficient
background in these concepts may proceed to Section 1.3, which provides an
overview of the RUSTIC linked modeling system, including major features and
limitations. The documentation consists of two volumes, and Section 1.4
gives a synopsis of the contents of each. This introduction is common to
both volumes.
1.1 BACKGROUND AND OBJECTIVES
The U.S. Environmental Protection Agency is continually faced with issues
concerning the registration and restriction of pesticides used for
agricultural purposes. Each of these regulatory processes requires that the
potential risk to human health resulting from the introduction or continued
use of these chemicals be evaluated. Recently, much of this attention has
been focused on exposure through leaching of pesticides to groundwater and
subsequent ingestion of contaminated water.
The capability to simulate the potential exposure to pesticides via this
pathway has two major facets:
• Prediction of the fate of the chemical, after it is applied, as it
is transported by water through the crop root zone, the vadose
zone, and saturated zone to a drinking water well
• Evaluation of the probability of the occurrence of concentrations
of various magnitudes at the drinking water well
There are a number of models which are capable of simulating the fate and
transport of chemicals in the subsurface and in the root zone of
agricultural crops. However, none of these models have been linked together
in such a way that a complete simulation package, which takes into account
the effects of agricultural management practices on fate and transport, is
available for use either by the Agency or the agricultural chemical industry
to address potential groundwater contamination problems. Without such a
package, the decision maker must rely on modeling scenarios that are either
incomplete or potentially incorrect. Each time a new scenario arises,
recurring questions must be answered:
-------
• What models should be used?
• How should mass transfer between models be handled?
The resolution of these issues for each scenario is both expensive and time
consuming. Furthermore, it precludes consistency of approach to evaluation
of contamination potential for various scenarios.
The modeling package described in this report seeks to overcome these
problems by providing a consistent set of linked models which have the
flexibility to handle a wide variety of hydrogeological, soils, climate, and
pesticide scenarios. However, the formulation of the risk analysis problem
requires more than a simple, deterministic evaluation of potential exposure
concentrations. The inherent variability of force, capacitance and
resistance in natural systems, combined with the inability to exactly
describe these attributes of the system, suggests that exposure
concentrations cannot be predicted with certainty. Therefore, the
uncertainty associated with the predictions must be quantified.
Consequently, this simulation package also seeks to provide this capability
by utilizing Monte Carlo simulation techniques.
Stated more concisely, the objectives of this model development effort were
to provide a simulation package which can:
• Simulate the fate and transport of field-applied pesticides in the
crop root zone, the unsaturated zone, and the saturated zone to a
drinking water well, taking into account the effects of
agricultural management practices
• Provide probabilistic estimates of exposure concentrations by
taking into account the variability in the natural systems and
uncertainty in system properties and processes
Furthermore, it was desirable that the simulation package be easy to use and
parameterize, and execute on the Agency's DEC/VAX machines. As a result,
considerable effort has gone into providing parameter guidance for both
deterministic and probabilistic applications of the model and software
development for facile model implementation.
1.2 CONCEPT OF RISK AND EXPOSURE ASSESSMENT
Exposure assessment, as defined in the Federal Register (1984) for human
impacts, is the estimation of the magnitude, frequency, and duration at
which a quantity of a toxicant is available at certain exchange boundaries
(i.e., lungs, gut, or skin) of a subject population over a specified time
interval. Exposure assessment is an element of the larger problems of risk
assessment and risk management, as demonstrated in Figure 1.1. The
concentration estimates generated during an exposure assessment are combined
with demographic and toxicological information to evaluate risk to a
population--which can be used, in turn, to make policy decisions regarding
the use or disposal of the chemical.
-------
REGULATORY CONCERN
SCIENTIFIC DATA Population
Exposure
Product Life Cycle
General Information Gathering
Preliminary Exposure
Assessment
Most Probable Areas of Exposure
Preliminary Exposure Awetsmem
ntification ToxiCHy
etc
*
Prelimm»'Y Ruk An»lym
D*cinon
IfvDepth Exposure
Assessment
Exposure A«es$mem
No Need for Future
Exposure A»«jsmeni
Mul1i-DitCiplin»ry
Pe«' Review
Design Atwtirnenl Sludy
Compreheniive D«t» G»th*ring
Conduct Retried Expoture Mo<*«Hini)
Regulatory Response
Dec t linn
nel Review
Hazard Input
al Risk
Deciston
Regulatory Proposal
Examined Exposures
Pre»ent No Unreasonable Ri»k
Figure 1.1. Decision path for risk assessment
-------
Major components of risk assessment are indicated below. Of these, the
first three constitute the important steps for exposure assessment and are
discussed in detail here.
• Characterization and quantification of chemical sources
• Identification of exposure routes
• Quantification of contaminant movement through the exposure routes
to the receptor population/location
• Characterization of the exposed population
• Integration of quantified environmental concentrations with the
characteristics of the exposed populations to yield exposure
profiles
Characterization of sources(s) requires in a broad sense the estimation of
the loading of a chemical into various environmental media. For the
groundwater contamination problem, on a regional scale, this requires data
on chemical uses and distribution of those uses (spatially and temporally).
It also requires information on the crops being grown, registered or
proposed chemical uses of those crops, and regional management practices.
For a specific field-scale area, similar data would be needed to support an
assessment; however, greater detail may be necessary.
The identification of exposure pathways involves a qualitative (or
semiquantitative) assessment of how the chemical is thought to move from the
source to the exposed population. Important fate processes which may serve
to reduce the concentration of the chemical(s) along various pathways in
different environmental media are also identified. For the case; of
groundwater exposure, important pathways and processes are predefined to a
large extent in the models to be used.
The quantification of concentrations in a medium, given the source strength,
pathways, and attenuation mechanisms along each pathway, is the next step,
and is the major benefit of using models such as RUSTIC. The guidelines are
very specific in the requirement that concentrations be characterized by
duration and frequency as well as magnitude. These characteristics can be
determined through the analysis of time series exposure data generated by
the model.
The model produces a time series of toxicant concentrations in groundwater
such as appears in Figure 1.2. The time series can be compared to a
critical value .of the concentration y. (This might be, for instance, the
ADI [average daily intake level].) This type of analysis easily shows
whether the criterion is exceeded and gives a qualitative feel for the
severity of the exceedance state. A frequency distribution of the values of
y (Figure 1.3) can be created by determining how often y is at a particular
level or within a specific range. By choosing any value of y in Figure 1.2
and determining the area under the curve to the right of that va.lue, Figure
1.4, which is a cumulative frequency distribution of the toxicant
-------
o
O
u
o
O
Time (t)
Figure 1.2. Time series plot of toxicant concentrations
c
0)
._>
"5
*•
CD
(0
ft
c
0
£
O)
«
•a
o
o
u
K
O
Concentration (y)
Concentration (y)
Figure 1.3.
Frequency distribution of
toxicant concentrations.
Figure 1.4.
Cumulative frequency
distribution of toxicant
concentrations.
-------
concentration can be plotted. The cumulative frequency distribution shows
the chance that any given value y will be exceeded. If the example time
series is long enough, then the "chance" approaches the true "probability"
that y will be exceeded.
Thus far, only the concentration to which the organism will be exposed has
been discussed and nothing has been said concerning the duration of the
event. If a window of length "t" is imposed on the time series at level y
(Figure 1.5), and moved incrementally forward in time, a statement can be
made concerning the toxicant concentration within the duration window.
Normally, the average concentration within the window is used. The
resulting cumulative frequency distribution shows the chance that the moving
average of duration tc will exceed the critical value of y, y The moving
average window should be the same length as that specified for y . For
instance, in the case of cancer risk, a 70-year (lifetime) window is
normally used to average the data in the simulated time series. The use of
the moving window for averaging the time series allows us to compare both
the concentration and duration against the standard. The chance or
probability that the moving average concentration exceeds the standard is
the essence of the exposure assessment. This type of information provides a
precursor to the estimates of risk taken in using this chemical under the
conditions of the model simulation. The use of models like RUSTIC which
provide data of environmental concentrations and probability of occurrence
ends here.
The next step in exposure assessment involves the characterization of the
exposed population. Such factors as habits, age, sex, and location with
respect to the source are of importance. The integration of concentration
estimates and population characteristics makes possible the counting of the
conditional events of concentration in an environmental medium and the
opportunity for the population to be exposed to these concentrations. The
exposure assessment ends at this point. The actual intake of chemicals,
their fate within the human body (e.g., pharmacokinetics, toxicology), and
their effects on the exposed population are not considered. These, however,
are elements of the risk assessment.
Although the concepts underlying an exposure assessment are relatively
simple, the actual application of these concepts is complicated because of
large variations in source-specific and environment-specific characteristics
and the necessity to integrate specialized knowledge from a number of
different fields. This variability underscores the need to use a model such
as RUSTIC in the evaluation of exposure concentrations.
1.3 OVERVIEW OF RUSTIC
This section gives an overview of the RUSTIC model highlighting the features
and limitations of the simulation package as a whole and the component
models PRZM, VADOFT, and SAFTMOD. The RUSTIC code was designed to provide
state-of-the-art deterministic simulation of the fate and transport of
pesticides, applied for agricultural purposes, in the crop root zone, the
vadose zone, and the saturated zone. The model is capable of simulating
multiple pesticides or parent/daughter relationships. The model is also
-------
c
o
0>
u
c
o
u
Time (t)
Figure 1.5.
Time series of toxicant concentrations with
moving average window of duration tc•
capable of estimating probabilities of concentrations or fluxes in or
fromthese various media for the purpose of performing exposure assessments.
To avoid writing an entirely new computer code, it was decided to make use
of existing codes and software to the extent possible. Thus, due to its
comprehensive treatment of important processes, its dynamic nature, and its
widespread use and acceptability to the Agency and the agricultural chemical
industry, the Pesticide Root Zone model (PRZM) (Carsel et al., 1984) was
selected to simulate the crop root zone.
Having selected PRZM, several options were evaluated for developing the
RUSTIC linked model to meet the objectives stated in Section 1.1. The first
involved use of PRZM and a saturated zone model only. In this
configuration, PRZM would be used to simulate both the root zone and the
vadose zone. This option was eliminated because the assumptions of the
elementary soil hydraulics in PRZM (i.e., drainage of the entire soil column
to field capacity in one day) were considered inadequate for simulating flow
in a thick vadose zone. The second involved the use of PRZM and a 2-D or 3-
D variably saturated flow model. In this configuration, the latter model
would represent both the vadose and saturated zones. Although the interface
between the unsaturated and saturated zones would be more rigorously
simulated, it was felt that the resulting code would be computationally too
intensive for the intended application scenarios, and especially for use in
a Monte Carlo (probabilistic) mode. Therefore, this option was also
dropped.
-------
The option finally selected is depicted in Figure 1.6. In this
configuration, an enhanced version of PRZM is linked to a two-dimensional
(single or two-layer) saturated zone model either directly or through a one-
dimensional vadose zone flow and transport model. Both the vadose and
saturated zone models simulate water flow and solute transport. A new code
(VADOFT) was written to perform the flow and transport simulation in
the vadose zone, and an existing code was modified to produce the two-
dimensional X-Y, X-Z or axisymmetric simulation model for the saturated zone
(SAFTMOD).
PRZM
(1—D Flow and Transport)
VADOSE
ZONE
MODEL
m
VADOFT
(1—D Flow and Transport)
SATURATED ZONE MODEL
SAFTMOD
(2-D Flow and Transport)
Figure 1.6. Linked modeling system configuration,
-------
A significant problem associated with this type of linkage is the difference
in time scales of the individual models. While the vadose zone models are
required to operate on a daily or less-than-daily time step, the time step
of the saturated zone model could vary from days to months. A second
problem is the correct interfacing of the vadose and saturated zone models,
especially for the case of a fluctuating water table. This requires special
handling of the spatial discretizations in the two models. The solution to .
these linkage problems is discussed in detail in Section 5.
1.3.1 Overview of PRZM
1.3.1.1 Features--
The Pesticide Root Zone Model (PRZM) is a one-dimensional, dynamic,
compartmental model that can be used to simulate chemical movement in
unsaturated soil systems within and immediately below the plant root zone.
It has two major components: hydrology (and hydraulics) and chemical
transport. The hydrologic component for calculating runoff and erosion is
based on the Soil Conservation Service curve number technique and the
Universal Soil Loss Equation. Evapotranspiration is estimated either
directly from pan evaporation data, or based on an empirical formula.
Evapotranspiration is divided .among evaporation from crop interception,
evaporation from soil, and transpiration by the crop. Water movement is
simulated by the use of generalized soil parameters, including field
capacity, wilting point, and saturation water content. With a newly added
feature, irrigation may also be considered.
The chemical transport component can simulate pesticide application on the
soil or on the plant foliage. Dissolved, adsorbed, and vapor-phase
concentrations in the soil are estimated by simultaneously considering the
processes of pesticide uptake by plants, surface runoff, erosion, decay,
volatilization, foliar washoff, advection, dispersion, and retardation. Two
options are available to solve the transport equations: (1) the original
backwards-difference implicit scheme that may be affected by excessive
numerical dispersion at high Peclet numbers, or (2) the 'method of
characteristics' algorithm which eliminates numerical dispersion while
increasing model execution time.
Predictions are made on a daily basis. Output can be summarized for a
daily, monthly, or annual period. Daily time series values of various
fluxes or storages can be written to sequential files during program
execution for subsequent analysis. In addition, a 'Special Action' option
allows the user to output soil profile pesticide concentrations at user-
specified times. With the Special Action option, the user can also change
the values of certain parameters during the simulation period.
1.3.1.2 Limitations--
There were significant limitations in the original (Release I) version of
PRZM. A few were obvious to the developers, while others were pointed out
subsequently by model users. These are broken into four categories:
• Hydrology
• Soil hydraulics
• Method of solution of the transport equation
• Deterministic nature of the model
9
-------
The Release II version of PRZM has been suitably modified to overcome many
of these limitations.
Hydrologic and hydraulic computations are still performed in PRZM on a daily
time step even though for some of the processes involved (evaporation,
runoff, erosion) finer time step might be used to ensure greater accuracy
and realism. For instance, simulation of erosion by runoff depends upon the
peak runoff rate, which is in turn dependent upon the time base of the
runoff hydrograph. This depends to some extent upon the duration of the
precipitation event. PRZM retains its daily time step primarily due to the
relative availability of daily versus shorter time step meteorological data.
This limitation has in part been mitigated by enhanced parameter guidance.
In PRZM, Release I, the soil hydraulics were simple--all drainage to field
capacity water content was assumed to occur within 1 day. (An option to
make drainage time dependent was also included, but there is not much
evidence to suggest that it was utilized by model users to any great
extent.) This had the effect, especially in deeper soils, of inducing a
greater-than-anticipated movement of chemical through the profile. While
this representation of soil hydraulics has been retained in PRZM, the user
has the option of coupling PRZM to VADOFT. PRZM is then used to represent
the root zone, while VADOFT, with a more rigorous representation of
unsaturated flow, is used to simulate the thicker vadose zone. The code
VADOFT is discussed in more detail in a subsequent section. For short
distances from soil surface to the water table, PRZM can be used to
represent the entire vadose zone without invoking the use of VADOFT as long
as no layers which would restrict drainage are present.
The addition of algorithms to simulate volatilization has brought into focus
another limitation of the soil hydraulics representation. PRZM simulates
only advective, downward movement of water and does not account for
diffusive movement due to soil water gradients. This means that PRZM is
unable to simulate the upward movement of water in response to gradients
induced by evapotranspiration. This process has been identified by Jury et
al. (1984) as an important one for simulating the effects of volatilization.
However, the process would seem less likely to impact the movement of
chemicals with high vapor pressures. For these chemicals, vapor diffusion
would be a major process for renewing the chemical concentration in the
surface soil.
Another limitation of the Release I model was the apparent inadequacy of the
solution to the transport equation in advection-dominated systems. The
backward difference formulation of the advection term tends to produce a
high degree of numerical dispersion in such systems. This results in
overprediction of downward movement due to smearing of the peak and
subsequent overestimation of loadings to groundwater. In this new release,
a new formulation is available for advection-dominated systems. The
advective terms are decoupled from the rest of the transport equation and
solved separately using the Method of Characteristics (MOC). The remainder
of the transport equation is then solved as before, using the fully implicit
scheme. This approach effectively eliminates numerical dispersion with,
10
-------
however, an increase in the computation time. In low-advection systems, the
MOC approach reduces to the original PRZM solution scheme, which becomes
exact as velocities approach zero.
The final limitation is the use of field-averaged water and chemical
transport parameters to represent spatially heterogeneous soils. Several
researchers have shown that this approach produces slower breakthrough times
than are observed using stochastic approaches. This concern has been
addressed by adding the capability to run PRZM in a Monte Carlo framework.
Thus, distributional, rather than field-averaged, values can be utilized as
inputs which will produce distributional outputs of the relevant variables
(e.g., flux to the water table).
1.3.2 Overview of the Vadose Zone Flow and Transport Model (VADOFT)
VADOFT is a finite-element code for simulating moisture movement and solute
transport in the vadose zone. It is the second part of the three-component
RUSTIC model for predicting the movement of pesticides within and below the
plant root zone and assessing subsequent groundwater contamination. The
VADOFT code simulates one-dimensional, single-phase moisture and solute
transport in unconfined, variably saturated porous media. Transport
processes include hydrodynamic dispersion, advection, linear equilibrium
sorption, and first-order decay. The code predicts infiltration or recharge
rate and solute mass flux entering the saturated zone. Parent/daughter
chemical relationships may be simulated. The following description of
VADOFT is adapted from Huyakorn et al. (1988a).
1.3.2.1 Features--
The code employs the Galerkin finite-element technique to approximate the
governing equations for flow and transport and allows for a wide range of
nonlinear flow conditions. Boundary conditions of the variably saturated
flow problems may be specified in terms of prescribed pressure head or
prescribed volumetric water flux per unit area. Boundary conditions of the
solute transport problem may be specified in terms of prescribed
concentration or prescribed solute mass flux per unit area. All boundary
conditions may be time dependent. An important feature of the algorithm is
the use of constitutive relationships for soil water characteristic curves
based on soil texture.
1.3.2.2 Limitations- -
Major assumptions of the flow model are that the flow of the fluid phase is
one-dimensional, isothermal and governed by Darcy's law and that the fluid
is slightly compressible and homogeneous. Hysteresis effects in the
constitutive relationships of relative permeability versus water saturation,
and water saturation versus capillary pressure head, are assumed to be
negligible.
Major assumptions of the solute transport model are that advection and
dispersion are one-dimensional, and fluid properties are independent of
contaminant concentrations. Diffusive/dispersive transport in the porous-
medium system is governed by Pick's law. The hydrodynamic dispersion
coefficient is defined as the sum of the coefficients of mechanical
11
-------
dispersion and molecular diffusion. Adsorption and decay of the solute is
described by a linear equilibrium isotherm and a lumped first-order decay
constant. Steady-state transport can not be simulated when decay is
considered.
The code handles only single-phase flow (i.e., water) and ignores the
presence of a second phase--i.e., air. The code does not take into account
sorption nonlinearity or kinetic sorption effects which, in some instances,
can be important. The code considers only single-porosity (granular) soil
media. It does not simulate flow or transport in fractured porous media or
structured soils.
1.3.3 Overview of the Saturated Zone Flow and Transport Model (SAFTMOD)
SAFTMOD is a finite element code for simulating groundwater flow and solute
transport in the saturated zone. It is the third part of the three-
component RUSTIC model for predicting the movement of pesticides within and
below the plant root zone and assessing subsequent groundwater
contamination. SAFTMOD performs two-dimensional simulations in an areal
plane or a vertical cross section. In addition, the code can also perform
axisymmetric simulations. Both single (unconfined or confined) and leaky
two-aquifer systems can be handled. Transport of dissolved contaminants may
also be simulated within the same domain. Transport processes accounted for
include hydrodynamic dispersion, advection, linear equilibrium sorption, and
first-order decay. Parent/daughter chemical relationships may be simulated.
The following description of SAFTMOD is adapted from Huyakorn et al.
(1988b).
1.3.3.1 Features--
The two dimensional analyses can be done in a transient or steady-state mode
using SAFTMOD. The code employs the Galerkin finite-element technique to
approximate the governing equations for flow and transport. For groundwater
flow simulations the code accommodates water table conditions, recharge by
infiltration or precipitation, and well pumping or injection. Boundary
conditions of the saturated flow problem are specified in terms of
prescribed hydraulic head (defined as the sum of pressure head and
elevation), or prescribed volumetric water flux. Boundary conditions of the
solute transport problem are specified in terms of prescribed concentration
or prescribed solute mass flux. All boundary conditions can be time
dependent.
1.3.3.2 Limitations--
The SAFTMOD code contains both saturated flow and solute transport models.
Major assumptions of the flow model are that: Darcy's law is valid and
hydraulic head gradients are the only significant driving mechanism for
fluid flow, the porosity and hydraulic conductivity are constant with time,
and gradients of fluid density, viscosity, and temperature do not affect the
velocity distribution. For areal flow simulation in two-aquifer systems,
vertical leakage is assumed in aquitards.
Major assumptions of the transport model are that fluid properties are
independent of concentrations of contaminants. Contaminants are miscible
12
-------
with in-place fluids. For areal transport simulations, advection in
aquitards is assumed to be negligible. Diffusive/dispersive transport in
the porous medium system is governed by Pick's law. The hydrodynamic
dispersion coefficient is defined as the sum of the coefficients of
mechanical dispersion and molecular diffusion. Adsorption and decay of the
solute may be described by a linear equilibrium isotherm and a first-order
decay constant. Steady-state transport can not be simulated when decay is
considered.
Other limitations of the SAFTMOD code are that it simulates only single -
phase water flow and solute transport in saturated porous media. It does
not consider unsaturated or fractured media. Non-Darcy flow that may occur
near pumping wells is neglected. The code does not take into account
sorption nonlinearity or kinetic sorption effects which, in some instances,
can be important.
1.3.4 Model Linkage
One of the more challenging problems in this model development effort was
the temporal and spatial linkage of the component models. In the section
which follows, these linkages are discussed.
1.3.4.1 Temporal Model Linkage--
The resolution of the temporal aspects of the three models was
straightforward. PRZM runs on a daily time step. The time step in VADOFT
is dependent upon the properties of soils and the magnitude of the water
flux introduced at the top of the column. In order for the nonlinear
Richards' equation to converge, VADOFT may sometimes require time steps on
the order of minutes. The SAFTMOD time step, on the other hand, would
normally be much longer than one day.
For simplicity, it was decided that the time step of SAFTMOD would always be
an integer multiple of PRZM's daily time step. This makes the direct
linkage of PRZM to SAFTMOD in a temporal sense very straightforward. PRZM
is simply run for the number of days of the SAFTMOD time step, SAFTMOD is
then run, and so forth. The water and pesticide fluxes from PRZM are summed
and averaged over the length of the SAFTMOD time step. SAFTMOD receives the
time-averaged fluxes as input.
For the linkage of PRZM, through VADOFT, to SAFTMOD, the resolution of time
scales is also straightforward. VADOFT is prescribed to simulate to a
"marker" time value, specifically to the end of a day. The last
computational time step taken by VADOFT is adjusted so that it coincides
with the end of the day. PRZM's daily water fluxes are used as input to
VADOFT. VADOFT utilizes this flux as a constant over the day and adjusts
its internal computational time step in order to converge. PRZM and VADOFT
execute for the number of days corresponding to the SAFTMOD time step. The
output from VADOFT is then time averaged and used as input to SAFTMOD.
1.3.4.2 Spatial Linkages--
The spatial linkages utilized for the models are more complex. The
principal problem is the presence of a fluctuating water table, which
complicates the interfacing of the vadose (or root zone) and saturated zone
13
-------
codes. A second problem is that of the incompatibility between the
hydraulics in PRZM and VADOFT. Of course, any linking scheme utilized must
provide a realistic simulation of the flow of water and transport of solutes
at the interfaces arid must ensure mass balance.
PRZM and VADOFT--The major problem with the interfacing of these two models
is that while VADOFT solves the Richards' equation for water flow in a
variably saturated medium, PRZM uses simple "drainage rules" to move water
through the soil profile. Because of this incompatibility, there may be
times when PRZM produces too much water for VADOFT to accommodate within one
day. This is very likely to happen in agricultural soils, where subsoils
are typically of lower permeability than those of the root zone, which have
been tilled and perforated by plant roots and soil biota. The result of
this would be water ponded at the interface which would belong neither to
PRZM or VADOFT.
The solution was to prescribe the flux from PRZM into VADOFT so that VADOFT
accommodates all the water output by PRZM each day. This eliminates the
problem of ponding at the interface. However, it does force more water into
the vadose zone than might actually occur in a real system, given the same
set of soil properties and meteorological conditions. The consequence is
that water and solute are forced to move at higher velocities in the upper
portions of the vadose zone. If the vadose zone is deep, then this
condition probably has little impact on the solution. If it is shallow,
however, it could overestimate loadings to groundwater, especially if
chemical degradation rates are lower in the vadose zone than in the root
zone.
VADOFT and SAFTMOD--The principal problem here is that of the presence of a
dynamic boundary (rising and falling water table) between these models. Two
approaches were considered to handle the problem. The first was to expand
or contract the spatial domains of the models to accommodate the moving
boundary. This is not a particularly insurmountable problem for the vadose
zone model; however, for the saturated zone model it would require the
addition of nodes at the upper boundary. This would require a constant
evaluation and switching of the set of nodes receiving fluxes from the
vadose zone. This appeared undesirable.
The.second option was to overlap the spatial domains of the models and
interpolate values for fluxes based upon the position of the water table.
This latter approach was ultimately utilized. It has the additional feature
of eliminating the effects of the bottom boundary conditions prescribed for
the vadose zone model on the simulation of solute transport just above the
water table. A detailed discussion of these spatial linkages is given in
Section 5.
1.3.5 Monte Carlo Processor
RUSTIC can be run in a Monte Carlo mode so that probabilistic estimates of
pesticide loadings to the saturated zone or concentrations in a well
downgradient from the source area can be made. The input preprocessor
allows the user to select distributions for key parameters from a variety of
14
-------
distributions; the Johnson family (which includes the normal and lognormal),
uniform, exponential and empirical. If the user selects distributions from
the Johnson family, he may also specify correlations between the input
parameters. The Monte Carlo processor reads the standard deterministic
input data sets for each model, then reads a Monte Carlo input file which
specifies which parameters are to be allowed to vary, their distributions,
the distribution parameters, and correlation matrix. The model then
executes for a prespecified number of runs.
The output processor is capable of preparing statistics of the specified
output variables including mean, maximum values and quantiles of the output
distribution. The output processor also can tabulate cumulative frequency
histograms of the output variables and send them to a line printer for
plotting.
1.3.6 Overview Summary
A modeling system (RUSTIC) has been developed for the U.S. Environmental
Protection Agency which is capable of simulating the fate and transport of
pesticides, following application, through the crop root zone, vadose zone
and saturated zone and into drinking water wells. The model is envisioned
for use in simulating farms or small communities, and was designed to handle
a variety of geometries likely to be encountered in performing evaluations
for pesticide registration or special reviews. A major objective was also
to keep the model simple and efficient, yet allow use in a Monte Carlo mode
to generate probabilistic estimates of pesticide loadings or well water
concentrations. The model consists of three major computational modules;
PRZM, which performs fate and transport calculations for the crop root zone
and is capable of incorporating the effects of management practices; VADOFT,
which simulates one-dimensional flow and transport within the vadose zone;
and SAFTMOD, which simulates two dimensional flow and transport in the
saturated zone. SAFTMOD is capable of simulating a variety of aquifer
geometries (water table, confined and leaky two-aquifer systems) and can
perform areal, cross-sectional and axisymmetric simulations. Linkage of
these models is accomplished through the use of simple bridging algorithms
which conserve water and solute mass.
1.4 OVERVIEW OF VOLUMES I AND II
Documentation for RUSTIC has been produced in two volumes. The subject of
Volume I is model theory and code verification or testing. It contains a
description of the theory underlying the PRZM, VADOFT, and SAFTMOD codes.
The description of each code includes a brief overview highlighting the
features and limitations of each code. This is followed by detailed
descriptions of the algorithms involved in each code and how they are solved
numerically. The description of each of these models concludes with a
section on algorithm testing. The fifth section of Volume I contains a
description of the theory behind the linkage of the three codes to provide a
cohesive simulation of the movement of pesticides following application,
through the root zone, the vadose zone, and the saturated zone to a drinking
water well. Section 6 of Volume I covers the theory behind the uncertainty
preprocessor. Volume I concludes with model development references and
appendices.
15
-------
Volume II is a model user's guide. It opens with an introduction and a
section on the installation of the code on the target computer systems.
Section 3 has a user-directed overview of the model software, simulation
modules, and a description of data bases for simulation support and
parameter estimation. The fourth section takes the user through problem
definition and model setup, and module input sequence building.. Section 5
covers parameter estimation for the execution supervisor and each
computational module. Section 6 takes the user through an example problem.
Following the references (Section 7), Section 8 contains appendices which
include a listing of error messages and warnings, a variable glossary, and
sample input and output for the example problem of Section 6.
16
-------
SECTION 2
PESTICIDE ROOT ZONE MODEL
(PRZM RELEASE II)
2:i INTRODUCTION AND BACKGROUND
This section describes the theoretical background for a mathematical
simulation model (PRZM) that has been developed and partially tested to
evaluate pesticide leaching from the crop root zone under field crop
conditions.
Following this short introduction, Section 2.2 describes the features and
limitations of the model. A description of the theory, including a detailed
description of the equations solved, is provided in Section 2.3. An outline
of the numerical implementation techniques used by the model to apply the
theory to the simulation of physical problems follows. This section
concludes with a discussion of testing results for new algorithms which have
been added in this release.
2.1.1 Introduction
Pesticide leaching from agricultural fields as nonpoint source loads can lead
to groundwater contamination. Nonpoint source contamination is characterized
by highly variable loadings, with rainfall and irrigation events dominating
the timing and magnitude of the loading of pesticides leaching below the root
zone. The potentially widespread, areal nature of resulting contamination
makes remedial actions difficult because there is no single plume emanating
from a "point source" (the more common groundwater problem) that can be
isolated and controlled. In any case, a more prudent approach to prevention
or reduction of pesticide groundwater contamination must be based on
understanding the relationships among chemical properties, soil system
properties, and the climatic and agronomic variables that combine to induce
leaching. Knowledge of these relationships can allow a priori investigation
of conditions that lead to problems, and appropriate actions can be taken to
prevent widespread contamination.
Many investigators have studied the factors contributing to pesticide
leaching. These investigations have shown that chemical solubility in water,
sorptive properties, volatility, formulation, and soil persistence determine
the tendency of pesticides to leach through soil. Similarly, the important
environmental and agronomic factors include soil properties, climatic
conditions, crop type, and cropping practices. In short, the hydrologic
cycle interacts with the chemical characteristics to transform and transport
pesticides within and out of the root zone. Vertical movement out of the
17
-------
root zone can result in groundwater contamination and is the problem which
the model is designed to investigate.
Numerical models for the movement of solutes in porous media for steady-
state, transient, homogenous, and multi-layered conditions have been
previously developed. Included in such studies have been linear and
nonlinear sorption, ion exchange, and other chemical-specific reactions.
These investigations have proven valuable in interpreting laboratory data,
investigating basic transport processes, and identifying controlling factors
in transport and transformation. As noted in a recent review of models for
simulating the movement of contaminants through groundwater flow systems,
however, the successful use of such models requires a great deal of detailed
field data. This unfortunate conclusion arises from the scaling problems
associated with laboratory experiments and the traditional solution of the
appropriate partial differential equations at points or nodes in a finite-
difference or finite-element grid network. Each spatial segment modeled must
be properly characterized--a most expensive, if not impossible, task for many
modeling problems.
Such problems in modeling pesticide leaching with existing procedures are
discouraging when one considers the need to evaluate future problems arising
from pesticides not yet widely distributed or used. Models used to perform
such evaluations should conform to the maximum possible extent to known
theory, but must be structured to enable efficient analysis of field
situations with minimal requirements for specialized field data. In short,
the goal is to integrate the essential chemical-specific processes for
leaching with reasonable estimates of water movement through soil systems.
Data input requirements must be reasonable in spatial and temporal
requirements and generally available from existing data bases. This model
attempts to meet these objectives.
2.1.2 Background
The Pesticide Root Zone Model (PRZM) (Carsel et al., 1984; Carsel et al.,
1985) was selected as the code to provide the capability to simulate the fate
and transport of agriculturally applied pesticides in the crop root zone.
PRZM was initially designed for this purpose and has attained a degree of
acceptability in both the regulatory community and in the agricultural
chemical industry. Therefore, its utility in accomplishing the objective of
this model development effort is obvious.
2.2 FEATURES AND LIMITATIONS
2.2.1 Features
The Pesticide Root Zone Model (PRZM Release II) is a one-dimensional,
dynamic, compartmental model for use in simulating chemical movement in
unsaturated soil systems within and immediately below the plant root zone
(see Figure 2.1). PRZM allows the user to perform simulations of potentially
toxic chemicals, particularly pesticides, that are applied to the soil or to
plant foliage. Dynamic simulations allow the consideration of pulse loads,
the prediction of peak events, and the estimation of time-varying mass
.18
-------
:
Figure 2.1. Pesticide root zone model
19
-------
emission or concentration profiles, thus overcoming limitations of the more
commonly used steady-state models. Time-varying transport by both advection
and dispersion in the dissolved phase or diffusion in the gas phase are
represented in the program.
PRZM has two major components: hydrology and chemical transport. The
hydrologic component for calculating runoff and erosion is based on the Soil
Conservation Service curve number technique and the Universal Soil Loss
Equation. Evapotranspiration is estimated from pan evaporation data, or by
an empirical formula if input pan data are unavailable. Evapotranspiration
is divided among evaporation from crop interception, evaporation from soil,
and transpiration by the crop. Water movement is simulated by the use of
generalized soil parameters, including field capacity, wilting point, and
saturation water content. Irrigation may also be considered.
Dissolved, adsorbed, and vapor-phase concentrations in the soil are estimated
by simultaneously considering the processes of pesticide uptake by plants,
surface runoff, erosion, decay, volatilization, foliar washoff, advection,
dispersion, and retardation. The user may elect to solve the transport
equations using one of two finite-difference numerical solutions, the
original backwards-difference implicit scheme featured in the first release,
or a Method of Characteristics algorithm which greatly reduces numerical
dispersion, but increases model execution time.
The hydrologic components of pesticide transport equations (i.e., moisture
content and soil-water velocities) are decoupled, solved separately, and used
to numerically integrate the equation in succeeding time steps. Predictions
are made on a daily basis. Output can be summarized on a dailjr, monthly, or
annual period. A daily time series value for various fluxes or storages can
be written to sequential files during program execution.
2.2.2 Limitations
There were severe limitations of the PRZM Release I Code, some which were
obvious to the developers and some which were pointed out subsequently by
model users. These are broken into four categories:
• Hydrology
• Soil hydraulics
• Method of solution of the transport equation
• Deterministic nature of the model
In Release II, many of these limitations have, to an extent, been overcome.
Hydraulic computations are performed in PRZM on a daily time step; however,
some of the processes involved (evaporation, runoff, erosion) are clearly
among those which might be simulated on a finer time step to ensure greater
accuracy and realism. For instance, simulation of erosion by runoff depends
upon the peak runoff rate, which is in turn dependent upon the time base of
the runoff hydrograph. This depends to some extent upon the duration of the
precipitation event. PRZM retains its daily time step in this release
primarily due to the relative availability of daily versus shorter time step
20
-------
meteorological data. Hopefully, a portion of this limitation has been
mitigated by enhanced parameter guidance.
The method of computing potential evapotranspiration using Hamon's formula,
in the absence of some evaporation data, has also been retained.
Evapotranspiration from irrigated citrus in Florida was found to be
substantially underpredicted when using this method to estimate potential
evapotranspiration (Dean and Atwood 1985). Users should carefully check the
model's hydrologic simulation when using this option.
The capability to simulate soil temperature has been added to this release of
PRZM in order to correct Henry's constant for the temperature occurring in
various depths in the soil when performing vapor-phase calculations. Removal
of water by evaporation versus transpiration from the profile may have a
pronounced effect on soil temperature. This is due to the fact that more
heat is removed during the process of evaporation because the energy
necessary to vaporize water leaves the system, producing a cooling effect.
No differentiation is made between evaporation and transpiration in PRZM at
this time.
In PRZM, Release I, the soil hydraulics were simple--all drainage to field
capacity water content was assumed to occur within 1 day. (An option to make
drainage time dependent was also included, but there is not much evidence to
suggest that it was utilized by model users to any great extent). This had
the effect, especially in larger soil cores, of inducing a greater-than-
anticipated movement of chemical through the profile. While this
representation of soil hydraulics has been retained in PRZM, the user has the
option, with the linked modeling system, of coupling PRZM to VADOFT. PRZM is
then used to represent the root zone, while VADOFT, with a more rigorous
representation of unsaturated flow, is used to simulate the thicker vadose
zone. The difficulties in parameterizing the Richards equation for
unsaturated flow in VADOFT is overcome by using the technique of van
Genuchten to generate soil water characteristic curves using soil textural
information. For short soil cores, PRZM can obviously be used to represent
the entire vadose zone.
The addition of algorithms to simulate volatilization has brought into focus
another limitation of the soil hydraulics representation. PRZM simulates
only advective, downward movement of water and does not account for diffusive
movement due to soil water gradients. This means that PRZM is unable to
simulate the upward movement of water in response, to gradients induced by
evapotranspiration. This process has been identified by Jury et al. (1984)
as an important one for simulating the effects of volatilization. However,
the process would seem less likely to impact the movement of chemicals with
high vapor pressures. For these chemicals, vapor diffusion would be a major
process for renewing the chemical concentration in the surface soil.
Another limitation of the Release I model was the inadequacy of the solution
to the transport equation in advection-dominated systems. The backward
difference formulation of the advection term tends to produce a high degree
of numerical dispersion in such systems. This results in overprediction of
downward movement due to smearing of the peak and subsequent overestimation
21
-------
of loadings to groundwater. In this new release, a new formulation is
available for advection-dominated systems. The advective terms are decoupled
from the rest of the transport equation and solved separately using a Method
of Characteristics (MOC) formulation. The remainder of the transport
equation is then solved as before, using the fully implicit scheme. This
approach effectively eliminates numerical dispersion, but with some
additional overhead expense in computation time. In low-advection systems,
the MOC approach reduces to the original PRZM solution scheme, which is exact
for velocities approaching zero.
The final limitation is the use of field-averaged water and chemical
transport parameters to represent spatially heterogeneous soils. Several
researchers have shown that this approach produces slower breakthrough times
than are observed using stochastic approaches. This concern has been
addressed by adding the capability to run PRZM in a Monte Carlo framework.
Thus, distributional, rather than field-averaged, values can be utilized as
inputs which will produce distributional outputs of the relevant variables
(e.g., flux to the water table).
2.3 DESCRIPTION OF THE EQUATIONS
The mathematical description of the processes simulated by PRZM are broken
down in the following discussion into five categories:
• Transport in Soil
• Water Movement
• Soil Erosion
• Volatilization
• Irrigation
The first three categories were simulation options previously available in
PRZM Release I. Since the capability to simulate ponding is new, the
mathematical basis of the ponding algorithms is described in detail. The
final process, volatilization, was not available in previous releases of
PRZM, and its theoretical basis is also described in detail.
2.3.1 Transport in Soil
The PRZM model was derived from the conceptual, compartmentalized
representation of the soil profile as shown in Figure 2.2. From
consideration of Figure 2.2, it is possible to write mass balance equations
for both the surface zone and the subsurface zones. Addition of the vapor
phase and ponded' water compartments in PRZM require the consideration of
additional terms. The surface zone expressions for each of the dissolved,
adsorbed, and vapor phases can be written as:
A Az 3(CW0)
JD " JV " JDW " JU " JQR + JAPP + JFOF ± JTRN (2-1)
22
-------
(JApp, JFOF)
(Surface Lc
Erosion)
V.OU.L j_cn;t
Runoff)
JER
SOLI
C5
PI
Adsorj
Desor]
-------
A Az a(Csps)
" JDS " JER (2-2)
at
A Az 5(Cga)
= JGD " JDG (2-3)
at
where
o
A - cross-sectional area of soil column (cm )
Az - depth dimension of compartment (cm)
GW — dissolved concentration of pesticide (g cm" )
Cs — sorbed concentration of pesticide (g g )
C - gaseous concentration of pesticide (g cm" )
00
6 - volumetric water content of soil (cm cm"J)
o o
a - volumetric air content of the soil (cm cm" )
o
ps - soil bulk density (g cm"0)
t - time (d)
JQ - represents the effect of dispersion and diffusion of dissolved
phase (g day)
Jv - represents the effect of advection of dissolved phase
(g day"1)
JQQ = represents the effect of dispersion and diffusion in vapor
phase (g day" )
Jjjy = mass loss due to degradation in the dissolved phase (g day" )
JJJQ — mass loss due to degradation in the vapor phase (g day )
Jjj = mass loss by plant uptake of dissolved phase (g day )
JQP — mass loss by removal in runoff (g day )
= mass ga*-n due to pesticide deposition on the soil surface
(g day"1)
= mass gain due to washoff from plants to soil (g day" )
= mass loss due to degradation of sorbed phase chemical
(g day'1)
= mass loss by removal on eroded sediments (g day )
= mass gai-n or loss due to parent/daughter transformations
24
-------
Equations for the subsurface zones are identical to Equations (2-1), (2-2),
and (2-3) except that JQR, JpoF' anci JER are not included. J^PP aPPlies to
subsurface zones only when pesticides are incorporated into the soil. For
subsurface layers below the root zone, the term J^ is also not utilized.
Note that terms representing phase transfers (e.g., volatilization) are
neglected in Equations (2-1) through (2-3) because they cancel when the
equations are added (see Equation (2-19) below).
Each term in Equations (2-1) through (2-3) are now further defined.
Dispersion and diffusion in the dissolved phase are combined and are
described using Pick's law as
A Az Dw d2(Cw0)
JD 5 (2'
dz2
where
DW = diffusion-dispersion coefficient for the dissolved phase,
f\ I
assumed constant (cm day )
O
GW = dissolved concentration of pesticide (g cm" )
o o
8 = volumetric soil water content (cm cm)
x — soil depth dimension (cm)
In a similar manner, dispersion and diffusion in the vapor phase are
described by Pick's law as
A Az D 32(C a)
JGD --- — - —
dz
where
D = Molecular diffusivity of the pesticide in the air-filled pore
p _ -I
space (cm day" )
Q
C = vapor-phase concentration of pesticide (g cm" )
&
o o
a = volumetric air content (cm cm" )
The dependence of the molecular diffusivity of the pesticide in air-filled
pore space on the volumetric air content is described by the Millington-Quirk
expression (Jury et al., 1983a)
25
-------
where
o o
a = the air-filled porosity (cm cm" )
o o
<£ = total porosity (cm cm~J)
Dfl = molecular diffusivity of the chemical in air, assumed constant
9 1
(cm day )
The mathematical theory underlying the diffusive and dispersive flux of
pesticide in the vapor phase within the soil and into the overlying air can
be found in the section describing volatilization.
The advective term for the dissolved phase, Jy, describes the movement of
pesticide in the bulk flow field and is written as
A Az V 3(CW0)
dz
where
V = velocity of water movement (cm day )
Vapor-phase advection has not been included as a flux in the transport
equation. A number of researchers have indicated a consensus that vapor -
phase advection is not likely to be significant for agricultural situations
(W. Jury, W. Spencer, W. Farmer, L. Thibodeaux - personal communications,
1987) . Early studies of water vapor movement suggested that the fluctuation
of barometric pressure at the soil surface could act as a pumping mechanism
for vapor-phase advective transport (Fukuda 1955; Farrell et al . , 1966;
Scotter and Raats 1970). However, using models for vapor emissions from
landfills, Thibodeaux et al. (1982) found that atmospheric pressure
fluctuations increased the total emission rate for benzene by only 15%,
compared to constant pressure conditions. Therefore, it appears to be a
reasonable assumption at this time to neglect vapor-phase advection in
modeling chemical migration for agricultural situations.
Degradation of a pesticide in or on soil may be due to such processes as
hydrolysis, photolysis, and microbial decay. If these processes follow
pseudo first-order kinetics, the rate coefficients may be combined into a
single decay coefficient. Assuming the same rate constants for the solid and
dissolved phases, the rate of change of chemical out of each phase due to
decomposition may be written as:
JDW - Ks Cw e A Az
26
-------
JDS - Ks Cs »s A Az <2-
JDG = Kg Cg a A Az (2'
where
K0 = lumped, first-order decay constant for solid and dissolved
s
phases (day )
K = lumped, first-order decay constant for vapor phase (day )
s i
C = solid-phase concentration of pesticide (g g )
Plant uptake of pesticides is modeled by assuming that uptake of a pesticide
by a plant is directly related to transpiration rate. The uptake is given
by:
Jn = f Cw 6 e A Az (2-10)
where
Ju - uptake of pesticide (g day )
f = the fraction of total water in the zone used for transpiration
(day'1)
e -an uptake efficiency factor or reflectance coefficient
(dimensionless)
Erosion and runoff losses as well as inputs to the surface zone from foliar
washoff are considered in the surface layer. The loss of pesticide due to
runoff is
JQR - A
in which
^OR = pesticide loss due to runoff (g day )
O 1
Q — the daily runoff volume (cm day )
n
A^ = watershed area (cm )
and the loss of pesticide due to erosion is
27
-------
JER
Xe rom Cs A
where
^ER " tne pesticide loss due to erosion (g day )
Xe = the erosion sediment loss (metric tons day )
rQm = the enrichment ratio for organic matter (g g )
p - a units conversion factor (g tons )
Soil erosion is discussed in more detail in Section 2.3.3.
Pesticides can be applied to either bare soil if pre-plant conditions prevail
or to a full or developing crop canopy if post-plant treatments are desired.
The pesticide application is an input mass rate which is calculated by one of
the application/deposition models discussed in Section 3.1. It is
partitioned between the plant canopy and the soil surface, and the rate at
which it reaches the soil surface is designated JAPP-
Pesticides applied to the plant canopy can be transported to the soil surface
as a result of rainfall washoff. This term, J „_..,, is defined as:
rur
JFOF - E pr M A <2-13>
where
E - foliar extraction coefficient (cm )
P = daily rainfall depth (cm day )
M = mass of the pesticide on the plant surface projected area basis
(g cm"2)
The foliar pesticide mass, M, is further subject to degradation and losses
through volatilization. Its rate of change is given by
- - KfMA - JFOF + AF b A (2-14)
where
= lumped first-order foliar degradation constant (day )
= application rate to the plant (g ha day )
= a units conversion factor (ha)
28
-------
Adsorption and desorption in Equations (2-1) through (2-3) are treated as
instantaneous, linear, and reversible processes. Using this assumption, the
sorbed phase concentration can be related to the dissolved-phase
concentration by:
Cs - Kd Cw (2-15)
where
KJ - partition coefficient between the dissolved and solid phases
(cm3 g'1)
A similar expression can be developed to express the vapor phase
concentration in terms of the dissolved-phase concentration as follows
Cg - KHCW (2-16)
where
Ku = Henry's constant, i.e., distribution-coefficient between liquid
phase and vapor phase (cm cm" )
The transformation of parent to daughter is assumed to be first order and
takes place according to
JTRN = - KTRN Cw A Az * (2-17)
where
&J.RN - the transformation rate constant (day )
When simulating an end-of-chain daughter, J-jRN may a^-so be a source term
equal to the sum of the first-order transfers from any and all parents.
JTRN - 2 KTRN cw A Az e <2
in which the superscript k denotes a parent compound. For intermediate
products, the solute transport equation may contain terms such as those shown
in both Equations (2-17) and (2-18). The transformation of parent to
daughter compounds is discussed in detail in Section 2.5.4. The section
includes a description of the equations used to simulate this process.
29
-------
Summing Equations (2-1), (2-2), and (2-3) and utilizing (2-15) and (2-16),
the following expression results for the mass balance of pesticide in the
uppermost soil layer:
K
at
aKR)]
w
3(aCwKH) 3CW<
- c.
w
Kdps) + KgaKH
JApp EPrM
- K
A Az Az
TRN
AwAz
KTRN Cw
PXeromKd
Equation (2-19) is solved in PRZM for the surface layer with f.6 =0, and an
upper boundary condition which allows vapor phase flux upward from the soil
surface to the overlying air. This upper boundary condition is described
more fully in the section on volatilization. The lower boundary condition is
one which allows advection, but no diffusion, out of the bottom of the soil
profile.
2.3.2 Water Movement
Because V and 8 are not generally known and not generally measured as part of
routine monitoring programs, it is necessary to develop additional equations
for these variables. In the general case, Darcy's law can be combined with
the continuity equation to yield the Richards equation (Richards 1931) :
where
K(0) - hydraulic conductivity at various heads (cm sec '
9 - soil water content (cnr cm" )
and
V = - K(0)
(2-21)
or, in simpler terms
at
(2-22)
30
-------
where
•3 O
d = soil water content (cm cm" )
V = soil water velocity (cm day )
Writing Equation (2-22) in an integrated backwards finite difference form
yields
Az (0t+- 0fc) = (V£ - V^) At (2-23)
or
6t+1 Az - (Vi - V1.1)At + 0cAz (2-24)
In these equations, t and t+1 denote the beginning and end of time step
values, respectively, and i is the soil layer index. These equations can be
further simplified by substituting the nomenclature SW for 6hz so that
SWt+1 - SWfc + (V - V) At (2-25)
where
SW = soil water content (cm)
The velocities in Equation (2-25) are a function of inputs to the soil
(precipitation, infiltration) and outflows from the soil (evapotranspiration,
runoff) .
Water balance equations are separately developed for (a) the surface zone,
(b) horizons comprising the active root zones, and (c) the remaining lower
horizons within the unsaturated zone. The equations are:
Surface Zone
(SW)£+1 = (SW){ -f INF - !]_ - Ej_ - Ux (2-26)
Root Zone
+ I^-L - Vj_ - I± (2-27)
Below Root Zone
= (SW)J + I^ - It (2-28)
31
-------
where
(SW)£ = soil water in layer "i" on day "t" (cm)
E^ = evaporation (cm day )
U^ = transpiration (cm day )
1^ = percolation out of zone i (cm day )
INF = infiltration into layer 1 (cm day )
Daily updating of soil moisture in the soil profile using the above equations
requires the additional calculations for infiltration, evaporation,
transpiration, and percolation.
Infiltration is calculated as
INF = P + SM - Q - E (2-29)
f\
where, assuming a unit area of 1 cm ,
P = precipitation as rainfall, minus crop interception (cm day )
SM = snowmelt (cm day )
Q = runoff (cm day )
E = evaporation (cm day )
The calculations of precipitation, snowmelt, and runoff on a daily time step
are described below. The disaggregation of these values and the calculation
of the change in the depth of ponding on a finer time step is included in
Sections 2.3.5.4 and 2.4.4 describing the simulation of furrow irrigation and
ponded surface water.
Input precipitation is read in and pan evaporation and/or air temperature are
inputs from which potential evapotranspiration (PET) is estimated. Incoming
precipitation is first partitioned between snow or rain, depending upon
temperature. Air temperatures below 0°C produce snow and may result in the
accumulation of a snowpack. Precipitation first encounters the plant canopy
and once the interception storage is depleted, the remaining depth is
available for the runoff or infiltration.
The runoff calculation partitions the precipitation between infiltrating
water and surface runoff. Infiltrating water may be ponded on the soil
surface for a period of time before it infiltrates, but this ephemeral
process is described in a following section. Runoff is calculated by a
modification of the USDA Soil Conservation Service curve number approach
(Haith et al. 1979). Snowmelt is estimated on days in which a snowpack
exists and above freezing temperatures occur as
SM = CMT (2-30)
32
-------
where
Cj. = degree-day snowmelt factor (cm °C day" )
T = average daily temperature (°C)
The precipitation and/or snowmelt are inputs to the SCS runoff equation
written as
(P+ SM - 0.2S)2
^ P + SM + 0.8S U ;
where S, the watershed retention parameter, is estimated by
S =* 1000/RCN - 10 (2-32)
where
RCN = SCS runoff curve number
Curve numbers are a function of soil type, soil drainage properties, crop
type, and management practice. Typically, specific curve numbers for a given
rainfall event are determined by the sum of the rainfall totals for the
previous 5 days, known as the 5-day antecedent moisture condition. In this
release of PRZM, as in the original version, the curve numbers are
continuously adjusted each day as a function of the soil water status in the
upper soil layers. These algorithms were developed and reported by Haith and
Loehr (1979).
The daily evapotranspiration demand is divided among evaporation from canopy,
ponded surface water, soil evaporation, and crop transpiration. Total demand
is first estimated and then extracted sequentially from crop canopy storage,
ponded surface water, and then from each layer until wilting point is reached
in each layer or until total demand is met. Evaporation occurs down to a
user-specified depth. The remaining demand, crop transpiration, is met from
the active root zone. The root zone growth function is activated at crop
emergence and increases stepwise until maximum rooting depth is achieved at
crop maturity.
Actual evapotranspiration from a soil layer is estimated as:
i-1
E1± - MIN [
-------
f^ - depth factor for layer 'i'
WP^ - wilting point water content in layer 'i' (cm)
ET - potential evapotranspiration (cm)
This equation states that the transpiration from any layer 'i' is the minimum
of the available water in layer 'i' or the demand remaining after extracting
available water from layers above 'i' in the profile.
The depth factor, f^, is internally set in the code. It linearly weights
the extraction of ET from the root zone with depth. A triangular root
distribution is assumed from the surface zone to the maximum depth of
rooting, with the maximum root density assumed to be near the surface. This
algorithm essentially views the plant as a pump and assumes that it will
expend the minimum energy possible in pumping. As long as the soil water is
equally available, water closest to the surface meets this criterion.
Evapotranspiration may also be limited by soil moisture availability. The
potential rate may not be met if sufficient soil water is not available to
meet the demand. In that case, PRZM modifies the potential rate by the
following equations:
ET = ET ; if SW > 0.6 FC (2-34)
ETp = SMFAC ETp; if WP < SW < 0.6 FC
ET = 0; if SW < WP
where
FC = soil moisture content at field capacity (cm)
WP = soil moisture content at wilting point (cm)
SMFAC - soil moisture factor
The SMFAC concept has been used in other similar water balance models (Haith
et al. 1979; Stewart et al. 1976) and is internally set in the code to
linearly reduce ET when soil water becomes limited. Finally, if pan
evaporation input data are available, ET is related to the input values as
ETp = Cp PE (2-35)
where
PE = pan evaporation (cm day )
C = pan factor (dimensionless)
34
-------
The pan factor is constant for a given location and is a function of the
average daily relative humidity, average daily wind speed, and location of
the pan with respect to an actively transpiring crop.
In the absence of pan evaporation data, ET is estimated by
ET = 14000 Ld2 (SVD) (2-36)
where
LJ = possible hours of sunshine per day, in 12-hour units
SVD = saturated vapor density at the mean air temperature (g cm" )
SVD = 0.622 SVP/(Rg Tabg)
where
SVP = saturated vapor pressure at the mean absolute air
temperature (mb)
R = dry-air gas constant
&
= absolute mean air temperature (°K)
The final term in the water balance equations that must be defined is the
percolation value, I. Because the Richards equation is not solved in PRZM
utilizing soil water characteristic curves to predict water movement, PRZM
resorts to "drainage rules" keyed to soil moisture storages and the time
available for drainage. Two options are included. Although these options
are admittedly simplistic representations of soil moisture redistribution,
they are consistent with the objectives of PRZM and its intended uses.
2.3.2.1 Option 1--
Percolation, I, in this option is defined in the context of two bulk soil
moisture holding characteristics commonly reported for agricultural soils:
field capacity and wilting point. Field capacity is a somewhat imprecise
measure of soil water holding properties and is usually reported as the
moisture content that field soils attain after all excess water is drained
from the system under influence of gravity, usually at tensions of about 0.3
bar. The difficulty with this concept is the fact that some soils will
continue to drain for long periods of time, and thus field capacity is not a
constant. Admitting the lack of theoretical and physical rigor, the concept
remains as a useful measure of soil moisture capacity and has been
successfully used in a number of water balance models (Haith et al., 1979;
Stewart et al., 1976). Wilting point is a function of both the soil and
plants growing in the soil. It is defined as the soil moisture content below
which plants are unable to extract water, usually at tensions of about 15
bar.
Field capacity and wilting point are used operationally to define two
reference states in each soil layer for predicting percolation. If the soil
35
-------
water, SW, is calculated to be in excess of field capacity, then percolation
is allowed to remove the excess water to a lower zone. The entire soil
profile excess is assumed to drain within 1 day. The lower limit of soil
water permitted is the wilting point. One outcome of these assumed "drainage
rules" is that the soil layers below the root zone tend to quickly reach
field capacity and remain at that value. When this condition is reached, all
water percolated below the root zone will displace the water within the lower
soil layer simulated, and so on. There is no allowance for lateral water
movement. Water balance accounting in this manner should be most accurate
for sandy soils in which water movement is relatively unimpeded and is least
accurate for clay soils (Stewart et al., 1976).
2.3.2.2 Option 2--
The second option is provided to accommodate soils having low permeability
layers that restrict the "free drainage" assumed in Option 1. In the context
of the field capacity reference condition, two things may occur. First,
conditions may prevail that raise the soil moisture levels above field
capacity for periods of time because the water is "backed up" above a
relatively impermeable layer. Second, the excess water may not drain during
the 1-day period assumed in Option 1. To accommodate these conditions, two
additional parameters are needed. Maximum soil moisture storage, 6 is added
to represent moisture contents under saturated conditions. The drainage rate
also must be modified to allow drainage to field capacity over periods in
excess of 1 day (one time step). The drainage rate is assumed to be a first-
order function of the water content above field capacity and is modeled by
d (8 - 0fc)
^ a (6 - *fc) (2-37)
which has the solution
exp(.aAt) + 0 (2-38)
where
o o
6 = soil layer water content (cm cm)
Q O
9fc= water content at field capacity (cm cm)
a - drainage rate parameter (day"1)
In this equation, t and t+1 denote beginning and end of time step values,
respectively, and i is the soil layer index. The value t* denotes a value of
time between beginning and end of time step. The variable 0? here denotes
current storage plus any percolation from the next layer above, before the
occurrence of any drainage from the current layer. Because Equation (2-38)
is solved independently for each layer in the profile, there is a possibility
of exceeding the storage capability (saturation water content, #s) of a low-
36
-------
permeability layer in the profile if a more permeable layer overlies it. At
each time step, once redistribution is complete, the model searches the
profile for any 6^ > 9S.. If this condition is found, the model
redistributes water back into overlying layers, as if the percolation of
additional water beyond that necessary to saturate the low-permeability layer
had not occurred. This adjustment is necessary due to the nature of Equation
(2-38) and the fact that these equations for each layer are not easily
coupled. The difficulty in coupling the equations for the entire profile
arises from the dichotomy that one of two factors limits percolation from a
stratum in the profile: either the rate at which that stratum can transmit
water, or the ability of the stratum below it to store or transmit water.
This dichotomy leads to an iterative (or at least corrective) approach to the
explicit solution of a system of equations for 8^ represented by Equation (2-
38). It should be noted, however, that the value of a selected by this
approach is only relevant if the permeability of the soil materials, and not
storage considerations in the profile (i.e., the presence of a water table),
is the limiting factor for percolation of water.
2.3.3 Soil Erosion
Removal of sorbed pesticides on eroded sediments requires estimates for soil
erosion. The Modified Universal Soil Loss Equation (MUSLE) as developed by
Williams (1975) is used to calculate soil loss:
Xe - a (Vr qp)°-56K LS C P (2-39)
where
Xg - the event soil loss (metric tons day )
Vr - volume of event (daily) runoff (m )
•3 1
qp - peak storm runoff (m sec'^)
K - soil erodability factor
LS = length-slope factor
C - soil cover factor
P - conservation practice factor
a = units conversion factor
Most of the parameters in Equation (2-39) are easily determined from other
calculations within PRZM (e.g., Vf), and others are familiar terms readily
available from handbooks. However, the peak storm runoff value, q , can vary
widely depending upon rainfall and runoff characteristics. A trapezoidal
hydrograph is assumed in PRZM. From the assumed hydrograph shape and the
storm duration, a peak runoff rate is calculated.
The enrichment ratio, rom, is the remaining term which needs to be defined to
estimate the removal of sorbed pesticides by erosion. Because erosion is a.
37
-------
selective process during runoff events, eroded sediments become "enriched" in
smaller particles. The sediment transport theory available to describe this
process requires substantially more hydraulic spatial and temporal resolution
than used in PRZM, leading to the adoption of an empirical approach (Mockus
1972) . The enrichment ratio for organic matter is calculated from
ln(rom) =2 + 0.2 In^/A^ (2-40)
2.3.4 Volatilization
As volatilization was not available in the previous releases of PRZM, its
theoretical basis is discussed in detail here. The following key processes
have been identified as important in volatilization algorithms to simulate
vapor-phase pesticide transport within the soil/plant compartments:
• Vapor-phase movement of the pesticide in the soil profile
• Boundary layer transfer at the soil -air interface
• Vertical diffusion of pesticide vapor within the plant canopy
• Pesticide mass transfer between the plant (leaves) and the
surrounding atmosphere
• Soil temperature effects on pesticide volatilization
The discussion of the volatilization algorithms is presented in four parts:
influence of vapor phase pesticide in soil and volatilization flux,
volatilization flux through the plant canopy, volatilization flux from plant
surfaces, and soil temperature modeling and effects. Figure 2,3 is a
schematic of the pesticide vapor and volatilization processes considered in
soil and plant compartments.
2.3.4.1 Soil Vapor Phase and Volatilization Flux--
The governing equations for chemical transport in the vapor phase were
introduced previously in the description of transport in the soil. Fluxes
from the soil column in the vapor phase are summarized in that discussion by
Equations (2-3), (2-5), and (2-9). The terms in these equations are summed
with the other flux terms to produce the transport equation (2-19). In
addition to these enhancements, the upper boundary of PRZM was changed from a
zero -concentration boundary to a stagnant -layer boundary to allow diffusive
transport upward from the soil to the overlying atmosphere. This enhancement
is discussed in detail below.
Surface boundary condition- -When a pesticide is incorporated into the soil,
the initial volatilization rate is a function of the vapor pressure of the
chemical at the surface as modified by adsorptive interactions with the soil.
As the concentration at the surface of the soil changes, volatilization may
become more dependent on the rate of movement of the pesticide to the soil
surface (Jury et al., 1983b) .
38
-------
___ — — — REFERENCE
PLANT COMPARTMENT HEIGHT
VOLATILIZATION FLUX
CANOPY HEIGHT
HORIZON 1
HORIZON 2
HORIZON 3
b
VOLATILIZATION
PLANTS
VOLATILIZATION
FROM SOIL
PLANT
COMPARTMENT
-^-U-jM-^J
( V%W-,W'> • • W)«'Mfo
AIR/SOIL
OUNDARY LAYER
i
VAPOR DIFFUSION'
SOIL LAYER
Figure 2.3. Schematic of pesticide vapor and
volatilization processes.
39
-------
The soil surface layer can be visualized as a membrane which only allows
water to pass through and keeps the solute behind. Experimental results show
that within the top centimeter of the soil surface, the pesticide
concentration can increase as much as 10-fold due to the accumulation of
chemical at the surface layer, resulting in higher vapor density. In order
to describe these phenomena, Jury et al. (1983a, 1983b) proposed a boundary
layer model which states that the controlling mechanism for pesticide
volatilization is molecular diffusion through the stagnant surface boundary
layer.
The layer of stagnant air may or may not form a significant barrier to
volatilization loss for a given pesticide, depending on a variety of factors.
In general, if the diffusion rate through the air layer is able to match the
upward flux to the soil surface without having the surface concentration
build up, then the stagnant layer is not acting as a barrier to loss and the
volatilization flux will not depend strongly on the thickness of the boundary
layer. Conversely, if the diffusion rate through the air is less than the
flow to the surface by diffusion or mass flow, then the concentration at the
soil surface will not be close to zero, and the thickness of the air layer
will regulate the loss by volatilization. In other words, the significance
of the boundary layer model depends on the ratio of the magnitudes between
the upward soil pesticide flux and the boundary layer diffusion flux. Only
downward, advective movement of water is treated in PRZM Release I. In this
case, the sources that contribute to the upward soil pesticide flux are only
the diffusion processes in the vapor and dissolved phases, but not upward
water advection.
The zero chemical concentration upper boundary condition in the first release
was modified in accordance with Jury's boundary layer model. The pesticide
volatilization flux from the soil profile can be estimated as follows:
(Cg)1 - C*jd) (2-41)
where
J-i - volatilization flux from soil (g day )
O I
Da = molecular diffusivity of the chemical in air (cm day )
r\
A = cross-sectional area of soil column (cm )
d = thickness of stagnant air boundary layer (cm)
C I = vapor-phase concentration in the surface soil layer (g cm" )
&> *•
C j = vapor-phase concentration above the stagnant air boundary layer
(g cm'3)
The thickness of the stagnant boundary layer can be estimated using a water
vapor transport approach (Jury et al., 1983a). However, Wagenet and Biggar
40
-------
(1987) assumed a constant value of 5 nun for this thickness, which is
consistent with the values estimated by Jury. Consequently, the same
assumption of a 5-mm thickness for the stagnant layer has been used here
pending the results of further sensitivity analyses. The value of C ^ can
take on a value of zero if the soil surface is bare or can be positive if a
plant canopy exists.
2.3.4.2 Volatilization Flux through the Plant Canopy--
In pioneering work on this topic, Parmele et al. (1972) discuss a number of
micrometeorological techniques for calculating pesticide volatilization flux
from observed aerial pesticide concentrations. Their procedures are based on
the assumption that the vertical diffusivity coefficient (Kz) for pesticide
vapor is analogous to the vertical diffusivity for water vapor, energy, or
momentum. The pesticide volatilization flux can be computed by Pick's first
law of diffusion, as follows:
JZ(Z) - -KZ(Z) (dP/dZ) (2-42)
where
JZ(Z) - pesticide flux at height Z (g m"2 s"1)
o
(dP/dZ) - pesticide concentration gradient (g m )
KZ(Z) - the vertical diffusivity at the height Z (m2 s"1)
The value of KZ depends on the turbulent flow of the atmosphere into which
the pesticide vapor is dissipated. Therefore, it is a function of the
prevailing meteorological conditions and not of any physical or chemical
property of the pesticide.
In order to apply these concepts, pesticide concentrations at two or more
heights are required to estimate the pesticide gradient and the subsequent
flux. For the estimation of vertical diffusivity, more extensive
meteorological information is also required. All of these data requirements
pose significant limitations for a predictive modeling approach consistent
with expected and current PRZM users.
In developing this PRZM module, the following approaches are proposed to
circumvent the intensive data requirements. First, a relationship for K is
derived as a function of height within the canopy. Then one need only
consider the pesticide concentration gradient (or a suitable surrogate) in
order to compute the pesticide volatilization flux.
Estimation of KZ(Z)--Mehlenbacher and Whitfield (1977) present the following
formula to compute KZ at various heights within the plant canopy:
KZ(Z) = KZ(ZCH) exp[4.0( A. . 1.0)1 (2-43)
L CH J
41
-------
KZ(ZCH) = U* k (ZCH - D
U*
k U
CH
ln[(ZCH - D/Z0)
(2-44)
(2-45)
where
KZ(ZCH>
ZCH
ID
k
U*
m
U
CH
thermal eddy diffusivity at height Z (m^ s'l)
9 -i
thermal eddy diffusivity at canopy height (m s )
canopy height (m)
roughness length (m)
zero plane displacement height (m)
von Karman's constant, 0.41
friction velocity (m s )
stability function for sensible heat
integrated momentum stability parameter as a function of ^
stability function for momentum
wind velocity at the canopy height (m s )
For agricultural applications, the canopy height is used as a reference
height for calculating U*. The user is required to input the wind speed and
the height where the measurement was made. The wind speed at the canopy
height (UQJ) is computed based on the logarithm law. The relationship is as
follows:
U,
CH
U
measured
In
In
"ZCH - 2
- Zo -
~7 D
measured —
Zo
(2-46)
The friction velocity U* can be visualized as a characteristic of the flow
regime in the plant canopy compartment in which the logarithmic velocity
distribution law holds. As shown in Equation (2-43), U* is calculated as a
function of U^jj, Z^^, ZQ, D, and V'm- Rosenberg (1974) describes ZQ + D as
the total height at which the velocity profile above the canopy extrapolates
to zero wind velocity. The values for both ZQ and D can be estimated with
the following equations presented by Thibodeaux (1979). For very short crops
(lawns, for example), ZQ adequately describes the total roughness length, and
little adjustment of the zero plane is necessary (i.e., D - 0). D is assumed
42
-------
to be zero in the current code when Z^ is less than 5 cm. For tall crops,
Z is related to canopy height (Z^jj) by
log ZQ - 0.997 log ZCH - 0.883 (2-47)
In tall crops Z_ is no longer adequate to describe the total roughness
length, and a value of D, the zero plane displacement, is needed. For a wide
range of crops and heights, 0.02 m < Z^ < 25 m, the following equation for D
has been presented by Stanhill (1969):
log D = 0.9793 log ZCR - 0.1536 (2-48)
This equation results from a linear regression analysis based on the
published data for 19 different crops with limited data measured for the same
crop at different growth stages. Strictly speaking, both ZQ and D should be
evaluated from experimental observations . In the calculation of KZ , the
module uses these two equations for estimation of ZQ and D, since there is no
method available to justify any variations for crop type, row spacing, or
canopy density.
With estimates of Z and D, U* (friction velocity) can be estimated if the
values of the stability parameters (^m and ^) are known. These two
variables are closely related to Ri, the Richardson number, which is the
measure of the rate of conversion of convective turbulence to mechanical
turbulence. It is defined as follows (Wark and Warner 1976):
(g/T) (3T/3Z)
Ri (2-49)
(SU/<9Z)2
where
e\
g = acceleration of gravity (m sec )
T - potential temperature (°K)
Z - elevation (m)
U — wind velocity (m s )
Potential temperature is defined as the temperature which a parcel of dry air
would acquire if brought adiabatically from its initial pressure to a
saturated pressure of 1000 millibars (Perkins 1974). In application of the
model, the measured temperature is used in the Richardson number estimation
as suggested by Rosenberg (1974).
The sign of Ri indicates the atmospheric condition, and its magnitude
reflects the degree of the influence. There are several different formulas
for relating Ri to the atmospheric stability parameters; for these purposes,
the sign of Ri is of greater concern than its magnitude. When Ri is larger
than 0.003, the atmosphere exhibits little vertical mixing, reflecting stable
43
-------
conditions; when the absolute value of Ri, |Ri|, is less than 0.003, neutral
stability conditions exist (Oliver 1971); and when Ri is less than -0.003,
convective mixing becomes dominant and atmospheric conditions are unstable.
To relate the atmospheric stability parameters to the Richardson number, Thorn
et al. (1975) proposed the following formulas based on the work by Dyer
(1974) and Dyer and Hicks (1970):
For stable conditions -
5'2Ri
For unstable conditions -
For neutral conditions -
*h - *m = 1 (2-52)
The integrated momentum stability parameter, ^m, can be evaluated based on
the following equation as derived by Lo (1977) :
ir/2 - 1 + ln(8) + ^ + 3 ln(^m) - 2
) - 2 tan-1(^m) (2-53)
Under neutral conditions, Vm = 0 and the equation is not used.
In the application of these procedures, the calculations are performed as
follows :
1) Evaluate Richardson number from temperature and wind velocity
gradients .
2) Determine stability condition based on calculated Ri.
3) Calculate 0n and ^m, based on the stability condition and associated
Equations (2-50), (2-51), or (2-52).
4) Calculate Vm, from Equation (2-53).
5) Calculate ZQ and D from canopy height using Equations (2-47) and (2-
48).
6) Estimate KZ(Z) by applying Equations (2-45), (2-44), and (2-43).
44
-------
The resistance approach for the estimation of volatilization flux from soil-
The calculation of the volatilization flux from the soil is based on a
resistance -type approach. For pre-plant pesticides, and time periods just
after emergence and post-harvest, transport by volatilization from plant
surfaces is much less than vapor phase transport by other mechanisms. For
those conditions in which the plant leaves do not act as significant sources
or sinks for pesticide vapor, the resistances of the air for the whole plant
compartment can be estimated as follows (Mehlenbacher and Whitfield 1977) :
Rbd + Rpc (2'54)
ZCR dz
RpC - - (2-56)
d K.(Z)
where
ZR — total vertical transfer resistance (day cm)
R^d = boundary layer resistance (day cm"1)
!LjC — plant canopy resistance (day cm" )
The flux is calculated as follows :
Pc
where
f\ -I
J — volatilization flux from plant canopy (g cm day" )
C -^ — pesticide vapor concentration in top soil layer (g cm)
For those conditions in which plants can act as significant pesticide sources
or sinks, another approach must be taken. The influences of plant canopy
require the formulation for the surface boundary condition as described in
the following two sections.
2.3.4.3 Volatilization Flux from Plant Surfaces--
A detailed description of the controlling factors for volatilization from
plant surfaces has been presented by Taylor (1978) . He indicated that the
distribution of the pesticide residues over the plant surface appeared to be
the dominant factor. This, together with the influence of the microscale
climate at the plant surface, makes accurate simulation of plant
volatilization processes very difficult.
45
-------
For organophosphate insecticides, Stamper et al. (1979) has shown that the
disappearance rates from leaf surfaces can be estimated by a logarithmic or a
first-order kinetics approach. Similar observations for first-order kinetics
were found for volatilization of 2,4-D iso-octyl ester from leaf surfaces by
Grover et al. (1985). Thus, a simple rate constant approach is possible,
which requires the user to input the first-order rate constant for
volatilization. The plant leaf volatilization flux can be estimated as
follows:
Jpl - M Kf (2-58)
where
o 1
J i = volatilization flux from the leaf (g cm day )
iy
M - foliar pesticide mass (g cm" )
KJ = first-order volatilization rate (day )
A resistance type approach is also applicable for volatilization flux
estimation from plant leaves. The current code employs the first-order
kinetics approach to calculate volatilization flux from plant leaf surfaces
described above. This approach, which requires the user to specify the
first-order rates constant for plant leaf volatilization, was selected
because it is consistent with the foliar fate model in PRZM Release I.
Average pesticide concentration in plant canopy--Volatilization flux from
plant leaves (Jp^) will exist only after pesticide application to the plant
foliage has been specified in the model input. When a plant canopy exists,
the average concentration in the air within the plant canopy can be estimated
as follows:
C* = (Jpc + Jpl) S R0.5 (2-59)
where
C — average concentration in the air between the ground surface and
g o
the plant canopy height (g cm )
ERQ c. = canopy resistance from half canopy height to the top of the
canopy
ZCR dz
= / (2-60)
0.5ZCH KZ(Z)
Equation (2-59) then calculates the mean plant compartment pesticide
concentration as the concentration at one-half of the canopy height. This
46
-------
approach assumes a linear concentration gradient from ground surface to
canopy height.
2.3.4.4 Soil Temperature Simulation--
Soil temperature is modeled in order to correct the Henry's law constant, K^,
for temperature effects. The interaction of microclimate with the soil
surface which results in a given soil temperature regime is complex and
dynamic. Soil surface configuration and plant residue cover, both affected
by tillage, have significant impact on soil heat flux and, therefore, soil
temperature. Studies of tillage and residue effects on soil temperature have
been dominated by qualitative observations and site-specific measurements.
The lack of mathematical evaluation and supporting field data has limited the
ability of researchers to predict, beyond qualitative terms, the tillage and
residue effect on soil temperature for soil and climatic conditions other
than those under which data have been collected.
The objective of the soil temperature model is to provide a scientifically
sound and usable approach: (i) to predict with reasonable accuracy the daily
average soil temperatures at the soil surface and in and below the root zone,
utilizing basic soil physical and thermal properties, and daily climatic
measurements taken at weather stations; and (ii) to allow consideration of
the residue, canopy, and tillage effects on soil temperature.
Several models are available to predict soil temperature under various soil
surface conditions, but there are restrictions to the general use of these
models because either they need large data bases which are not available at
many places, or they are site specific. Existing soil temperature models
form two general groups: (1) process-oriented models, which require detailed
information on soil and surface characteristics, initial and boundary
conditions, and inputs, and (2) semi- or non-process-oriented models, which
often utilize weather station information and soil temperature information at
one depth to develop empirical relationships.
Table 2-1 summarizes the key characteristics of the soil temperature models
reviewed in this work. For both the process and semi-process oriented
models, the two primary components are estimation of soil surface (or upper
boundary) temperatures and soil profile temperature utilizing the calculated
or estimated surface temperature as the upper boundary condition. A number
of the models utilize the same procedure for calculating temperature in the
soil profile (Gupta et al., 1981; Wagenet and Hutson 1987) and differ only in
the procedures for specifying the surface boundary condition.
Van Bavel and Hillel (1975, 1976) developed a dynamic numerical procedure to
link together the process-oriented simulations of heat movement in the soil
and the partition of heat and energy at the soil surface. Soil surface
temperature, TQ, is calculated as a factor in predicting evaporation from a
bare soil. Their technique utilized simultaneous solutions of seven
equations with seven unknowns: net radiative flux, evaporation rate, air
sensible heat flux, soil sensible heat flux, surface soil temperature,
Richardson's number, and the saturation humidity at the surface soil
temperature. Heat and water (liquid) flows are each coupled at the soil
surface. An iterative procedure was used at each update to find the proper
47
-------
CJ
h-
1/1
(-4
o;
UJ
u
|
3:
o
LU
8
X
UJ
or
ID
1 —
LJJ
Q_
X
Ul
h-
_J
O
CO
u.
O
SUMMARY
7
r\i
_j
m
i—
— • ro
C (OtO
CJ 01 v-'
o ^
0) 4-£
Ul~
J
(~
u
SS
a. *-^
*
— «
(0 ^
•MT-^CO
"S-.
CO •— C\J
ST?
3
3
(0 '•N
01 O
>«"• ^*
.c
t—
— > ft)
0)— • <-»
>— -in
03 — S-
ED 3= O
CTJ ^
<0 C
> (0
0)
o
-C
4-1
^
"cj
"2
6
X
X XX
X X
X XX
X
X X
X X
X X
X X
X XX
ft) ~O
4-* 0)
C 4J
0) C
TV — 0)
0) l_ —
4-> O (- OI
C t o ol
0) 01 i 01
"••- 0) 0) O
«< f_ 0) O) O C C
OIOUOI L.OOC
•Q i O U a..— .— o
O 01 L. O 4J 4-"^-
x 01 a. i- 3uo<-'
ajio. o 3 0) • —
a o £ c u. c CTJ
(_ 01 O O O <0
0)O.CO3E 4-iOOQ£
Q. (0
i- (n_Q o n: ra X) o
*— OJ
t-
UJ X
•*
UJ O
X O
X X
*
X
X X
UJ
X X
X
x
a.
0)1
— o o
C"- • r~.
ai o 4-> c o
t_ ••- (0 •• CT ^~
3 4-i — ' 0) UJ ^^ ^^
4->.~ 0) — •
L. (_ *4- O O «— "-
01 CO — ' O — ' S. 3
Q.Q- (0 t— u_ 0) CH
(C >»-^- 4J 3 ^"O
)— CD C_ 0) CO -Q O)
[_..- L. OJ 01 •'O C
XOIQ. 3 3: O — ' C —
I-CE 4-" O(OC04->
(0 UJ Ul (0 a L. 4-1
•b L. I Q. 4-1 (0 •—
cxx o>«— ajcou.
D JD_a Q. 01 C
O G OI^Z O) 0) 0)
m • • 0) c *-• ^£ L. >
4-» 4-> 1— >*- C 0> (—
(-O1O) >CD(0"-3
(U UJ UJ — ^ ^ C Z 3E CJ
=) ro JD co^unj-QU
ro ^
X
X X X X X O X
X X X X X
X XX XX XXX X
XX X
E
u
x xxx uixxx xx x
a
XX XX
X
XX X XXXX XXX X
XXXXXX X X
XXXXXXXXXX X >< X
.
5 oi c 0)
1— • 0) O 4->
— ' CO C Q-^> "- 4-i
u— i_o s w a
— O 3"- WH- — X
4-> i— o o> c
C, CO CO 4-1 Q U O
C C 3 I. — 4J.C>- g.— — C
3 c c Q.I— i-$4J3ocjoj— 'O oa)-n
p* co fo E co ^ *^- 0) °u a ra cj > c *>^
U Oil ' — • U O <0 L. L.O) OOtX
QC X X 1— •- O < O -~-,-C OJ -^ 01 CO OJ CJ CJ O
CO CO < CO — ' XCO 4-> — ' C CD U C
COXXOI 0)0)4-1 CD 3 •— 4-> CDOIOIC3
4J u >.>> o>— >-Dtm x c >4- 3 cocj ax
ro>>>»co — . — j co "OO- 0) L.T3CO OJ— 'X
Q — — H-t-t-S-Tj— 0-. — — .0 3—- o^ t— 1 3 Oj:j=a
>-i COJD OT3 OIH- o)j=.— — ..*— • co co^J cj t- ra J3
LA <1 N-
0)
L.
3
4-1
CO
(_
TJ E*
0) Q)
O) 4-*
-3 ^
01 L.
i- ra
raia
01 3
3 ^
(0 (1)
fl» &
Q.D
S CO •
B
>*•!_•*-(.
0 3T3 4)
01 4J ra *4-
(Q L. H-
ft) 0) C. O
C- Q. CO
3 e — ' o>
O) d) O 4->
(0 4J O)"-
s 5 —
ft) CO 4-1
O14J T3"-
CO C 0
C_ 0> — .f-
ft) .— CO — ^
^ Q I ' O.
< S o x
: < I— UJ
1 1 1 1
> I— X X ^ en 6
3 CO O C
(T «1 — TJ
0) Q. > C
g^4J _ CO
3 OJ O O
0 4J — > C
— • TJ C' —
H .fi •— E
• i- L.
4-> O Q. D-X
(0 O> " '
OJ o .£:••-
j: t_ 4-» 01 ra
o • c ~a
Q>*--O 01 o
i cj j= — -n
CsJ TJ Ol 4-< 4J C
o> 3 a. ro a
CO 4-» 0) — •
4-» Ol T3 CJ «
TJ — — t- 0
4) 4- *-* T3
01 C. C.— $
3 Ol O 0) co JJ
— 4-> <- O — •
^r c £ H- 4-*
O— c_.^ E o
^ ro Q. ^: H-
3 4-* 4-1 L,
• CT-C (0 g 3
CO Q.O)
..g-S^'Sc
0) .^- 4-> ••- O
W O> « 4--^
C «) C <- '— 4-*
o Q) •«— 4> — • ca
+j t. Q. a. Q.—
t- o> fc p E "O
O 4* CO 4) ••— to
Z tK O •»-• 00 U
III 1
* Q LU
* * Q X
48
-------
soil surface temperature. Soil temperatures were then estimated (Wierenga
and de Wit 1970) by using these estimates of T as the surface boundary
condition. Inputs required for this model include solar radiation, air and
dewpoint temperature, wind speed, initial soil temperature profile, and the
surface roughness evaluated by its effect on the aerodynamic roughness
parameter. No comparisons were made between predicted and measured soil
temperatures. Thibodeaux (1979) describes a similar energy-balance procedure
for calculating soil surface temperatures.
For modeling soil profile temperatures, Hanks et al. (1971) used a numerical
approximation for the one-dimensional soil-heat flow equation. This method
requires the input of initial and boundary conditions, as well as the soil
thermal conductivity and heat capacity as a function of depth and time.
Predicted root zone soil temperature profiles were within 1°C of observed
values for a 3-day period, but this model needs estimated or measured soil
surface temperatures as upper boundary condition.
Using the Hanks et al. (1971) procedure for the root zone, Gupta et al.
(1981-1984) developed a model for estimating hourly soil temperature by depth
from meteorologic data. Inputs needed for this model include hourly air
temperature at the 2-m height; daily maximum and minimum soil temperatures;
initial soil temperature with depth; and soil thermal diffusivity, which may
be estimated from soil mineral composition, organic matter percentage, bulk
density, and soil water content. The upper boundary temperatures are
estimated by a sine function. The amplitude of the function is equal to the
difference between daily maximum temperatures of air and soil surface or
daily minimum temperatures of air and soil surface. Empirical curves
relating daily maximum air temperature to daily maximum soil surface
temperature, and daily minimum air temperature to daily minimum soil surface
temperature, were developed for different residue and tillage conditions for
the specific application site. These relationships provided a means of
accounting for residue and tillage effects on soil temperature, but require
site-specific data.
The soil temperature model in PRZM is derived from a combination of the work
by van Bavel and Hillel (1976) and Thibodeaux (1979) for estimating the soil
surface/upper boundary temperature. The soil profile temperature procedures
were developed by Hanks et al. (1971) and applied by Gupta et al. (1981,
1982, 1983) and Wagenet and Hutson (1987)
Estimating upper boundary temperature--An energy balance procedure is used in
PRZM to estimate soil surface temperature (Thibodeaux 1979; van Bavel and
Hillel, 1976). The same procedure is used in the POSSM model (Brown and
Boutwell, 1986), which employs PRZM as a framework for PCB fate simulation.
-2 -1
The basic energy-balance equation with terms having units of cal cm day
at the air/soil interface may be described as:
Rn - Hs - LEs - Gs - ATH (2-6D
49
-------
where
Rn = net radiation (positive downward)
Hg = sensible air heat flux (positive upward)
LEg = latent heat flux (positive upward)
GS = soil heat flux (positive downward)
ATH = change in thermal energy storage in the thin soil layer (cal
9 1
cm day )
The term ATH can be evaluated as :
ATH - (/>bd) s (Ti+1 - T!) (2-62)
where
o
Pfo = bulk density of soil (g cm" )
d = thickness of a thin, surface soil layer (cm)
s = the specific heat capacity of soil (cal g °C )
Ti'Ti+l = t*ie rePresentati-ve temperature for the surface layer at two
consecutive time steps and can be represented as the average
of temperatures at the top and bottom of the soil layers.
For evaluating the heat exchange across the air/soil interface, the
thickness, d, can be set to a small value so that ATH may be neglected. As a
result, the right side of Equation (2-61) is set equal to zero.
Net radiation flux at any surface can be represented as :
Rn -
-------
Rsr = a Rs (2-64)
where
a = the albedo of the surface (dimensionless)
Therefore, the short-wave radiation component of the energy balance is
Rs - Rgr = Rs(l - a) (2-65)
The incident short-wave radiation can either be measured directly using
pyranometers or else calculated using a variety of available empirical
relationships or nomographs. The model requires input of a radiation time
series, whether measured or calculated, in order to simulate soil
temperature.
The albedo of a canopy-covered land surface can be estimated as:
o(t) - ac C(t) + as (1 - C(t)) (2-66)
where
a(t) = albedo on day t
a = albedo of canopy cover (0.23 for vegetation)
C(t) - canopy cover on day t (fraction)
ag = albedo of soil surface (dimensionless)
Since the albedo of soil surface changes with the soil surface condition, it
is defined by the user as 12 monthly values corresponding to the first day of
each month; the albedo value for each day is interpolated between the
neighboring monthly values. For snow cover less than 0.5 cm, the surface
albedo is estimated using Equation (2-66), and for snow cover above 0.5 cm,
the surface albedo is set equal to the snow albedo value (0.80).
The incident long-wave atmospheric radiation, R , is represented as
.Lcl
Rla - ea ° Ta4 (2-67)
where
ea = emissivity of the atmosphere [dimensionless]
a = the Stefan-Boltzmann constant [11.7*10~8 cal cm"2 °K"4 day"1]
Ta = the air temperature [°K]
51
-------
Wunderlich (1972) has proposed a correction to Equation (2.67) for the
effects of cloud cover, which could increase R, by up to 25 percent under
overcast conditions. However, this correction is not included in the model
because it would require input of a cloud cover timeseries, and the effect on
the calculated soil surface temperature would be small.
The emissivity of the atmosphere varies from a low of 0.7 to almost unity.
Numerous empirical relationships for estimating e have been proposed
(Salhotra 1986). A simple reliable method is the use of Swinbank's formula:
ea = 0.936*10'5 Ta2 (2-68)
The reflected long-wave radiation, R. , can be expressed as:
Rlar = Rla d ' 7> <2-69>
where
7 = the reflectivity of the surface for long-wave radiation
[dimensionless ]
The resulting net atmospheric long-wave radiation component becomes:
Rla - Rlar = ^a'1 ~ ?> = 0.936*10'5 TJ> * (1 - 7) (2-70)
The long-wave radiation component emitted by the soil surface is represented
in an analogous equation to the atmospheric component, as follows:
Rls = es * Ts <2
where
es = infrared emissivity of soil (dimensionless)
T = soil surface temperature (°K)
Since the soil emissivity and reflectivity are related as eg-l-7, we can
replace (1 - 7) in Equation (2-70) with eg.
Combining the radiation components from Equations (2-65), (2-70), and (2-71),
the net radiation flux is calculated as follows:
R - R (1 - a) + 0.936*10'5 a Ta6 eg - es a Ts4 (2-72)
52
-------
The evaporative heat flux, LE , is estimated by:
LES - /i E pw (2-73)
where
H = latent heat of vaporization/unit quantity of water
(580.0 cal g-1)
E - evaporation rate (cm day )
pv — density of water (1.0 g cm" )
The evaporation rate is obtained from the evapotranspiration (EVPOTR)
subroutine of PRZM. It is assumed that the calculated evapotranspiration
from the top 5 cm depth of soil represents the potential evaporation energy
loss at the air/soil interface. However, only a fraction of the
evapotranspiration loss calculated by PRZM contributes to this heat flux.
This fraction is estimated as the portion of the land surface not covered by
vegetation, (i.e., 1.0 - canopy cover).
The sensible air heat flux, H , is given by:
Hs - Pa Cpa h (Ts - Ta> (
where
•3
pa = air density (g cm" )
- (-0.0042 Ta + 1.292)10"3 (2-75)
C_a= specific heat of air at constant pressure
(0.2402 cal g'1 'IT1)
h = heat transfer coefficient at air-soil interface (cm day"-'-)
Ta - the air temperature (°C)
The air density is computed based on the daily air temperature using a simple
linear correlation (2-75) developed from data in Thibodeaux (1979). The heat
transfer coefficient is given by:
h - K-,
where
K-^ - Von Karman's number (0.41)
(2-76)
53
-------
Vz — wind velocity (cm day )
ZRJJ= reference height at which V is measured (m)
D = zero plane displacement (m)
ZQ = roughness height (m)
Equation (2-76) is valid only when the air temperature does not: vary greatly
with height, as is often the case near sunrise or sunset or under cloudy
skies or when canopy heights are relatively small. It appears to be a
reasonable approximation for most agricultural crops. Correlations have been
developed relating D and ZQ to the canopy height as described previously in
this section by Equations (2-47) and (2-48).
From the fundamental equation of heat conduction, the soil heat flux, Gs, is
given by:
Gs - (Ts - Tx) A1/D1 (2-77)
where
T~L — temperature of the soil at bottom of layer 1 (°K)
TS = soil surface temperature (°K)
AI = thermal conductivity of layer 1 (cal cm day" °K )
DI = thickness of layer 1 (cm)
Substituting Equations (2-70), (2-71), (2-72), and (2-75) into Equation (2-
59), the following fourth-order equation in terms of Tg results:
T - [(l-a)R_ + 0.936*10"5 a Tfl6 e_
o o d o
= 0 (2-78)
The value of T at each time step is estimated by solving the above equation
using an iterative solution based on the Newton-Raphson method. The initial
estimate of soil surface temperature is taken to equal measured air
temperature, and TcL., LES, Hg, and GS are calculated as explained above. The
value for Ti is obtained from the previous time step. These calculations are
repeated until the difference between two consecutive estimates for soil
surface temperature is less than the convergence criteria (set to 0.1°C).
Simulation of heat flow through soil profile--The soil profile temperature
model is based on the one-dimensional partial differential equation
describing heat flow in soils:
H - fc (« 0)
54
-------
where
d = the thermal diffusivity.
The thermal diffusivity is equal to the ratio of thermal conductivity and
heat capacity of the soil. The procedures used to estimate soil thermal
conductivity and heat capacity are taken from de Vries (1963). They are
calculated from basic soil properties; soil water content, mineral
composition, texture, and thermal conductivity of the individual soil
particles. These parameters are either input or supplied by the model in the
simulation. The thermal diffusivity is given by:
d = A/C (2-80)
where
f\ -I
d = thermal diffusivity of the soil layer (cm day )
A = thermal conductivity of the soil layer
(cal cm'1 day'1 "(T1)
C - heat capacity per unit volume of the soil layer
(cal cm'3 °C'1)
Temperature effect--A detailed discussion of the temperature effect on the
volatilization behavior of pesticides is presented by Streile (1984). Two
parameters that influence the vapor-phase transport in the soil profile are
Henry's constant and the vapor diffusion coefficient.
The equation used to correct Henry's constant for temperature effects is
(Streile 1984):
KR(T) = KR -L exp
R
(2-81)
where
^H 1 = Henry's constant at the reference temperature T-^
AHV - partial molar enthalpy of vaporization from solution
(J mole"1)
The temperature effect on the vapor phase diffusion coefficient can be
estimated from the Fuller correlation as presented in Liley and Gambill
(1973). However, it is not implemented in the code due to the general lack
of information required to use it.
55
-------
2.3.5 Irrigation Equations
PRZM irrigation algorithms determine depths of irrigation water to be applied
at the soil surface. These depths are computed from the soil water deficit
and are added as infiltration to the first PRZM soil compartment. Above and
below canopy sprinklers, flooding, and furrow irrigation can be simulated.
Methods for computing water application depths for each type of irrigation
are described in the following paragraphs .
2.3.5.1 Soil Moisture Deficit- -
Irrigation is triggered when the average root-zone soil moisture volume falls
below a level fc defined by the user as a fraction of the available water
capacity. The soil moisture deficit, D, is then given by:
D - (fc - 0Z) Zr (2-82)
where
D - soil moisture deficit (cm)
__ Q 0
Sz - average root-zone soil moisture content (cnrcm )
0fc- average root-zone soil moisture content at field capacity
o o
(cnr cm"J)
Zr - root zone depth (cm)
D is the depth of water over the unit area that must be added to the soil by
irrigation to bring the soil water content up to field capacity.
2.3.5.2 Sprinkler Irrigation- -
Irrigation water from sprinklers may be applied either above or below the
crop canopy. When applied above the crop canopy, irrigation water is
intercepted by the canopy and may run off when it reaches the soil surface.
The depth of water applied during a daily PRZM time step by overcanopy
sprinklers is estimated from the soil moisture deficit:
Da = D (1 + LF) + If (2-83)
where
Da = depth of irrigation water applied to the field (cm)
If = crop canopy interception capacity (cm)
LF = a factor specified by the user to allow for the practice in
saline soils of adding water to leach salts out of the root zone
(fraction of Da)
The water depth Da is applied as precipitation above the crop canopy, and
canopy interception is computed for the current crop in the PRZM crop growth
56
-------
subroutines. Sprinkler runoff from the soil surface is estimated using the
SCS curve number approach, assuming that runoff characteristics of sprinkler
water are similar to those of precipitation. Water that does not run off
infiltrates into the first PRZM soil compartment.
Irrigation water applied below the crop canopy is not subject to canopy
interception losses. The depth of water applied by undercanopy sprinklers is
therefore given by:
Da - D (1 + LF) (2-84)
The irrigation water depth APDEP is applied as throughfall to the soil
surface, and sprinkler runoff is estimated using the SCS curve number
approach.
In some instances the sprinkler system may be unable, due to hydraulic
limitations, to deliver water at the rate needed to meet the required daily
application depth. In these cases the sprinkler application depth Da is set
equal to the maximum depth that the system can deliver. The user is 1
therefore required to input the maximum water application rate R (cm hr )
for the sprinkler system.
2.3.5.3 Flood Irrigation--
Flood irrigation in this case refers to the practice of flooding entire
fields with irrigation water. Flood-irrigated fields are diked around the
edges to allow water to pond and infiltrate into the soil. In the PRZM
irrigation algorithm, it is assumed that this water ponds uniformly over the
entire field. The amount of water applied to the soil surface is then:
Da = D (1 + LF) (2-85)
Since the field is assumed to be diked around the edges, no water is allowed
to run off from the field.
2.3.5.4 Furrow Irrigation--
Furrow irrigation involves the release of water into numerous small channels
which cut across the planted field. Infiltration depths within furrows vary
due to differences in times at which water reaches various locations down the
furrow, with less water infiltrating at the downstream end (Figure 2.4).
Hydraulic characteristics of the furrow determine how quickly water moves
down the channel, while soil characteristics determine the rate of
infiltration once water reaches a location in the furrow.
The PRZM furrow irrigation model computes daily infiltration depths at
various locations down the length of the furrow. This requires solution of
the open channel flow equations of motion coupled with a soil infiltration
model. Numerous attempts have been made in the literature to solve the
furrow-irrigation advance problem, ranging in complexity from empirical
57
-------
CL
UJ
O
g
i-
DISTANCE ALONG THE FURROW
Inflow
End
Outflow
End
Figure 2.4. Variability of infiltration depths within
an irrigation furrow.
58
-------
volume-balance solutions (Wilke and Smerdon 1965; Fok and Bishop 1965) to
numerical solutions of the full open channel flow equations of motion
(Bassett and Fitzsimmons 1974). In general, solutions of the full equations
of motion are too computationally intensive for this application, while
simpler empirical models involve infiltration parameters which are not easily
related to physical soil characteristics.
The PRZM furrow advance model uses the kinematic wave simplification of the
equations of motion coupled with the Green-Ampt infiltration model to
determine furrow infiltration depths. Kinematic-wave theory neglects
inertial accelerations and assumes that the water surface slope is equal to
the ground slope. The equations of motion then reduce to:
3Q 3A -3q ,„ Rfi.
+~ = - (2'86)
where
O 1
Q = flow rate in the channel (m s )
9
A — cross-sectional area of flow (m )
x — distance down the furrow (m)
3 -1
q = volume infiltrated per unit length of channel (m m )
The flow area A is related to the flow rate Q by Manning's equation:
Q — AR2/3 S1/2 (2-87)
where
n — Manning's roughness coefficient
R = the hydraulic radius of flow (m)
S = the channel slope (vertical/horizontal)
Section 2.4.4 explains how the solution of the horizontal furrow irrigation
equation is applied to PRZM.
2.4 NUMERICAL SOLUTION TECHNIQUES
This section describes the numerical techniques which are used to solve the
differential equations introduced in the preceding section. Section 2.4.1
discusses the two numerical techniques available to solve the chemical
transport equations, a backwards-difference implicit scheme, and a Method of
Characteristics algorithm. The additional terms and the adjustment in the
upper soil boundary which are added into these transport equations to
simulate volatilization are described in Section 2.4.2. The numerical
approximations used to calculate soil temperature are presented in Section
59
-------
2.4.3 and the numerical solution for furrow infiltration depths are presented
in Section 2.4.4.
Although the equations in this section are shown with a constant compartment
depth (i.e. Az) for simplicity and clarity, the PRZM module allows the use of
variable compartment sizes so that each soil horizon can use a different
compartment depth. For example, this capability allows use of smaller depths
in the surface horizon for better simulation of surface processes, and larger
depths in the lower horizons to minimize computational time. This is
discussed in greater detail in Section 2.5.2.1.
2.4.1 Chemical Transport Equations
The second-order partial differential equation outlined in Section 2.3 must
be solved with appropriate boundary conditions. The calculations for
moisture contents, air contents, pore velocities, erosion, and runoff are
decoupled from, and solved in advance of, the transport equation. The
resulting values, treated as constant for each specific time step, are then
used as coefficients in a discretized numerical approximation of the chemical
transport equation.
Two techniques are currently available to solve the discretized chemical
transport equation for the new dissolved pesticide concentration at the end
of the time step. The available techniques are:
• A backward-difference, implicit scheme to simulate all chemical
transport processes
• A method of characteristics (HOC) algorithm which simulates
diffusion, decay, erosion, runoff, and uptake by the backward-
difference technique, but uses the method of characteristics to
simulate advective transport
The user is allowed to select the desired solution technique in the input
sequence. Details of these techniques are provided below. Results from test
simulations are provided in Section 2.5.1.
Identical discretizations and initial and boundary conditions are used with
both numerical simulation techniques. A spatial and time discretization step
is used equal to those applied in the water balance equations. For boundary
conditions at the base of the soil column, the numerical technique uses
- CWi «1
0 (2-88)
Az
in which the subscripts "i" refer to soil layer numbers.
This condition corresponds to a zero concentration gradient at the bottom of
the soil profile. The upper boundary condition is discussed in more detail
in Section 2.4.2.
60
-------
A backwards-difference solution algorithm was the only solution option
available in the original PRZM model. In this method, the first derivative
in space, the advection term, is written as a backward difference (i.e.,
involves the difference C[i,j]-C[i-l,j]). The second spatial derivative, the
diffusion term, is centered in space (i.e., based on the terms C[i-
1,j]+C[i+l,j]-2C[i,j]). The time derivative is also calculated as a backward
difference in the original code, (C[i,j]-C[i,j-l]). The equations are then
made implicit by writing each concentration for the (j+l)th time step. The
advantage of this numerical scheme is that it is unconditionally stable and
convergent. However, the terms truncated in the Taylor's series expansion
from which the finite difference expression are formulated lead to errors
which, in the advection terms, appear identical to the expressions for
hydrodynamic dispersion. In the simulation results these terms manifest
themselves as "numerical dispersion," which is difficult to separate from the
physical dispersion which is intentionally simulated. In systems with
significant advection (i.e., high Peclet number), the artificial numerical
diffusion may dominate the physical dispersion. It can be up to orders of
magnitude larger, leading to difficulty in the interpretation of simulation
results.
To minimize the effects of numerical dispersion in systems having high Peclet
numbers, a method of characteristics solution was added as an option to PRZM.
This solution method avoids the backwards-difference approximation for the
advection term and the associated numerical dispersion by decomposing the
governing transport equation. In advection-dominated systems, as the
dispersion term becomes small with respect to the advection term, the
advection-dispersion equation approaches a hyperbolic equation. According to
the MOC theory, advection of the solute can be simulated separately from the
other processes governing the fate and transport of that advected solute. M.
Baptista et al. (1984) state that no error is introduced by this
decomposition provided that the advection equation is solved first by an
explicit procedure, and the diffusion equation is solved next by an implicit
technique. This order was preserved in the PRZM model by utilizing a new
explicit algorithm for advection which is always called first, and is
immediately followed by execution of a modified version of the existing
implicit algorithm for simulation of other processes. The advection
algorithm employed was adapted from those described by Khalell and Reddell
(1986) and Konikow and Bredehoeft (1978) . These techniques were modified to
allow simulation of changes in saturation and adsorption of the pesticide and
variable compartment size.
In the new explicit advection algorithm, in addition to the fixed grid
system, a set of moving points is introduced. These points can be visualized
as carrying the chemical mass contained within a small region in space
surrounding the point. Initially, these points are uniformly distributed
throughout the flow domain. At each time interval, these moving points are
redistributed according to the local solute velocity in each compartment.
New points may enter the top of the flow domain, while old points may move
out the bottom. When the moving points are transported in horizons where the
compartment size is larger and numerical resolution is less, the points may
be consolidated to conserve computational effort. After the new locations
have been assigned to each point, the average concentration in each
61
-------
compartment is computed based on the number and mass carried by the points
contained within the compartment at that time. This temporary average
concentration is returned to the main program, and a subroutine which
assembles the terms in the transport equation (without advection) is called.
Changes in concentration due to all other fate and transport processes
(diffusion, decay, sources, etc.) are calculated for each compartment exactly
as in the original version of PRZM. These values are then returned to the
main program, and one transport step is complete.
When the MOC algorithm is called during the next time step, the exact
location of each moving point has been saved. The first task is to update
the masses carried by each moving point using the changes calculated during
the last time step. Increases in mass are simply added equally to each point
in the compartment, while decreases are weighted by the actual value at each
point before subtraction to avoid simulating negative masses. The updated
moving points are then relocated and the two-step process is repeated again
until the end of the simulation.
Under some extreme conditions (e.g. high rainfall on soils with very low
field capacity values), mass balance problems can occur when advective
transport is simulated by means of the MOC option. The mass balance errors
are related to the generation of extremely high velocities. Because a daily
time step is used, high velocities can cause the moving points used in the
MOC option to pass through several soil horizons without the properties of
those horizons being considered. If significant chemical mass is in the soil
column when the high velocities occur, mass balance problems can result. A
statement printed in the PRZM output file warns users who are using the MOC
option when high velocities are generated. If unacceptable mass balance
errors are noted, the simulation should be reproduced using the backward-
difference option instead of the MOC scheme.
2.4.2 Volatilization
The numerical techniques discussed in the previous section (2.4.1) are the
basis of the simulation of chemical transport in all phases. However, the
following modifications have been made to the upper boundary condition in
order to model volatilization of chemical from the soil surface.
In order to simulate vapor-phase pesticide movement past the soil surface,
the zero concentration upper boundary conditions used in the original PRZM
code has to be modified. Jury's boundary layer model (1983a, 1983b) has been
incorporated into the PRZM code. The model states that the controlling
mechanism for pesticide volatilization is molecular diffusion through the
stagnant surface boundary layer. The volatilization flux from soil profile
can be estimated as follows:
Jl=¥(Cg,l - Cg,d) (2-89)
where
volatilization flux from soil (g day )
62
-------
9 1
Da = molecular diffusivity of the chemical in air (cm day )
C -i = vapor-phase concentration in the surface soil layer (g cm" )
&> ••-
j*
C ^ — vapor-phase concentration above the stagnant air boundary layer
O
(= 0, for the no -canopy field condition) (g cm" )
d — thickness of stagnant air boundary layer (cm)
This equation defines the new flux- type boundary condition for the
volatilization simulation. In order to incorporate the new flux-type
boundary condition into the PRZM code , new mass balance equations were
derived for the surface soil and stagnant air layers. Figure 2.5(a) is a
schematic of the top two soil layers and the stagnant surface boundary layer
when no plant canopy exists. Zero concentration is assumed for C ^ under
the no -canopy field condition.
A mass balance equation for the uppermost soil compartment is
A(a C
v
AC Da
AT - A d- cg,i - v a
g,i
<2-90)
where
D =
V
A
a
K
O
molecular diffusivity of pesticide in air filled pore space
(cm2 day"1)
= volume of the compartment (cm )
r\
= area of the compartment (cm )
o o
— volumetric air content (cm cm" )
first-order reaction rate constant (day )
The first term of the right-hand side of Equation (2-90) represents the gas
diffusive flux into the surface soil layer, and the second term denotes the
gas diffusive output as governed by the stagnant boundary layer above the
soil surface. By using backward implicit finite differencing, the following
is derived:
K
H
At
Az
- D [2,n] KR Cw[2,n]
&-
At
Az'
K
H
K
R
Kg)
At D
Az d
KH Cw[l,n] (2-91)
where
n = time index
63
-------
Cg,d* —
//\\
Ax
= KHC1
Ax
= KHC2
(a) without plant canopy
C* *~ 0
-------
By substituting Equation (2-91) into the overall (i.e., all phases) mass
balance equation for the uppermost soil layer, a flux-type upper boundary
condition is obtained. Figure 2.5(b) reflects the field situation when a
plant canopy exists. Zero concentration is now assumed to exist above the
top of the canopy compartment. The volatilization flux from the plant canopy
is defined as follows:
'PC ,^n . ™, (C, 1 - C*) (2-92)
where
9 1
J = volatilization flux through the plant canopy (g cm day )
SR = vertical transfer resistance (day cm, described in
Section 2.3.4.3)
C = concentration above the plant canopy (assumed to be zero)
By carrying out a similar mass balance using finite differences, the boundary
condition which describes the field with canopy existing is obtained.
2.4.3 Soil Temperature
Soil temperature is solved for numerically. Section 2.3.4.4 describes the
theoretical basis for the simulation of soil temperature. The distribution
of temperature within the soil profile is summarized by Equation (2-79).
This equation is solved numerically for soil temperature, T, as a function of
depth, Z, and time, t, based on the input thermal diffusivity, d, for each
soil compartment, and the following initial and boundary conditions:
Initial Condition:
Tz>0 = T(z) (2-93)
Boundary Conditions:
To,t
TL,t
where
T(z) — initial soil temperature in each soil compartment (°C)
Ts(t) - calculated soil surface temperature for each time step (°C)
65
-------
= lower boundary temperature condition at the bottom of the soil
core (°C)
The lower boundary temperature is defined by the user as 12 monthly values
corresponding to the first day of each month; the value for each day is
interpolated between the neighboring monthly values.
The following numerical approximation used in the model is taken from Hanks
et al. (1971):
At Az2
Equation (2-96) is solved using a modified numerical solution procedure of
Hanks et al. (1971) which involves the same finite difference technique and
tridiagonal matrix solver (Thomas algorithm) used in PRZM (Carsel et al.
1984).
2.4.4 Furrow Irrigation
To simplify the algebra required to calculate the furrow infiltration volume
as Manning's equation is substituted into the kinematic wave model (equation
2-86), Manning's equation is approximated as follows:
A = a Qm (2-97)
a and m are constants which are estimated by the model from the parameters of
Manning's equation as follows:
ln(A2) -
m -- (2-98)
ln(Q2) - ln(Qx)
(2-99)
where
o
T , Ao= cross -sectional areas (m ) at depths y^ and y2
O I
!» ^2= fl°w rates (nr s"1) computed from Manning's equation [(Equation
2-75)] at depths y-^ and y2
- = 1 cm
66
-------
j2 = 10 cm
The depths y-^ and y>^ were chosen to represent the range of depths likely to
occur in furrows.
Substituting Equation (2-97) into Equation (2-86):
No closed- form solution to the above equation is known when infiltration is
time -variable. Equation (2-100) is therefore solved for Q using the
backwards -space, backwards -time, finite -difference solution described by Li
et al. (1975). Writing Equation (2-100) in finite -difference form:
where
Q^ = flow rate at time k, station i
Az = spatial step
At = time step
Infiltration volumes are computed using the Green-Ampt model:
aij f (H + H )«1
3£ = Ks [l + - i - J
(2-102)
where
1^ = infiltration depth at time k (m) , station i
KS = saturated hydraulic conductivity of the soil (m s )
H — ponded water depth (m)
Hg - suction parameter (m)
6 — available porosity (fraction)
I = total volume of infiltrated water (m)
The Green-Ampt model has long been accepted as a model of the advance of the
wetting front through the soil column, and involves parameters which can be
related to well-known soil properties. The volume of infiltration is
67
-------
computed assuming l£ is an average infiltration depth for the channel at
location i:
q = w i (2-103)
where
q£ = volume infiltrated at location i (m m" )
W^ = current flow width at location i (m)
Furrow channels are assumed to be trapezoidal in shape. Equation (2-87) is
solved at each station at the end of each time step for the new flow rate
Qn -• . Because the equation is non-linear with respect to Q, the new value of
flow is found using second-order Taylor series iteration. Given the flow
rate in the furrow, infiltration depths at each location are then computed
using the Green-Ampt model [Equation (2-102)].
The PRZM furrow irrigation model determines infiltration depths at various
locations in the furrow. Irrigation continues until the depth of water
infiltrated at the downstream end of the furrow is sufficient to meet the
soil moisture deficit SMDEF. The depth of water applied as irrigation to the
first PRZM soil compartment is then set equal to either the average furrow
infiltration depth or the infiltration depth at a specific location in the
furrow, depending upon options selected by the user. This depth of water
then infiltrates through the root zone as determined by the PRZM soil
hydraulic algorithms.
2.5 RESULTS OF PRZM TESTING SIMULATIONS
This section includes the results of testing the two solute transport
solution techniques and the volatilization algorithm. Simulated results are
compared with those from analytic solutions. Sensitivity analyses were also
performed to evaluate the effects of key model parameters on the prediction
of volatization rates. A test comparison of the model with field data from
Georgia (soybeans) concludes the section.
The PRZM model has undergone additional performance testing with field data
in New York and Wisconsin (potatoes) , Florida (citrus) , and Georgia (corn)
(Carsel et al., 1985; Jones 1983; Jones et al., 1983). The results of these
tests demonstrate that PRZM is a useful tool for evaluating groundwater
threats from pesticide use. Please refer to these references for information
regarding the further testing of PRZM under field conditions.
2.5.1 Transport Equation Solution Options
Currently, two numerical solution options are available to the PRZM user for
the chemical transport equation. As discussed in Section 2.4.1, the finite
difference option (utilizing subroutine SLPSTO) is unconditionally stable and
convergent, but may result in excessive numerical dispersion in high Peclet
68
-------
number systems. The HOC method of characteristics algorithm (utilizing
subroutines HOC and SLPST1) eliminates or reduces that numerical dispersion.
Two examples are provided which compare the alternate solutions methods at
high Peclet number (greater than 5.0) and at low Peclet number (less than
.5).
2.5.1.1 High Peclet Number--
Figure 2.6 presents the analytical solution (Hunt 1978) together with the „
SLPSTO and MOC/SLPST1 solutions at 6 days for the transport of a 69 mg cm
pesticide application in the uppermost compartment. The physical parameters
are as presented in the figure--notably the Peclet number is 5.1. The
following table details pertinent features of the simulation:
Location Value
of of Peak % Error Runtime
Method Peak (mg/cm3^ at Peak (sec)
Analytical 5.8 11.2
SLPSTO 4.5 5.07 -54 88.5
MOC/SLPST1 5.5 12.09 +7 112.4
At this relatively high Peclet number, the SLPSTO algorithm shows excessive
numerical dispersion, capturing only about half the amplitude of the peak
concentration, while showing excessive mass in both tails. In addition, the
SLPSTO algorithm does not predict the location of the peak precisely (it is
lagged behind the location of the peak given by the analytical solution and
the MOC/SLPST1 solution). In this example, the MOC/SLPST1 algorithm requires
27% more runtime, but errs by only 7% in the peak and shows good agreement in
the tails. Users should note that the increase of runtime could be
significantly high for longer simulation periods and deeper soil profiles.
2.5.1.2 Low Peclet Number--
Figure 2.7 illustrates the results of a SLPSTO and MOC/SLPST1 simulation 8
days after an incorporation of 69 mg/cm in the sixth compartment using the
parameters listed. The predicted concentrations at this lower Peclet number,
0.46, are very similar in the peaks and the tails, and apparently little
additional resolution is gained from utilizing the HOC algorithm. However,
the additional computational burden associated with the HOC algorithm is only
7% in this example.
2.5.2 Testing Results of Volatilization Subroutines
To test and validate the operation of the volatization algorithms, model
results were compared with Jury's analytical solution (Jury et al., 1983a),
and against field data for trifluralin from Watkinsville, GA. Sensitivity
analyses were also performed to evaluate effects of key parameters on model
predictions. The intent of this preliminary model testing was to evaluate
model operation by comparing the results for the volatization flux from a
soil surface application.
69
-------
7 9 11 13
COMPARTMENT NUMBER
Velocity =11.1 cm/day
Diff coef = 2.0 crn2/day
Retardation Coef = 11.74
Decay Constant = .10 /day
Core length = 20 cm
Delta x
Delta t
1 cm
1 day
Peclet = 5.1
Figure 2.6. Comparison of simulation results at high Peclet number.
70
-------
2
u
z
o
\
7 9 11 13
COMPARTMENT NUMBER
15
17
19
Velocity =1.82 cm/day
2
Diff Coef = 4.0 cm /day
Retardation Coef = 11.74
Decay = .1 /day
Delta x = 1 cm
Delta t = 1 day
Core length = 20 cm
Peclet = .46
Figure 2.7. Comparison of simulation results at low Peclet number.
71
-------
2.5.2.1 Comparison with Analytical Solution--
Jury et al. (1983a) presented a mathematical model for describing volatile
loss and movement of soil-applied organic chemicals. By making the following
assumptions, they derived an analytical solution for evaluating the chemical
concentration profile within the soil and the volatilization flux at the soil
surface:
1) Uniform soil properties consisting of a constant water content, bulk
density, liquid water flux (either upward, downward, or zero), and a.
constant organic carbon fraction
2) Linear equilibrium adsorption isotherm
3) Linear equilibrium liquid-vapor partitioning (Henry's law)
4) Uniform incorporation of a quantity of chemical to a specified depth
below the surface
5) Pesticide loss by volatilization through a stagnant air boundary
layer at the soil surface
6) Infinite depth of uniform soil below the depth of incorporation
Assumptions 2 to 5 are satisfied by the current PRZM code. Assumption 6
defines zero concentration for the bottom layer, which is somewhat different
from PRZM's zero gradient bottom boundary condition. However, as long as no
chemical reaches the bottom layer, these two types of boundary conditions
produce identical results. Our test runs for volatilization were designed to
satisfy this requirement. In order to comply with assumption 1, the
hydrological computation subroutines in PRZM were bypassed and replaced with
a constant value for water flux. A positive flux value indicates a leaching
condition, while a negative flux value indicates an evaporating condition.
The hydrological subroutines in PRZM are based on a moisture-routing method
in which daily accounting of water inflow and outflow is recorded. One
limitation of the moisture-routing method is that it is unable to properly
describe the upward movement of evaporating water. Evaporation loss is
removed from specific surface soil layers without accounting for movement
between layers.
The pesticide 2,4-D was chosen as the test compound for our simulation; the
input parameters are listed in Table 2-2 and were obtained from Jury et al.
(1983a). The test run results for daily volatilization flux are presented in
Figures 2.8(a), 2.8(b), 2.9(a), and 2.9(b), corresponding to the four test
cases listed at the bottom of Table 2-2. Two different soil compartment
depths (DELX) of 1.0 and 0.1 cm were used to investigate the sensitivity of
the volatilization algorithms to the spatial discretization in the surface
soil horizon.
Figure 2.8(a) shows the steady state situation (i.e., no evaporation and no
leaching) without any advective movement. The daily volatilization flux
values predicted by the two different DELXs are almost identical. In this
case, the magnitude of DELX is relatively unimportant. The simulation
72
-------
TABLE 2-2. INPUT PARAMETERS FOR THE TEST CASES - ANALYTICAL SOLUTION
DG Air diffusion coefficient 0.43 (m2 day"1)
DL Water diffusion coefficient 4.3 x 10" (m2 day"1)
Porosity 0.5
p Bulk density 1.35 (kg 111.3)
T Temperature 25°C
foc Organic carbon fraction 0.0125
8 Water content 0.3
a Air content 0.2
M Pesticide applied 1 (kg ha"1)
L Depth of incorporation O.lm
KH Henry's constant for 2,4-D 5.5 x 10'9
KQC Organic carbon partition
coefficient for 2,4-D 0.02 (m3 kg"1)
-2 i
u Decay coefficient for 2,4-D 4.62 x 10 (day"1)
1 Total depth of soil column 0.3 m
t Simulation period 30 days
Jw Water flux
E Evaporation flux
Test case #1: no evaporation and no leaching (J = E = 0)
Test case #2: with leaching (Jw = 0.01 cm day"1)
Test case #3: with evaporation (E - 0.01 cm day"1)
Test case #4: with evaporation (E = 0.25 cm day )
73
-------
c
0
•H
•p
0
N
0
o
CE-6 Kg/ho-doy)
15
10
— Analytic Sol'n
PRZM Results
D DELX -= O. 1 cm
* DELX - 1 cm
a.
TO 12 14 16 18 2O 22 24 26 28 3O
4 68
No Evapor-at i on S. No Leaching
0
a
Figure 2.8.
CE-6 Kg/ha-day)
15
0
N
•H
•-)
•H 10
4J
0
— Analytic Sol'n
PRZM Results
O DELX - O. 1 cm
* DELX - 1 cm
b.
6 8 10 12 14 16
LQaching Rats
18 2O 22 24 26 28 3O
• Q. 01 cm/day
Comparison of volatilization flux predicted by PRZM
and Jury's analytical solution: Test cases #1 and #2.
74
-------
C 15
0
•*4
4J
0
N
•H 10
+J
0
r-l
0
0
D
-6 K/no-doy)
— Analytic Sol'r
PRZM Results
O DELX - O. 1 cm
* DELX - 1 cm
a.
B1O 12 T4 16 18 2O2224 26 28 3O
Cday>
4 e
Evaporation Rata •• O. 01 cmXday
CE-6 KgXho-doy?
C 60
0
.H 4O
0
ft
0
0
O
— Analytic Sol'n
PRZM Results
O DELX - O. 1 cm
•«• DELX - 1 cm
b.
4 B
Evapor-at 1 on Rate «• O. 25 cm/day
B TO T2 14 IB TB SD 22 24 2B 28 3O
Cday)
Figure 2.9. Comparison of volatilization flux predicted by PRZM
and Jury's analytical solution: Test cases #3 and #4.
75
-------
results with a leaching rate of 0.01 cm/day are shown in Figure 2.8(b).
Because of the leaching influence, the predicted daily flux is smaller than
the corresponding daily value shown in Figure 2.8(a). The differences beween
the analytical solution and the PRZM predictions are due to the finite
difference solution technique and the occurrence of advective movement by
leaching. The simulation results using the smaller DELX (0.1 cm) more
closely match the analytical solution results, and an even smaller DELX would
have improved the agreement further. The slope of both DELX curves is the
same as the analytical solution, and the maximum differences (for the 1.0 cm
DELX) from the analytical solution are 10% or less.
Figure 2.9 shows the simulation results under evaporating conditions with the
upward advective velocity at 0.01 (Figure 2.9(a)) and 0.25 (Figure 2.9(b))
cm/day. The "wick effect" phenomenon (described in Section 2.3.4) leading to
enhanced upward movement of the pesticide can be observed in these two
figures. The maximum daily flux occurs on the first day for the leaching
conditions. Depending on the magnitude of the evaporating water velocity,
the maximum daily flux no longer occurs on the first day of the pesticide
application. Also the magnitude of the maximum daily flux is enhanced by the
magnitude of the evaporating water velocity. The effect of DELX becomes more
critical as the influence of advective movement increases. For simulations
using a 1.0-cm DELX, Figure 2.9(a) shows stable numerical behavior with a
small discrepancy when compared to the analytical solution result. As the
advective movement becomes larger, the numerical behavior becomes more
unstable, as shown in Figure 2.9(b). The smaller 0.1-cm DELX showed good
agreement with the analytical solution for both test cases shown in Figure
2.9.
Based on these test cases, it appears that a finer DELX, in the range of 0.1
to 0.5 cm, is needed for top soil layers when volatilization processes are
simulated with PRZM. However, this finer DELX requirements poses an
additional computational burden for PRZM applications due to the increase in
the number of soil compartments. To circumvent this burden, the PRZM code
was modified to allow a variable compartment depth, which allows the user to
select a smaller DELX for the top horizon (or any other horizon) and a bigger
DELX for the rest of the soil profile. By selecting this variable
compartment depth capability, a significant saving in CPU time may be
achieved while a better representation is provided for calculation of the
surface volatilization flux. In conjunction with field data comparisons
(presented below), the results of model runs and CPU time are presented for
both uniform and variable compartment depth simulation runs.
2.5.2.2 Comparison with Field Data--
Preliminary model testing with field observations was also performed to
assess the ability to predict the general magnitude of volatilization losses
and daily fluxes under field conditions. Based on a review of available
volatilization field data sets, a USDA experimental watershed site in north-
central Georgia was selected because of its use of a volatile pesticide
(trifluralin), surface-applied to a major crop (soybeans), with a
comprehensive micrometeorological and soil sampling plan.
76
-------
The study site was located at Watkinsville, GA, on a 1.26-ha watershed
comprised of Cecil soil (63.9% sand, 23.6% silt, and 12.5% clay) with 0.55%
organic carbon, a pH of 6.5, and a slope of 3.0%. Harper et al. (1976)
present a detailed description of the site, the equipment, and the
installation procedures required for collecting microclimate data. They also
summarize the method, assumptions, and calculations used for determining
pesticide volatilization flux rates. Trifluralin was surface-applied as a
spray to a bare soil surface, using a ground sprayer equipped with flat-fan
nozzles, at a rate of 1.12 kg/ha between 1220 and 1247 eastern daylight time
(EDT) on 15 June 1973.
The field results shown in Table 2-3 were obtained from White et al. (1977).
The values in columns 2, 4 and 5 of Table 2-3 provide the cumulative
volatilization flux, remaining pesticide in soil, and total cumulative decay
losses, respectively. A discrepancy is noted for the data in column 4 of
Table 2-3; the pesticide remaining in soil at the 35th day is smaller than
that at the 49th day. This discrepancy is most likely due to sampling
variations, although data were not available to establish accuracy limits on
the data points. Meteorological data required for applying PRZM to the site,
which include daily precipitation and pan evaporation, were obtained from
Smith et al. (1978).
The PRZM input parameters for trifluralin and the Watkinsville site are
listed in Table 2-4. Two additional key parameters which influence the
volatilization results are the decay rate and the adsorption partition
coefficient. The magnitude of the decay rate can be estimated from the data
in column 5 of Table 2-3, assuming that decay accounts for all losses from
the soil other than volatilization. A value of 0.0206 per day for the first-
order decay rate constant obtained from these data points is consistent with
the value of 0.0198 per day used by Donigian et al. (1986) after reviewing
the literature. An initial value for K^ was obtained from the organic carbon
content of 0.55% and an organic-carbon partition coefficient (Koc) value of
13,700, resulting in a Kd of 75 ml/g. Figure 2.10 shows the results of
sensitivity analyses runs for K^ and the decay rate; the observed data for
trifluralin from Table 2-3 is also included for comparison. Figure 2.10(a)
shows a good representation of the observed cumulative volatilization curve.
Figure 2.10(b) shows that a value of 40 for K^, and a decay rate of 0.02 per
day provides the best representation of the decay rate values analyzed.
The simulation results for cumulative volatilization flux and cumulative
pesticide decay are shown in Figure 2.11 for four different DELX
combinations. For these simulations, DELX values of 1.0, 0.5, 0.25, and 0.1
cm were chosen for the first horizon and 5-cm DELX for the rest of the
profile. The field data are also included in the figures for comparison.
Table 2-5 shows the total volatilization flux for each of the 4 using
variable DELX, as well as for a simulation using simulations, a constant 1.0-
cm DELX throughout the whole soil profile. The CPU requirements for each run
are also included in Table 2-5. The predicted total volatilization flux
using the smallest DELX of 0.1 cm is closest to the field-measured value; the
values for DELX of 0.25 cm and 0.50 cm are also quite close to the field
value. The saving of CPU time can be observed from Table 2-5. One hundred
twenty-nine seconds are required for the simulation using 1.0-cm DELX for the
77
-------
TABLE 2-3. TRIFLURALIN VOLATILIZATION LOSSES, AMOUNTS REMAINING IN SOIL, AND
ESTIMATED LOSSES VIA OTHER PATHWAYS FOR THE 120-DAY FIELD TEST
Time, (day)
Application
1
2
6
18
35
49
63
76
120
Source:
Cumulative
% of Total
Applied
3.5
3.8
5.3
10.9
20.5
23.4
24.4
25.1
25.4
25.9
White et al. ,
Volatilized
% of Total
Volatilized
13.3
14.8
20.3
42.2
79.1
90.2
94.1
96.9
98.2
100.0
1977.
...
Remaining
in Soil,
% Applied
89
72
64
51
33
35
23
20
11
j_ _ f\ .*_ _ ~I C
Estimated
Other Losses,
% of Applied
7.2
22.7
25.1
28.5
43.6
40.6
48.9
54.6
63.1
compared with an initial 1.0 /ig/g level at application (rate was
1.12 kg/ha).
TABLE 2-4. INPUT PARAMETERS FOR THE TEST CASES - WATKINSVILLE SITE
Simulation start date
Simulation end date
Trifluralin: Henry's constant
Diffusion coefficient in air
Application date
Amount applied
Incorporation depth
14 June 1973
31 December 1973
6.7 x 10
-3
0.43 m2 day'1
15 June 1973
1.12 kg ha'1
5 cm
Horizon
1
2
3
4
Thickness
(cm)
5
10
15
60
DELX
(cm)
0.1
5.0
5.0
5.0
Field
Capacity
.207
.207
.339
.320
Wilting
Point
.095
.095
.239
.239
Initial
Water Content
0.166
0.217
0.318
0.394
78
-------
0
>
01
D
E
3
L)
4O
20
10
of oppliad)
Field Data
PRZM Results
KD - 3O
KD - 4O
KD - 50
KD - 75
a.
2436 48SO 72 84
Sensitivity of KD
96
108 12O
0
rH
0
D
E
D
(J
40
30
2O
1O
of applied)
Fiald Date
PRZM Results
K - O. Ol
K - O. 02
K - O. O3
12
b.
"24 §S 48 63 72 84
Sensitivity oF Decay Rata
•95—ros—rso
Cday>
Figure 2.10. Sensitivity of cumulative volatilization flux to
K and decay rate.
d
79
-------
c
0
l-t
+J
0
N
4J
0
r-t
3
E
3
U
X of o
pplisd)
Fiald Data
PRZM Results
DEUX - 1.D a 5
_ DEUX - 0. 5 & 5 cm
DELX - 0. 25 & S o
DELX - O. 1 a 5 cm
12 24 36 48 BD 72 ff4 5E TDB 120
a. Effact oF DELX on Volatilization Fluxd
X
D
0
01
D
(V
T3
+J
I/)
01
D.
D
E
D
U
Figure 2.11.
7O
6O
so
4O
01 30
20
10
of applied?
Field Data
PRZM Results
DELX - O. 1 a S cm
DELX - O. 25 a 5 cr
DELX - O. 5 a 5 cm
DELX - l.O a 5 cm
b.
12 24 SB 4£T
EFfect of DELX
6O 72 B4
n Pest i cide
96 1OB
Decay
12O
(day)
Effects of DELX on volatilization flux and pesticide
decay.
80
-------
whole soil profile, compared with only 39 seconds for the simulation using
1.0 cm for the top horizon and 5.0 cm for the rest of the profile. The
results in Table 2-5 indicate that a DELX of 0.25-0.50 cm for the top horizon
may be a reasonable compromise between simulation accuracy and CPU costs.
TABLE 2-5. SIMULATION RESULTS OF USING DIFFERENT COMPARTMENT DEPTH (DELX)
Horizon
1
2
3
4
Total
Constant DELX
Depth DELX
(cm) (cm)
5
10
15
60
1.0
1.0
1.0
1.0
0.393
Variable DELX
DELX
(cm)
1.0
5.0
5.0
5.0
0.398
DELX
(cm)
0.5
5.0
5.0
5.0
0.338
DELX
(cm)
0.25
5.0
5.0
5.0
0.317
DELX
(cm)
0.1
5.0
5.0
5.0
0.316
Field
Value
0.290
Volatilization
Flux (kg/ha)
CPU (Sec)
129
39
46
67
106
Figure 2.12(a) demonstrates significant differences between the observed
pesticide decay and the simulated values during the first few weeks following
application. In fact, the observed data appear to indicate a much higher
attenuation rate during the first few days following application, with a
lower rate for the remaining period. To better match the decay
characteristics, and evaluate the potential impact on the volatilization
simulation, a two-step decay procedure was used with a rate of 0.1 per day
for 5 days following application and a rate of 0.01 per day for the remaining
period. The results of these simulations in terms of pesticide remaining in
the soil, shown in Figure 2.12, indicate a much better agreement with the
observed field values in Figure 2.12(b). The impact of the two-step decay on
both cumulative decay and volatilization flux is shown in Figure 2.13. The
cumulative pesticide decay shown in Figure 2.13(a) improves considerably
(compared to Figure 2.11(b)), while the results for cumulative volatilization
flux (Figure 2.13(b)) are slightly better than those in Figure 2.11(a).
2.5.2.3 Conclusions from Volatilization Model Testing--
The primary conclusions derived from this preliminary model testing are as
follows:
1) Comparisons with Jury's analytical solution indicate that the
volatilization algorithms are operating correctly, and that with a
very small DELX (0.1 cm or less), the results are in excellent
agreement.
81
-------
0C% of appl led>
^3l_J
" 80
0
c 70
•H
01 BO
TJ
•H
i 50
0)
Q) 4O
QL
0> 30
c
'lUo
0
E
0) 10
cr
n
\
~\V. F,.= 1H Hnl-n
yA PRZM Results
DELX = O. 1 & 5 cm
DELX - 0. 25 & 5 cm
\ \v DELX - O. 5 & 5 cm
- \ VY — — DELX - 1. O 8 S em
\\\
- \\
\\
V^x
^r\
s-^*>^ \
»'**>» \
^o*. \
v^s-N^—_
^^s**"v>» —
"""^ ;^^ -^___^
~;5?^-^^_^ ' ' ~-
^=5r— -..___
~
12 24 36 48 6O 72 84 96 IDS 1 2O
9O
7. of applied)
80
0
in
c 7D
0) 6O
U 50
B
Q) 4O
Q-
^ 30
D
E
01 10
a
Field Data
PRZM Results
DELX - O. 1 S. 5 cm
DELX - O. 25 & 5 cr
DELX - O. 5 8. 5 cm
DELX - l.O & S cm
12 24 36 48 6O 72 84
b. Simulations with Two—Step Decay Rates
96 JOB 12O
Cday>
Figure 2.12. Comparison of constant and two-step decay rates.
82
-------
X
D
0
01
D
01
T)
4}
(I)
01
Q.
•P
D
E
3
U
BO
7O
BO
50
4O
3D
2O
10
(7. of opplaecD
Field Data
PRZM Results
DEUX - O. 1 S 5 cm
DEUX - O. 25 & 5 em|
DEUX - O. 5 8. S cm
DEUX - l.O & 5 cm
a.
12 24 35 48 6O72§4SB 1 OB
S i mi_i 1 at i or-is with Two —Step Decay Rates
120
Cday)
of app 1 i ed5
U.
c
0
4J
0
4J
0
r-M
D
E
U
PRZM Results
DEUX - 1. O S 5 cm
DEUX = O. 5 a 5 cm
DEUX - O. 25 S 5 or
DEUX - O. 1 a 5 cm
b.
TS 24 36 48 6O 72
Simulations with Two — S
84 96 1O8 12O
(day)
Rates '
Figure 2.13. Effects of two-step decay rates on volatilization
flux and pesticide decay.
83
-------
2) The preliminary field testing results with trifluralin in
Watkinsville, GA, indicate good agreement between measured and
predicted volatilization flux when measured decay rates and adjusted
K^ values are used.
3) Small soil layer depths in the range of 0.25 and 0.50 cm are needed
to provide the best presentation of volatilization flux at
reasonable CPU times, based on the Watkinsville testing.
4) A two-step decay rate best represents the attenuation behavior of
trifluralin, using a higher rate for the period immediately
following application, and a lower rate for the remaining period.
Further testing of the volatilization model should be performed to evaluate
its capabilities for different compounds, different regions, and other crops.
In addition, the plant compartment vapor transport and concentration
calculations should be tested with the additional data available from the
Watkinsville site and from other field data sets (e.g., Grover et al., 1985;
Willis et al., 1983).
2.5.3 Testing Results of Soil Temperature Simulation Subroutine
Preliminary testing of the soil profile temperature simulation subroutine was
performed by comparing predicted values with values obtained by an analytical
solution to the governing heat flow equation. These testing results are
discussed in this section. Field testing of the soil surface/upper boundary
temperature simulation, estimated by the energy balance procedure in the
model, was not performed due to problems in obtaining observed meteorological
and soil temperature data for the Watkinsville, GA, test site.
An analytical solution presented in Kreysig (1972) for the classical one-
dimensional heat flow partial differential equation (described in Section
2.3.4.4) was used to calculate changes in the soil temperature profile with
time, due to a change in the upper boundary temperature. In order to develop
a valid comparison between the analytical and finite difference methods, the
following assumptions were made:
a) Uniform properties throughout the soil profile
b) Constant lower-boundary temperature
c) Uniform initial temperatures throughout the profile
To compare the results of the analytical solution with the finite difference
solution from the soil temperature model, the following parameters were used:
Depth of the soil profile - 100 cm
Compartment thickness (DELX) = 1.0 cm
Diffusivity of the soil profile = 864 cm2 day "1
Upper-boundary temperature, T/Q t\ = 30°C
Lower-boundary temperature, T/^'J-N - 20°C
Initial temperature, T/x o\ ' - 20°C
84
-------
Figures 2.14 and 2.15 show the comparison of soil temperature profiles
predicted by both the analytical solution and the finite difference soil
temperature model after 1 day and 5 days of simulation. In Figure 2.14 the
finite difference solution is obtained by using a 1-hour time step, while in
Figure 2.15 a 1-day time step is used. The following observations are
evident from these testing results:
1) Comparison of the soil temperature profiles predicted by both
methods indicate excellent agreement when the smaller, 1-hour time
step is used in the finite difference procedure, as shown in Figure
2.14.
2) The finite difference solution obtained by using the daily time
steps deviates from the analytical solution by about 1°C, in the
upper and middle portions of the soil profile (Figure 2.15). This
deviation is due to the assumption of a constant initial temperature
profile and the abrupt change in the upper-boundary temperature from
20°C to 30°C for the first daily time step.
3) As the steady-state condition is approached, irrespective of the
time step used in the finite difference solution, the soil
temperature profiles predicted by both methods are in good agreement
(Figures 2.14(b) and 2.15(b)).
Table 2-6 shows that reducing the depth of the compartment from 1 cm to 0.1
cm does not produce any significant change in the finite difference solution.
These depths bracket the range of values for DELX (i.e., compartment
thickness) likely to be used for the surface soil horizon.
These test results show that for smaller time steps the finite difference
solution will be in complete agreement with the analytical solution. For a
daily time step as used in PRZM, under expected environmental conditions,
with a non-uniform initial temperature profile, non-uniform soil
characteristics, and smaller daily changes in the upper-boundary temperature,
the soil temperature profile estimated by the finite difference method used
in the model is expected to be capable of providing close agreement with
observed temperature profile data. In addition to further testing of the
soil profile temperature model with field data, the procedure to estimate the
upper-boundary temperature should be tested to evaluate and demonstrate the
validity of the entire soil temperature simulation model.
2.5.4 Testing of Daughter Products Simulation
The fate of pesticides in soils is a complex issue. There are many processes
(i.e., volatilization, degradation, etc.) which must be considered in order
to adequately address this issue. One of these processes, which has been
largely neglected in pesticide leaching models, is that of the transformation
of the parent compound to various toxic daughter products. The tendency has
been to lump all the toxic family into a "total toxic residue" and model the
fate of this composite as a single chemical. This assumption may not be
acceptable, especially if the daughters have very different decay rates or
adsorption partition coefficients from the parent or from each other.
85
-------
(a) . Soil Temperature Profile After 1 Day
Profile Initiol Temp" 20 C
Upper Boundory Temp « 30 C
Lower Boundary Temp » 20 C
Time Step • 1 hr
Analytical Soln.
Finite Diff. Soln.
100
Depth of Soil Profile in cm
(b) . Soil Temperature Profile After 5 Days
Profile Initial Temp« 20 C
Upper Boundory Temp * 3D C
Lower Boundary Temp • 20 C
Time Step « 1 hr
Analytical Soln.
Finite Diff. Soln.
Depth of Soil Profile in cm
Figure 2.14
Comparison of soil temperature profiles
predicted by analytical and finite
difference solutions.
86
-------
(a) . Soil Temperature Profile After 1
30-*
27--
18--
15
Profile Initial Temp- 20 C
Upper Boundary Temp - 30 C
Lower Boundary Temp • 20 C
Time Step - Idcy
Analytical Soln.
Finite Diff. Soln.
10
20
30
50
60
70
80
90
100
Depth of Soil Profile in cm
(b). Soil Temperature Profile After 5 Davs
30
27-
01
3
?21
01
IB-
^v
Profile Initial Temp- 20 C
Upper Boundary Temp » 30 C
Lower Boundary Temp « 20 C
Time Step • Iday
Analytical Soln.
Finite Diff. Soln.
10 20 30 40 50 60 70
Depth of Soil Profile in em
60
90
100
Figure 2.15. Comparison of soil temperature profiles
predicted by analytical and finite
difference solutions.
87
-------
Table 2-6. SIMULATED SOIL TEMPERATURE PROFILE AFTER ONE DAY FOR DIFFERENT
COMPARTMENT THICKNESSES (TIME STEP = I DAY)
Depth (cm)
0.0
1.0
2.0
3.0
4.0
5.0
10.0
20.0
30.0
40.0
50.0
60.0
75.0
99.0
100.0
DELX = 1 cm
30.000
29.665
29.341
29.028
28.725
28.432
27.109
25.048
23.577
22.524
21.766
21.215
20.638
20.023
20.000
DELX - 0.1 cm
30.000
29.664
29.340
29.026
28.723
28.431
27.106
25.045
23.574
22.520
21.760
21.206
20.627
20.020
20.000
Algorithms have been included in PRZM to simulate parent/daughter
relationships. An analytical solution to the decay and transformation model
was derived to check the numerical model.
The system which was modeled is shown in Figure 2.16. The C^ are dissolved
concentrations and the C^ are adsorbed concentrations. The K^ are adsorption
partition coefficients, the k^ are decay and transformation rates in the
dissolved species, the kf are adsorbed phase decay coefficients and
8 and p are the water content and soil bulk densities, respectively. Notice
that only the dissolved forms may be transformed from one toxic form to
another. A system of first order differential equations describing this
system can be written as:
C^ (2-104)
- (k3 + k4) C20 + k2 C^e (2-105)
88
_
d C
-------
f ki* f V f V
: *
Kl ft K2
: •.*-
i tl
c * r * '
C2 C3 ,
ft K3 ft
k,, j
C> ^ =3 J
T k3 I k5
ADSORBED PHASE
DISSOLVED PHASE
Figure 2.16. Schematic of a system of parent and
daughter pesticides.
d C-)6
-^ k5 C28 + k4 C26 (2-106)
d Ci*fl
(2-107)
C2* ^ (2-108)
(2-109)
Making use of C^ K£ = C^ the six equations above can be reduced to three
equations in three unknowns, namely:
d C1
d C2
d C3
Cj_ (2-110)
Cx + a3 C2 (2-111)
C2 + a5 C3 (2-112)
89
-------
in which
*
k2) - kx7
^ v „ (2-113)
k4) -
(2-114)
*
- k3* K3p
These ordinary differential equations with constant coefficients can be
solved analytically for C^, C2 and C3 using the initial conditions Ci —
Ci when t = 0 and C2 = C3 = 0 at t = 0 . The solutions as given in Dean and
Atwood (1985) are:
Cl = Cl e*l (2-118)
f a2 } , a,t f
- C, e ! -
Ui - a3J 1 Ul
C2 " llTT^rl Cj. a-" - I-—— Ci e^ (2-119)
and
p _ | I /-«
3 I fa-i - a.r-} (a-i - a^)J
a2
a4 1 c, ea3t
(ax - a3)J 1
a2 a4 a2 a4 , act
+ T7 ^77 -T - T- -^T- —A C, e 5 (2-120)
In PRZM, the equations are solved numerically as part of the general
advection-dispersion equation for a solute in a porous medium by using an
implicit scheme. A new subroutine was added to set up the transformation
90
-------
(source and sink) terms for the system. The relationship C-^ -* Co -* C-j may be
modeled or the system can be configured for C-i -* Co and Ci -* Co or for
independent C-, , £^ and C-j simply by selecting zero or positive values for the
appropriate transformation rate constants.
Figures 2.17 through 2.19 show the results of a series of tests performed on
the numerical model and checked by the analytical model. In these figures,
the solid line represents the "true" or analytical solution, while the dashed
line represents the approximate numerical solution. ' In Figure 2.17, there
was no decay of the dissolved phase chemicals and no adsorption of any
species. The rate of transformation from Ci to Cj was 0.2 day" and that
from C2 to C-j, 0.5 day"-1-. After 20 days nearly all the chemical is in form
Cj. The numerical model traces the decay and formation of each constituent
closely, being poorer in those regions where the rate of change of the
concentrations are more rapid. Figure 2.18 shows the same system with a
decay rate of 0.01 day" in the dissolved phase. Figure 2.19 shows the
solution for the same system except that the adsorption coefficients have
been given values (K-^ = 0.5, 1^ = 1.0 and K-j = 5.0) as have the adsorbed
phase decay coefficients (^ = 0.01). Notice that the transformation rate is
retarded when adsorption is introduced.
Using the analytical model, the assumption of modeling the "total toxic
residue" decay as a first-order process was tested. Adsorption coefficients
for aldicarb, aldicarb sulfoxide and aldicarb sulfone in a Woburn sandy loam
(Ki = 0.55, Ko — 0.16 and Ko = 0.185) and decay and transformation rate
constants (^ - 0.07, k2 = 0.55, k3 - 0.01, k4 = 0.031 and k5 = 0.0152) were
taken from Bromilow et al. (1980). A soil bulk density of 1.45, a water
content of 0.27 cm cm and an initial aldicarb parent mass of 100 rag were
also used. The model was run for 90 days and the results are shown in Figure
2.20.
The results show that the decay of the sum of the dissolved aldicarb
concentrations does not follow first-order kinetics. The reason for this is
the conversion of aldicarb parent to aldicarb sulfoxide. Because the
sulfoxide has a lower partition coefficient, the dissolved concentration
increases until most of this conversion is complete. Once this happens,
however, the sum of the sulfoxide and the sulfone concentrations does follow
a first-order decay curve.
91
-------
o
c
ro
O
"8
o f
gg
• H -H
'(I)
>
C
O
U
CM
0)
S-l
3
XN30U3W Nl 'NOUVWJ.N3DNOO
92
-------
100
i 1 1 1 1 1 1 1 r
i
NUMERICAL
ANALYTICAL
• 10 12
TIME. IN DAYS
Figure 2.18
Conversion of Ci to C2 to C3 with decay
but no adsorption.
NUMERICAL
ANALYTICAL
8 10
TIME, IN DAYS
Figure 2.19,
Conversion of C
and adsorption.
93
to C2 to C3 with decay
-------
0)
•o
•H
X
O
CO
U
-H
•a
rH
ro
O
"s
•o S
OJ
CO
o
3
c
O o
o i;
o
CN
•
CM
-H
•J/fcr/
'NOI1YHXN30NOD
94
-------
SECTION 3
VADOSE ZONE FLOW AND TRANSPORT MODEL (VADOFT)
3.1 INTRODUCTION
VADOFT is a finite-element code for simulating moisture movement and solute
transport in the vadose zone. It is the second part of the three-component
RUSTIC model for predicting the movement of pesticides within and below the
plant root zone and assessing consequent groundwater contamination. The
VADOFT code simulates one-dimensional, single-phase moisture movement in
unconfined, variably saturated porous media. The code considers only single-
porosity media and also ignores the effects of hysteresis. Transport of
dissolved contaminants may also be simulated within the same domain.
Transport processes accounted for include hydrodynamic dispersion, advection,
linear equilibrium sorption, and first-order decay. VADOFT also simulates
solute transformations in order to account for parent/daughter relationships.
3.2 OVERVIEW OF VADOFT
3.2.1 Features
3.2.1.1 General Description--
The VADOFT code can be used to perform one-dimensional modeling of water flow
and transport of dissolved contaminants in variably or fully saturated
soil/aquifer systems. For demonstrative purposes, we consider a scenario
involving the migration of pesticide through the root zone, the vadose zone,
and the saturated zone of an unconfined aquifer (Figure 3.1). VADOFT can be
operated as a stand-alone code or operated in conjunction with the root zone
model, PRZM, and/or the saturated-zone model, SAFTMOD. In the latter cases,
boundary conditions at the interfaces of the modeled domains are established
via model linkage procedures.
3.2.1.2 Process and Geometry--
VADOFT performs one-dimensional transient or steady-state simulations of
water flow and solute transport in variably saturated porous media. The code
employs the Galerkin finite-element technique to approximate the governing
equations for flow and transport. It allows for a wide range of nonlinear
flow conditions, and handles various transport processes, including
hydrodynamic dispersion, advection, linear equilibrium sorption, and first-
order decay. Steady-state transport can not be simulated when decay is
considered. Boundary conditions of the variably saturated flow problems are
specified in terms of prescribed pressure head or prescribed volumetric water
flux per unit area. Boundary conditions of the solute transport problem are
95
-------
SATURATED ZONE
/////>
Figure 3.1.
////// /////////////////// //////
Cross-section view showing the solute migration
path in the vadose and saturated zones of an
unconfined aquifer.
specified in terms of prescribed concentration or prescribed solute mass flux
per unit area. All boundary conditions may be time dependent.
3.2.1.3 Assumptions--
The VADOFT code contains both flow and solute transport models. Major
assumptions of the flow model are as follows:
• Flow of the fluid phase is one-dimensional and considered isothermal
and governed by Darcy's law.
• The fluid considered is slightly compressible and homogeneous.
• Hysteresis effects in the constitutive relationships of relative
permeability versus water saturation, and water saturation versus
capillary pressure head, are assumed to be negligible.
Major assumptions of the solute transport model are as follows:
• Advection and dispersion are one-dimensional.
• Fluid properties are independent of concentrations of contaminants.
• Diffusive/dispersive transport in the porous-medium system is
governed by Pick's law. The hydrodynamic dispersion coefficient is
defined as the sum of the coefficients of mechanical dispersion and
molecular diffusion.
96
-------
• Adsorption and decay of the solute may be described by a linear
equilibrium isotherm and a first-order decay constant.
• Vapor transport can be neglected.
3.2.1.4 Data Requirements--
Data required for the simulation of variably saturated flow include values of
the saturated hydraulic conductivity and specific storage of the porous
media, the geometry and configuration of the flow region, as well as initial
and boundary conditions associated with the flow equation. Soil moisture
relationships are also required. These include relative permeability versus
water phase saturation and capillary head-versus water phase saturation.
These relationships may be supplied to the code using tabulated data or
functional parameters.
Data required for the simulation of solute transport in variably saturated
soil include dispersivity and porosity values, retardation and decay
constants, values of Darcy velocity and water saturation, as well as initial
and boundary conditions associated with the transport- equation.
3.2.2 Limitations
Major limitations of the VADOFT code are summarized as follows:
• In performing a variably saturated flow analysis, the code handles
only single-phase flow (i.e., water) and ignores the flow of a
second phase (i.e., air) which, in some instances, can be
significant.
• The code ignores the effects of hysteresis on the soil moisture
constitutive relations.
• The code does not take into account sorption nonlinearity or kinetic
sorption effects which, in some instances, can be important.
• The code considers only single-porosity (granular) soil media. It
cannot handle fractured porous media or structured soils.
• The code does not take into account transverse dispersion, which can
be important for layered media.
3.3 DESCRIPTION OF FLOW MODULE
3.3.1 Flow Equation
VADOFT considers the problem of variably saturated flow in a soil column in
the vadose zone of an unconfined aquifer. The code solves the Richards'
equation, the governing equation for infiltration of water in the vadose
zone:
97
-------
where ^ is the pressure head (L), K is the saturated hydraulic conductivity
(LT"1), krw is the relative permeability, z is the vertical coordinate
pointing in the downward direction (L) , t is time (T), and ri is an effective
water storage capacity (I/-*-) defined as
j C
AJ? (3-2a)
where Sg is specific storage (L"1), Sw is water saturation, and is the
effective porosity. Specific storage is defined by
Ss = pg[ct + (I - 4) cg] (3-2b)
where Cf is the fluid compressibility (LT M" ) , cg is the solid skeleton
o 1
M" ) g
compressibility (LT2!!"1), p is the fluid density ML"3), and g the
gravitational acceleration
The initial and boundary conditions of the one -dimensional infiltration
problem may be expressed as:
) = t± (3-3)
either
V(0,t) = I (3-4a)
or
*(0,t) = ^0 (3-4b)
either
tf(L,t) = ^ (3-5a)
or
V(L,t) = 0 (3-5b)
where V^ i-s the initial pressure head value (L) , ^o is the pressure head at
the upper boundary (L) , VT i-s tne pressure head at the lower boundary (L) , I
is the rate of infiltration at the soil surface (LT"1) , L is the thickness of
the vadose zone (L) , and V the vertical Darcy velocity (LT"1) (defined by
Equation 3-8). The boundary condition in equation 3-5b is valid because the
bottom boundary of VADOFT coincides with the bottom of the SAFTMOD aquifer
(see the section on linkage) .
98
-------
To solve the variably saturated infiltration problem, it is also necessary to
specify the relationships of relative permeability versus water saturation,
and pressure head versus water saturation. Two alternative function
expressions are used to describe the relationship of relative permeability
versus water saturation. These functions are given by (Brooks and Corey
1966; van Genuchten 1976):
krw - Sen (3-6a)
and
krw = SeV2 [1 - (1 - Se1/T)7]2 (3_6b)
where n and j are empirical parameters and Se is the effective water
saturation defined as S& = (Sw - Swr)/(l - Swr), with Swr being referred to
as the residual water saturation.
The relationship of pressure head versus water saturation is described by the
following function (van Genuchten 1976; Mualum 1976):
swr
for rj> < ^a
(3-7)
1 for ^ ^ ^_
where a and ft are empirical parameters, i/>a is the air entry pressure head
value (L), and S „ is the residual water phase saturation. The parameters ft
and 7 are related by 7 - 1 - I./ft.
Descriptive statistical values for a, ft, and 7 have been determined by Carsel
and Parrish (1987) for 12 soil classifications (see Volume II, Section 5).
Using the mean parameter values, the relationships of effective saturation
versus capillary head and relative permeability versus effective saturation
are plotted. Logarithmic plots are shown in Figures 3.2-3.4. To show more
vividly the high degree of nonlinearities, the relationships of relative
permeability versus effective saturation are also plotted on arithmetic
scales and presented in Figures 3.5-3.7. It is important that the finite
element flow module be capable of handling such high nonlinearities to be
successful in performing a Monte Carlo study of infiltration in the
unsaturated zone.
Equation 3-1 is solved using the Galerkin finite element subject to the
initial and boundary conditions given in Equations 3-3 through 3-5. After
the distributions of T/I and S have been determined, the Darcy velocity is
computed from
V - - Kkrw( - 1)
-------
Id1 10? I03
CAPILLARY HEAD, cm
(a)
SATURATION
(b)
Figure 3.2.
Logarithmic plot of constitutive relations for
clay, clay loam, loam, and loamy sand: (a)
saturation vs. capillary head and (b) relative
permeability vs. saturation.
100
-------
2
QC
(0
SILTY CLAY
SILTY CLAY LOAM
SILT
SILT LOAM
JO'-
102
103
CAPILLARY HEAD, cm
(a)
s
i
UJ
tu
OC
to-? -
10-3 -1
Figure 3.3.
10'- 10'^ JO"1 ! 0
SATURATION
(b)
Logarithmic plot of constitutive relations for silt,
silty clay loam, silty clay and silty loam: (a)
saturation vs. capillary head and (b) relative
permeability vs. saturation.
101
-------
LU
OC 10
101 10J 103
CAPILLARY HEAD, cm
(a)
SATURATION
(b)
Figure 3.4. Logarithmic plot of constitutive relations for
sandy clay, sandy clay loam, sandy loam and sand:
(a) saturation vs. capillary head and (b) relative
permeability vs. saturation.
102
-------
CO
LU
5
LT
LU
LU
_J
LU
tr
Figure 3.5.
0.4 0.6
SATURATION
Standard plot of relative permeability
vs . saturation for clay/ clay loam,
loam and loamy sand.
3.3.2 Numerical Solution
3.3.2.1 Numerical Approximation of the Flow Equation--
A numerical approximation of the one-dimensional flow equation in the vadose
zone is obtained using a Galerkin finite-element formulation with spatial
discretization performed using linear elements. Time integration is
performed using a backward finite difference approximation. This leads to a
system of nonlinear algebraic equations. For a typical node i in the finite-
element grid (see Figure 3.8), the equation may be expressed as
k+l
.k+1
.
k+l
(3-9)
where k+1 is the current time level, and a,
and di are given by
-i-l
AZJ
(3-10a)
103
-------
{Kkrw>i-l { {Kkrw>i
Az^_^ Az^
-{Kkrw)i
>ifAzi-l + Azil
^l 5 J
(3-10b)
(3-10c)
Azi
{Kkrw>i-l -
-------
1.0
0.8 -
5 o.e •-
m
(X
LLJ
°-
LU
LU
0.4 ••
0.2 ••
0.0
SAND
SANDY LOAM
SANDY CLAY LOAM
SANDY CLAY
0.0
0.2
0.4 0.6
SATURATION
0.8
Figure 3.7.
Standard plot of relative permeability
vs. saturation for sandy clay, sandy
clay loam, sandy loam and sand.
In the Picard scheme, the matrix coefficients, a^, ;?£, 7^, and d^, are first
evaluated using an initial estimate of pressure head values, i/y^. The
resulting system of linearized equations is then solved for rf& using the
Thomas algorithm. Updating of the matrix coefficient is performed by
recomputing values of nonlinear soil parameters. Iterations are performed
until the successive change in pressure head values is within a prescribed
tolerance.
In the Newton-Raphson scheme, the nonlinear system of equations is treated by
applying the Newton-Raphson technique (see Huyakorn and Finder 1983, pp. 159-
162) to Equation 3-9. This leads to the following system of linearized
algebraic equations:
,o
081
(3-11)
where superscript r is used to denote the r-th iterate; a^, 0^,
are as defined previously; and a^, f)^, 7^, and d^ are given by
and
105
-------
-------
The initial solution and subsequent iterations of the Newton-Raphson scheme
are performed in the same manner as that described for the Picard scheme.
3.3.2,2 General Guidance on Selection of Grid Spacings and Time Steps, and
the Use of Solution Algorithms--
In designing a finite-element grid for variably saturated flow simulations,
one should select nodal spacings that will yield reasonable approximations to
the expected moisture profiles.
In the analysis of the given variably saturated flow problem, small nodal
spacings should be used in the zones where head gradients or moisture fronts
are steep. The nodal spacings may be gradually increased in the zone where
no abrupt changes in hydraulic conductivities occur and the head gradients
are gradually sloping. The variably saturated flow simulation can be
performed using either the Picard or either Newton-Raphson solution
algorithm. For one-dimensional cases where convergence difficulties are not
expected, the efficiency of these algorithms have been found to be similar.
For certain steady-state cases involving highly nonlinear soil moisture
characteristics, the use of either Newton-Raphson algorithm is preferable,
particularly when the Picard algorithm fails to converge within a reasonable
number of iterations (say between 10 and 20).
3.4 DESCRIPTION OF THE TRANSPORT MODULE
3.4.1 Transport Equation
The governing equation for one-dimensional transport of a nonconservative
solute species in a variably saturated soil takes the form
Si
9 1
where D is the apparent dispersion coefficient (L T ) , c is the solute
concentration (ML ) , 9 is the volumetric water content (0 = ^Sw) , R is the
retardation coefficient, and A is the first-order decay constant (T ) . Note
that the apparent dispersion coefficient is defined as D = a-^V + <^D , where
a-r is the longitudinal dispersivity , and D is the effective molecular
diffusion coefficient.
The initial and boundary conditions of the one -dimensional transport problem
may be expressed as:
c(z,0) = GI (3-14)
either
-D (0,t) - V (CQ - c) (3-15a)
107
-------
or
c(0,t) - CQ (3-15b)
(L,t) - 0 (3-16)
where c^ is the initial concentration (ML~), and CQ is the leachate
concentration at the source (ML ) .
*
3.4.2 Numerical Solution of the Transport Equation
3.4.2.1 Numerical Approximation of the Transport Equation- -
A numerical approximation of the one -dimensional transport equation is
obtained using an upstream-weighted finite -element formulation with spatial
discretization performed using linear elements. Time integration is
performed using a central finite -difference approximation. This leads to a
system of linear algebraic equations. The equation corresponding to node i
takes the form:
k+1 . 0 k+1 . k+1 . ,,
- di <
where a^, fi^, 7^, and d^ are given by
*
at - ra + (flR}^ Azi.1/(6Atk) (3-18a)
) (3-18b)
L AZi/(6Atk) (3-18c)
- (r - l)(a c.! + j9 c + -y
JT (3-18e)
A
•o I l?K.Ji_l ^Zi-l +
(1 -_»). MTI. -L . C1 + "> {VJj^ (3-18f)
L-^fi (3-18e)
o
108
-------
with T and w denoting the time weighting factor and the upstream weighting
factor, respectively.
To obtain a second-order temporal approximation, the value of T is set equal
to 1/2. This corresponds to using the Crank-Nicholson central difference
time stepping scheme. The upstream weighting factor co is introduced in the
above numerical approximation to curb numerical oscillations that may occur
when the selected finite-element grid is not sufficiently refined for a given
value of longitudinal dispersivity. For each time step, the linear system of
algebraic equations is solved using the Thomas algorithm.
Transport of a daughter species in a decay chain can also be handled by the
VADOFT code. In this case, the right hand side of the governing equation for
single species transport (3-13) is modified by adding a source term
accounting for transformation of parent components. This source term is
given by
m = -S Sw e^ Xji R^ c^ (3-19)
where subscript S. is used to denote the parent species, n^ is the number of
parent species, and e^ is the mass fraction of parent component that is
transformed into the daughter species under consideration.
The numerical solution of the modified transport equation can be performed in
the same manner as that described previously for a single species. The
source term from Equation 3-19 is incorporated into the finite element matrix
equation by adding d^ to the right hand side. The term d^ is given by
* 1 "p k+1
dj = - [ S e, A, Rj eg? L] (Az, + AZii)^ S )± (3-20)
2 je-1 1
In performing the solute transport analysis, the selection of nodal spacing
(Az) and time step value (At) should follow the so-called Peclet number and
Courant number criteria where possible. These two criteria are given as
follows:
<4 (3-21)
VSQl At/Az :< 1 (3-22)
Vsol - V/*R (3-23)
where «L is the longitudinal dispersivity, Vgo^ is the solute velocity, V
Darcy velocity, S - water content, and R = retardation coefficient.
109
-------
The VADOFT code also provides the user with the option of using upstream
weighting to curb numerical oscillations that may occur in solving the
advective- dispersive transport equation. The recommended value of w, the
weighing factor, is determined by using the following formulae:
w - 1 - 4aL/^, £ > 4aL (3-24)
w - 0, i. < 4aL (3-25)
where c*L is the longitudinal dispersivity, and £ is the length of the
element .
3.5 RESULTS OF VADOFT TESTING SIMULATIONS
Three sets of benchmark problems were used to test the VADOFT code. The
first set consists of two steady and transient problems designed to test the
variably saturated flow component of the code. The second set consists of
four transient one -dimensional transport problems. The third set consists of
two coupled flow- transport problems. Numerical results obtained from VADOFT
are compared with analytical solutions and results obtained using two other
finite -element codes, UNSAT2 and SATURN. These test problems were simulated
using VADOFT before it was linked to the other RUSTIC modules.
3.5.1 Flow Module (Variably Saturated Flow Problems)
3.5.1.1 Transient Upward Flow in a Soil Column- -
This problem concerns transient, vertically upward moisture movement in a
soil column 20 cm in length. The soil column is subject to zero pressure
head at the base and zero flux at the top. The initial distribution of
pressure head is hydrostatic: (t = 0) = -90 + z cm, where z is the depth
below the top of the soil column. Soil properties and discretization data
used in the simulation are presented in Table 3-1. The simulation was
performed for 15 time steps with constant time step value of t = 0.01 d.
Numerical results given by the Picard and the Newton -Raphs on schemes are
virtually identical. Both schemes require between 2 and 3 iterations per
time step to converge to a head tolerance of 0.01 cm. The simulation results
obtained from VADOFT are compared with those obtained from UNSAT2 and SATURN
(the two-dimensional finite -element codes described by Davis and Neuman
[1983], and Huyakorn et al. [1984]) respectively. Shown in Figures 3.9 and
3.10 are plots of distributions of pressure head and water saturation,
respectively. As can be seen, the results of VADOFT are in good agreement
with the results of the other two codes.
3.5.1.2 Steady Infiltration in a Soil Column- -
This problem concerns steady-state infiltration in a soil column. The column
is 550 cm in length and is subject to an infiltration rate of 4.07 cm/d at
the top and zero pressure head at the bottom. Soil properties used in the
simulation are presented in Table 3-2. Five cases of varying degree of
nonlinearity of relative permeability function (krw = Sen) were simulated.
Both the Picard and the Newton- Raphson schemes were used in conjunction with
110
-------
Table 3-1. SOIL PROPERTIES AND DISCRETIZATION DATA USED IN SIMULATING
TRANSIENT FLOW IN A SOIL COLUMN
Parameter Value
Length of soil column, L 20 cm
Saturated hydraulic conductivity, K 10 cm d
Porosity, 0.45
Residual water phase saturation, Swr 0.333
Air entry value, \f>a 0.0 cm
Constitutive relations:
krw -
<* - *a)/(*r - ^a> -
where ^r — -100 cm.
Az - 0.5 cm
At - 0.01 d
a finite-element grid having constant nodal spacing, z - 10 cm. The
performance of the two iterative schemes are illustrated in Table 3-3. Note
that the Newton-Raphson scheme converges for all cases, whereas the Picard
scheme fails to converge when the nonlinear exponent n exceeds 4. Simulated
distributions of pressure head and water saturation are shown in Figure 3.11
and 3.12, respectively. These results of'the VADOFT code are virtually
identical to corresponding results obtained using the SATURN code.
3.5.2 Transport Module
3.5.2.1 Transport in a Semi-Infinite Soil Column--
This problem concerns one-dimensional transport of a conservative solute
species in a saturated soil column of infinite length. The solute is
introduced into the column at the inlet section where z - 0. The initial
concentration is assumed to be zero, and the dimensionless constant inlet
concentration is prescribed as 1. Values of physical parameters and
discretization data used in the numerical simulation are given in Table 3-4.
The finite-element grid representing the soil column was 400 cm in length.
The simulation was performed for 20 time steps. Thus the duration of the
simulation time of transport in the soil column was 50 hours. For this
duration, the selected grid length is sufficient to avoid the end boundary
effect. The numerical solution obtained from the VADOFT code was checked
against the analytical solution of Ogata and Banks (1961). Shown in Figure
3.13 and Table 3-5 are concentration values at t = 25 hours and t = 50 hours.
As can be seen, the numerical and analytical solutions are in excellent
agreement.
Ill
-------
Table 3-2. SOIL PROPERTIES USED IN SIMULATING STEADY-STATE INFILTRATION
Parameter Value
Length of soil column, L 550 cm
Saturated hydraulic conductivity, K 25 cm d
Porosity, 0.331
Residual water saturation, ST7T. 0.0
W L
Air entry value, TJ>a 0 .0 cm
Constitutive relations:
krw = Sen, n = 3, 4, 6, 8, 10
se -_ 1 —
where Se - (Sw - Swr)/(l - Swr), a - 0.014 cm'1, ^a - 0 cm,
/3 - 1.51, 7 - 0.338
3.5.2.2 Transport in a Finite Soil Column--
In this problem, downward vertical transport of dissolved contaminants in a
soil column above the water table of an unconfined aquifer is considered.
The length of the soil column is 20 m and the Darcy velocity and water
content are assumed to be constant and equal to 0.25 m/d and 0.25,
respectively. The initial concentration is zero, and water with
dimensionless solute concentration of 1 enters the soil surface at a rate of
0.25 m/d. At the water table, a zero dispersive-flux boundary condition is
assumed. A list of physical parameter values and discretization data used in
the simulation is provided in Table 3-6. Two cases involving conservative
and nonconservative species were simulated. Results obtained from the VADOFT
code are compared in Figure 3.14 and Table 3-7 with the analytical solution
given by van Genuchten and Alves (1982). There is excellent agreement
between the numerical and analytical solutions for both cases.
3.5.2.3 Transport in a Layered Soil Column--
This problem concerns one-dimensional transport of a conservative solute
species in a soil column consisting of thr.ee layers. The initial
concentration in the soil column is assumed to be zero, and the two boundary
conditions prescribed are a unit concentration at the top and a zero
dispersive flux boundary condition at the bottom. A list of physical
parameter values and discretization data used in the simulation is provided
112
-------
20
16
12
E
u
UNSAT2
© VADOFT
X SATURN
-20 -40 -00
PRESSURE HEAD.cm
-BO
-100
Figure 3.9. Simulated pressure head profiles for the problem of
transient upward flow in a soil column. (Adapted
from Battelle and GeoTrans, 1988).
in Table 3-8. Two cases corresponding to those considered by Shamir and
Harleman (1967) were simulated. Both cases have contrasting longitudinal
dispersivity values among the three layers. The dispersivity values of the
second case are ten times those of the first case for the same layers. The
intention here is to test the numerical scheme used in the VADOFT code, as
well as to check the validity of an approximate analytical solution presented
by Shamir and Harleman (1967) and Hadermann (1980). It should be noted here
that the approximate solutions by Shamir and Harleman (1967) and Hadermann
(1980) are valid only for relatively small values of dispersivity.
Therefore, for a small dispersivity value, the solutions can be employed to
verify the VADOFT code. Then with appropriate discretization, the VADOFT
code could be used to determine the validity of the analytical solutions at
large dispersivity values.
113
-------
t>
UNSAT2
O VADOFT
X SATURN
0.4 0.6
WATER PHASE SATURATION
1.0
Figure 3.10. Simulated profile of water saturation for the problem
of transient upward flow in a soil column. (Adapted
from Battelle and GeoTrans, 1988).
Table 3-3. ITERATIVE PROCEDURE PERFORMANCE COMPARISON
Case
Number of Nonlinear Iterations
Newton-
Raphson Picard
n
n
n
n
n
3
4
6
8
10
12
13
19
27
31
33
56
n.c.*
n.c.
n.c.
No convergence. Head tolerance = 0.0001 cm. Grid spacing z - 10 cm.
114
-------
Table 3-4. VALUES OF PHYSICAL PARAMETERS AND DISCRETIZATION DATA USED IN
SIMULATING ONE-DIMENSIONAL TRANSPORT IN A SEMI-INFINITE SOIL
COLUMN
Parameter Value
Darcy velocity, V 1 cm hr
Porosity, 0.25
Longitudinal dispersivity, a^ 5 cm
Concentration at the source, CQ 1
Az = 10 era
At = 2.5 hr
Using the discretization data given in Table 3-8, the VADOFT code was run for
180 time steps. Simulated breakthrough curves at the bottom end of the
column (z = 86.1 cm) are presented in Figures 3.15 and 3.16 and in Tables 3-9
and 3-10. As can be seen, the'numerical solution of the VADOFT code compares
very well with the analytical solution for case 1: The small dispersivity
case, where the analytical assumption of infinite ratio of layer thickness to
layer dispersivity--i. e. , each layer extends to infinity—is fairly accurate.
There is a slight discrepancy of the analytical solution from the numerical
solution for case 2, where the analytical assumption is less accurate.
3.5.3 Combined Nonlinear Flow and Transport Modules
3.5.3.1 Transport During Absorption of Water in a Soil Tube--
This problem is selected to provide simultaneous testing of the flow and the
transport modules of VADOFT. The problem is depicted schematically in Figure
3.17. A conservative solute species has a uniform initial concentration and
moisture content. The initial concentration is assumed to be zero, and the
inlet concentration CQ is assumed to be 1 ppm. The solute is transported by
dispersion and advection. Note that the solute front and the wetting front
advance at different rates. The solute velocity, VSQ-p was previously
defined as Equation 3-23. The velocity of the wetting front is dependent
upon the rate of water sorption into the soil, which is dependent on moisture
diffusivity; thus, calculation of the wetting front velocity requires
integration of the mass balance equation. For the sake of convenience, all
physical data pertaining to the geometry of the soil tube and the physical
parameter values are kept the same as those used in the paper by Huyakorn et
al. (1985). The complete set of data is listed in Table 3-11. The
simulation was performed in two stages. In the first stage, the transient
water flow problem was analyzed to determine the distributions of Darcy
velocity and water saturation for each time level. These results are written
on an output file. In the second stage, the transient solute transport
115
-------
problem was analyzed to determine concentration distributions using the
velocity and water saturation data file obtained from the flow simulation.
The spatial and temporal discretization data used in running the VADOFT code
are also given in Table 3-11. Both the flow and the transport analyses were
performed for 50 time steps. Results of the flow analysis are plotted in
Figure 3.18. The water saturation profiles given by VADOFT compare well with
those obtained using the semi-analytical solution of Phillip (1955) and the
UNSAT2 finite-element flow code. Results of the transport analysis are
ICC -
too
E
u
400 -
too -
•oo
too -
-•c
PRESSURE HEAD (cm)
Figure 3.11.
Simulated pressure head profiles for five cases
of the problem of steady infiltration in a soil
column. (Adapted from Springer and Fuentes, 1987)
116
-------
plotted in Figure 3.19. The concentration distributions given by VADOFT also
compare well with those obtained using the semi-analytical solution of Smiles
et al. (1978) and the FEMWASTE finite-element transport code documented by
Yeh and Ward (1981).
3.5.3.2 Transient Infiltration and Contaminant Transport in the Vadose Zone-
This problem, schematically depicted in Figure 3.20, involves variable
infiltration and contaminant transport in a layered system in which layer
100
too
E
I
fr-
it*
o
900
400
500
A n
•f n
x n
v n
e n
- 10
• 8
• 6
• 4
• 3
o.o
.2
1.0
SATURATION
Figure 3.12.
Simulated profiles of water saturation for five cases
of the problem of steady infiltration in a soil column.
(Adapted from Springer and Fuentes, 1987).
117
-------
permeabilities differ by more than two orders of magnitude. The problem was
chosen to demonstrate the capability of VADOFT to handle a higher nonlinear
situation involving soil materials with sharp contrast in drainage
properties. Shown in Table 3-12 are values of physical parameters and
discretization data used in the flow and transport simulations. For the
unsaturated flow simulation, the transient infiltration rates illustrated in
Figure 3.21 were used. It was assumed that the initial condition
corresponded to a hydrostatic pressure head distribution in the soil with
pressure head values at the water table and the top of the soil equal to 0
and -420 cm, respectively. The simulation was performed for 20 time steps
1.0
.£ -
c
o
•^
c '6
0}
o
c
o
OJ
o:
.2 --
0.0
Analytic Soln
O VADOFT
150. 200. 250.
Distance, cm
400.
Figure 3.13. Simulated concentration profiles for the problem of
solute transport in a semi-infinite soil column.
118
-------
Table 3-5. CONCENTRATION PROFILE CURVES AT t = 25 hr AND t = 50 hr
SHOWING COMPARISON OF THE ANALYTICAL SOLUTION AND RESULTS
FROM VADOFT
Concentration Values
z
Distance
(cm)
00.0
10.0
20.0
30.0
40.0
50.0
60.0
70.0
80.0
90.0
100.0
110.0
120.0
130.0
140.0
150.0
160.0
170.0
180.0
190.0
200.0
210.0
220.0
230.0
240.0
250.0
260.0
270.0
280.0
290.0
300.0
310.0
320.0
330.0
t =
Analytical
1.0000
0.9997
0.9983
0.9945
0.9854
0.9662
0.9313
0.8745
0.7924
0.6858
0.5619
0.4321
0.3099
0.2060
0.1264
0.0713
0.0369
0.0175
0.0075
0.0030
0.0011
0.0003
0.0000
25 hr
VADOFT
1.0000
0.9998
0.9987
0.9954
0.9870
0.9688
0.9346
0.8781
0.7956
0.6889
0.5660
0.4394
0.3222
0.2235
0.1474
0.0928
0.0560
0.0327
0.0184
0.0101
0.0054
0.0029
0.0015
t =
Analytical
1.0000
1.0000
1 . 0000
1 . 0000
0.9999
0.9999
0.9996
0.9991
0.9981
0.9960
0.9921
0.9854
0.9743
0.9570
0.9313
0.8953
0.8475
0.7872
0.7151
0.6331
0.5447
0.4541
0.3660
0.2845
0.2129
0.1532
0.1058
0.0701
0 . 0444
0.0270
0.0157
0.0087
0.0046
0.0000
50 hr
VADOFT
1.0000
1.0000
1 . 0000
1 . 0000
1.0000
0.9999
0.9997
0.9994
0.9985
0.9967
0.9933
0.9871
0.9767
0.9599
0.9348
0.8991
0.8513
0.7908
0.7186
0.6368
0.5491
0.4598
0.3736
0.2942
0.2246
0.1662
0.1193
0.0831
0.0563
0.0371
0.0239
0.0150
0.0092
0.0055
119
-------
e.e
4.
8. 12.
Distance, m
(a) X • 0 d'1
e. 12.
Distance, in
(b) X - 0.25 d
-1
2C.
— Analytic Soln.
e VADOFT
2E.
Figure 3.14.
Simulated concentration profiles for two cases of the
problem of solute transport in a soil column of finite
length, (a) X = 0 d"1, and (b) X = Q.25 d"1.
120
-------
Table 3-6. VALUES OF PHYSICAL PARAMETERS AND DISCRETIZATION DATA USED
IN SIMULATING ONE-DIMENSIONAL TRANSPORT IN A FINITE SOIL COLUMN
Parameter Value
Thickness of soil column, L 20 m
Darcy velocity, V 0.25 m d
Water content, 6 0.25
Retardation coefficient, R 1
Longitudinal dispersivity, a^ 4m
Source leachate concentration, CQ 1
Case 1:
Decay constant, A 0 d
Case 2:
Decay constant, A 0.25 d"^
Az = 1.0 m
At - 0.5 d
using At - 1 d. Shown in Figures 3.22 through 3.24 are simulated profiles of
water saturation, pressure head, and vertical Darcy velocity, respectively.
As expected, the two sand layers exhibit fast drainage response whereas the
intervening clay-loam layer exhibits slow drainage response. This behavior
is seen in Figure 3.22. The pressure head and velocity profiles depicted in
Figures 3.23 and 3.24 directly reflect the effect of temporal change in the
infiltration rate. Note that the values of Darcy velocity at the soil
surface (Figure 3.24) are equal to the values of infiltration rate for the
same time values. Following the unsaturated flow simulation, the transport
simulation was performed using the Darcy velocity file from the flow
computation as an input file for the transport computation. Concentration
profiles determined by the code are plotted in Figure 3.25. As illustrated,
the contaminant front exhibits slow movement through the clay loam layer.
121
-------
Table 3-7. CONCENTRATION PROFILE CURVES SHOWING COMPARISON OF THE
ANALYTICAL SOLUTION AND VADOFT
Distance
z, (m)
0.0
2.0
4.0
6.0
8.0
10.0
12.0
14.0
16.0
18.0
20.0
t =
Analytical
0.764
0.638
0.502
0.371
0.256
0.164
0.097
0 053
0.027
0.013
0.009
5 d
VADOFT
0.751
0.624
0.489
0.360
0.247
0.158
0.094
0.052
0.027
0.014
0.009
Case 1: A =
t = 10
Analytical
0.884
0.820
0.742
0.655
0.561
0.466
0.375
0.293
0.224
0.176
0.157
• 0 cT1
d
VADOFT
0.878
0.812
0.733
0.645
0.552
0.457
0.367
0.286
0.219
0.171
0.152
t = 20
Analytical
0.963
0.942
0.914
0.881
0.841
0.796
0.748
0.698
0.652
0.617
0.602
d
VADOFT
0.961
0.939
0.911
0.877
0.837
0.791
0.742
0.692
0.646
0.610
0.595
Case 2: A = 0.25 d'1
Distance
z, (m)
0.0
2.0
4.0
6.0
8.0
10.0
12.0
14.0
16.0
18.0
20.0
t =
Analytical
0.593
0.416
0.283
0.186
0.116
0.069
0.038
0.020
0.009
0.004
0.002
5 d
VADOFT
0.588
0.411
0.279
0.182
0.113
0.067
0.037
0.019
0.009
0.004
0.002
t = 10
Analytical
0.615
0.449
0.326
0.236
0.169
0.119
0.083
0.057
0.039
0.028
0.024
d
VADOFT
0.613
0.447
0.325
0.234
0.167
0.118
0.083
0.057
0.039
0.028
0.024
t - 20
Analytical
0.618
0.453
0.333
0.244
0.179
0.131
0.096
0.071
0.053
0.042
0.038
d
VADOFT
0.617
0.452
0.332
0.243
0.178
0.131
0.096
0.071
0.053
0.042
0.038
122
-------
1.0
c
o
c
QJ
o
c
o
QJ
(X
CASE 1
—Analytic Soln
O VADOFT
0.0
600.
650.
700. 750.
Time, s
Figure 3.15. Simulated outflow breakthrough curve for case 1 of the
problem of solute transport in a layered soil column.
123
-------
1.0
.S -•
.6 ••
c
o
fO
s_
c
O)
c
o
a>
a: .2 •
0.0 <•»
CASE 2
A Analytic Solution
O VADOFT
400. 500. 600. 702.
Time, s
900.
1000.
Figure 3.16. Simulated outflow breakthrough curve for case 2 of the
problem of solute transport in a layered soil column.
Table 3-8. VALUES OF PHYSICAL PARAMETERS USED IN THE SIMULATION OF
TRANSPORT IN A LAYERED SOIL COLUMN
Parameter
Layer thickness, t^
Seepage velocity, u^
Retardation coeff . , R£
Decay constant, A^
Source concentration, CQ
Case 1:
Dispersivity, OL^
Case 2:
Dispersivity, a^
Az - 0.6888 cm
At - 5 s
Layer 1
25. 48
0.127
1.0
0
1.0
0.076
0.76
Value for Laver 1
Layer 2
30.31
0.123
1.0
0
0.174
1.74
Layer 3
30.31 cm
0.121 cm s'1
1.0
0 s'1
0.436 cm
4.36 cm
124
-------
Table 3-9. BREAKTHROUGH CURVES (at z - 86.1 cm) COMPUTED USING
THE ANALYTICAL SOLUTION AND VADOFT: CASE 1
Time, t (s)
600
610
620
630
640
650
660
670
680
690
700
710
720
730
740
750
760
770
780
790
800
810
820
830
840
850
Concentration Values
Analytical Solution
0.0204
0.0361
0.0596
0.0923
0.1354
0.1887
0.2514
0.3217
0.3971
0.4748
0.5518
0.6255
0.6935
0.7544
0.8072
0.8517
0.8881
0.9172
0.9400
0.9573
0.9704
0.9800
0.9870
0.9919
0.9950
0.9970
for Case 1
Numerical VADOFT
0.0262
0.0427
0.0665
0.0989
0.1410
0.1930
0.2543
0.3234
0.3981
0.4755
0.5526
0.6266
0.6951
0.7564
0.8096
0.8542
0.8907
0.9197
0.9421
0.9590
0.9715
0.9805
0.9869
0.9913
0.9943
0.9964
125
-------
Table 3-10. BREAKTHROUGH CURVES (at z - 86.1 cm) COMPUTED USING THE
ANALYTICAL SOLUTION AND VADOFT= CASE 2
Time, t (s)
600
610
620
630
640
650
660
670
680
690
700
710
720
730
740
750
760
770
780
790
800
810
820
830
840
850
900
Concentration Values
Analytical Solution
0.303
0.330
0.357
0.384
0.412
0.439
0.466
0.493
0.519
0.544
0.569
0.593
0.617
0.639
0.661
0.681
0.701
0.720
0.738
0.755
0.771
0.787
0.801
0.815
0.828
0.840
0.889
for Case 2
Numerical VADOFT
0.310
0.337
0.365
0.394
0.422
0.450
0.478
0.505
0.532
0.558
0.584
0.608
0.632
0.655
0.677
0.698
0.718
0.737
0.755
0.772
0.788
0.804
0.818
0.831
0.844
0.856
0.904
SOLUTE
FRONT
I j
xxx
WETTING
FRONT
c = 1 ppm
1
1
1
1
'/////
FLOW
-- ^
=-83.33 cm
c or 8
Figure 3.17.
One-dimensional solute transport during absorption
of water in a soil tube. (Adapted from Huyakorn
et al., 1985).
126
-------
Table 3-11. VALUES OF PHYSICAL PARAMETERS AND DISCRETIZATION DATA USED
IN SIMULATING TRANSPORT IN A VARIABLY SATURATED SOIL TUBE
Parameter
Value
Length of soil column, L
Saturated hydraulic conductivity, K
Initial pressure head, i/>^
Remaining flow parameters
Initial concentration, c^
Longitudinal dispersivity, a^
Molecular diffusion, D
Decay constant, A
Retardation coefficient, R
20 cm
1 cm d'1
83.33 cm
See Table 3-2
0 ppm
0 cm
2 '1
1
0
1
cm d
'1
Az = 0.25 cm
At = 0.0025 d
Table 3-12. VALUES OF PHYSICAL PARAMETERS AND DISCRETIZATION DATA USED IN
SIMULATING TRANSIENT INFILTRATION AND CONTAMINANT TRANSPORT IN
THE VADOSE ZONE
Property
Saturated conductivity, K
Porosity, 4>
Residual Water Saturation, S
Air entry value , i/>a
Soil moisture parameter, a
Soil moisture parameter, /?
Soil moisture parameter, 7
Longitudinal dispersivity, a-r
Retardation coefficient, R
Decay coefficient, A
Az = 4 cm
At - 1 d
Material 1
(Sand)
713
0.43
0.105
0.0
0.145
2.68
0.63
1.0
1.1
0.00274
Material 2
(clay loam)
6.24 cm d'1
0.41
0.232
0.0 cm
0.019 cm"1
1.31
0.24
1.0 cm
1.5
0.00274 d'1
127
-------
1.0
VADOFT
UNSAT2
SEMI-ANALYTIC SOLUTION
0.2
Figure 3.18. Simulated profiles of water saturation during
absorption of water in a soil tube. (Adapted
from Huyakorn et al., 1984a).
128
-------
u
(J
z"
o
tr
t-
r
u
u
z
o
u
lu
UJ
c
1.00
0.75 -
0.50
0.25 -
VADOFT
—*•— FEMWASTE
SEMI-ANALYTIC
SOLUTION
Figure 3.19.
Simulated concentration profiles for the problem of
one-dimensional solute transport during absorption
of water in a soil tube. (Adapted from Huyakorn,
et al., 1985).
129
-------
20 cm.
120 cm
280 cm
l,cm/d
I I I I I I I
SAND Ksat= 713 cm/d
CLAY LOAM
Ksat= 6-24 cm/d
SAND
Ksat= 713 cm/d
WATER TABLE „
Figure 3.20,
Problem description for transient infiltration
and transport in the vadose zone.
5
4 •
5 3-
E
o
c
£ 2
(0
L.
4^
»^-
C
1
0
L_T
4 8 12 16 20
time , days
Figure 3.21
Infiltration rate vs. time relationship
used in numerical simulation.
130
-------
420
0.4 0.6
SATURATION
Figure 3.22. Simulated water saturation profiles.
131
-------
420f
0
-100 -200 -300
PRESSURE HEAD, cm
-400
Figure 3.23. Simulated pressure head profiles.
132
-------
420
2 3
VELOCITY, cm/d
Figure 3.24. Simulated vertical Darcy velocity profiles.
133
-------
E
o
Q.
LLJ
Q
100
0.2
0.4 0.6
CONCENTRATION
0.8
1.0
Figure 3.25. Simulated solute concentration profiles.
134
-------
SECTION 4
SATURATED ZONE FLOW AND TRANSPORT MODEL (SAFTMOD)
4.1 INTRODUCTION
SAFTMOD is a finite-element code for simulating groundwater flow and solute
transport in the saturated zone of an unconfined aquifer system. It is the
third part of the three-component RUSTIC model for predicting the movement of
pesticides within and below the plant root zone and assessing consequent
groundwater contamination. SAFTMOD performs two-dimensional simulations in
an areal plane or a vertical cross section. In addition, the code also
performs axisymmetric simulations. Both single and leaky two-aquifer systems
can be handled. The code considers only single-porosity media and ignores
the effects of fracture flow in the aquifer system. Transport of dissolved
contaminants may also be simulated within the same domain. Transport
processes accounted for include hydrodynamic dispersion, advection, linear
equilibrium sorption, and first-order decay. Parent/daughter transformations
are also simulated.
4.2 OVERVIEW OF SAFTMOD
The SAFTMOD code may be used in conjunction with the root zone model, PRZM,
and/or the vadose-zone model, VADOFT, to simulate groundwater flow and
transport of dissolved contaminants in the saturated zone. SAFTMOD predicts
groundwater velocity field, flow rates, concentration plumes, and temporal
distributions of contaminant concentration at downstream locations.
4.2.1 Features
4.2.1.1 General Description--
The SAFTMOD code can be used to perform two-dimensional or axisymmetric
modeling of groundwater flow and solute transport in the saturated zone of an
unconfined aquifer system. For demonstrative purposes, we consider a
situation involving the migration of pesticide in an aquifer subjected to
well pumping (Figure 4.1). SAFTMOD can be used in three ways. First, the
code can be used to perform either an axisymmetric or vertical cross-section
analysis to assess the local well flow or migration of contaminants. Second,
the code can be used to perform an areal analysis to assess flow and the
migration of contaminants on a regional scale. SAFTMOD can be operated as a
stand-alone code or operated in conjunction with PRZM and/or VADOFT. In the
latter case, PRZM and/or VADOFT need to be run to obtain the necessary
information pertaining to flux boundary conditions on the water table beneath
the unsaturated soil column.
135
-------
o
t.,
LU
Q
o
h-
co
LU
Q.
2C
D
J
_i
l.
T_
*^m
<
LU /
i o /
^ N /
^ /
Q /
LU /
h- /
< /
oc /
=>
H-
< 1
CO
-7
_JA
—— — ™ ' ~—1 "" """"
t>
/
/
/
s
1
LU
Z
o
N
Q
LU
h-
<
CC
D
<
CO
/
/
/
/
/
\
\
\
\
\
\
\
\
\
\
\
\
\
\
c
•H
C CD
O C
•H -H
4J CX
nj E
u p
•^i
E O
-P
(U
T) -P
•H U
U 0)
-P JD
W D
0) Ul
U-l 0)
O -P
V)
c >,
o w
•H
CUO)
•H
-------
4.2.1.2 Processes and Geometry--
SAFTMOD performs two-dimensional and axisymmetric analysis of saturated
groundwater flow and associated solute transport problems. The analyses can
be done in a transient or steady-state mode. The code employs the Galerkin
finite-element technique to approximate the governing equations for flow and
transport. It handles single nonleaky aquifer systems as well as leaky
aquifer systems having two aquifers separated by semi-permeable layers or
aquitards. For groundwater flow simulations, the code accommodates water
table conditions, recharge by infiltration or precipitation, and well pumping
or injection. For solute transport simulations, the code accounts for
various transport processes, including hydrodynamic dispersion, advection,
linear equilibrium sorption, and first-order decay. Steady-state transport
can not be simulated when decay is considered. Boundary conditions of the
saturated flow problem are specified in terms of prescribed hydraulic head
(defined as the sum of pressure head and elevation), or prescribed volumetric
water flux. Boundary conditions of the solute transport problem are
specified in terms of prescribed concentration or prescribed solute mass
flux. All boundary conditions can be time dependent.
4.2.1.3 Assumptions•-
The SAFTMOD code contains both saturated flow and solute transport models.
Major assumptions of the flow model are as follows:
• Darcy's law is valid and hydraulic head gradients are the only
significant driving mechanism for fluid flow.
• The porosity and hydraulic conductivity are constant with time.
• Gradients of fluid density, viscosity, and temperature do not affect
the velocity distribution.
• For areal flow simulation in two-aquifer systems, vertical leakage
in aquitards is assumed.
Major assumptions of the transport model are as follows:
• Fluid properties are independent of concentrations of contaminants.
• Contaminants are miscible with in-place fluids.
• For areal transport simulations, advection in aquitards is assumed
to be negligible.
• Diffusive/dispersive transport in the porous-medium system is
governed by Pick's law. The hydrodynamic dispersion coefficient is
defined as the sum of the coefficients of mechanical dispersion and
molecular diffusion.
• Adsorption and decay of the solute may be described by a linear
equilibrium isotherm and a first-order decay constant.
137
-------
4.2.2 Limitations
Major limitations of the SAFTMOD are summarized as follows:
• The code simulates only single-phase water flow and solute transport
in saturated porous media. It does not consider unsaturated or
fractured media.
• Non-Darcy flow that may occur near pumping wells is neglected.
• Only two-dimensional or quasi -three -dimensional problems can be
simulated.
• The code does not take into account sorption nonlinearity or kinetic
sorption effects which, in some instances, can be important.
4.3 DESCRIPTION OF FLOW MODULE
4.3.1 Flow Equation
The form of the flow equation used in SAFTMOD depends on the type of
analysis. In a vertical cross -sectional analysis, the following equation is
used:
ah a ( ah ah
where x and z are the horizontal and vertical coordinates, h is the hydraulic
head (L), KXX and KZZ are principal components of hydraulic conductivity
tensor in the x and z directions, respectively (LT ), Sg is the specific
storage (L ), t is time (T), and q is the volumetric flow rate of sources or
sinks per unit volume of the porous medium (L T L ).
In an axisymmetric analysis, the following equation is used:
ahl j. d (v 3h
^} + ^ [Kzz 3z
c ah ,/ ON
Ss at ' q (4'2)
where r and z are the radial and vertical coordinates of the cylindrical
coordinate system, and Krr and KZZ are the radial and vertical components of
the hydraulic conductivity tensor (LT ) .
In an areal analysis, flow in a horizontal plane (x-y) of each aquifer is
described by the following vertically averaged equation:
„ . a fv „ ah"] ah * // ^
H + ay" [*yy H a7J = * a^ ' q ' 7 (4'3)
138
-------
where q is the volumetric flow rate of sources or sinks per unit area (L T"
L"^), 7 is the vertical leakage flux from aquitard to aquifer (LT"1), H is
the saturated thickness (L) , and ij is a storage coefficient defined as
S H for a confined aquifer
(4-4)
S + Ss H for an unconfined aquifer, where S » SSH
in which Sy is the specific yield of the aquifer.
The use of equations 4-1, 4-2, and 4-3 implies that the principal directions
of anisotropy coincide with the grid coordinate axes.
In the areal analysis, flow in aquitards is described approximately by the
following one -dimensional vertical leakage equation:
Ji *• - *
-- ~~i S ~~ ~ -* ~ -* * ~ *•
storage (L ) of the aquitard, and h' is the hydraulic head in the aquitard.
where K' and Sg are the vertical hydraulic conductivity (LT ) and specific
The initial and boundary conditions of the groundwater flow problem may be
expressed as
h(xifO) - h0(Xi) (4-6)
h(xi,t) = h on B1 (4-7a)
Vini = Vn on B2 (4-7b)
where hQ is the initial head (L), h is the prescribed head (L) on boundary
portion B-^, n^ is the outward unit normal vector on boundary portion ~S>^,
where the fluid flux is prescribed as -Vn (LT ). The sign convention for V
is positive for inflow and negative for outflow, x^ (i=l,2) are the
coordinates that correspond to x and z for cross-sectional, x and r for
axisymmetric, and x and y for areal simulation.
The SAFTMOD code performs a finite element solution of the stated flow
problem, and determines values of the hydraulic head at nodal points of a
two-dimensional finite-element grid. After the head distribution has been
determined, values of Darcy velocity components at element centroids are
evaluated using Darcy's law, which may be expressed as
139
-------
vx • -KXX
Vy - -Kyy (4-8a)
Equations (4-8) and (4-8a) are used in the areal analysis. They are also
applicable in the vertical cross -sectional and axisynunetric analyses if the
coordinate system (x,y) is replaced by (x,z) and (r,z), respectively.
4.3.2 Numerical Solution of Flow Problems
Numerical approximations of the groundwater flow equations describing cross -
sectional, axisymmetric, and areal problems are obtained using the Galerkin
finite element technique. In the Galerkin finite element approximation
procedure, the flow region is first discretized into a network of finite
elements , and an interpolating trial function is used to represent the
unknown function (hydraulic head) over the discretized region. An integral
approximation of the flow equation is then obtained using the Galerkin
weighted residual criterion. Spatial integration is performed piecewise over
each element. Upon assembly of the elements and incorporation of boundary
conditions, a system of nodal equations is obtained. For a steady-state
simulation, the nodal equations are algebraic equations. For a transient
simulation, the nodal equations are first-order time differential equations
that can be integrated using a finite difference approximation. For each
time step, this gives rise to a system of algebraic equations that are solved
using a direct Gaussian elimination procedure.
4.3.2.1 Cross -sectional and Axisymmetric Flow --
For a cross -sectional or an axisynunetric simulation, the flow system is
modeled as one aquifer unit discretized into a number of elements. (Simple
rectangular plane elements or axisymetric rectangular ring elements are used
in SAFTMOD.) Heterogeneity of the flow system is treated by assigning
different sets of hydraulic properties to individual elements. (A detailed
description of the finite element technique can be found in standard finite
element textbooks, e.g., Finder and Gray 1977; Huyakorn and Finder 1983.) In
each case, the governing equation (4-1) or (4-2) is written in a general
tensorial form as follows:
It
where x* (i-1,2) are coordinates as defined before, and K^ is component ij
of the hydraulic conductivity tensor.
A trial function h^. is used to approximate h. The expression for h-j. is given
by
hT(x,z,t) - Nj(x,z) hj(t), J - 1,2, ...,n (4-10)
140
-------
where Nj(x,z) and hj(t) are interpolating basis functions and nodal values of
h-j. at time t, respectively, n is the number of nodes in the finite element
network, and repeated indicies imply nodal summation.
An integral approximation of (4-9) may be obtained using the Galerkin
weighted residual criterion, which requires that the following condition be
satisfied:
/ Nj. L(h) dR - 0, 1 = 1, 2 ..... n (4-11)
R
where L(h) corresponds to the left-hand-side of (4-9) with h replaced by
equation (4-10), R is the solution domain, and dR - dx-^ dX2 for the plane
case and dR - 2jrx^dxidx2 for the axisymmetric case.
Combining equations (4-9) - (4-11) and using Green's theorem to transform the
second derivative terms, one obtains
5Nj dhj
Sg NI Nj
R I dx± axj dt
dR - / Kj, n, dB
B J ax.
- / N! q dR - 0, I = 1, 2, ... n (4-12)
R
where B is the boundary of the solution domain, and n^ is the outward unit
vector normal to the boundary.
Equation (4-12) can be written more concisely as
AIJ hJ + BIJ
dhj
where
AIJ - I / Kij -- ^ (4-13a)
Re
BIJ = I S ss NI NJ ^ (4-13b)
Re
FI " I < / Vn Nj dB + / Nj q dR) (4-13c)
S Be Re
141
-------
in which Re is the element subdomain with boundary Be, and V denotes the
normal velocity at the boundary. The sign convention for V is the same as
for q. That is, Vn is positive for inward flow and negative for outward
flow.
Equation (4-13) represents a system of n linear ordinary differential
equations. Time integration is performed as follows:
Au hj+ A: (hJ+1
where w is a time-weighting factor, superscripts k and k+1 are used to denote
the previous and current time level, respectively, and At^ = t^+-i - t^. The
value of w is set to 0.5 for a centered-in-time scheme. A value of 1.0 for w
results in a fully implicit (backward difference) time scheme.
Equation (4-14) is a set of linear algebraic equations that may be rearranged
in the form:
In the SAFTMOD code, the value of u> is normally set equal to 1. The
resulting set of algebraic equations is solved using a direct Gaussian solver
that takes advantage of the symmetric and banded feature of the global
coefficient matrix.
4.3.2.2 Areal Flow --
Consider an areal flow problem involving a general situation where the system
consists of two aquifers separated by an intervening aquitard. SAFTMOD
represents the domain of each aquifer by a two-dimensional areal grid with
rectangular elements. Both aquifers are discretized in an identical manner.
Vertical one-dimensional column grids are used to represent the aquitard flow
domain. Each aquitard column has the two end nodes that coincide with
aquifer nodes.
For each aquifer, the governing areal flow equation is written in a tensorial
form as follows:
ah
where T^ is the i-j component of the transmissivity tensor defined as T^.= -
K^H with Kji and H being the component of the hydraulic conductivity tensor
ana the saturated thickness, respectively.
142
-------
Applying the Galerkin procedure to (4-16) and performing a temporal
approximation using finite differences yield a set of algebraic equations
that can also be represented by equation (4-15) with Ajj, BJJ, and Fj being
defined as follows:
Au - S Ti
Re 8*1 ax,
BIJ - S 1 N Nj dR (4-17a)
Re
FI - I ( / Vn H Nj dB + J Nj q* dR + / Nj 7 dR) (4-17b)
e Be Re Re
Because the right-hand side of (4-17b) contains the unknown aquitard leakage
flux, 7, it is necessary to consider the one -dimensional governing equation
for vertical leakage in the aquitard column intersecting aquifer node I. The
equation corresponds to equation (4-5). To approximate (4-5), the aquitard
column is discretized into a network of n'-l one -dimensional linear elements.
Next the Galerkin procedure and the finite difference temporal approximation
technique are applied. This leads to a tridiagonal set of algebraic
equations for a current time level k+1. The tridiagonal system can be
represented by:
dl
dlf i - 2, 3, .... n'-l
where a^, /3^, and £^ are tridiagonal matrix elements, and d^(i = 1, 2, ...,
n' ) are elements of the right -hand- side vector.
For the general case of areal flow in a two-aquifer system, the solution
scheme implemented in the SAFTMOD code is a sequential solution scheme that
consists of three stages performed for each aquifer. Stage 1 involves using
the Thomas algorithm to perform a forward elimination of the tridiagonal
matrix equation (4-18), thereby allowing the vertical leakage flux -yk+l to be
expressed in terms of head value at node I of the aquifer domain. Stage 2
involves the direct solution of the aquifer matrix equation (with the
incorporation of the aquitard leakage flux) to determine the head values in
the aquifer. At this stage, Picard iterations are performed if the aquifer
concerned is an unconfined aquifer. Stage 3 of the overall scheme involves
performing back- substitution in the Thomas algorithm to determine nodal
values of head in each aquitard column.
143
-------
4.4 DESCRIPTION OF THE TRANSPORT MODULE
4.4.1 Transport Equation
The form of the advective-dispersive transport equation used in SAFTMOD also
depends on the type of analysis. In a vertical cross -sectional analysis, the
following equation is used:
n
3x
-------
o
where c' is the solute concentration in the aquitard (ML ) , D' , $' , and R'
are the vertical diffusion coefficient (L T ) , effective porosity, and the
retardation coefficient for the aquitard, respectively.
For the sake of completeness, components of the hydrodynamic dispersion
tensor also need to be defined. For the aquifer, the components of D-n ,
where i and j denote coordinate subscripts, are defined as
Dll - <*L V12/| V | + aT V22/| V | + DQ (4-23)
D12 - D21 ~ Vl V2/l V I (4-23a)
where QT and a™ are the longitudinal and transverse dispersivities,
respectively, and D is the effective molecular diffusion coefficient.
The initial and boundary conditions of the transport problem are as follows :
c(xi,0) - c0 (4-24)
c(xi,t) = c on Bj_' (4-24a)
Dij ffr ni = q? °n V (4-24b)
ni - qc on B3' (
where B-^' is the portion of the boundary where concentration is prescribed
as c, 62' is the portion of the boundary where dispersive mass flux of solute
is prescribed as q^, and B-j' is the portion of the boundary where the total
solute mass flux is prescribed as q£. Note that the sign convention
for qj; and q£ are positive for inward mass flux and negative for outward mass
flux, x^ (i-1,2) and x^ (j-1,2) are coordinates defined earlier.
4.4.2 Numerical Solution of the Transport Problems
Numerical approximations of the solute transport equations are performed in a
similar manner to that of the groundwater flow equations. For a cross-
sectional or an axisymmetric flow simulation, the following set of linear
algebraic equations is derived:
BIJ
n
145
-------
where w is the time weighting factor (set equal to 4 if the Crank-Nicholson
time stepping scheme is required, and set equal to 1 if a fully-implicit
scheme is desired), and
3N
dR (4-26)
(4-26a)
dR) (4-26b)
AIJ - j
BIJ = j
BJJ = ]
r r D
I J Pij
3 Re [ J aXi
r / * R % N
5 Re
^ ^ J n I
0 O
1 T7 M
axj dx±
j
-------
that is generated by transformation of parent species. The source term
required is given by
m - -X 6 ft R0 *o cf (4-28)
where subscript £ is used to identify the parent species, ru is the number of
parent species, and £j> is the mass fraction of parent component £ that
transforms into the daughter species under consideration.
The numerical solution of the transport equation for a daughter species of a
decay chain can be performed in the same manner as that described previously
for a single species. The source term from equation (4-28) is incorporated
into the finite element matrix equation by modification of its right-hand
side. The modified right-hand term of the finite element equation is given
by
F* - FT + T (f NT m dR) (4-29)
£ Re
*}? •
where Fj and Fj are the original and modified right-hand terms, and m is
given by equation (4-28).
For areal transport in a system with aquitards, it is also necessary to treat
the aquitard diffusion equation (4-22) by adding the following source term to
the right-hand side:
m' = -S^p RJ, \'f c'p (4-30)
In the finite element solution, m' is treated in a similar manner to m.
4.4.4 Spatial and Temporal Discretizations
Spatial discretization of the solution region is performed automatically in
the SAFTMOD code. A rectangular grid (Figure 4,2) is used to represent the
modeled two-dimensional region. In a cross-sectional or an axisymmetric
analysis, the modeled region corresponds to the entire system, which is
considered as one unit where different materials are accommodated by
assigning different sets of material properties for different elements. In
an areal analysis, a distinction is made between aquifer and aquitard
regions. Each aquifer region is discretized using a network of rectangular
elements. The same nodal arrangement and element configuration are used for
different aquifers. Nodes and elements are numbered sequentially as shown in
147
-------
YO
8
12
16
,0
,
0
XO
20
19
18
18
LEGEND
ELEMENT
NUMBER
NODE
NUMBER
Y
(a)
o.
0
0
7
o
0
0
0
0
0
0
0
0
0
(b)
Figure 4.2.
Discretization of aquifer regions
region, and (b) irregular region.t
(a) regular
148
-------
Figure 4.2a. Irregular boundaries are accommodated by assigning zero
material number to elements outside the region of interest (Figure 4.2b). In
the areal simulation of a system comprising two aquifers separated by an
aquitard, numbering of nodes and elements is performed aquifer by aquifer
starting from the bottom aquifer. The confining aquitard region is treated
as a series of vertical columns each of which is bounded at both ends by two
opposite aquifer nodes. Each aquitard column is represented by a one-
dimensional linear finite element grid. The one-dimensional grid is
generated automatically using a local dimensionless vertical coordinate, £,
defined as £»z'/b', where z' is the height measured from the bottom of the
aquitard, and b' is the thickness of the aquitard column. Aquitard nodes in
the column are numbered sequentially from the bottom to top of the column.
In SAFTMOD, temporal discretization can also be performed automatically if
required. Time steps and time values at the end of the time steps are
generated using a simple algorithm: At^ - fAt^.^ < Atmax, and t^+]_- t^+At^,
where f is the time multiplier and Atmax is the maximum value of time step
allowed. As an alternative to automatic time step generation, the time
discretization data required for transient analysis can be supplied to the
code by the user.
4.5 RESULTS OF SAFTMOD TESTING SIMULATIONS
Three sets of benchmark problems were used to test the SAFTMOD code. The
first set consists of three steady and transient problems designed to test
the flow component of the code. The second set consists of four transient
two-dimensional transport problems. The third set consists of a coupled
flow-transport problem designed to illustrate a field application of SAFTMOD.
Numerical results obtained from SAFTMOD are compared with analytical
solutions and results obtained using another finite element code, VAM3D
(Huyakorn et al., 1986b). The SAFTMOD testing simulations were conducted
before being linked to the other RUSTIC modules.
4.5.1 Flow Module
Three flow problems are presented which test various formulations of the flow
module. The first problem is a well flow problem designed to check the areal
and axisymmetric analyses as applied to two-aquifer systems. The second
problem is a drain flow problem used to compare the areal and cross-sectional
analyses of steady and transient flow in unconfined aquifers. The third
problem is a confined flow problem used to test the incorporation of flux and
head boundary conditions in the flow module.
1) Well Flow in a Two-Aquifer System with a Confining Aquitard
This problem concerns transient flow to a pumping well fully penetrating a
system consisting of two aquifers separated by a confining aquitard (Figure
4.3). For such a system, the analytical solution of the problem can be found
in Neuman and Witherspoon (1969). The case used to test SAFTMOD had been
solved numerically by Huyakorn et al. (1986). Two alternative simulations
were performed using the SAFTMOD code. In the first simulation, the
axisymmetric grid shown in Figure 4.4a was used. This grid consists of 120
elements and 143 nodes. The elements that lie in the two identical aquifers
149
-------
are assigned material number 1, and the elements that lie in the aquitard are
assigned material number 2. In the second simulation, the areal grid shown
in Figure 4.4b was used to represent the first (bottom) aquifer, and a second
areal grid with identical nodal spacings to the first one was used to
represent the second (top) aquifer. Node numbers of this second grid start
from 170 and finish at 338, and element numbers start from 145 and finish at
288. Both the areal and the axisymmetric simulations were performed for 40
time steps using an identical set of time step values. The time step values
were generated as follows: At^ = 36 s, Atlc=1.414Atk.1 < Atmax = 3.6*106 s.
Numerical results obtained from SAFTMOD are compared with the analytical
solution. Figure 4.5 shows the plot of time versus drawdown in the pumped
aquifer. Agreement between the numerical and the analytical solutions is
fairly good despite the fact that the numerical results were obtained using a
rather coarse grid and large time steps. Note that the same degree of
accuracy was obtained with the axisymmetric and areal numerical simulations.
2) Flow to Parallel Drains in an Unconfined Aquifer
This problem concerns flow to parallel drains in an unconfined aquifer
receiving vertical recharge due to precipitation (Figure 4.6). A steady-
state analytical solution can be found in Bear (1979). Both areal and cross-
sectional simulations were performed using SAFTMOD. The simulations were
transient and based on the areal and cross-sectional grids shown in Figure
4.7. The initial hydraulic head in the aquifer was assumed to be 50 ft. The
hydraulic conductivity and effective porosity of the aquifer were assumed as
1 ft/day and 0.2, respectively. The aquifer was subjected to vertical
Q= 1.58 x 10~3 m3/B
////s/ss ////////
AQUIFER K= 2 x 1O^5 m/s
S § — 1 0 m
K = 10"8 m/s
AQUITARD _ _4 _-i
Ss=8x 10 * m '
K = 2 x 10~5 m/S
' AQUIFER .. -
i r Ss = 1 0~^ m~'
i _ — _^~ — . — . — . — . . — . — ^ — . — . — . — . — , — , — • — • — i — i — ,
1
5 m
<
1
10
5
m
m
i
1,000 m
Figure 4.3. Problem description of flow to a well in a two-
aquifer system with a confining aquitard.
150
-------
ftl11
1
A >
13
/
~* — I —
(,
a)
14.
13
3
Material 1
Material 2
Material 1
3
X = T
169
fc.
157
(b)
Figure 4.4.
Finite element grids used in simulating well flow
in two aquifer system: (a) axisymmetric grid, and
(b) areal grid.
151
-------
10
O
D
cr
Q
0.1 -:
rj.01
— ANALYTICAL SOLN.
* F.E. (Axsym. grid)
0 F.E. (Areal grid)
1 — i i 1 1 1 n| - 1 — i i i i ni
i i 1 1 1 I
0.01
0.1
1.0 10
TIME, hours
100
1000
Figure 4.5. Simulated drawdown-time relationships for the
pumped aquifer.
recharge at a rate of 0.01 ft/day. Note that it was sufficient to consider
only half of the saturated flow region because of symmetry about the vertical
axes midway between the drains. A constant head of 50 ft was prescribed at
the drain where x = 0. Zero normal velocity was specified on the aquifer
base and at x = 500 ft (midway between the two drains). Both simulations
were performed for 35 time steps using the following time step values: At-^ =
1 day, At^ = l^At^-^ < 100 days, for k - 2 35. Results obtained from
the simulation are plotted in Figure 4.8 for five typical time values. As
expected, the plot shows a continuous water table mounding, and a steady-
state condition at t > 2000 days. There is excellent agreement between the
steady-state analytical and numerical solutions.
3) Flow to a Gallery
This problem concerns flow to a gallery (horizontal well) penetrating a
shallow confined aquifer system with two partially penetrating parallel
drains (Figure 4.9). The aquifer is anisotropic with a 10 to 1 contrast in
the horizontal versus vertical hydraulic conductivity. The gallery is
152
-------
LAND SURFACE
M = 0.01 tl/d
60 ft
"Vy/rr
K = 1 tt/d
p= 0.2
/////// /\7 /////////////////
\. 1OOO tt ^|
Figure 4.6. Problem description for flow to parallel drains in an
unconfined aquifer.
pumping at a constant rate of 0.5 m /hr. The water levels in both drains are
assumed to remain constant with time. No general analytical solution is
available for the case considered. For this reason, numerical results from
SAFTMOD were checked against corresponding results from another published
numerical code, SATURN (Huyakorn et al., 1984). The transient flow problem
was simulated with both codes using the cross-sectional grid shown in Figure
4.10. The simulations were performed for 20 time steps. Time step values
were computed using the following algorithm: At-^ = 0.1 hr, At^ = 1.414At^_-^
< 10 hr, k = 2 20. Numerical results obtained from both codes are
compared in Figure 4.11, where hydraulic head profiles along the horizontal
line passing through the gallery are plotted. As can be seen, the two
numerical solutions are in excellent agreement.
4.5.2 Transport Module
Four transient solute transport problems were selected. The first problem is
a one-dimensional problem used to check the cross-sectional formulation of
the Galerkin finite element technique. The second and third problems are
two-dimensional dispersion problems used to check the areal formulation and
to assess the accuracy of the numerical approximations for cases involving
conservative and nonconservative contaminants and source boundary conditions
of continuous and instantaneous releases. The fourth problem is an injection
well transport problem used to test the axisymmetric formulation applied to a
situation involving aquifer-aquitard systems.
1) One-Dimensional Dispersion in Uniform Flow
This problem involves one-dimensional advection-dispersion of a conservative
solute species through a semi-infinite porous medium (Figure 4.12a). As
illustrated, the contaminant is released from a channel fully penetrating a
shallow confined aquifer. The analytical solution of the problem was
developed by Ogata and Banks (1961). This solution can also be found in Bear
(1979). For the case simulated by SAFTMOD, the source concentration at the
inlet was prescribed as 1 mg/m . The values of Darcy velocity, effective
153
-------
100
1 ^ *•
t
... O
. ^
50 ft
1
1
* 41 *"
mnnft ^ ^ '
1 U U U I I ' • ' ^
(a)
, Y = Z
42
41 n r\ n ft fc-
• IUUU Tl ^
(b)
X
Figure 4.7. Finite element grids used in the simulation of flow
through parallel drains: (a) areal grid, and (b) cross-
sectional grid.
154
-------
60
70
e>
t_
c
r 60
2
LU
_J
UJ
• ANALYTIC
FINITE ELEMENT
t= 0
100
200 300
DISTANCE, meters
400
600
Figure 4.8. Simulated steady-state and transient water table profiles.
(Adapted from Battelle and GeoTrans, 1988).
porosity, and longitudinal dispersivity were specified as 1 ra/d, 0.25, and
5m, respectively. The modeled length of the aquifer domain was specified as
400 m. A uniform rectangular grid consisting of 82 nodes and 40 elements,
each with Ax - 10 m, was used. The simulation was performed for 25 time
steps with At kept constant and equal to 2.5 d. Simulated concentration
distributions are plotted in Figure 4.12b for two typical timevalues. As can
be seen, the numerical and analytical solutions are in excellent agreement.
2) Areal Dispersion in Uniform Flow
This problem concerns two-dimensional dispersion of solute in a uniform and
steady groundwater flow field. In practice, the situation may correspond to
that involving the areal migration of contaminants continually released into
an extensive aquifer from an injection well treated as a point source in the
areal plane (Figure 4.13a). Assuming that the injection rate is small so
that the natural groundwater flow is undisturbed, the analytical solution of
the problem can be found in Wilson and Miller (1978). Two cases involving
conservative and non-conservative solute species were simulated using the
SAFTMOD code with the parameter values shown in Table 4-1. The selected
parameter values for Case 1 were based on data on the field study of the
hexavalent chromium contamination problem reported by Perlmutter and Lieber
155
-------
15 m
5m
10 m
K2 = 0.1 m/hr
Ss = 0.0005 m~'
10 m
Figure 4.9. Problem description of flow to a gallery..
b
4
3
2
1
_^ or\n ^, .___..... k^
IUO 1
104
103
i
*
15 m
102
101 -
r
Figure 4.10. Finite element grid used in simulation of flow to a
gallery.
156
-------
30
26 -•
D
-------
Channel
Piezometric Head
0
(a)
UJ
a. - 5m
ANALYTIC
• SAFTMOD
100
ISO 200 750
DISTANCE, x (m)
350
400
(b)
Figure 4.12.
One dimensional dispersion in uniform flow: (a)
problem description, and (b) simulated concentration
profiles. (Adapted from Huyakorn et al. , 1984b) .
158
-------
(a)
ac/ay
300 m
^source
v ./X..- 360 m
=0
ac/ay = o
360 m
1980 m
(b)
Figure 4.13.
Areal dispersion in uniform flow: (a) problem
description, and (b) modeled region and boundary
conditions. (Figure 4.13(b) is a half-system
representation and is not drawn to scale.)
159
-------
TABLE 4-1. VALUES OF PHYSICAL PARAMETERS FOR PROBLEM 2
Parameter Value
Darcy velocity, V 0.161 m d
Porosity, <£ 0.35
Longitudinal dispersivity, a^ 21.30 m
Transverse dispersivity, QJ 4.30 m
Aquifer saturated thickness, b 33.50 m
Contaminant mass flux, QcQ 704.0 g m d
(per unit thickness of aquifer)
Case 1:
Retardation coefficient, R 1.0
Decay constant, A 0.0 d
Case 2:
Retardation coefficient, R 2.0
Decay constant, A 0.00019 d"1
Selected grid
Fine: Ax = 60 m, Ay = 15 m
Medium: Ax = 60 m, Ay = 30 m
Time steps
At = 100 d (Use 28 time steps)
element grid has 49 rows of nodes having variable spacing, and 11 columns of
nodes having constant spacing of Ay = 5 m. The grid consisted of a total of
539 nodes and 480 elements. Values of the physical parameters used for the
simulation and coordinates of the grid lines are provided in Table 4-2. One
node located at x = y = 0 was used to represent the contaminant source.
o
Initially the concentration was specified as 400 mg/nr at this node and zero
elsewhere. Considering that the total nodal area of the source node ^
(counting the upper and the lower halves of the xy plane) is equal to 25 m
and the porosity is 0.35, this represents placing a finite contaminant slug
at x = y = 0. The mass of this slug per unit aquifer thickness is equal to
3500 mg/m.
In order to accurately simulate movement of the dispersing slug, temporal
discretization was performed using an automatic scheme (also presented in
Table 4-2). The initial time step value was assigned as 1 day. Values of
the remaining 29 time steps were generated using a multiplier of 1.2 and a
maximum allowable time step value of 2 days. The problem was solved for 30
time steps and was subject to zero normal concentration gradient boundary
conditions.
160
-------
1.2 r
c"
o
c
Q)
O
C
o
O 0.6
V)
c
O 0.4
C
Q)
E
QO.Z
Case 1: X = 0, R = 1
Analytic Solution
SAFTMOD (medium-mesh)
SAFTMOD (fine-mesh)
- xs) / 2 OL
VC/QCO
6 12 16 20
Dimensionless Distance, x,
2«
28
Figure-4.14.
Simulated concentration profiles along the x-axis for case 1.
(Adapted from Huyakorn et al., 1984b).
The numerical solution is compared in Figure 4.16b with the analytical
solution computed by using equations (4.13) and (4.14). It can be seen that
good predictions of the concentration distributions were achieved. At the
early time value of 3.96 days, the numerical result deviates slightly from
the analytical solution, particularly near the peak of the concentration
curve.
4) Transport from an Injection Well in an Aquifer Confined by Aquitards
This problem concerns contaminant transport from a well in an aquifer
confined by two very low permeability aquitards (Figure 4.17). The well
injects a fluid containing dissolved contaminant of prescribed concentration,
CQ. To make the problem amenable to analytical solution, it is assumed that
the dispersion coefficient of the aquifer is constant and advection in the
aquitard is insignificant. Based on these assumptions, the analytical
solution for the case of infinite radial extent and infinite aquitard
thickness can be obtained by converting the solution derived by Avdonin
(1964) for the equivalent heat transport problem.
161
-------
Case 2: X = 1.9 x ICT" d~'
R = 2
Analytic solution
--~ SAFTMOD (fine-mesh)
XD = (x - xs) / 2 QL
CD = 2n(aLaT)w VC/Qc0
Q
< B 12 16 20 24 28 32
Dimensionless Distance, x
Figure 4.15. Simulated concentration profiles along the x-axis for case 2.
(Adapted from Huyakorn et al., 1984b).
Two cases were simulated: (1) in which the aquitard is leaky, and (2) in a
perfectly confined system. Values of the physical parameters used for both
these cases are presented in Table 4-3. Both simulations were conducted in
axisymmetric mode using a vertically-oriented, rectangular, finite element
grid. Except for the first node, the other 40 nodal spacings in the radial
direction were kept constant with Ar = 0.5 m. The 15 nodal grid lines in the
vertical (z) direction were variable and equal to 0.0, 1.0, 2.0, 2.1, 2.25,
2.75, 3.0, 3.5, 4.0, 6.0, 10.0, 12.0, 20.0, 30.0, and 52.0 m, respectively.
The grid thus consisted of 615 nodes and 560 elements. The simulations were
run for 50 time steps using a constant At of 0.5 d.
The concentration breakthrough curves for case 1 (leaky system) and for case
2 (nonleaky system) are compared for t = 10 d and t = 25 d in Figure 4.18.
In both cases there is excellent agreement between the numerical and
analytical solutions.
The breakthrough curves for case 1 are steeper and more retarded than these
for case 2 due to the influence of dispersion into the aquitard.
162
-------
Initial
injectionA
V
(a)
— Analytic
SAFTMOD
t - 16.59 d
Figure 4.16. Areal dispersion of a solute slug in uniform groundwater
flow: (a) problem description, and (b) simulated
concentration profiles. (Adapted from Huyakorn et al.,
1984b).
163
-------
TABLE 4-2. VALUES OF PHYSICAL PARAMETERS AND DISCRETIZATION DATA FOR THE
PROBLEM OF AREAL DISPERSION OF A SOLUTE SLUG IN UNIFORM
GROUNDWATER FLOW
Parameter
Value
Darcy velocity, V
Porosity, $
Longitudinal dispersivity, a^
Transverse dispersivity, aj
Retardation coefficient, R
Decay constant, A
Solute mass per unit aquifer thickness, m
2 m d'1
0.35
4 m
1 m
1
0.0 d'1
3500 mg m"1
x,.' (m) ,
j. N ' *
i-l 49
-50.,
5.,
50.,
95.,
180.,
270.,
-40.
10.
55.
100.
190.
280.
-30
15.
60.
110.
200.
290.
•25.,
20.,
65.,
120.,
210.,
300.
-20.
25.
70.
130.
220.
-15.,
30.,
75.,
140.,
230.,
-10
35.
80.
150
240
, -5. ,
. 40.,
, 85.,
, 160.,
, 250.,
o.,
45.,
90.,
170.,
260.,
j-i,.!. ,ii
Atj^ - 0.1 d
Ati, -1.2 Ati, -i < 2
0., 5., 10.-, 15., 20., 25., 30., 35., 40.,
45., 50.
d, k - 2 30
//
AQUIFER v >'.•.-: •"•-• • '•'•••'.;• -V.-'.~.
•;-;r; •_.-/•.•.•.•.-«'•.'••'•'•'*.•'•:'•':'.•••'.•.'.•'•.'!•.'•.••.
b = 4 m
J_
b'= 50 m
r0= 20 m
H
Figure 4.17.
Solute transport from an injection well fully
penetrating an aquifer overlain and underlain
by aquitards.
164
-------
Analytic
0 SAFTMOD
5 8 10 12
RADIAL DISTANCE,
(a)
10 12 14 16 IB
RADIAL DISTANCE, m
Analytic
0 SAFTMOD
(b)
Figure 4.18. Simulated concentration profiles in the aquifer: (a) case 1,
and (b) case 2.
165
-------
TABLE 4-3. VALUES OF PHYSICAL PARAMETERS FOR THE PROBLEM OF SOLUTE TRANSPORT
FROM AN INJECTION WELL
Parameter
Fluid injection rate, Q
Solute concentration at the well, CQ
Well radius, rw
External radius , rQ
Porosity of aquifer, ^
Radial dispersion coefficient for aquifer, DQ
Aquifer thickness, b
Assigned
Value
10.0
1.0
0.1
20.0
0.25
0.152
4.0
Unit
m3 d'
ppm
m
m
--
m2 d~
m
1
1
Vertical dispersion coefficient of
aquitard, DQ 0.152 m2 d"1
Aquitard porosity, «£' 0.025
Case 2: Nonleaky System
No aquitard
166
-------
SECTION 5
LINKAGE OF VADOSE AND SATURATED ZONE MODELS
As stated earlier in this volume, the objective of RUSTIC is to simulate the
fate and transport of chemicals applied in agricultural settings, from the
source area, through the crop root zone, the vadose and saturated zones, and
into drinking water wells. Instead of writing an entirely new computer code
to meet this objective, it was decided early on in the development effort
that existing software would be used where it was available and met, or could
be modified to meet, the Agency's standards for computer codes. PRZM was
selected to provide the simulation of the crop root zone, due to its degree
of recognition and acceptability both to the Agency and the agricultural
chemical industry.
Several options were studied for linkage of models to PRZM to provide the
balance of the required simulation capabilities. The first involved use of
PRZM and a saturated zone model only. In this configuration, PRZM would be
used to represent both the root zone and the vadose zone. This option was
eliminated because the assumptions of the elementary soil hydraulics in PRZM
(i.e., drainage of the entire soil column to field capacity in one day) were
not considered to be adequate for representing a thick vadose zone. The
second involved the use of PRZM and a 2-D or 3-D variably saturated flow
model. In this configuration, the latter model would represent both the
vadose and saturated zones. Although the interface between the unsaturated
and saturated zones would be more rigorously simulated, it was felt that the
resulting code would be too computationally intensive for the intended
application scenarios, especially when Monte Carlo simulation was involved.
Therefore, this option was also dropped.
The combination of models finally selected is depicted in Figure 5.1. In
this configuration, PRZM is linked to a two-dimensional (single or two-layer)
saturated zone model through a one-dimensional vadose zone model. Both the
vadose and saturated zone models simulate water flow and solute transport.
Two new codes were written to provide the simulation of the vadose zone
(VADOFT) and the two-dimensional X-Y, X-Z or axisymmetric simulation of the
saturated zone (SAFTMOD). A significant problem associated with this type of
linkage is difference between time domains of the individual models. While
the vadose zone models are required to operate on a daily or less-than-daily
time step, the time step of the saturated zone model could be from days to
months. Another significant problem is the interfacing of the vadose and
saturated zone models, where a fluctuating water table may require special
handling of model spatial domains.
167
-------
3«U
/fv> <^
m
i
VADOSE
ZONE
MODEL
m
PR2M
(1-D Flow and Transport)
VADOFT
(1-D Flow and Transport)
SATURATED ZONE MODEL
SAFTMOD
(2-D Flow and Transport)
Figure 5.1. Linked modeling system configuration
168
-------
5.1 TEMPORAL MODEL LINKAGE
The resolution of the temporal aspects of the three models was straight-
forward. PRZM runs on a daily time step. The time stepping in VADOFT is
dependent upon the properties of soils and the magnitude of the water flux
introduced at the top of the column. In order for the non-linear Richard's
equation to converge, VADOFT may sometimes require time steps on the order of
minutes. The SAFTMOD time step, on the other hand, would normally be much
longer than one day.
For simplicity, it was decided that the time step of SAFTMOD would always be
an integer multiple of PRZM's daily time step. This makes the direct linkage
of PRZM to SAFTMOD in a temporal sense very straightforward. PRZM is simply
run for the number of days of the SAFTMOD time step, SAFTMOD is then run, and
so forth. The water and pesticide fluxes from PRZM are summed and averaged
over the length of the SAFTMOD time step. SAFTMOD receives the time-averaged
fluxes as input.
For the linkage of PRZM, through VADOFT, to SAFTMOD, the resolution is not
much more complicated. VADOFT is prescribed to simulate to a "marker" time
value, specifically to the end of a day. The last computational time step
taken by VADOFT is adjusted so that it coincides with the end of the day.
This is depicted in Figure 5.2. PRZM's daily water fluxes are used as input
to VADOFT. VADOFT utilizes this flux as a constant over the day and adjusts
its internal computational time step in order to converge. PRZM and VADOFT
execute for the number of days of the SAFTMOD time step. The output from
VADOFT is then time averaged and used as input to SAFTMOD.
5.2 SPATIAL LINKAGES
The spatial linkages utilized for the models are more complex. The principal
problem is that of the presence of a rising and falling water table, which
complicates the interfacing of the vadose (or root zone) and saturated zone
codes. A second problem is that of the incompatibility between the
hydraulics in PRZM and VADOFT. Of course, any linking scheme utilized must
provide a realistic simulation of the flow of water and transport of solutes
at the interfaces and must provide for continuity (mass balance).
5.2.1 PRZM and VADOFT
The major problem with the interfacing of these two models is that while
VADOFT solves the Richard's equation for water flow in a variably saturated
medium, PRZM uses simple "drainage rules" to move water through the soil
profile. Because of this incompatibility, there may be times when PRZM
produces too much water for VADOFT to accommodate within one day. This is
very likely to happen in agricultural soils, where subsoils are typically of
lower permeability than those of the root zone, which have been tilled and
perforated by plant roots and soil biota. The result of this would be water
ponded at the interface which would belong neither to PRZM nor VADOFT.
Several solutions were proposed to handle this situation, including:
169
-------
Marker Time Value-
PRZM TIMESTEP, Atd
At.
Last Aty Corrected to Coincide
with Marker Time Value
Atv
h
VADOFT TIMESTEP. At.
SAFTMOD TIMESTEP, AT
Figure 5.2. Resolution of component model timesteps,
170
-------
• Automatically switching from a flux (unponded) condition to a head
(ponded) condition in VADOFT at the top node if this latter
condition occurs
• Changing PRZM hydraulics to a Green-Ampt formulation, which has a
parameter for the saturated hydraulic-conductivity
• Applying the VADOFT flow and transport solutions to the root zone as
well
• Checking the velocities coming out of PRZM to make sure that they
are within an acceptable range and if not, using PRZM's restricted
drainage option to correct them
For a variety of reasons, all of these proposals proved to be unattractive.
In the end it was decided that the solution was to prescribe the flux from
PRZM into VADOFT so that VADOFT accommodates all the water output by PRZM
each day. This eliminates the problem of ponding at the interface. However,
it does force more water into the vadose zone than might actually occur in a
real system, given the same set of soil properties and meteorological
conditions. The consequence is that water and solute are forced to move at
higher velocities in the upper portions of the vadose zone. If the vadose
zone is deep, then this condition probably has little impact on the solution.
If it is shallow, however, it could overestimate loadings to groundwater,
especially if chemical degradation rates are lower in the vadose zone than in
the root zone.
5.2.2 VADOFT and SAFTMOD
The principal problem here is that of the presence of a dynamic boundary
(rising and falling water table) between these models. There are two
approaches available to handle the problem. The first is to expand or
contract the spatial domains of the models to accommodate the moving
boundary. This is not a particularly insurmountable problem for the vadose
zone model; however, for the saturated zone model it would require the
addition of nodes at the upper boundary. This would require a constant
evaluation and switching of the set of nodes receiving fluxes from the vadose
zone. This appeared undesirable.
The second option is to overlap the spatial domains of the models and
interpolate values for fluxes based upon the position of the water table.
This latter approach was ultimately utilized. It has the additional feature
of eliminating the effects of the bottom boundary conditions prescribed for
the vadose zone model on the simulation of solute transport just above the
water table.
5.2.2.1 VADOFT/SAFTMOD (2-D, X-Y) Linkage
The following paragraphs describe the interface of VADOFT and SAFTMOD where a
2-D (X-Y) configuration is used for SAFTMOD. The linkage is depicted in
Figure 5.3. The sequence in which the models are run is VADOFT flow, then
SAFTMOD flow, VADOFT transport, then SAFTMOD transport.
171
-------
SAFTMOD TIMESTEPS
VADOFT TIMESTEPS
(VADOFT)
Flow
Boundary
Area Ab —
Control Volume
Boundary Nodes
.4-
Flux from
Boundary Water F|ux
Node,mb to SAFTMOD, w
VADOFT
Nodes
- VADOFT Nodes
below Water Table
Common Aquifer
Base Node
Figure 5.3.
Definition sketch for a VADOFT/SAFTMOD
control volume (x-y) configuration.
172
-------
VADOFT and SAFTMOD are linked through the use of a common node at the base of
the aquifer, thus ensuring continuity of pressure head and mass flux. The
VADOFT flow model is run using At length time steps (7 in Figure 5.3), and
the sum of the water flow predicted by VADOFT from the base of the aquifer is
the water flux input to SAFTMOD. Next, SAFTMOD flow is run and a new
pressure head at the aquifer base ($0) and a new water table elevation (Z^)
are predicted (by SAFTMOD) at time T2 (see Figure 5.3). The new water table
elevation (Z 2) is used, along with (Z 1), to interpolate water table
elevations(Zj) for each VADOFT time step (At). Once these water table
elevations and (V^) are known, VADOFT transport is run to compute c£, the
contaminant concentrations at the water table. These results (V™ and G£ )
can be used to compute the contaminant mass flux (M) to the water table over
the interval AT. Once M is known, SAFTMOD transport can be run. Because
SAFTMOD predicts a new mass in the control volume beneath the water table
different from VADOFT, the concentrations at each VADOFT node 'i' beneath the
water table, C.2, must be corrected. This is done in the next AT time step
by multiplying concentrations at the VA.DOFT nodes beneath the water table by
a correction factor which exactly conserves mass.
The steps which the model takes to achieve this linkage are outlined below:
(1) VADOFT flow is run for AT (T^ to T2) taking small time steps At with
i>s (VADOFT) - \J>Q (SAFTMOD) as the initial condition.
P r*
(2) The water fluxes from the base node of the aquifer are summed and
divided by AT. This scalar quantity, W, is the water flux to
appropriate SAFTMOD nodes for AT. Each SAFTMOD node receives a total
flux in each AT defined as the water flux times the area assigned to
each node.
T
(3) SAFTMOD flow is run to obtain Z^2.
T T
(4) Zj is computed for each At by interpolating between Z-^- and Z ^.
(5) VADOFT transport is then run.
(6) M, the mass flux of contaminant to the water table, is determined by
12 ^ r Ac
M = S (C| V£ At + 8 D — At)/(T2 - T-.) (5-1)
TI AZ
where
G.JJ = the concentration at the water table at time t (ML )
V.p = the vertical Darcy velocity at the water table at time t, as
computed by VADOFT flow in Step (1) above (LT"1)
173
-------
6 - the volumetric water content of the vadose zone just above the
water table (L3L~3)
D - the dispersion/diffusion coefficient in the vadose
zone (L2!'1)
Ac - the concentration difference in the vadose zone between the
node at the water table and the node above (ML )
AZ = is the distance between the two nodes (L)
(7) SAFTMOD transport is run using M as the mass flux boundary condition at
the appropriate SAFTMOD nodes. The actual flux to each node during AT
is determined by multiplying M by the appropriate nodal area.
(8) The contaminant mass in the saturated control volume, ¥, at the end of
AT is determined from:
To N To
V = 2 (C 2 R)i Vt (5-2)
where
T
C 2 - the SAFTMOD concentration values at node i at time value (ML"3)
¥j - the portion of the saturated control volume assigned to
node i (L3)
<£^ - porosity of node i
Ri — retardation coefficient of node i
N - the number of SAFTMOD nodes in the control volume
(9) A mass correction is made to VADOFT nodes below the water table at the
end of AT. The basis for the correction is
T2 s
To To
-------
? V
(M. ) - the total mass in the control volume below the water table
given by VADOFT at T2 (M)
5.2.2.2 PRZM/SAFTMOD (2-D, X-Y) Linkage
The procedure for linking PRZM to SAFTMOD is very similar to that for linking
VADOFT to SAFTMOD. The difference is that PRZM does not compute head values
and, therefore, head continuity between PRZM and SAFTMOD can not be used as a
common point. PRZM flow and transport are solved using a single-pass
solution (i.e., unlike in VADOFT, flow and transport equations are solved in
a single application of the model). The water flux to the water table is
computed over time step AT. PRZM is initialized so that the water contents
below the root zone are at field capacity (#fc)• Therefore, the water flux
below the root zone can be taken as the water flux to the water table.
SAFTMOD is run using this water flux to obtain Z 2 and values of the
water table are interpolated for each At between T-^ and T2. Values of mass
flux to the water table are then interpolated from the midpoints of PRZM
compartments on either side of the water table. From this point on, the
procedure is exactly the same as for the VADOFT/SAFTMOD linkage.
5.2.2.3 VADOFT/SAFTMOD (2-D, X-Z) Linkage
This linkage is virtually the same as for the 2-D (X-Y) case. The difference
is, of course, that SAFTMOD is capable of simulating vertical variations in
pressures, velocities, and contaminant concentrations. The linkage is most
easily accomplished if VADOFT and SAFTMOD nodes below the water table
coincide.
First, VADOFT flow is run. The initial VADOFT pressure heads at nodes below
the water table are equal to the average pressure heads across the SAFTMOD
nodal rows within the two-dimensional control slice (see Figure 5.4). Then
SAFTMOD flow is run to obtain Z_,2. Following interpolation to
obtain Zt, VADOFT transport to the water table is computed, just as in the 2-
D (X-Y) linkage. SAFTMOD transport is then run and the mass flux from the
boundary of the control slice is computed. The VADOFT contaminant profile is
adjusted exactly as in the 2-D (X-Y) case, over the next AT. The new average
pressure heads for each row within the control slice become the new pressure
heads for VADOFT nodes below the water table for the next VADOFT flow
simulation.
5.2.2.4 PRZM/SAFTMOD (2-D, X-Z) Linkage
This link is accomplished in exactly the same way as that described for
VADOFT and SAFTMOD except that, as in the 2-D, X-Y PRZM to SAFTMOD linkage,
head continuity at the aquifer base is not maintained.
175
-------
VADOFT
Nodes
Control Slice
Boundary
SAFTMOD
Node Rows
VADOFT
Nodes
Figure 5.4.
Schematic of VADOFT /SAFTMOD
x-z configuration
176
-------
SECTION 6
UNCERTAINTY PREPROCESSOR
6.1 INTRODUCTION
In recent years the use of quantitative models to assess the fate and
transport of contaminants in the environment has increased significantly.
Typically these models include a set of algorithms that simulate the fate and
transport of a contaminant within a media (e.g., unsaturated zone, saturated
porous media, air or a surface water body) based on a number of user-
specified parameters. These parameters describe the properties of the
chemical, the transport medium, and the effects that man has on the system.
Unfortunately, the values of these parameters are not known exactly due to
measurement errors and/or inherent spatial and temporal variability.
Therefore, it is often more appropriate to express their value in terms of a
probability distribution rather than a single deterministic value and to use
an uncertainty propagation model to assess the effect of this variability on
the fate and transport of the contaminant.
This section describes the Monte Carlo method of uncertainty propagation and
a Monte Carlo shell which is coupled with the RUSTIC model (subsequently
referred to as the deterministic code in this report). The composite code
(i.e., the uncertainty shell coupled with the deterministic code) can be used
for the quantitative estimate of the uncertainty in the concentrations at the
monitoring point due to uncertainty in the (fate and transport) model input
parameters.
6.2 OVERVIEW OF THE PREPROCESSOR
The objective of the uncertainty analysis/propagation method is to estimate
the uncertainty in model output (e.g., the concentration at a monitoring
point) given the uncertainty in the input parameters and the fate and
transport model. Alternatively stated, the objective is to estimate the
cumulative probability distribution of the concentration at a receptor
location given the probability distribution of the input parameters. If C
represents the concentration at the receptor, then
Cw - g(X) (6-1)
where the function g represents the fate and transport model and X represents
the vector of all model inputs. Note that some or all of the components of X
177
-------
may vary in an uncertain way, i.e. they are random variables defined by
cumulative probability distribution functions. Thus the goal of an
uncertainty propagation method is to calculate the cumulative distribution
function F^ (Cw) given a probabalistic characterization of X. Note that
Fc (Cw) is defined as:
FCw(Cw> - Probability (Cw < GW) (6-2)
where GW is a given output concentration.
6.2.1 Description of the Method
Given a set of deterministic values for each of the input parameters, X-^, X2,
. . . Xn, the composite model computes the output variable (e.g., a
downgradient receptor well concentration C ) as:
Cw = g (X1( X2, X3 . . . Xn) (6-3)
Application of the Monte Carlo simulation procedure requires that at least
one of the input variables, X-^ . . . Xn, be uncertain and the uncertainty
represented by a cumulative probability distribution. The method involves
the repeated generation of pseudo-random number values of the uncertain input
variable(s) (drawn from the known distribution and within the range of any
imposed bounds) and the application of the model using these values to
generate a series of model responses i.e. values of Cw. These responses are
then statistically analyzed to yield the cumulative probability distribution
of the model response. Thus, the various steps involved in the application
of the Monte Carlo simulation technique involve:
i) Selection of representative cumulative probability distribution
functions for describing uncertainty in the relevant input
variables.
ii) Generation of pseudo-random numbers from the distributions
selected in (i). These values represent a possible set of values
for the input variables.
iii) Application of the model to compute the derived inputs and
output(s).
iv) Repeated application of steps (ii) and (iii).
v) Presentation of the series of output (random) values generated in
step (iii) as a cumulative probability distribution function
(CDF).
vi) Analysis and application of the cumulative probability
distribution of the output as a tool for decision making.
178
-------
6.2.2 Uncertainty in the Input Variables
The parameters required by a fate and transport model can be broadly
classified into two different sets that exhibit different uncertainty
characteristics. These are:
• Chemical parameters. Examples of these variables include the
octanol-water partition coefficient, acid, neutral, and base
catalysed hydrolysis rate, soil-adsorption coefficient, Henry's Law
Constant, etc.
• Media parameters. Examples of these variables include the
groundwater velocity, soil porosity, organic carbon content,
dispersivity values, etc.
• Meteorological parameters. Examples include precipitation,
evaporation, solar radiation.
• Management parameters. Examples include irrigation timing,
pesticide application timing, well pumping rates, etc.
Uncertainty in chemical parameters primarily arises due to laboratory
measurement errors or theoretical methods used to estimate the numerical
values. In addition to experimental precision and accuracy, errors may arise
due to extrapolations from controlled (laboratory) measurement conditions to
uncontrolled environmental (field) conditions. Further, for some variables,
semi-empirical methods are used to estimate the values. In this case, errors
in using the empirical relationships also contribute to errors/uncertainty in
the model outputs.
Uncertainty in the second and third sets of parameters, identified above, may
include both measurement and extrapolation errors. However, the dominant
source of uncertainty in these is the inherent natural (spatial and temporal)
variability. This variability can be interpreted as site-specific or within-
site variation in the event that the fate and transport model is used to
analyze exposure due to the use and/or the disposal of a contaminant at a
particular site. Alternatively it can represent a larger scale
(regional/national) uncertainty if the model is used to conduct exposure
analysis for a specific chemical or specific disposal technology on a
generic, nation-wide or regional basis. Note that the distributional
properties of the variables may change significantly depending upon the
nature of the application. Uncertainty in the fourth set of parameters may
arise from a complex variety of factors including climate, sociology,
economics, and human error.
Whatever the source of uncertainty, the uncertainty preprocessor developed
here requires that the uncertainty be quantified by the user. This implies
that for each input parameter deemed to be uncertain, the user select a
distribution and specify the parameters that describe the distribution.
The current version of the preprocessor allows the user to select one of the
following distributions:
179
-------
i) Uniform
ii) Normal
iii) Log-normal
iv) Exponential
v) Johnson SB distribution
vi) Johnson SU distribution
vii) Empirical
viii) Triangular
Depending on the distribution selected, the user is required to input
relevant parameters of the distribution. The first requires minimum and
maximum values. The second and third distributions require the user to
specify the mean and the variance. The fourth distribution requires only one
parameter - the mean of the distribution. For the empirical distribution,
the user is required to input the coordinates of the cumulative probability
distribution function (minimum 2 pairs, maximum 20 pairs) which is
subsequently treated as a piece-wise linear curve. For the triangular
distribution the user is required to input the minimum, maximum and the most
likely value. Finally, the Johnson SB and SU distribution requires four
parameters - mean, variance, and the lower and upper bounds.
In addition to the parameters of the distribution, the user is required to
input the bounds of each model parameter. These bounds may be based on
available data or simply physical considerations e.g. to avoid the generation
of negative values. Values generated outside these bounds are rejected.
Of the above eight distributions, the characteristics of the majority are
easily available in literature (Benjamin and Cornell, 1970). The triangular
distribution has been discussed in Megill (1977). Details of the Johnson
system of distributions are presented in McGrath and Irving (1973) and
Johnson and Kotz (1970). Additional details for each of these distributions
are presented in the following discussion.
In some cases, it may be desirable to include correlations among the
variables. For example there may be correlation between hydraulic
conductivity and particle size or between adsorption and degradation
coefficients. The uncertainty processor allows the generation of (linearly)
correlated variables for cases where the underlying distribution of the
variables is either normal and/or lognormal.
6.3 DESCRIPTION OF AVAILABLE PARAMETER DISTRIBUTIONS
The Monte Carlo Shell has the ability to generate data from a number of
probability distributions listed above. A description of each of these
distributions is provided in the following paragraphs, including parameters
of the distributions, equations for the probability and cumulative density
functions, and a brief discussion of the properties of each distribution.
6.3.1 Uniform Distribution
A uniform distribution is a symmetrical probability distribution in which all
values within a given range have an equal chance of occurrence. A uniform
180
-------
distribution is completely described by two parameters: 1) the minimum value
(lower bound) A, and 2) the maximum value (upper bound) B. The equation for
the uniform probability density distribution of variable x is given by:
fu(x) = 1/(B - A)
(6-4)
where fu(x) is the value of the probability density function for x. The
cumulative distribution F(x) is obtained by integrating Equation (6-4). This
yields the probability distribution:
Fu(x) - (x - A)/(B - A)
(6-5)
where F(x) is the probability that a value less than or equal to x will
occur.
6.3.2 Normal Distribution
The term normal distribution refers to the well known bell-shaped probability
distribution. Normal distributions are symmetrical about the mean value and
are unbounded, although values further from the mean occur less frequently.
The spread of the distribution is generally described by the standard
deviation. The normal distribution has only two parameters: the mean and
the standard deviation. The probability density function of x is given by:
fn(x) — exp
x -
-0.5
(6-6)
where Sx is the standard deviation, and n^ is the mean of x. The cumulative
distribution is the integral of the probability density function:
Fn(x) - Jf fn(x)dx
(6-7)
The above integration must be performed numerically, but tables of
numerically-integrated values of Fn(x) are widely available in the
statistical literature.
6.3.3 Log-Normal Distribution
The log-normal distribution is a skewed distribution in which the natural log
of variable x is normally distributed. Thus, if y is the natural log of x,
181
-------
then the probability distribution of y is normal with mean n^, and standard
deviation S and a probability density function similar to Equation (6-10).
The mean and standard deviation of x (n^ and Sx) are related to the log-
normal parameters IIL^ and S as follows:
n^ - exptniy + 0.5(Sy)2] (6-8)
Sx2 = n^2 [exp(Sy2) - 1] (6-9)
To preserve the observed mean and standard deviation of x, the parameters of
the log-normal distribution (IIL^ and S ) are therefore selected such that the
above relationships are satisfied. Note that m^ and S do not equal the
natural logs of n^ and Sx respectively. Log-normal distributions have a
lower bound of 0.0 and no upper bound, and are often used to describe
positive data with skewed observed probability distributions.
6.3.4 Exponential Distribution
The probability density function for an exponential distribution is described
by an exponential equation:
exp(-x/mx)
fe
where n^ is the mean of x. The cumulative distribution is given by:
Fe(x) = 1 - exp(-x/mx) (6-11)
The exponential distribution is bounded by zero; the probability density
function peaks at zero and decreases exponentially as x increases in
magnitude .
6.3.5 The Johnson System of Distributions
The Johnson system involves two main distribution types: SB (Log-ratio or
bounded) and SU (unbounded or hyperbolic arcsine) . These two distribution
types basically represent two different transformations applied to the random
variable such that the transformed variable is normally distributed. The
specific transformations are:
SB: Y = In \^-3-\ (6-12)
182
-------
where
In — natural logarithm transformation
x = untransformed variable with limits of variation from A to B.
Y = the transformed variable with a normal distribution
Selection of a particular Johnson distribution for sample data set is
accomplished by plotting the skewness and kurtosis of the sample data. The
location of the sample point indicates the distribution for the sample data.
For additional details of the Johnson system of distributions , the reader is
referred to McGrath and Irving (1973) and Johnson and Kotz (1970) .
6.3.6 Triangular Distribution
A triangular distribution is a relatively simple probability distribution
defined by the minimum value, the maximum value, and the most frequent value
(i.e., the mode). Figure 6.1 shows an example triangular probability density
function. The cumulative distribution for values of x less than the most
frequent value, x^, is given by:
(x - xx)
F(x) - - - r-, - r- (6-14)
' xl><*2 - xl>
where x-^ is the minimum value and X2 is the maximum value . For values of x
greater than the most frequent value, the cumulative distribution is:
F(x) - 1 (6-15)
6.3.7 Empirical Distribution
At times it may be difficult to fit a standard statistical distribution to
observed data. In these cases it is more appropriate to use an empirical
piecewise-linear description of the observed cumulative distribution for the
variable of interest.
Cumulative probabilities can be estimated from observed data by ranking the
data from lowest (rank = 1) to highest (rank = number of samples) value. The
183
-------
Figure 6.1. Triangular Probability Distribution
184
-------
cumulative probability associated with a value of x is then calculated as a
function of the rank of x and the total number of samples. The cumulative
probabilities of values between observed data can be estimated by linear
interpolation.
6.3.8 Uncertainty in Correlated Variables
In many cases model input variables are correlated due to various physical
mechanisms. Monte Carlo simulation of such variables requires not only that
parameters be generated from the appropriate univariate distributions, but
also that the appropriate correlations be preserved in the generated input
sequences. The Monte Carlo module currently has the ability to generate
correlated normal, log-normal, Johnson SB, and Johnson SU numbers; the
procedures used are described in the following paragraphs.
The correlation coefficient is a measure of the linear dependence between two
random variables and is defined as:
x y
x-y PX py
where
px v - the correlation coefficient between random variables x any y
cov(x,y) — the covariance of x and y as defined below
pK, pv — the standard deviation for x and y.
The covariance of x and y is defined as:
cov(x.y) - E [(x-mx)(y-n^)]
-HO
- J J (x-ra^ (y-my) fXjy(x,y) dx dy (6-17)
where
E - the expected value
n^, IIL^ - the mean of the random variables x and y
fx y(x,y) - the joint probability distribution of x and y.
Note that the linear correlation coefficient between x and y can be computed
using
185
-------
n
I xi Yi - n x y
- -
X (x^ - nx2) X (yi2 - ny2)
i-1 1=1
To generate correlated random variables three steps are required. First
uncorrelated normally distributed random numbers are generated. This vector
is then transformed to a vector of normally distributed numbers with the
desired correlation. Finally, the normally distributed numbers are
transformed to numbers with the desired distribution.
The transformation of uncorrelated to correlated normal numbers consists of
multiplying the uncorrelated vector of numbers with a matrix B:
Y' = B e (6-19)
where
e - the vector of uncorrelated, normally distributed random numbers.
B - and N by N matrix
Y' — a vector of standard normal deviates of mean zero and standard
deviation of unity.
The matrix B is related to the variance-covariance matrix S as follows:
S - B BT (6-20)
rri
where B is the transpose of the B matrix. Since the normal variables Y'
have means of zero and unit variances, the variance-covariance matrix is
equivalent to the correlation matrix.
Thus, if the correlation matrix S is known, B can be found from Equation (6-
20) by using a Choleski decomposition algorithm. This algorithm will
decomposes a symmetric positive definite matrix, such as S, into a triangular
matrix such as B (de Marsily, 1986, p. 381).
Having generated a vector of correlated normally distributed random numbers,
the vector Y' can be converted, through appropriate transformations, to the
distribution of choice. Thus for parameters X^ that have a normal
distribution, the Y' numbers are transformed as follows:
Xi = "be + o^i) (6-21)
186
-------
For parameters that follow the lognormal distribution, the following
transformation applies:
Xi- exp[(Y!) (1] (6-22)
where
£ - the log mean of the i parameter.
-i £ — the log standard deviation of the i parameter
For parameters with Johnson SB and SU distributions, the Y' are first
transformed to normally distributed variables Y with mean H^ and standard
deviation a :
Yj_ = My + ay Y[ (6-23)
Johnson SB numbers are then computed from Y^ as follows :
Xj_ = (B exp(Yi) - A)/(l + exp(Yi)) (6-24)
Johnson SU numbers are computed by:
Xj_ - A + (B-A)[exp(Yi) - exp(-Yi)]/2 (6-25)
Other distributions can be easily incorporated into the analyses at a later
time when suitable transformations from the normal distribution can be found.
It is important to note that in using this technique, the correlations are
maintained in normal space so if these correlations are estimated using
actual data, the data should be transformed to a normal distribution before
correlation coefficients are estimated.
For two correlated variables, one with a normal distribution (X2) and the
other with a log normal distribution (x^) , the following equation is used to
transform correlations to normal space (Meija and Rodriguez -Iturbe, 1974):
(6'26)
187
-------
where
a = the correlation coefficient between the two variables in the
normal space
°v, v~= the correlation coefficient between the two variables in the
X1'X2
arithmetic space
2
CT — the variance of y-, derived from Equation (6-9)
If both x^ and X£ are log-normally distribution then the correlation
coefficient is transformed using Meija and Rodriguez-Iturbe (1974):
^ 1 + *Xl.Xo L „ I (6-27)
where the relationships between S (Sx ) and S (S ) are given by Equations
(6-8) and (6-9). ±2 Jl J2
Thus, for log-normal variables the user enters the values of the correlation
coefficients in log-normal space; Equations 6-26 and 6-27 are then used to
transform the correlation coefficients into normal space.
No direct transformation of Johnson SB or SU correlations to normal
correlations is currently known. For these distributions the user must
therefore supply the correlation coefficients between normal-transformed
numbers. This may be accomplished by first transforming Johnson SB and SU
data to normal data using Equations 6-12 and 6-13. The covariance matrix S
is then derived using only normal, log-normal, and normal-transformed SB and
SU data.
6.3.9 Generation of Random Numbers
Having selected the distribution for the various input parameters, the next
step is the generation of random values of these parameters. This requires
the use of pseudo-random number generating algorithms for Normal and Uniform
numbers. There exist numerous proprietary as well as non-proprietary
subroutines that can be used to generate random numbers. A number of these
are comparable in terms of their computational efficiency, accuracy, and
precision. The performance of the algorithms included in this preprocessor
has been checked to ensure that they accurately reproduce the parameters of
the distributions that are being sampled (Woodward-Clyde Consultants, 1988).
6.4 ANALYSIS OF OUTPUT AND ESTIMATION OF DISTRIBUTION QUANTILES
Model output generally will consist of a volume of data that represents a
sample of outcomes. Given the natural variability and the uncertainty of
various model components, there will be variability in the output. All of
the factors that were allowed to vary within the model contribute to
variability in model predictions. Taken as a whole, the model output depicts
188
-------
possible events in terms of their relative frequency of occurrence. Values
produced by the model generally are treated as if they were observations of
real field events. In interpreting these values, it is important to maintain
the perspective dictated by the design and scope of the study.
Model output can be analyzed in various ways depending upon current
objectives. Many features of the distribution may be characterized. Quite
often, for example, it is of interest to estimate certain quantiles or
percentiles of the distribution. Since the model output is treated as a
sample from an unknown parent population, the methods of statistical
inference normally are used to estimate distribution parameters and to
associate with them measures of uncertainty.
One of the most frequently asked questions concerns the number of samples
required for some given purpose. In modelling, this translates into the
number of model runs needed. For the most part, since methods of basic
inference are being applied in a Monte Carlo framework, resulting model
output values are treated as observations forming a random sample. The
sample size required to estimate a given parameter depends on a number of
factors. These include the nature of the parameter that is being estimated,
the form of the underlying distribution, the variability in the observations,
the degree of precision and/or accuracy desired, the level of confidence to
be associated with the estimate, and the actual statistical estimator used to
provide the estimate.
Generally, if the output distribution is to be accurately characterized with
respect to its many features, the number of model runs needed will be higher
than if only a few parameters are to be estimated. The simulation strategy
should be determined by the issues addressed by the modelling effort. It may
be important, for example, to estimate the extreme upper percentiles of the
output distribution. In this case, the choice of simulation design should
account for the relative difficulty of obtaining such estimates. If it is
not known exactly how the data will be utilized, then the problem becomes one
of establishing a distributional representation that is as good as possible
under the most extreme usage or estimation scenario. For example, if only a
distribution mean were to be estimated, the sample size required could be
determined without concern for estimating, say, the 99th percentile.
6.4.1 Estimating Distribution Quantiles
In the following section, a summary is given for statistical techniques used
to estimate distribution quantiles. There are many such methods available to
estimate a given percentile of an unknown distribution on the basis of sample
data. In the RUSTIC code four such methods can be used. Among these are
distribution-free or nonparametric techniques as described below. Others
include methods specific to certain distributions that assume a knowledge of
the distributional form. First, the point estimators are given, then the
method for constructing a confidence interval is briefly described.
The order statistics of a sample are merely the ordered values denoted by
x(l)' x(2) x(n)' wnere n represents the sample size. The empirical
cdf can be defined simply as
139
-------
g(x) =
0, if x^ < x,
1/n, if x^ < x < X(i+1), for i-1 n-1 (6-28)
1, if x > x(n).
Mathematically, g(x) is a step function, discontinuous at each value
By definition, the lOOp-th percentile (i.e., the p-level quantile) is given
by XL where
p = Pr{ X < Up } (6-29)
If F(x) denotes the cumulative distribution function,
p - F(Up) and Up = F'^p) (6-30)
When only sample information is available, u^ is unknown, but it can be
estimated by forming an appropriate function of the observations .
Nonparametric point estimates of UL can be constructed as linear combinations
of the order statistics. In particular, each of Y^ through Yo below is an
estimator of vu. Let [z] denote the largest integer less than or equal to z.
Define
j = [np], g = np - j (6-31)
i - [np + 0.5], r = (np + 0.5) - i (6-32)
k - [(n+l)p], h = (n+l)p - k (6-33)
Then,
Yl -
Y2 -- - — , if g=0 (6-35)
' if
Y3 - (0.5+i-np) X(i) + (0.5-i+np) X(i+1) (6-36)
= (1 - r) X(i) + r X(1+1)
In each of these definitions, only the values of n and p determine which
order statistics are used in forming an estimate of u^. Thus, the estimators
190
-------
do not depend on the underlying distributions. However, the relative
performance of these estimators is dependent upon several criteria involving
the level p, the sample size n, the type of parent distribution from which
samples are drawn, estimator bias, and the mean squared error. If the sample
size is very large, the differences among the estimates are not very great.
Of the estimators available, the three shown above exhibit the best
performance in relatively small samples (n<50) from normal and lognormal
distributions.
Another simple estimator used in the model is calculated by constructing the
cdf of the output
F(x) = i/n (6-37)
in which i is the rank of the outcome in the sample. The specific quantile
of interest is then determined by interpolation.
6.4.2 Confidence of u^
Approximate confidence statements can be placed on XL by selecting
appropriate order statistics to serve as the upper and lower confidence
bounds. The rationale is given as follows.
For a given distribution, the value u_ is such that exactly 100p% of all
values of this distribution are less than u^, and 100(1-p)% exceed this
value. An individual value selected randomly from the distribution has
probability p of being less than u^. In a random sample of size n from this
distribution, the probability of not exceeding LL remains constant for each
individual element of the sample. Thus, the number of values in the sample
that are less than or equal to XJL is distributed binomially. The probability
that the random interval (X/^\, X/j.i\) will contain u^ is equivalent to the
\J/ \ J * "^ ^ P
probability that exactly i of the n elements of the sample will be less
than VL. Hence, this probability is
p1 (l-p)"'1 (6-38)
which is a simple binomial probability.
This expression can be calculated for each pair of consecutive order
statistics X/j\, X/^+^s, for i— 1, ..., n-1. However, it is more convenient
to deal with these several intervals by calculating cumulative probabilities
of the form
i0 ( j ) PJ
<6-39>
191
-------
For practical convenience, the normal approximation
F {[(i+0.5)-np]/y[np(l-p)]} (6-40)
can be used, where F represents the cdf of the standard normal distribution.
All of this is utilized for determining two order statistics, denoted below
with subscripts i and j, with the property
Pr{X(i) < Up < X(j)) = 1 - a (6-41)
where 1-a is the predetermined confidence coefficient; typically, 1-a = 0.95.
Computationally, i and j can be determined by solving the equations
a/2 = F{[(i+0.5)-np]//[jipTl-p)]} (6-42)
and
1 - a/2 = F{[(j+0.5)-np]/y[np(l-p)]} (6-43)
This results in
i - (np-0.5) + y[np(l-p)] F'1(a/2) (6-44)
j = (np-0.5) + 7[np(l-p)] F'1(l-a/2) (6-45)
where F denotes the inverse cdf of the standard normal distribution (e.g.,
for 1-a - 0.90, F'1(l-a/2) = 1.645). For example, with n-100, p-0.95, and
l-a=0.90, i=90 and j=98, so that (X/90), X/ggO forms the approximate 90%
confidence interval on u^.
Although the expressions for the confidence interval do not depend in any way
on the underlying distribution, the expected width of the interval does. In
particular, it depends on the expected values of the order statistics
involved. In the example above, if the sample is from a standard normal
distribution, IL = 1.645 and the expected half-width of the interval is
0.349. If the sample is from a lognormal distribution based on a standard
normal, u^ - 5.180 and the expected half-width is 1.858. Also, note that in
normal sampling the expected confidence interval half-width for n-500 is
0.192 for the same estimate.
192
-------
SECTION 7
REFERENCES
Avdonin, N.A. 1964. Some Formulas for Calculating the Temperature Field of
a Stratom during Termal Injection. Izvestiya Vysshikh Uchebnykl
Zavedenii Neft. Gaz. 7(3) 37-41.
Bassett, D.L., and D.W. Fitzsimmons. 1974. A Dynamic Model of Overland Flow
in Border Irrigation. Am. Soc. Ag. Eng. Paper No. 74-2529, St. Joseph,
MI.
Battelle, Pacific Northwest Laboratories, and GeoTrans. 1988. FASTCHEM
Package Volume 2: User's Guide to the EFLOW Groundwater Flow Code.
Electrical Power Rsearch Institute, EA-5870-CCM, Volume 2, Research
Project 2485-2.
Bear, J. 1979. Hydraulics of Groundwater. McGraw-Hill, New York, NY.
Benjamin, J.R. and C.A. Cornell. 1970. Probability. Statistics, and
Decision for Civil Engineers. McGraw Hill Book Company.
Bromilow, R.H., M. Richard and M. Leistra. 1980. Pesticide Sci. 11:389-395.
Brooks, R.H. and A.T. Corey. 1966. Properties of Porous Media Affecting
Fluid Flow. ASCE J. Irrig. Drain Div. 92 (IR2):61-68.
Brown, S.M., and S.H. Boutwell. 1986. Chemical Spill Exposure Assessment
Methodology. RP 2634-1. Prepared for Electric Power Research Insti-
tute, Palo Alto, CA.
Carsel, R.F., andR.S. Parrish. 1988. Developing Joint Probability
Distributions of Soil-Water Retention Characteristics. Water Resour.
Res. 24(5):755-769.
Carsel, R.F., L.A. Mulkey, M.N. Lorber, and L.B. Baskin. 1985. The
Pesticide Root Zone Model (PRZM): A Procedure for Evaluating Pesticide
Leaching Threats to Ground Water. Ecological Modeling 30:49-69.
Carsel, R.F., C.N. Smith, L.A. Mulkey, J.D. Dean, and P. Jowise. 1984.
User's Manual for the Pesticide Root Zone Model (PRZM): Release 1.
U.S. Environmental Protection Agency, Environmental Research Laboratory,
Athens, GA. EPA-600/3-84-109.
193
-------
Chen, C.W., S.A. Cheriai, R.J.M. Hudson, andJ.D. Dean. 1983. The
Integrated Lake-Watershed Acidification Study. Vol. I. Model
Principles and Application Procedures. Electric Power Research Inst.
Rep. No. EA-3221, Vol. 1. Palo Alto, CA.
Cruse, R.M., D.R. Linden, J.K. Radke, W.E. Larson, and K. Larntz. 1980. A
Model to Predict Tillage Effects on Soil Temperature. Soil Sci. Soc.
Am. J. 44:378-383.
Davis, L.A., and S.P. Neuman. 1983. Documentation and User's Guide: UNSAT2
- Variably Saturated Flow Model. U.S. Nuclear Regulatory Commission
Report, NUREG/CR-3390, Washington, D.C.
de Marsily, G. 1986. Quantitative Hydrogeology: Groundwater Hydrology for
Engineers. Academic Press, San Diego, CA.
de Vries, D.A. 1963. Thermal Properties of Soils. In: Physics of the
Plant Environment, ed. W.R. van Wijk, pp 210-235. New York: John Wiley
and Sons, Inc.
Dean, J.D., and D.F. Atwood. 1985. Exposure Assessment Modeling for
Aldicarb in Florida. EPA Contract No. 68-03-3116. Work Assignment No.
23, Final Report.
Donigian, A.S., Jr., C.S. Raju, andR.F. Carsel. 1986. Impact: of
Conservation Tillage on Environmental Pesticide Concentrations in Three
Agricultural Regions. Prepared for U.S. EPA, Integrated Environmental
Management Division, Office of Policy Analysis, Washington, D.C.
Dyer, A.J. 1974. A Review of Flux-Profile Relationships. Boundary-Layer
Meteorology 7:363-372.
Dyer, A.J., and B.B. Hicks. 1970. Flux-Gradient Relationships in the
Constant Flux Layer. Quarterly Journal of the Royal Meteorological
Society 96:715-721.
Farrell, D.A., E.L. Greacen, and C.G. Gurr. 1966. Vapor Tramjfer in Soil
Due to Air Turbulence. Soil Sci. 102:305-313.
Federal Register. 1984. Proposed Guidelines for Exposure Assessments:
Request for Comments. U.S. Environmental Protection Agency. Vol. 49
(227), November 23.
Fok, Y.S., and A.A. Bishop. 1965. Analysis of Water Advance In Surface
Irrigation. ASCE of Irrig. and Drainage 91(IR1):99-116.
Fukuda, H. 1955. Air and Vapor Movement in Soil Due to Wind Gustiness.
Soil Science 79:249-258.
Grover, R., S.R. Shewchuk, A.J. Cessna, A.E. Smith, andJ.H. Hunter. 1985.
Fate of 2,4-D Iso-Octyl Ester after Application to a Wheat Field. J.
Environ. Qual. 14:203-210.
194
-------
Gupta, S.C., W.E. Larson, and R.R. Allmaras. 1984. Predicting Soil
Temperature and Soil Heat Flux under Different Tillage-Surface Residue
Conditions. Soil Sci. Soc. Am. J. 48:223-232.
Gupta, S.C., W.E. Larson, and D.R. Linden. 1983. Tillage and Surface
Residue Effects on Soil Upper Boundary Temperatures. Soil Sci. Soc. Am.
J. 47:1212-1218.
Gupta, S.C., J.K. Radke, W.E. Larson, andM.J. Shaffer. 1982. Predicting
Temperature of Bare and Residue-Covered Soils from Daily Maximum and
Minimum Air Temperatures. Soil Sci. Soc. Am. J. 46:372-376.
Gupta, S.C., J.K. Radke, and W.E. Larson. 1981. Predicting Temperature of
Bare and Residue-Covered Soils with and without a Corn Crop. Soil Sci.
Soc. Am. J. 45:405-412.
Hadermann, J. 1980. Radionuclide Transport Through Heterogeneous Media.
Nuclear Technology 47:311-323.
Haith, D.A., and R.C. Loehr, eds. 1979. Effectiveness of Soil and Water
Conservation Practices for Pollution Control. U.S. Environmental
Protection Agency, Athens, GA. EPA-600/3-79-106.
Hanks, R.J., D.D. Austin, and W.T. Ondrechen. 1971. Soil Temperature
Estimation by a Numerical Method. Soil Sci. Soc. Am. Proc. 35:665-667.
Harper, L.A., A.W. White, Jr., R.R. Bruce, A.W. Thomas, and R.A. Leonard.
1976. Soil and Microclimate Effects on Trifluralin Volatilization.
J. Environ. Qual. 5:236-242.
Hasfurther, V.R., and R.D. Burman. 1974. Soil Temperature Modeling Using
Air Temperature as a Driving Mechanism. Trans. ASAE 17:78-81.
Hunt, B. 1972. Dispersive Sources in Uniform Ground-Water Flow. ASCE
Journal of the Hydraulics Division. 104(HYI).
Huyakorn, P.S., H.O. White, Jr., J.E. Buckley, and T.D. Wadsworth. 1988a.
VADOFT: Finite Element Code for Simulating One-Dimensional Flow and
Solute Transport in the Vadose Zone. Project Report to Woodward-Clyde
Consultants. February.
Huyakorn, P.S., and J.E. Buckley. 1988b. SAFTMOD: Saturated Zone Flow and
Transport Two-Dimensional Finite Element Model. Project Report to
Woodward-Clyde Consultants. March.
Huyakorn, P.S., B.C. Jones, and P.F. Andersen. 1986a. Finite Element
Algorithms for Simulating Three-Dimensional Groundwater Flow and Solute
Transport in Multilayer Systems. Water Resour. Res. 22(3)361-374.
Huyakorn, P.S., E.P. Springer, V. Guvanasen, andT.P. Wadsworth. 1986b. A
Three-Dimensional Finite-Element Model for Simulating Water flow in
Variably Saturated Porous Media. Water Resour. Res. 22(13)1790-1808.
195
-------
Huyakorn, P.S., J.W. Mercer, and D.S. Ward. 1985. Finite Element Matric and
Mass Balance Computational Schemes for Transport in Variably Saturated
Porous Media. Water Resour. Res. 21(3)346-358.
Huyakorn, P.S., S.D. Thomas, and B.M. Thompson. 1984a. Techniques for
Making Finite Elements Competitive in Modeling Flow in Variably
Saturated Porous Media. Water Resour. Res. 20(8):1099-1115.
Huyakorn, P.S., A.G. Kretschek, R.W. Broome, J.W. Mercer, and B.H. Lester.
1984b. Testing and Validation of Models for Simulating Solute
Transport in Ground-Water. International Ground Water Modeling Center,
Holcomb Research Inst., HRI No. 35.
Huyakorn, P.S., andG.F. Finder. 1983. Computational Methods in Subsurface
Flow. Academic Press. 473 pp.
Johnson, N.L., and S. Kotz. 1970. Distributions in Statistics: Continuous
Univariate Distributions. Boston: Houghton Mifflin Company.
Jones, R.L. 1983. Movement and Degradation of Aldicarb Residues in Soil and
Ground Water. Presented at the CETAC Conference on Multidisciplinary
Approaches to Environmental Problems, November 6-9, Crystal City, VA.
Jones, R.L., P.S.C. Rao, and A.G. Hornsby. 1983. Fate of Aldicarb in
Florida Citrus Soil 2. Model Evaluation. Presented at the Conference
on Characterization and Monitoring of Vadose (Unsaturated) Zone,
December 8-10, Las Vegas, NV.
Jury, W.A., W.F. Spencer, and W.J. Farmer. 1984. Behavior Assessment Model
for Trace Organics in Soil: III. Application of Screening Model. 13:
573-579.
Jury, W.A., R. Grover, W.F. Spencer, and W.J. Farmer. 1983a. Behavior
Assessment Model for Trace Organics in Soil: I. Model Description.
J. Environ. Qual. 12:558-564.
Jury, W.A., R. Grover, W.F. Spencer, and W.J. Farmer. 1983b. Use of Models
for Assessing Relative Volatility, Mobility, and Persistence of
Pesticides and Other Trace Organics in Soil Systems. In: Harvard
Assessment of Chemicals, ed. J. Saxena, pp. 1-43. New York: Academic
Press.
Khalell, R., and D. Reddell. 1986. MOC Solutions of Convective-Dispersion
Problems. Ground Water 24(6):798-807.
Konikow, L.F., and J.D. Bredehoeft. 1978. Computer Model of Two-Dimensional
Solute Transport and Dispersion in Ground Water. Techniques of Water-
Resources Investigations of the United States Geological Survey, Book 7,
Chapter C2.
Kreysig, E. 1972. Advanced Engineering Mathematics, pp. 430-434. John
Wiley and Sons Inc., New York, N.Y.
196
-------
Lester, B.H., P.S. Huyakorn, H.O. White, Jr., T.D. Wadsworth, andJ.E.
Buckley. 1986. Analytical Models for Evaluating Leachate Migration in
Groundwater Systems. Technical report prepared by GeoTrans, Inc., for
U.S. Environmental Protection Agency, Office of Solid Waste, Washington,
D.C. August.
Li, R. , D.B. Simons, and M.A. Stevens. 1975. Non-linear Kinematic Wave
Approximation for Water Routing. Water Res. Res. 11(2):245-252.
Liley, P.E., and W.R. Gambill. 1973. Section 3: Physical and Chemical
Data. In: Chemical Engineers Handbook, eds. P.H. Perry and C.H.
Chilton. New York: McGraw Hill.
Lo, A.K. 1977. An Analytical-Empirical Method for Determining the Roughness
Length and Zero-Plane Displacement. Boundary-Layer Meteorology 12:141-
151.
M. Baptista, A.E., E.E. Adams, and K.D. Stolzenbach. 1984. Eularian-
Lagrangian Analysis of Pollution Transport in Shallow Water. MIT Report
296.
McGrath, E.J. and D.C. Irving. 1973. Techniques for Efficient Monte Carlo
Simulation, Volume II. Random Number Generation for Selected
Probability Distributions. Report prepared for Office of Naval
Research. Project No. NR 366-076/1-5-72, Code 462.
Megill R. E. 1977. An Introduction to Risk Analysis. Tulsa, Oklahoma:
Petroleum Publishing Company.
Mehlenbacher, L.A., and D.W.A. Whitfield. 1977. Modeling Thermal Eddy
Diffusivity at Canopy Height. Boundary-Layer Meteorology 12:153-170.
Meija, J.M., and I. Rodriguez-Iturbe. 1974. On the Synthesis of Random
Field Sampling From the Spectrum: An Application to the Generation of
Hydraulic Spatial Processes. Water Resour. Res. 10(4), 705-712.
Mockus, V. 1972. Estimation of Direct Surface Runoff from Storm Rainfall.
In: National Engineering Handbook. Section IV, Hydrology. U. S. Soil
Conservation Report NEH-Notice 4-102. August.
Mualem, Y. 1976. A New Model for Predicting the Hydraulic Conductivity of
Unsaturated Porous Media. Water Resour. Res. 12(3):513-522.
Neuman, S.P., R.A. Feddes, and E. Bresler. 1974. Finite Element Simulation
of Flow in Saturated-Unsaturated Soils Considering Water Uptake by
Plants. Hydrodynamics and Hydraulic Engineering Laboratory Report for
Project No. ALO-SWC-77, Haifa, Israel.
Neumen, S.P., and P.A. Witherspoon. 1969a. Theory of Flow in a Confined
Two-Aquifer System. Water Resource Res. 5(3), 803-816.
197
-------
Neuman, S.P., and P.A. Witherspoon. 1969b. Applicability of Current
Theories of Flow in Leaky Aquifers. Water Resour. Res. 5(4), 817-829.
Ogata, A., and R.B. Banks. 1961. A Solution of the Differential Equation of
Longitudinal Dispersion in Porous Media. U.S. Geological Survey
Professional Paper No. 411-A.
Oliver, H.R. 1971. Wind Profiles in and above a Forest Canopy. Quarterly
Journal of the Royal Meteorological Society 97:548-553.
Papadopulos, I.S. 1965. Nonsteady Flow to a Well in an Infinite Anisstropic
Aquifer. International Association of Scientific Hydrology Symposium of
of Dubrovnik.
Parmele, L.H., E.R. Lemon, and A.W. Taylor. 1972. Micrometeorological
Measurement of Pesticide Vapor Flux from Bare Soil and Corn under Field
Conditions. Water, Air, and Soil Pollution 1:433-451.
Parton, W.J. 1984. Predicting Soil Temperature in a Shortgrass Steppe.
Soil Sci. 138(2):93-101.
Perkins, H.C. 1974. Air Pollution. New York: McGraw-Hill, Inc.
Phillip, J.R. 1955. Numerical Solution of Equations of the Diffusion Type
with Diffusivity Concentration-Dependent. Trans. Faraday Soc.
51:885-892.
Pinder, G.F., and W.G. Gray. 1977. Finite Element Simulation in Surface and
Subsurface Hydrology. Academic Press, London.
Richards, L.A. 1931. Capillary Conduction of Liquids through Porous
Mediums. Physics 1:318-333.
Rosenberg, N.J. 1974. Microclimate: The Biological Environment. New York:
Wiley Interscience.
Salhotra, A.M. 1986. A Coupled Heat, Salt and Water Balance Model of
Evaporation and Stratification in Saline Terminal Lakes: An Application
to the Dead Sea. Ph.D. Dissertation. Massachusetts Institute of
Technology, Department of Civil Engineering.
Sauty, J.P. 1980. An Analysis of Hydro Dispersive Transfer in Aquifers.
Water Resour. Res. 16(1)145-158.
Scotter, D.R., and P.A.C. Raats. 1970. Movement of Salt and Water near
Crystalline Salt in Relatively Dry Soil. Soil Sci. 109:170-178.
Shamir, U.Y., and D.R.F. Harleman. 1967. Dispersion in Layered Porous
Media. J. of the Hydraulics Division, Amer. Soc. Civ. Engr.
93(HY5):237-260.
198
-------
Smiles, D.E., J.R. Phillip, J.H. Knight, andD.E. Elrick. 1978.
Hydrodynamic Dispersion During Absorption of Water by Soil. Soil Sci.
Soc. Am. J. 42:229-234.
Smith, C.N., R.A. Leonard, G.W. Langdale, and G.W. Baily. 1978. Transport
of Agricultural Chemicals from Small Upland Piedmont Watersheds. U.S.
EPA, Environmental Research Laboratory, Athens, GA. EPA-600/3-78-056.
Springer, E.P., and H.R. Fuentes. 1987. Modeling Study of Solute Transport
in the Unsaturated Zone. U.S. Nuclear Regulatory Commission, NUREC/CR-
4615, LA-10730-MS, Volume 2.
Stamper, J.H., H.N. Nigg, andJ.C. Allen. 1979. Organophosphate Insecticide
Disappearance from Leaf Surfaces: An Alternative to First-Order
Kinetics. Environ. Sci. and Tech. 13:1402-1405.
Stanhill, G. 1969. A Simple Instrument for the Field Measurement of Tur-
bulent Diffusion Flux. J. Applied Meteorology 8:504-513.
Stewart, B.A., D.A. Woolhiser, W.H. Wischmeier, J.H. Caro, and M.H. Fere.
1976. Control of Water Pollution from Cropland, Volume II: An Over-
view. U.S. Environmental Protection Agency, Athens, GA. EPA-600/2-75-
026b.
Streile, G.P. 1984. The Effect of Temperature on Pesticide Phase Parti-
tioning, Transport, and Volatilization from Soil. Ph.D. dissertation,
University of California, Riverside, CA.
Szeicy, G., G. Endrodi and S. Tajchman. 1969. Aerodynamic and Surface
Factors in Evaporation. Water Resour. Res. 5(2) .'380-394.
Taylor, A.W. 1978. Post-Application Volatilization of Pesticides under
Field Conditions. Journal of the Air Pollution Control Association
28:922-927.
Thibodeaux, L.J., C. Spencer, and L.M. Riley. 1982. Models of Mechanisms
for the Vapor Phase Emission of Hazardous Chemicals from Landfills.
Journal of Hazardous Materials 7:63-74.
Thibodeaux, L.J. 1979. Chemodynamics: Environmental Movement of Chemicals
in Air. Water, and Soil. New York: John Wiley & Sons.
Thorn, A.S., J.B. Stewart, H.R. Oliver, and J.H.C. Gash. 1975. Comparison of
Aerodynamic and Energy Budget Estimates of Fluxes over a Pine Forest.
Quarterly Journal of the Royal Meteorological Society 101:93-105.
van Bavel, C.H.M., and D.I. Hillel. 1976. A Simulation Study of Soil Heat
and Moisture Dynamics as Affected by a Dry Mulch. In: Proceedings of
1975 Summer Computer Simulation Conference, San Francisco, California,
pp. 815-821. Simulation Councils, Inc., La Jolla, CA.
199
-------
van Bavel, C.H.M. and D.I. Hillel. 1975. Calculating Potential and Actual
Evaporation from a Bare Soil Surface by Simulation of Concurrent Flow of
Water and Heat. Agri. Meteorol. 17:453-476.
van Genuchten, M. T., and W.J. Alves. 1982. Analytical Solutions of the
One-Dimensional Convective-Dispersive Solute Transport Equation. U.S.
Technical Bulletin No. 1661. 151 pp.
van Genuchten, M. T. 1976. A Closed-form Equation for Predicting the
Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. J.
44:892-898.
Wagenet, R.J., and J.W. Biggar. 1987. Measurement and Interpretation of
DBCP Fate in Agricultural Soils. Journal of Environmental Quality (in
press).
Wagenet, R.J., and J.L. Hutson. 1987. Leaching Estimation and Chemistry
Model. Department of Agronomy, Cornell University, Ithaca, NY.
Wark, K., and C.F. Warner. 1976. Air Pollution: Its Origin and Control.
New York: Harper and Row Publishers.
White, Jr., A.W., L.A. Harper, R.A. Leonard, and J.W. Turnball. 1977.
Trifluralin Volatilization Losses from a Soybean Field. Journal of
Environmental Quality 6:105-110.
Wilke, 0., and E.T. Smerdon. 1965. A Solution of the Irrigation Advance
Problem. ASCE Journal of Irrigation and Drainage 91(IR3):23-24.
Williams, J.R., C.A. Jones, and P.T. Dyke. 1984. A Modeling Approach to
Determining the Relationship Between Erosion and Soil Productivity.
Trans. ASAE 27:129-144.
Williams, J.R. 1975. Sediment Yield Prediction with Universe Equation Using
Runoff Energy Factor In: Present and Prospective Technology for
Predicting Sediment Yields and Sources. U.S. Department of Agriculture.
AR5-5-40.
Willis, G.W., L.L. McDowell, L.A. Harper, L.M. Southwick, and S. Smith.
1983. Seasonal Disappearance and Volatilization of Toxaphene and DDT
from a Cotton Field. J. Environ. Qual., 12:80-85.
Wilson, J.L., and P.J. Miller. 1978. Two-Dimensional Plume in Uniform
Groundwater Flow. J. Hydraulics Div., ASCE 104,HY4.
Woodward-Clyde Consultants. 1988. Multimedia Exposure Assessment Model for
Evaluating the Land Disposal of Hazardous Wastes. EPA Contract No. 68-
03-6304. Work Assignment No. 1. Draft Report.
Wunderlich, W.O. 1972. Heat and Mass Transfer Between a Water Surface and
the Atmosphere. Laboratory Rpt. No. 14 T.V.A. Eng. Lab., Norris, TN.
Yeh, G.T., and D.S. Ward. 1981. FEMWASTE: A Finite-Element Model of Waste
Transport Through Saturated-Unsaturated Porous Media. Rep. ORNL-5601,
Oak Ridge National Laboratory, Oak Ridge, TN. 137 pp.
200
-------
SECTION 8
APPENDICES
8.1 QUALITY ASSURANCE
The quality assurance program for this model development program effort was
broken down into two phases:
• Algorithms
• Software
The quality assurance aspects for these entities are discussed in the
following sections.
8.1.1 Algorithms
The purpose of the quality assurance program for model algorithm development
is to reduce or eliminate "model error" caused by the use of algorithms that
are inappropriate, have assumptions which are frequently violated, have
extensive data requirements, or produce biased results. A second purpose is
to ensure that algorithms which are appropriate are implemented properly in
the computer. The final purpose is to ensure that the code developed is
capable of simulating results that are physically realistic and meaningful.
The quality assurance plan for algorithm development has three facets:
• Selection
• Implementation
• Testing
Algorithm selection was initially done by the individual or organization
responsible for its development. Normally, the selection process was
initiated by a literature review. The key criterion for initial selection
was that the algorithms be physically based with parameters which could be
easily estimated from readily available information. In addition, it was
desirable that the algorithms had been tested against field or laboratory
data and shown to be robust under a wide variety of circumstances. These
initially selected algorithms were then proposed in a memorandum or oral
presentation format to the model development team at one of a series of
technical meetings which were held every 6 to 8 weeks during project
execution. Changes suggested by reviewers were considered before committing
the algorithms to a FORTRAN format.
201
-------
Algorithm implementation was carried out by the same individual(s) who
selected the algorithm. Each algorithm implemented was tested in one of
several ways. First, if mass balance was a consideration, an acceptable
closure for mass had to be obtained. Second, the "realism" of simulation was
checked. This was done, under the most favorable circumstances, by checking
the results of the simulations against published results of other similar
models, analytical solutions, or data with which other models were developed.
If other results were not available, algorithms were checked simply to see if
they produced results which were reasonable based on the judgment of the
developer and the rest of the model development team. Results of
implementation testing were often reperformed or at least reviewed by other
members of the model development team, often on different computer systems.
The final phase of algorithm QA was testing. Once all the individual
pretested modules of code were assembled, the complete package was tested for
mass balance realism. "Realism" was judged by comparing the simulation
results to observed field data. In this case, the observed data were for the
pesticide aldicarb in a contaminated, unconfined aquifer in Long Island, New
York. The final testing was done by Woodward-Clyde Consultants, while model
algorithm and software development were done primarily by Woodward-Clyde,
HydroGeologic, and Aqua Terra Consultants. Thus, a high degree of
objectivity was maintained in the initial code testing phase, included in the
development effort. Continued model testing and evaluation is being
performed under the support of the EPA Center for Exposure Assessment
Modeling (CEAM) in Athens, GA by EPA scientists, contractors, and selected
'beta' test sites/organizations across the country.
8.1.2 Software
Given the large amount of new computer code required for this project,
combined with the fact that several different firms were involved in code
development, there was a great need for a systematic approach to software
development and testing. Our software quality assurance program consisted of
the following activities:
Software Development
• Use of structured code
• Use of standardized coding conventions
Software TestinE
• Review of code capabilities versus project specifications
• Code analysis using Maintainability Analysis Tool (MAT)
• Testing as many code options as possible
• Compilation and testing of code on different machines
8.1.2.1 Quality Assurance in Software Development
In developing the RUSTIC program, the use of structure in design,
implementation, and documentation was emphasized. The concepts of
"structured programming technology" (IBM 1974) have been used extensively.
It has been found that the use of these methods leads to computer programs
202
-------
and documentation which are easier to understand than those produced by
conventional methods. Also, the software has a unified underlying structure
and is easier to maintain.
Coding conventions are particularly useful when different groups are jointly
developing software and when groups that support the software are different
from the groups that develop the software. With the use of conventions, the
software can more easily be read, understood, debugged, and maintained. We
have used the USGS Office of Surface Water's FORTRAN coding conventions
(Kittle et al. 1986) as the basis for code development related to this
project.
8.1.2.2 Quality Assurance in Software Testing
Considerable effort was expended to ensure the functional capabilities of the
RUSTIC code prior to its release. Throughout the project, the code
capabilities were reviewed to ensure that the evolving code was capable of
performing according to the project goals and requirements. The quality of
new subroutines and progress was assessed by applying the Maintainability
Analysis Tool (Berns 1987) to completed code. The use of this FORTRAN static
analyzer allows the checking of new or modified code after it has been
successfully compiled. It supplements the compiler by analyzing the program
as a whole rather than as individual source modules. By having this
viewpoint, it identifies problems with common blocks and argument lists.
Additionally, it notes bad programming practices, possible logic errors
(i.e., variables referenced but not set), and deviations from standard
FORTRAN.
After improving the code based on the suggestions of the static analyzer, a
series of tests were performed which exercised as many program options as
possible. While this testing resulted in the correction of many problems, it
is anticipated that initial application studies using the new RUSTIC program
will uncover additional oversights. As with any computer program, the user
should carefully review input and results. Output review by the user should
always include analysis of mass balance summaries reported by the RUSTIC
process modules.
A final concern in software development is code portability. The RUSTIC code
has been compiled and run on the IBM personal computer and the Prime computer
running PRIMOS. Implementation of the current version on the IBM-PC requires
greater than 640 K RAM. Adaptation to other machines should require minimal
changes due to the coding conventions which were enforced during program
development.
8.1.3 REFERENCES
Berns, G.M. 1987. "Reliable FORTRAN with MAT," Computer Language, Vol. 4,
No. 6.
International Business Machines Inc. 1974. Structured Programming Textbook
& Workbook -- Independent Study Program,
Kittle, J.L., Jr., A.L. Lumb, and T.D. Barnwell. 1986. ANNIE/WDMS FORTRAN
Coding Convention and Documentation Software. U.S. Geological Survey
Open File Report XXXX (Draft). Reston, VA.
203
oU.S. GOVERNMENT PRINTING OFFICE: 1 989 -6 "» 8 • 1 6 V 0032*
------- |