United States
Environmental Protection
Agency
Office of Research and
Development
Washington DC 20460
EPA/600/R-01/074
November 2001
US EPA OifKe of Research and Development
A User's Guide to IPX, the
In-Place Pollutant Export
Water Quality Modeling
Framework
Version 2.7.4
-------
-------
EPA/600/R-01/074
November 2001
A USER'S GUIDE TO IPX, THE IN-PLACE POLLUTANT
EXPORT WATER QUALITY MODELING FRAMEWORK
VERSION 2.7.4
Mark Velleux
Stephen Westenbroek
James Ruppel
Wisconsin Department of Natural Resources
Madison, Wisconsin 53701
Michael Settles1
OAO Corporation
Douglas Endicott1
U.S. Environmental Protection Agency
Lakes Research Station
Grosse He, Michigan 48138
October 31,2000
U.S. ENVIRONMENTAL PROTECTION AGENCY
OFFICE OF RESEARCH AND DEVELOPMENT
NATIONAL HEALTH AND ENVIRONMENTAL EFFECTS RESEARCH LABORATORY
MID-CONTINENT ECOLOGY DIVISION
LARGE LAKES RESEARCH STATION
GROSSE ILE, MICHIGAN 48138
Recycled/Recyclable
Printed with vegetable-based ink on
paper that contains a minimum of
50% post-consumer fiber content
processed chlorine free.
-------
DISCLAIMER
Development of the IPX framework was funded wholly or in part by the United States
Environmental Protection Agency (USEPA). It has been subjected to the Agency's peer and
administrative review, and has been approved for publication as an USEPA document. Mention
of trade names, commercial products, or services, does not convey, and should not be interpreted
as conveying, official USEPA approval, endorsement, or recommendation.
11
-------
FOREWORD
The need to consider the environmental benefits and consequences of various
management options continues to intensify as (1) environmental problems become more
complex, (2) the means necessary to solve problems become more technical, time-consuming,
and expensive, and (3) the costs to implement remedial strategies increase. Decisions which
include remediation must be weighed carefully against the costs and effectiveness of the
solutions. One environmental problem receiving considerable attention is that of contaminated
sediments. The primary question is whether contaminated sediments should be remediated or
left in place. Contaminated sediments are typically a legacy issue; however, continuing loads
from other point and non-point sources confound clear solutions. Given the reality of limited
funds for clean-up, environmental managers need to define what degree of remediation would
achieve desired results and which sites are the priorities. In this regard, managers must justify
expenditures and need to understand the environmental improvements expected to result from
any remedial action on both local and lake-wide scales.
Mathematical models of environmental processes and effects offer an approach to help
answer these questions and guide decisions on management options. These models should
incorporate appropriate physical, chemical, and biological processes, as well as inputs of solids
and chemicals. Model results must be confirmed using field-collected data for reliability and
scientific credibility. Once the models are confirmed, the predictive capabilities for forecasting
the degree of remediation, the environmental benefits, and the time necessary to realize those
benefits can be used in the decision-making process. The IPX framework is one such tool and
was specifically designed to address the transport and fate of in-place sediment pollutants.
IPX was initially developed during a comprehensive analysis of the Lower Fox River as
part of the Lower Fox River/Green Bay Mass Balance Project and was used to predict future
loadings to Green Bay under alternative remedial actions. Subsequently, IPX was successfully
applied to the Buffalo and Oswego Rivers in New York. This document describes revisions and
updates to IPX which have occurred during the Lake Michigan Mass Balance Study through
application to rivers in the State of Wisconsin. Sufficient detail has been provided to allow other
investigators to apply the IPX framework to other sites.
iii
-------
Development and application of complex, dynamic models should be approached with
care, and requires the cooperation of a team of experts for any modeling endeavor. As is the case
with the IPX framework, cross-agency effort and coordination has been required for many years;
involving aspects such as sample collection, analyses, and database development; and including
federal, state, and private sector efforts. It is recommended that users be familiar with basic
water quality modeling techniques and associated theory regarding the transport and fate of
chemicals in the aquatic environment. Finally, model development should be conducted in
concert with monitoring and experimental research to provide data, estimates of model process
rates, and in turn, reduce model uncertainties.
This report represents a collaborative effort among scientists and engineers at the
USEPA, ORD, Mid-Continent Ecology Division, Large Lakes and Rivers Forecasting Research
Branch and the Wisconsin Department of Natural Resources (WDNR). It is our desire that these
efforts contribute to a more complete understanding of our environment and lead the way to
practical improvements that enhance environmental quality.
IV
-------
AUTHORSHIP AND ACKNOWLEDGMENTS
Version 2.7.4 of the IPX Water Quality Modeling Framework was the result of a
collaborative effort between the USEPA Large Lakes Research Station (LLRS) and the WDNR.
The principal contributors were: Mark Velleux (WDNR), Stephen Westenbroek (WDNR), James
Ruppel (WDNR), Michael Settles (OAO Corporation), and Douglas Endicott1 (USEPA).
Additional contributions in support of source code debugging were made by Michael Riley (S. S.
Papadopulos and Associates). The version number of this release (Version 2.7.4) corresponds to
the revision level assigned to the source code with the Revision Control System (RCS) for file
management.
IPX Version 2.7.4 User's Guide
The IPX Version 2.7.4 user's guide was prepared by M. Velleux. This document was
directly derived from the IPX "(Version 1.0) user's guide (Velleux et al. 1994b). Portions of the
document were adapted in whole or in part from the WASP4 User's Guide (Ambrose et al.
1988), the WASPS User's Guide (Ambrose et al. 1993), and other documentation prepared
during prior stages of IPX framework development. While it was not always possible to
determine complete authorship for each element of the documentation, every effort has been
made to appropriately credit the source authors for their contributions.
Chapter 1. IPX Model Theory: this chapter was adapted, in part, from the
WASP4/WASP5 User's Guide. The source authors of the text from the user's guides were
Robert B. Ambrose (USEPA), Tim Wool (Computer Sciences Corporation), John P. Connolly
(Manhattan College), Robert W. Schanz (Woodward-Clyde Consultants), and James Martin
(AScI Corporation). Text describing the settling and resuspension processes, photolysis, and
volatilization was developed by M. Velleux, J. Gailani, and D. Endicott. Additional text
describing the probability of deposition was summarized from QEA (1999). Text describing the
semi-Lagrangian frame of reference for the sediment bed was developed by M. Velleux.
Chapter 2. IPX Input Data File Structure: this chapter was also adapted from the
WASP4/WASP5 User's Guide. The source authors for much of this text were R. Ambrose, T.
1 Present affiliation: Great Lakes Environmental Center, 739 Hastings Road, Traverse City, Michigan 49686.
V
-------
Wool, J. Connolly, and R. Schanz, and J. Martin. Additional source authors were Kirk Freeman
(Computer Sciences Corporation), D. Endicott, M. Velleux, and J. Gailahi.
Chapter 3. IPX Output and Post-Processing. The source authors of this text were M.
Velleux and K. Freeman.
Chapter 4. IPX Input Data File Pre-Processing. The source authors of this text were J.
Gailani and M. Velleux.
Chapter 5. IPX Programmer's Guide. Significant portions of this chapter were adapted
from the WASP4/WASP5 User's Guide with additional text and editorial modifications by M.
Velleux. Text describing UNIX environment development tools for IPX was prepared by M.
Settles.
The foreword to this report was prepared by Russell Kreis (USEPA) based on text written
by William Richardson (USEPA). The preface was adapted from the WASP4/WASP5 User's
Guide with text and editorial modifications by M. Velleux.
IPX Framework
In addition to the documentation effort, many people also contributed to source code
development and programming efforts for the IPX framework. IPX Version 2.7.4 is an
outgrowth of WDNR and USEPA efforts in support of the USEPA Lake Michigan Mass Balance
Study. This version was developed from the first release of the IPX framework (Velleux et al.
1994b) and was derived from the WASP family of water quality modeling frameworks. Source
code management using RCS was implemented by M-.- Settles at the LLRS. Additional
programming support was also provided by M. Settles and Xiangsheng Xia (Computer Sciences
Corporation). Initial development of IPX Version 2.7.4 was performed by M. Velleux with the
assistance of S. Westenbroek and J. Ruppel and with contributions of D. Endicott. M. Riley
provided source code debugging support and contributed the concept for ghost element collapse.
Final programming and code development for Version 2.7.4 was performed by M. Velleux.
As noted above, IPX was derived from the WASP4 water quality modeling framework,
which was itself derived from earlier frameworks, such as WASP and TOXIWASP developed
with support from USEPA. During the earliest stages of the IPX development effort, WASP4
was significantly revised by K. Freeman and D. Endicott. Subsequent programming was
performed by M. Velleux, K. Freeman, and Frank Mitchell (Computer Sciences Corporation).
vi
-------
Final programming and code development was performed by Keqin Shen (Computer Sciences
Corporation), M. Velleux, J. Gailani, and D. Endicott.
Post-Processing Program
The W4DIS post-processing program was written by K. Freeman. Programming for
Version 2.7.4 (W4DIS274) was performed by M. Velleux.
Pre-Processing Programs
The REDUCE pre-processing program was written by F. Mitchell with modifications by
M. Velleux. The SETTLE and RESUSPND pre-processing programs were written by Anne
Dame (USEPA), M, Velleux, and J. Gailani. Additional programming was performed by S.
Westenbroek. The TIMESTEP pre-processing program was written by A. Dame and M.
Velleux.
Other Contributions and Acknowledgments
Beyond those contributions already noted, many other individuals also contributed to the
early development of IPX from WASP4. In particular, J. Martin and T. Wool (AScI
Corporation) were often consulted during the initial stages of WASP4 use and debugging efforts
at LLRS. W. Richardson, during his tenure as Chief of the LLRS, was very supportive of all IPX
framework development efforts at each stage of the process. Debra Caudill (Computer Sciences
Corporation) provided word processing support. All figures in this manual were revised and
prepared by Kay Morrison (OAO Corporation). Finally, the support of Russell Kreis (USEPA),
present Chief of the LLRS, is gratefully acknowledged.
vn
-------
PREFACE
Numerous frameworks are available to simulate the transport and fate of toxic chemicals
in surface waters. Water quality modeling frameworks can range in complexity from simple,
steady-state analytical solutions that describe a single first-order decay process to detailed, time-
variable numerical solutions that describe the physical, chemical, and biological transfer and
transformation processes affecting a chemical. As an outgrowth and extension of previous water
quality modeling frameworks such as WASP (Di Toro et al. 1981), WASP4 (Ambrose et al.
1988), TOXIWASP (Ambrose et al. 1983), and WASTOX (Connolly and Winfield, 1984), the
IPX framework was specifically designed to address the transport and fate of in-place pollutants
in rivers and streams while retaining the flexibility needed to analyze a variety of water quality
problems.
This manual describes IPX and is divided into five main sections. The first section, IPX
Model Theory, documents the equations and assumptions underlying the IPX framework and
components. Some guidance on using the framework is offered, along with typical input data
values, when appropriate. More general summaries of the equations and possible input data
values are provided in the "Rates Manual" (Bowie et al. 1985), the "Screening Manual" (Mills et
al. 1985), and the 'Toxicant Rates Manual" (Schnoor et al. 1987) as well as other texts cited in
the references.
The second section, IPX Input Data File Structure, documents the input data
specifications necessary to run IPX. Each data group is described, with input names, formats,
and definitions. Convenient tabular summaries" of environmental parameters, chemicals
constants, and kinetic time functions are also provided.
The third section, IPX Simulation Output and Post-Processing, documents the output files
created by IPX during a simulation as well as the W4DIS274 post-processing program. Short
descriptions of each file are presented. A description of the W4DIS274 post-processor and
examples of its use are also presented.
vin
-------
The fourth section, IPX Input Data File Pre-Processing. documents the pre-processing
programs that accompany IPX. The pre-processors include REDUCE, SETTLE, RESUSPND,
and TIMESTEP. Short overviews detailing the use of these programs and the structure of their
input data are provided.
The fifth section, IPX Programmer's Guide, documents the general structure of the IPX
program source code. First, a general warning is presented suggesting that users verify code
operation. Second is a short section that describes porting IPX to various computational
platforms. Finally, descriptions of the major subroutines and their calling sequence, common
blocks (including how to re-dimension IPX for custom applications), and input/output units are
presented.
IX
-------
TABLE OF CONTENTS
DISCLAIMER ii
FOREWORD iii
AUTHORSHIP AND ACKNOWLEDGEMENTS v
PREFACE viii
TABLE OF CONTENTS x
FIGURES .' xiv
TABLES xv
CHAPTER 1: IPX MODEL THEORY : 1
1.1 IPX FRAMEWORK OVERVIEW 1
1.2 THE GENERAL MASS BALANCE EQUATION 2
1.3 THE MODEL NETWORK/MODEL SEGMENTATION 4
1.4 Advective, Dispersive, and Evaporative Transport 9
1.4.1 Water Column Advection 9
1.4.2 Water Column Dispersion 13
1.4.3 Pore Water Advection 13
1.4.4 Pore Water Diffusion 14
1.4.5 Evaporation and Precipitation 15
1.5 Sediment Transport and the sediment bed 16
1.5.1 Settling 17
1.5.2 Resuspension, Sediment Aging, and the Surficial
Sediments 23
1.5.3 Burial and Scour (Unburial) 26
1.5.3.1 Eulerian Option '. 27
1.5.3.2 Semi-Lagrangian Option 28
1.5.4 The Sediment Bed 29
1.5.4.1 Surficial Sediments 29
1.5.4.2 Subsurface Sediments 30
1.5.5 Screening-Level IPX Sediment Transport Simulation
Guidelines , 33
1.6 Simulating the Transport and Fate of Toxic Chemicals using IPX 34
-------
1.6.1 Equilibrium Sorption 4 .-. 38
1.6.2 Kinetic Transformations 42
1.6.3 Hydrolysis , 44
1.6.4 Photolysis 47
1.6.5 Oxidation 51
1.6.6 Bacterial Degradation 52
1.6.7 Volatilization 56
1.6.8 Extra Reaction 65
1.7 Simulating the Transport and Fate of Heavy Metals using IPX 66
1.8 Simulation Complexity Levels 67
1.8.1 Simulating Solids 69
1.8.2 Simulating Equilibrium Partitioning '..: 69
1.8.3 Simulating Kinetic Reactions 71
1.9 Summary of Data Requirements 74
1.10 Summary of Model Equations 77
1.L1 Other Model Inputs 80
1.11.1 Model Identification and Control Parameters , 80
1.11.2 Boundary Parameters 82
1.11.3 Point Source and Non-Point Source Loads 83
1.11.4 Initial Conditions 84
1.12 Applying the Model 84
CHAPTER 2: IPX INPUT DATA FILE STRUCTURE 87
2.1 IPX Input Data File Structure Overview '...... 87
2.2 IPX Input File Data Group Descriptions 89
2.2.1 Data Group A: Model Identification and Simulation
Control 89
2.2.2 Data Group B: Exchange (Dispersion) Coefficients 94
2.2.3 Data Group C: Volumes 96
2.2.4 Data Group D: Water Column Flows; Pore Water Flows;
Settling Velocities; Resuspension Velocities; and
Precipitation/Evaporation 98
2.2.5. Data Group E: Boundary Concentrations 106
2.2.6 Data Group F: Point and Non-Point Source Loads (Forcing
Functions) 107
2.2.7' Data Group G: Parameters.... 110
2.2.8 Data Group H: Physicochemical Constants Ill
xi
-------
2.2.9 Data Group I: Kinetic Time Functions
2.2.10 Data Group J: Initial Concentrations
2.2.11 Data Group K: Surficial Sediment Age Layer Conditions
2.2.12 Data Group L: Initial Conditions for Deep Sediment Layers
(Ghost Stack) (Seml-Lagrangian Bed Option).
2.3 IPX MODEL PARAMETERS, CONSTANTS, AND TIME-
FUNCTIONS OVERVIEW
2.3.1 Environmental Parameters: Data Group G
2.3.2 Physicochemical Constants: Data Group H
2.3.3 Kinetic Time Functions: Data Group I
CHAPTER 3: IPX SIMULATION OUTPUT AND POST-PROCESSING
3.1 IPX OUTPUT OVERVIEW
3.2 The *.DMP File
3.3 The *.DMA File
3.4 The *.EXP File
3.5 The *.MSB File
3.6 The*.OUTFile
3.7 The *.RST File
3.8 W4DIS274: the *.DMP and *.DMA File Post-Processor ,
CHAPTER 4: IPX INPUT DATA FILE PRE-PROCESSING
4.1 IPX Input Data File Pre-Processing Program Overview....
4.2 REDUCE: The Simulation Hydrograph Pre-Processor
4.3 SETTLE: The Settling Velocity Pre-Processor
4.4 RESUSPND: The Resuspension Velocity Pre-Processor ,
4.5 TIMESTEP: The Simulation Time-Step Pre-Processor ,
CHAPTER 5: IPX PROGRAMMER'S GUIDE
5.1 IPX PROGRAMMING OVERVIEW
5.2 AVAILABILITY OF IPX PROGRAMS AND Porting to Other
Platforms
5.3 UNIX Environment Development Tools for IPX
5.4 IPX PROGRAM DESCRIPTION
5.4.1 Main Program :
5.4.2 Input Subroutines
5.4.3 Output Subroutines
5.4.4 Simulation Control and Process Subroutines to Compute
5.4.5 Utility Subroutines
113
114
116
118
119
120
122
133
136
136
136
136
137
137
137
137
138
142
142
142
144
145
146
148
148
148
149
150
151
151
153
154
159
xn
-------
5.4.6 Toxic Chemical Kinetic Subroutines 160
5.5 Input and Output Units 164
5.6 COMMON BLOCKS ....„ 165
REFERENCES 167
APPENDIX A: IMPLEMENTATION OF THE MASS BALANCE EQUATION 175
Xlll
-------
FIGURES
NUMBER PAGE
1.1. Coordinate system for mass balance equation 3
1.2. The model network: Model segmentation and segment types 4
1.3. Spatial scales used for analysis of Lake Ontario 6
1.4. Example observed and calculated frequency distributions. , 7
1.5. Mass transport, transfer and transformation processes in IPX 35
1.6. Hydrolysis reactions 45
1.7. pH dependence of example hydrolysis reactions • 46
1.8. Photolysis reactions •• 48
1.9. Example microbial transformations (from Alexander, 1980) 53
1.10. Volatilization : 57
5.1. Hierarchy of major subprogram units in IPX 151
5.2. Hierarchy of major kinetic process subroutines in IPX 161
Al. Definition sketch for finite difference equation .- 176
xiv
-------
TABLES
NUMBER
PAGE
1.1. Comparison of hydraulic exponents 12
1.2. Stokes' settling velocities for a range of particle sizes and densities
at 20°C 18
1.3. Generalized grain size distribution functions for regulated tributaries 22
1.4. Generalized grain size distribution functions for unregulated tributaries 22
1.5. Concentration related S3onbols used in mathematical equations. 36
1.6. IPX partitioning/binding data 42
1.7. IPX general kinetic data 44
1.8. IPX Hydrolysis Data 47
1.9. Average light intensity Lk of wavelength k at 40°N, mE/cm2/day 50
1.10. IPX photolysis data. 50
1.11. IPX oxidation data 52
1.12. IPX Biodegradation Data 54
1.13. Size of Typical Bacterial Populations in Natural Waters 55
1.14. IPX volatilization data. 65
1.15. Speciation of priority metals between dissolved and particulate phases as
a function of suspended solids concentrations in streams 68
1.16. Environmental properties affecting interphase transport and
transformation processes 75
1.17. Chemical properties affecting interphase transport and transformation
processes 76
1.18. Time variable environmental forcing functions (combinations of IPX
time functions and constants/parameters) 77
1.19. Numerical Dispersion (rn^/s) values resulting from u and At for a
Segment Length of L = 2000m 83
2.1. Conceptual levels of simulation complexity 120
2.2. Data Group G environmental parameters 121
2.3. Partition coefficients and rate constants for simple reactions 123
xv
-------
NUMBER PAGE
2.4. General chemical constants 124
2.5. Sorption constants 124
2.6. Hydrolysis constants , 125
2.7. Photolysis constants 126
2.8. Global constants (generally not used at this time) 126
2.9. Oxidation constants 127
2.10. Biodegradation constants 127
2.11. Volatilization constants 128
2.12. "Extra reaction" constants 129
2.13. Data Group I kinetic time functions 134
3.1. W4DIS274 display variables 139
3.2. W4DIS274 sediment archive stack display variables 141
5.1. Structure of the non-point source load file NPS.DAT. 153
xvi
-------
CHAPTER 1
IPX MODEL THEORY
1.1
IPX FRAMEWORK OVERVIEW
IPX, the In-place Pollutant eXport model, is a dynamic compartment (control volume)
water quality modeling framework specifically designed to simulate the transport and fate of
sediments and associated hydrophobic contaminants in tributaries. IPX includes process
descriptions for sediment aging, decreased sediment resuspendability with increasing age, and
resuspension of freshly deposited sediments as a function of water velocity (shear stress at the
sediment-water interface). These processes are needed to realistically simulate contaminant
transport and substantially improve the model framework for application to tributary systems
subject to significant deposition and resuspension events (where the shear stress is a function of
water velocity). IPX has been applied, in a number of case studies, for estimating contaminant
export from tributaries with contaminated sediments to receiving waterbodies (Velleux and
Endicott, 1994a; Velleux et al. 1995; Velleux et al. 1996).
IPX evolved from a project to simulate the transport of polychlorinated biphenyls (PCBs)
in the Lower Fox River, Wisconsin, as part of the EPA Green Bay Mass Balance Study (Velleux
and Endicott, 1994a). For that project, the WASP4 water quality modeling framework (Ambrose
et al. 1988) was used. In the course of applying WASP4 to model transport in the Fox River,
limitations were encountered which led to its modification and, eventually, its evolution into
IPX. In addition to sediment transport enhancements, IPX expands the contaminant simulation
capabilities of WASP4 to allow simultaneous simulation of an unlimited number of contaminants
(constrained only by available computer memory) while retaining most of the computation speed
and flexibility of the original WASP4 framework. IPX may be considered a hybrid of WASP4,
but because of significant changes in its capabilities and function, IPX is documented as an
independent modeling framework. The publication by Velleux et al. (1996) provides further
detail regarding the distinctions between IPX and WASP4.
-------
The equations solved by IPX are based on the key principle of the conservation of mass.
This principle requires that the mass of each water quality constituent (chemical or solid)
simulated must be accounted for in one way or another. IPX traces each water quality constituent
from each point of spatial and temporal input to all final points of export, thereby accounting for
(conserving) mass in space and time. To perform these mass balance computations, the user
must supply the framework with input data defining seven important characteristics:
• simulation and output control
• model segmentation
• advective, dispersive, and sediment transport
• boundary concentrations
• point and diffuse (non-point) source waste loads
• kinetic parameters, constants, and time functions
• initial concentrations
The input data define transport and transformation process parameters in the generalized
mass balance equations that uniquely define the spatial and temporal domain of the water quality
problem to be solved. These equations are numerically integrated using Euler's method. At
user-specified print intervals, the values of all display variables are saved for subsequent retrieval
by the W4DE3274 post-processor program. W4DIS274 allows the user to produce data tables of
model output that can be subsequently graphed in any number of software packages. The output
options available are discussed in Chapter 3. ,
1.2 THE GENERAL MASS BALANCE EQUATION
A mass balance equation for a water body accounts for all material entering and leaving
the system by direct and diffuse loading, advective and dispersive transport, settling and
resuspension, and physical, chemical, and biological transformations. Consider the coordinate
system shown in Figure 1.1, where the x- and y-coordinates are in the horizontal plane, and the z-
coordinate is in the vertical plane. The mass balance (partial differential) equation for an
infinitesimally small fluid volume is:
af p 3C
3z Z9z
(1.1)
+ SL+SB+SK
-------
= Coordinate Direction of flow
X direction >
Figure 1.1. Coordinate system for mass balance equation.
where: C =
t =
UX,Uy,UZ =
SL
SB
SK
concentration of the water quality constituent, mg/L (g/m3) [M/L3]
time, days [T]
longitudinal, lateral, and vertical advective velocities, m/day [L/T]
longitudinal, lateral, and vertical diffusion (dispersion) coefficients,
m2/day [L2/T]
direct and diffuse loading rate, g/m3-day [M/L3/t]
boundary loading rate (including upstream, downstream, sediment, and
atmospheric), g/m3/day [M/L3/T]
total kinetic transformation rate; positive indicates a source, negative a
sink, g/m3/day [M/L3/T]
By expanding the infinitesimally small control volumes into larger adjoining "segments"
and specifying transport, loading, and transformation parameters, a finite-difference form of
Equation 1.1 is derived. The finite-difference form of the general mass balance equation is
implemented in IPX. For brevity and clarity, however, the derivation of the finite-difference
form of the mass balance equation is for a one-dimensional reach. Assuming vertical and lateral
homogeneity, Equation 1.1 can be integrated over y and z to obtain:
_
dt
(1.2)
-------
X
= Coordinate
direction
Segment Types:
Surface water I i
Subsurface water !•'<•'. -•!
Surficial sediment ^^m
Subsurface sediment
Figure 1.2. The model network: Model segmentation and segment types.
where: A = cross-sectional area, m2 [L2]
This equation represents the three major classes of water quality processes: transport
(term 1), loading (term 2), and transformation (term 3). The finite-difference form of this
equation is derived in Appendix A. The model network and the major processes are discussed in
the following sections.
1.3
THE MODEL NETWORK/MODEL SEGMENTATION
The model network is a set of expanded control volumes or "segments" that together
represent the physical configuration of the water body. As Figure 1.2. illustrates, the network
may subdivide the water body laterally and vertically as well as longitudinally. The model
-------
network can include sediment segments as well as water column segments. Concentrations of
water quality constituents are calculated by solving the mass balance equation within each
segment. Mass transport rates of each constituent are then calculated across the interface of
adjoining segments.
Segments in IPX may be one of six types, as specified by the input variable 1TYPE. A
value of 1 indicates an epilimnion layer (surface water), 2 indicates a hypolimnion layer
(subsurface water), 3 indicates a Eulerian surficial sediment layer, 4 indicates a Eulerian
subsurface sediment layer; 5 indicates a semi-Lagrangian surface sediment layer, and 6 indicates
a semi-Lagrangian subsurface sediment layer. The segment type is used to define aspects of
sedimentation and transformation processes. The user must also specify the vertical organization
of the model network. The segments immediately above and below each segment are specified
by the input variables ITOPSG and IBOTSG. This alignment is important for determining mass
transport due to burial and scour.
Segments are explicitly specified active compartments of the model network. During
numerical integration of the mass balance equations, the model framework computes mass
transport, transfer, and transformation terms for all active compartments. In addition to explicitly
specified compartments, the semi-Lagrangian option for the sediment bed (Section 1.4.5.3)
allows the user to also specify inactive compartments of the model network for the sediment bed.
The network of inactive model segments is referred to as the "ghost" stack. Inactive
compartments can enter or leave the active model network (thereby entering or leaving the ghost
stack) in response to interactions within the sediment bed (decreases or increases of sediment bed
elevation). With the semi-Lagrangian sediment bed option, it is also possible for the model to
generate additional inactive compartments (beyond any initially specified) in response to
interactions that increase sediment bed elevation. During numerical integration, the model
framework does not perform mass balance computations on ghost stack elements.
In tributary models, segment volumes, river flows, and the simulation time-step are
directly related. Segment volume and flow determine the hydraulic residence time (HRT). As
HRT decreases, the simulation time-step must also decrease to insure stability and numerical
accuracy. Segment size can vary dramatically, as illustrated in Figure 1.3. For any simulation,
segment sizes may be dictated as much by the spatial and temporal scale of the simulation as by
the physical characteristics of the system. For example, analyzing a simulation involving the
upstream migration of a pollutant into a water supply near a river mouth due to flow reversals
might require a time-step on the order of minutes. By contrast, analyzing a problem involving
. ' • " 5
-------
Model
Designation
• direction of
water flow
Number of
Segments
Horizontal
Scale (km2)
Epilimnion
Segments
Lakel
13,000
LakeS
67
200-1000
Rochester
Embayment
72
10-100
Figure 1.3. Spatial scales used for analysis of Lake Ontario.
-------
TO
as
a
i
Observed
Time Scale 2
Time Scale 1
Steady - State
5 50 95
Cumulative Probability (%)
Figure 1.4. Example observed and calculated frequency distributions.
the total residence time of that pollutant in the same water body could allow a time-step of hours.
As an example from a lake analysis consider the networks shown in Figure 1.3. The first
network was used to study the general eutrophic status of Lake Ontario. The second network was
used to investigate the lake-wide spatial and seasonal variations in eutrophication. The third
network was used to predict changes in near-shore eutrophication of Rochester Embayment
resulting from specific pollution control plans.
As part of the problem definition, the user must determine the accuracy to which the
temporal and spatial variability of a contaminant's frequency distribution must be predicted. For
example, a daily-average dissolved oxygen concentration of 5 mg/L would not sufficiently
protect fish if fluctuations result in concentrations less than 2 mg/L for 10% of the time.
However, predicting extreme concentration values is generally more difficult than predicting
average values. Figure 1.4 illustrates typical frequency distributions predicted by three model
7
-------
time scales and a typical distribution observed by rather thorough sampling as they would be
plotted on probability paper. The straight lines imply normal distributions. While reducing
model segment sizes (and consequently the model time-step) allows more accurate simulation of
the frequency distribution, increased model accuracy requires increased temporal and spatial
resolution of the model inputs such as loads and flow patterns.
When determining the temporal resolution of a simulation, care must be taken to preserve
important non-linear interactions. When two or more important processes have different periods
of variation, the temporal resolution of the simulation will generally be controlled by the time
scale of the process with the greatest temporal variability. Temporal resolution may also be
affected by the physical characteristics of the model system. For example, discontinuous batch
discharges into a large lake might safely be averaged over a day or week, because large scale
transport variations are relatively infrequent. However, it may not be possible to average the
same batch discharge into a tidal estuary or river because of the semi-diurnal or diurnal tidal
variations or seiche action. A third example is salinity intrusion in estuaries. Tidal variations in
flow, volume, and dispersion can interact such that even long-term simulations may require
model inputs and simulation time-steps on the order of hours. In general, the simulation time-
step must be somewhat less than the temporal variability of the dominant mass transport/transfer
pathways. For advectively dominated tributaries (where advective transport is the greatest mass
transport pathway), a rough guide for determining the simulation time is to limit the time-step to
one third the hydraulic residence time of the smallest water column segment in the model
network.
Once the temporal variability has been determined, then the spatial variability of the
water body must be considered. Generally, spatial characteristics must be homogeneous within a
segment. In some cases, this restriction can be relaxed by judicious averaging over width, depth,
and/or length. For example, depth governs the impact of reaeration and sediment oxygen
demand in a column of water. Nevertheless, averaging the depth across a river is generally
acceptable for conventional waste load allocation, whereas averaging the depth across a lake
would not generally be acceptable. Other important spatial characteristics to consider (depending
upon the problem being analyzed) include temperature, light penetration, velocity, pH, and the
physical characteristics of the sediments, as well as sediment contaminant concentrations.
The expected spatial variability of the water quality concentrations also affects segment
sizing. The user must determine to what extent averaging concentration gradients is acceptable.
For example, water quality conditions change rapidly near a loading point but may stabilize
8
-------
downstream. Studying the effects of the discharge on a beach one quarter mile downstream of a
discharge point may require smaller segments than when studying the effects of the discharge on
a beach several miles downstream.
A final, general guideline for obtaining accurate simulations is that water column segment
volumes should be of roughly similar hydraulic retention times to minimize numerical
dispersion. If river flows increase significantly downstream in the model network, then segment
volumes should also increase proportionately. The user should define segment volumes for
critical reaches of the system and then determine appropriate sizes for upstream and downstream
segments accordingly.
The segment volumes specified must also be sized to best represent the actual spatial
variability, as discussed above. Nonetheless, this guideline will allow use of larger time-steps
and result in greater numerical accuracy over the entire model network, as explained in Section
1.11.
1.4 ADVECTIVE, DISPERSIVE, AND EVAPORATIVE TRANSPORT
Transport processes in IPX are divided into seven distinct types, or fields. Each field is
associated with a data group of the model input file. Transport fields 1, 2, and 7 are used to
specify advective, dispersive, and evaporative transport. Transport fields 3, 4, 5, and 6 are used
to specify sediment transport. This section describes advective, dispersive, and evaporative
transport and the operation of transport fields 1, 2, and 7. Sediment transport and the operation
of transport fields 3,4, 5, and 6 are described in Section 1.5.
i
The first transport field type moves all phases of each state variable (dissolved, bound,
and paniculate) and is used to specify advective flow and dispersive mixing in the water column.
The second transport field moves only the dissolved and bound phases of each state variable and
is used to specify the movement of pore water in the sediment bed. The seventh field moves
water across the air-water interface (and can also be used to move state variables into but not out
of the model) and is used to specify evaporation and precipitation.
1.4.1 Water Column Advectioii
Advective water column flows directly control the transport of dissolved and particulate
pollutants in many water bodies. Advective flow moves water quality constituents downstream
with the water and accounts for instream dilution. In addition, changes in velocity and depth
which result from Variable flows can affect such kinetic processes as reaeration, volatilization,
-------
and photolysis. An important early step when developing a simulation is to properly describe
water column advection. In IPX, water column flow is input via transport field 1 in Data Group
D. As a user option, water column segment volumes are held constant or are allowed to vary in
response to changes in flow as specified in Data Group C.
IPX tracks each separate inflow specified by the user from its point of origin and through
each segment until it exits the model network. For each inflow, the user must supply a continuity
or unit flow response function and a time function. The continuity function describes the spatial
extent of the inflow response as it varies throughout the model network. The time function
describes the temporal variability of the inflow. The actual flow between segments that results
from a given inflow is the product of the time function and the continuity function. If several
inflow functions are specified between any segment pair, then the total flow between segments is
computed as the sum of the individual flow functions. In this manner, the effect of several
tributaries joining, density currents, and wind-induced flow patterns can be described in a simple
manner. For flow option 1 (IQOPT = 1), IPX sums all the flows within a segment to determine
the direction of net flow and moves mass only in the direction of the net flow (one direction
only). For flow option 2 (IQOPT = 2) IPX moves mass independently of net flow (more than one
direction) to give the effect of large circulation patterns.
A good description of segment geometry as a function of flow conditions can be
important when using IPX to simulate rivers. For flow option 1, a set of user-specified hydraulic
discharge coefficients from Data Group C defines the relationship between velocity, depth, and
stream flow. This method, described below, follows the implementation in QUAL2E (Brown
and Bamwell, 1987).
Discharge coefficients giving depth and velocity from stream flow are based on empirical
observations of the velocity-depth-stream flow relationship (Leopold and Maddox, 1953). It is
important to note that these coefficients are only important when calculating reaeration. The
velocity calculations are not used in time travel. The equations relate velocity, channel width,
and depth to streamflow through power functions:
v = aQh (1-3)
D = cQd (1-4)
B = eQ/ (1-5)
10
-------
where: v =
D
B
Q
a,b,c,d,e,f =
average water velocity, m/s [L/T]
average depth, m [L]
average width, m [L]
flow, m3/s [L3/T]
empirical coefficients
Given that area is a function of channel width (B) and depth (D):
A = DB
From continuity, it is clear that:
Q = vA
= vDB
(1.6)
= (ace)Q(b+d+f)
Therefore, the following relationships hold:
ace = 1
b+d+f=l
(1.7)
(1.8)
IPX only requires that the relationships for velocity (Equation 1.3) and depth (Equation
1.4) be specified; coefficients for channel width (Equation 1.5) are implicitly specified by
Equations 1.7 and 1.8.
These options can be put into perspective by noting that, for a given specific channel
cross-section, the coefficients (a, c, e) and exponents (b, d, f) can be derived from the Manning
equation. For example, if a channel of rectangular cross-section is assumed, then channel width
(B) is not a function of streamflow (Q): the exponent (f) is zero (0.0) and the coefficient (e) is the
width of the rectangular channel (B). By noting that hydraulic radius (R) is approximately equal
to depth (D) for wide streams and that A = D B, the discharge coefficients for rectangular cross
sections can be shown to be 0.4 for velocity and 0.6 for width.
Leopold et al. (1964) noted that stream channels in humid regions tend towards a
rectangular cross-section because cohesive soils promote steep side slopes whereas non-cohesive
11
-------
Table 1.1. Comparison of hydraulic exponents.
Channel Cross-Section
Rectangular
Average of 158 U.S. Gaging
Stations
Average of 10 Gaging
Stations on Rhine River
Ephemeral Streams in
Semiarid U.S.
Exponent for
Velocity (b)
0.40
0.43
0.43
0.34
Exponent for
Depth (d)
0.60
0.45
0.41
0.36
Exponent for
Width (f)
0.00
0.12
0.13
0.29
soils encourage shallow sloped, almost undefined banks. Table 1.1 compares hydraulic
exponents for a rectangular channel with data reported by Leopold et al. (1964). Note that the
average velocity exponent is relatively constant for all channel cross sections. The major
variation occurs as a decrease in the depth exponent and concomitant increase in the width
exponent as channel cross-sections change from the steep side slopes characteristic of cohesive
soils to the shallow slopes of arid regions with non-cohesive soils.
For water bodies such as lakes and reservoirs, velocity and depth may not be functions of
flow. For these cases, the velocity and depth exponents (b and d) can be chosen to be zero (0.0).
Because Q to the zero power equals one (1.0), the coefficients a and c must be the velocity and
depth: if b = 0, then a = V; if d = 0, then c = D. When the depth exponent is zero, IPX adjusts
segment depth as a function of segment volume, assuming rectangular sides.
For site-specific river or stream simulations, hydraulic coefficients and exponents must be
estimated. Brown and Barnwell (1987) recommended estimating the exponents (b and d) and
then calibrating the coefficients (a and c) to observed velocity and depth. The exponents may be
chosen based on observations of channel shape noted in a reconnaissance survey. If cross
sections are largely rectangular with vertical banks, the first set of exponents shown in Table 1.1
should be useful. If channels have steep banks typical of areas with cohesive soils, then the
second set of exponents is appropriate. If the stream is in an arid region with typically non-
cohesive soils and shallow sloping banks, then the last set of exponents is recommended.
A key property of the channel that should be noted in a reconnaissance survey is the
condition of the bank slopes or the extent to which width may expand with increasing
12
-------
streamflow. Clearly the bank slopes and material in contact with the streamflow at the flow-
rate(s) of interest are the main characteristics to note in a reconnaissance.. Table 1.1 gives general
guidance but it should be noted that values are derived for bank-full flows. Even in streams with
vertical banks, the low flows may be in contact with a sand bed having shallow sloped, almost
nonexistent banks more representative of ephemeral streams in semi-arid areas.
1.4.2 Water Column Dispersion
Dispersive water column exchanges significantly influence the transport of dissolved and
particulate pollutants in such water bodies as lakes, reservoirs, and estuaries. Dispersion causes
mixing and dilution between regions of high concentrations and regions of low concentrations.
Even in advectively dominated systems such as rivers, longitudinal dispersion can be the most
important process in diluting peak concentrations that may result from dynamic (unsteady) loads
or spills. Natural or artificial tracers such as dyes, salinity, conductivity, or heat (temperature)
are often used to calibrate dispersion coefficients for a model network.
In IPX, water column dispersion is input via transport field 1 of Data Group B. Any
number of groups of exchanges may be defined by the user. For each group, the user must
supply a time function giving dispersion coefficient values (in m2/s) as they vary in time. For
each exchange in the group, the user must supply an interfacial area, a characteristic mixing
length, and the segments between which the exchange takes place. The characteristic mixing
length is typically the distance between segment midpoints. The interfacial area is the area
normal to the characteristic mixing length shared by the exchanging segments (cross-sectional
area for horizontal exchanges, or surface area for vertical exchanges). The dispersive exchange
between segments i and j at time t is given by:
(1.9)
3r
where: Mj
C- C-
\^j, \_j
Eij(t)
Aij
Lcij
Lcij
mass of chemical in segment i, g [M]
total chemical concentration in segments i and j, mg/L (g/m3) [M/L3]
dispersion coefficient time function for exchange "ij", m2/day [L2/T]
interf acial area shared by segments i and j,m2 [L2]
characteristic mixing length between segments i and j, m [L]
1.4.3 Pore Water Advection
Pore water "flows into or out of the sediment bed can significantly influence sediment
pollutant concentrations. Dissolved or bound water quality constituents are carried through the
13
-------
sediment bed by pore water flow. Depending on the direction of these flows, and the source of
the pollutants, pore water advection may be a source or sink of pollutants to overlying waters.
If sediment segments are included in the model network, the user may specify advective
transport of dissolved chemicals in the pore water. In IPX, pore water flows are input via
transport field 2. Pore water advection transports water and dissolved chemical; sediment and
particulate chemical are not transported. The mass derivative of chemical due to pore water flow
from segment j to segment i is given by:
(1.10)
where: MI = mass of chemical in segment i, g [M]
Cj = total chemical concentration in segment j, mg/L (g/m3) [M/L3]
n; = porosity of segment j, Lw/L (volume of water/volume total solution)
[dimensionless]
fdj = dissolved fraction of chemical in segments i and j [dimensionless]
Qy = pore water flow-rate from j to i, m3/day [L3/T]
The dissolved fraction of a contaminant (fd) may be input by the user in Data Group J or
computed from sorption kinetics.
IPX tracks each separate pore water inflow through the network of sediment segments.
For each inflow (or outflow), the user must supply a continuity function and a time function. The
actual flow through the sediment segments that results from each inflow is the product of the
time function and the continuity function. If a flow originates in or empties into a surface water
segment, then a corresponding surface water flow function must be described in flow field 1 that
matches the pore water function.
1.4.4 Pore Water Diffusion
Diffusive exchanges between pore waters and with the overlying water column can
significantly influence sediment pollutant concentrations, particularly for relatively soluble
chemicals and water bodies with little sediment loading. Dissolved and bound water quality
constituents can be exchanged between the sediment bed and the water column by pore water
diffusion. Depending on the dissolved chemical concentration gradient, pore water diffusion
may be a source or sink of pollutants to overlying waters.
14
-------
If sediment segments are included in the model network, the user may specify diffusive
transport of dissolved chemicals in the pore water. In IPX, pore water diffusion is input via
transport field 2 in Data Group B. Any number of groups of exchanges may be defined by the
user. For each group, the user must supply a time function giving dispersion coefficient values
(in m2/s) as they vary in time. For each exchange in the group, the user must supply an
interfacial area, a characteristic mixing length, and the segments between which exchange takes
place. For pore water diffusion between two sediment segments, the characteristic mixing length
is typically the distance between two sediment segment midpoints; this length is multiplied by a
tortuosity factor, which is internally computed in the framework as the inverse of porosity. For
pore water exchange with a water column segment, the characteristic mixing length is usually
taken to be the depth of the surficial sediment segment. The interfacial area is the surficial area
of the sediment segment and is internally multiplied in the framework internally by the sediment
porosity. The concentration of chemical diffusing is the dissolved fraction per unit volume of
pore water. The actual diffusive exchange between sediment segments i and j at time t is given
by:
n}
n.
(1.11)
where: fdt, fdj = dissolved fraction of chemical in i and j [dimensionless]
ny = average porosity at interface "ij", Lw/L (volume of water/volume total
solution) [dimensionless]
Ejj(t) = diffusion coefficient time function for exchange "ij", m2/day [L2/T]
Ay = interfacial area shared by segments i and j, m2 [L2]
Lcij = characteristic mixing length between segments i and j, m [L]
1.4.5 Evaporation and Precipitation
Evaporation and precipitation can play a significant role in the overall water balance for
some water bodies. These processes can also affect contaminant concentrations. Precipitation
may be a low level source of some pollutants. As water evaporates, pollutant concentrations may
increase. In IPX, evaporation and precipitation are input via transport field 7 of Data Group D.
Different groups of segments having the same evaporation or precipitation rates may be defined
by the user. For each group, the user must supply a time function giving evaporation or
precipitation velocities (in m/sec). For each segment in the group, the user must also supply the
surface area (in m2) and a segment number pair. For evaporation the pair is from "ISEG" to "0"
and for precipitation the pair is from "0" to "ISEG", where ISEG must be the number of a surface
water segment. The precipitation input to segment i (i = ISEG) at time t is given by:
15
-------
(1.12)
where: Pi(t) = precipitation velocity time function for segment i, m/day [L/T]
Aj = surface area of segment i, m2 [L2]
C0i(t) = concentration of pollutant in precipitation, mg/L [M/L3]
In addition, segment volumes are adjusted:
(1.13)
where: Ei(t) = evaporation velocity time function for segment i, m/day [L/T]
Although it is possible to simulate contaminant inputs to a water body by precipitation,
there are limitations to this approach. IPX only allows the user to specify a single boundary
condition for any given model segment. For example, consider a simple model network that is a
single water column segment. This segment could have boundaries with surrounding water and
sediments as well as with the atmosphere. However, the user can only input a single boundary
condition for this segment and this one boundary condition will be applied to all transport
processes across all boundaries. If contaminant concentrations at each boundary significantly
differ, the boundary load computed for any process may not accurately represent the contaminant
mass input to the segment. Consequently, precipitation inputs to the model are best treated as
point loads rather than boundary conditions (loads).
It should be noted that no chemical or solids mass is transported out of the model by
evaporation. Volatilization (Section 1.6.7) is used to describe the movement (loss or gain) of
chemical state variables across the air-water interface. Also note that concentration changes
attributable to changes in water segment volumes are automatically computed.
1.5
SEDIMENT TRANSPORT AND THE SEDIMENT BED
As noted in Section 1.4, transport processes in IPX are divided into seven distinct types,
or fields. Transport fields 1, 2, and 7 are used to specify advective, dispersive, and evaporative
transport. Transport fields 3, 4, 5, and 6 are used to specify sediment transport. This section
describes sediment transport and the operation of transport fields 3, 4, 5, and 6. Advective,
16
-------
dispersive and evaportative transport and the operation of transport fields 1, 2, and 7 are
described in Section 1.4.
The third, fourth, and fifth transport field types move particles and particulate phase
pollutants and are used to specify settling for each of three possible solids state variables that can
be simulated. Transport field 3 moves particle state variable 1 and all chemical state variables
associated with particle state variable 1. Similarly, transport field 4 moves particle state variable
2 and field 5 moves particle state variable 3, as well as all chemical state variables associated
with those particles. The sixth transport field type moves particles and particulate phase
pollutants and is used to specify resuspension for all possible particle state variables that can be
simulated.
Note that transport fields 3-6 describe interactions between the water column and
sediment bed. Interactions within the sediment bed, described by burial and scour (unburial), are
automatically computed from the difference between the gross settling and resuspension fluxes.
Two options are available to represent interactions within the sediment bed: 1) Eulerian, and 2)
semi-Lagrangian. These options are described in Sections 1.5.3 and 1.5.4. Using either option,
the user must specify the particulate fraction or partitioning parameters for each chemical. Initial
(at t = 0) dissolved chemical fractions are entered in Data Group J.
1.5.1 Settling
Particles (sediments) and particulate phase chemicals in the water column may settle
through water segments and deposit onto surficial sediment bed segments. Settling rates in IPX
are described by velocities and surface areas in transport fields 3, 4, and 5 of Data Group D.
Particle settling velocities may vary both in time and in space, and are multiplied by cross-
sectional areas to obtain transport rates for solids and the particulate fractions of chemicals.
An upper bound for particle settling velocities is given by Stokes' Law:
where: vs =
g =
u =
pp =
pw =
(1.14)
Stokes' velocity for a particle of diameter dp and density pp, m/day [L/T]
acceleration of gravity = 981 cm/sec2 [ITT2]
absolute viscosity of water = 0.01 poise (g/cm-sec) at 20°C [M/L/T]
density of particle, g/cm3 [M/L3]
density of water, g/cm3 [M/L3]
17
-------
Table 1.2. Stokes' settling velocities for a range of particle sizes and densities at 20°C.
Particle Diameter, mm
Fine Sand
0.3
0.05
Silt
0.05
0.02
0.01
0.005
0.002
Clay
0.002
0.001
Particle Density, g/cm3
1.80
2.00 | 2.50
2.70
Stokes' Law Settling Velocity (m/day)
300
94.0
94.0
15.0
3.80
0.94
0.15
0.15
0.04
400
120
120
19.0
4.70
1.20
0.19
0,19
0.05
710
180
180
28.0
7.10
1.80
0.28
0.28
0.07
800
200
200
32.0
8.00
2.00
0.32
0.32
0.08
dp = particle diameter, mm [L]
Stokes' Law settling velocities for a range of particle sizes and densities are presented in Table
1.2.
While Stokes' Law is a valid means to approximate the settling velocity of discrete,
spherical particles under quiescent conditions, the actual rate at which particles settle is often
substantially less than predicted by Equation 1.14. In addition, sediment particle sizes, and
therefore settling velocities, typically vary over several orders of magnitude (Lee et al. 1981;
Filkins et al. 1993). As a result, it is difficult to estimate settling velocities without detailed
consideration of the characteristics of the particles to be simulated. To accommodate this range
of characteristics, it is useful to divide the range of suspended solids into three particle classes:
coarse, medium, and fine. These classifications divide the particle grain size distribution into
more easily defined groups whose properties are less difficult to characterize.
Coarse particles (>62 urn) are generally non-cohesive and, compared to finer particles,
may settle at rates similar to those described by Stokes' Law under quiescent conditions.
However, since natural particles are rarely spherical Stokes' Law may overestimate settling
velocities. A number of empirical relationships describing settling velocities of non-cohesive
18
-------
particles such as sands are available. A summary of representative relationships is presented by
Yang (1996). For non-cohesive (fine sand) particles with diameters from 62 to 500 ^m, settling
velocity can be computed as (Cheng, 1997):
(1.15)
(1.16)
d
<*. =.
1/3
where: vs
v
dp
d*
S
g
quiescent settling velocity, m/s [L/T] .
kinematic viscosity of water, m2/s [L2/T] = 1.007 x 10~6 at 20 °C
average particle diameter, m [L]
non-dimensional particle diameter [dimensionless]
specific gravity of particle [dimensionless] = 2.65 for pure sands
gravitational acceleration, m/s2 [L/T2] = 9.81
Medium particles (<62 (am) are often cohesive and may flocculate. Floe size and settling speed
depend on the conditions under which the floe was formed (Burban et al. 1990). Settling
velocities of cohesive particles can be approximated by:
v =
where: vs =
BI =
G =
dm =
62=
Cd =
VA =
(1.17)
(1.18)
floe settling velocity cm/s [ITT]
experimentally determined constant = 9.6 x 1Q-4
internal fluid shear stress, dynes/cm2 [M/L/T2]
median floe diameter, cm [L]
experimentally determined constant = 7.5 x 1Q-4
coefficient of drag [dimensionless] = 0.0015 for internal fluid shear stress
advective velocity, cm/s [L/T]
Under conditions typically found in tributaries, floe settling speeds range from 60-160 |im/s
(Burban et al. 1990).
Fine particles (<10 jam) generally may not extensively flocculate and typically have the
smallest settling velocities as a result of their size, shape, density, and other physicochemical
19
-------
properties. For example, clay particles often have negative electrical charges which can inhibit
direct particle aggregation in dilute suspensions. Also, biotic materials such as algae often fall
into this size class of particles. Algae in particular generally possess a variety of mechanisms to
minimize their settling velocities (Wetzel, 1983). As a result of these attributes and other
conditions, fine particles may have near-zero settling velocities. For simplicity, the settling
velocity of fine particles is often represented by a small constant value.
Not all particles settling through the water column will necessarily reach the sediment-
water interface or otherwise be incorporated into the sediment bed. As a result, effective settling
velocities of particles in flowing water are usually less than quiescent settling velocities. The
effective settling velocity can be described as a reduction in the quiescent settling velocity by the
probability of deposition:
v
v
= P, v
fdep v i
(1.19)
where: vse = effective settling velocity [L/T]
pd = Probability of deposition [dimensionless]
The probability of deposition varies with shear stress near the sediment bed and particle size. As
particle size decreases or shear stress increases, the probability of deposition decreases.
For non-cohesive particles, the probability of deposition has been described as a function
of bottom shear stress (Gessler, 1967):
(1.20)
(1.21)
where: P
a
= probability integral for the Gaussian distribution
= experimentally determined constant = 0.57
= critical shear stress for deposition of non-cohesive particles, defined as the
shear stress at which 50% of the particles in the size class settle
The critical shear stress for non-cohesive deposition, Tcd, can be computed by determining the
upward force needed to balance the net downward force (gravity force minus the sum of the
20
-------
buoyant and drag forces) for a particle with a diameter equal to the mean diameter of the size
class (i.e. dp = d50).
For cohesive particles, the probability of deposition has been described as a function of
bottom shear stress (Patheniades, 1992):
P =i_p = i
"
cr
(1.22)
(1.23)
where: a
experimentally determined constant = 0.49
critical shear stress for deposition of cohesive particles, defined as the shear
stress at which 100% of the particles in the size class settle
The probability integral in Equations 1.20 and 1.22 can be approximated as (Abramowitz and
Stegun, 1972):
P = 1 - F(y)(o.4362X - 0.1202X 2 + 0.9373X3) for Y > 0
for7<0
(1.24)
= (l + 0.3327F)~1
(1.25)
(1.26)
It is not always necessary or desirable to have more than one state variable to represent
solids. When explicit representation of multiple particle types is not desired, it is possible to
conceptually represent multiple particle types as a single state variable by characterizing the
particle grain size distribution as a function of flow as described by Gailani et al. (1991). Under
low flow conditions, fine particles predominate with a small fraction of medium materials. At
moderate flows, medium materials predominate with a smaller amount of fine particles and the
potential for a minor fraction of coarse materials. At high flows, one of two conditions may
exist: 1) on the rising limb of the hydrograph (at the beginning of a flow event) most suspended
solids originating from the sediment bed will be medium particles until the flow peak is reached
(when shear stresses at the sediment water interface are highest) and significant erosion of coarse
21
-------
Table 1.3. Generalized grain size distribution functions for regulated tributaries.
Condition
Non-Event
Event:
Rising Limb
Event:
Falling Limb
River Flow
Q r7nr.)
U;>+ ^, \*L ysoj
650
Coarse (fc)
0.0
0.0
0 0
0.7
07 °'31fO OCA)
° 650 ^ G80J
Eq.
(1.27)
(1.28)
n 291)
v1-i='y
(1.30)
a<21-\
.ji;
Table 1.4. Generalized grain size distribution functions for unregulated tributaries.
Condition
Nfon-Event
Event:
Rising Limb
Event:
Falling Limb
River Flow
Q
-------
When simulating solids as a single state variable (TSS), overall particle settling rates are
computed as a function of the particle size distribution and the settling velocities of each particle
class: . '
(1.36)
where: vs =
Vsc =
vsm=
vsf =
fc =
fm =
ff =
total suspended solids (TSS) settling velocity [L/T]
coarse particle class settling velocity [L/T]
medium particle class settling velocity [L/T]
fine particle class settling velocity [L/T].
fraction of coarse class particles
fraction of medium class particles
fraction of fine class particles
When simulating solids as two or three independent state variables, the particle grain size
distribution functions can be used to divide total solids loads into loads for each state variable.
1.5.2 Resusoension. Sediment Aging, and the Surficial Sediments
Particles and particulate phase chemicals in the sediment bed may resuspend from the
surficial sediments and return to the water column. Resuspension rates in IPX are described by
velocities and surface areas in transport field 6 of Data Group D. Unlike settling, where the user
must input settling velocities for each particle type state variable simulated, only one transport
field is used to specify resuspension velocities even if more than particle state variable is
simulated. IPX is designed so that all particle types simulated resuspend in proportion to their
abundance. Particle resuspension velocities may vary both in time and in space, and are
multiplied by cross-sectional areas to obtain transport rates for solids and the particulate fractions
of chemicals.
Particles leaving the water column by settling enter the surficial layer of the sediment bed.
IPX allows the user to define two sediment bed types on the basis of vertical location: surficial
and subsurface sediments. Surficial sediments are comprised of age layers. Sediments settling to
the bed over any defined interval (generally one day) enter the surficial age layer of surficial
sediments. IPX allows the user to define an arbitrary number of sediment age layers. Sediments
in each layer undergo an aging process characterized by changes in the physical characteristics of
the layer. These changes alter the erosion potential of the sediment bed. Freshly deposited
sediments are more resuspendable than those that have undergone some degree of aging. With
respect to its impact on sediment erosion potentials, the sediment aging process is typically
23
-------
completed in seven days (Tsai and Lick, 1987; Xu, 1991). Sediments in the final age layer of
surficial sediment have the lowest resuspension potential. Not all age layers within surficial
sediment segments will necessarily hold sediment at all times during a simulation. Frequently,
the oldest layer will hold most of the sediment mass.
For any given resuspension event, the particle resuspension flux can be described as a
function of the shear stress at the sediment-water interface, which can in turn be approximated as
a function of flow (Ziegler et al. 1988; Gailani et al. 1991):
(1.37)
(1.38)
m
where: ST>TC = amount of sediment resuspended when the shear stress exerted at the
sediment-water interface (t) exceeds the critical shear for entrainment (tc),
g/cm2 [M/L2]
ao = empirical sediment yield constant
Z = empirical sediment age constant
T = shear stress exerted at the sediment-water interface, dynes/cm2 [M/L/T2]
tc = critical shear stress for entrainment, dynes/cm2 [M/L/T2]
empirical sediment entrainment exponent
the amount of sediment resuspended when t is less than or equal to the critical
shear for entrainment, g/cm2 [M/L2]
Cf = coefficient of friction [dimensionless] = typically 0.002-0.004
VA = advective water velocity, cm/sec [L/T]
The parameters tc, m, ao, and Z depend on the physical characteristics (age, water content,
cohesiveness, etc.) of a particular sediment. Once tc is exceeded, sediments are quickly
entrained. From the amount resuspended (e), a resuspension velocity can be computed:
(1.39)
Ptfe
where: vr = resuspension velocity [L/T]
Pb = bulk density of sediments [M/L3]
24
-------
te = time to entrain sediments [T]
The most significant factor controlling erosion is the shear stress exerted at the sediment-
water interface by water flowing over the sediments, x. When x exceeds the critical shear for
entrainment, significant resuspension can occur. IPX dynamically computes x for each surficial
sediment segment as a function of the river flow and channel cross-sectional area as described in
Equation 1.38.
When average shear stresses are below the critical value, resuspension may occur at a
background level that is significantly less than when shear stresses exceed the critical threshold.
This is defined as background resuspension. Background resuspension is a lumped parameter
used to represent any particle movements that can occur when the average shear stresses are less
than the average critical shear stress for erosion. Very fresh sediments can form "fluff layers
that have significantly different erosion properties (greater yields coefficients and much lower
critical shear stresses) than older sediments (Gailani et al. 1991). In addition, bed sediments are
comprised of a mixtures of particle types and sizes. Even within a single size class, particles
smaller than the average particle size (dso) may begin to resuspend before larger particles within
the class. Background resuspension is typically so small that it does not significantly impact
sediment bed morphometry of resultant suspended solids concentrations in the water column.
However, in tributaries with significant in-place pollutant reservoirs, hydrophobic contaminants
are typically present in the sediments at levels much higher than found in the water column
relative to paniculate materials. As a consequence, even minute resuspension (a little as several
mm/year) can result in significant increases in water column contaminant concentrations.
Therefore, although background resuspension is typically unimportant for accurate sediment
transport simulation, it is necessary to accurately simulate contaminant transport.
Other significant factors that influence erosion are sediment armoring and the extent of
sediment aging. These factors are described through the parameters xc, m, ao, and Z. Sediments
armor in response to the differential erosion of finer-grained, more easily resuspendable particles
from the sediment bed, which leaves behind less resuspendable materials. IPX tracks the
computed shear stress exposure history of the sediments. As a resuspension event progresses, the
shear stress to which the sediments have been exposed increases (TC increases), armoring the
sediment bed, and making further resuspension more difficult. IPX also automatically
resuspends fresh, readily erodible sediments (with a low xc) that settle over an armored layer
(with a high Tc) during periods such as the declining limb of an event when the shear stress
exerted at the sediment-water interface is less than the critical shear for resuspending the armored
25
-------
sediments but is greater than the shear needed to resuspend fresh (unarmored, less consolidated)
sediments (TCfreSh < T < ^caimored)- Tne parameters ao and m express how readily erodible a given
sediment is. Sediment resuspension is a highly nonlinear function of flow. Laboratory
experiments have determined that for many sediments m is in the range of 2 to 3 (Xu, 1991; Lick
et al. 1995). For a limited number of sediments, ao has been determined to be in the range of
0.27 x 10~3 to 8 x 10~3 (Xu, 1991, Lick et al. 1995). The parameter Z expresses the effect of
sediment age on resuspension. Z is an exponential function of the time after deposition of the
sediments in a given age layer and can vary from 0.1 to 50 (Tsai and Lick, 1987; Xu, 1991).
When computing resuspension velocities from Equations 1.37 - 1.39, the sediment age
factor (Z) of the first age layer should be used for Equation 1.37. Z for the first age layer may
differ from river reach to reach, however. In this manner, the resuspension potential of each age
layer can vary relative to the resuspension potential of the first age layer as resuspension proceeds
through each age layer. Users should also be aware that Equation 1.37 is applicable to cohesive
sediments. However, sediments largely comprised of non-cohesive particles do not armor as
extensively as cohesive sediments. As a consequence, during a resuspension event non-cohesive
sediments will continue to resuspend (beyond the amount £ that would be applicable to cohesive
sediments) until the shear stress at the sediment-water interface falls below the critical shear or
the supply of resuspendable sediments is exhausted.
1.5.3 Burial and Scour (Unburial)
Burial and scour (unburial) are the sediment transport processes by which particles and
contaminants interact within the sediment bed. Two options are available for representing these
interactions: 1) Eulerian, and 2) semi-Lagrangian. The sediment bed option is determined by
value of the input variable TTYPE (defined in Section 2.2.3). Sediment segments where ETYPE
is 3 or 4 use an Eulerian frame of reference. With an Eulerian frame of reference, an open
boundary exists at the bottom of the sediment stack. For each Eulerian sediment stack, the user
specifies a boundary condition at the bottom of the stack and material can move into or out of the
model network across the bottom boundary. Sediment segments where ITYPE is 5 or 6 use a
semi-Lagrangian frame of reference. With a semi-Lagrangian frame of reference, a barrier (such
as bedrock or hardpan materials) exists at the bottom of the sediment stack. For each semi-
Lagrangian sediment stack, no boundary condition at the bottom of the stack is specified (since
there is no open boundary) and no material moves into or out of the model network across the
bottom barrier. For any stack of sediments (from the sediment-water interface to the bottom-
most layer identified), consistent specification of the sediment column as Eulerian or semi-
Lagrangian is necessary. For example, specification of a sediment stack with an Eulerian surface
26 .
-------
layer (TTYPE = 3) and semi-Lagrangian subsurface layers (TTYPE = 6) directly below will cause
a simulation error condition.
1.5.3.1 Eulerian Option
With the Eulerian sediment bed option, the model frame, of reference is fixed at the
sediment-water interface. Particles and associated contaminants move between surficial and
subsurface. sediment segments by burial and scour. Burial (which may be positive or negative in
direction) is automatically computed in the framework from the difference between the particle
settling and resuspension fluxes; the user does not specify burial/scour velocities.
The solids concentration in all sediment segments is held constant. The volume (depth)
of all surficial sediment segments is allowed to vary from their initial values in response to
erosion and deposition. The volume of all subsurface sediment segments is held constant. At a
user defined interval, representative of the physical, biological, and chemical processes that
change sediment properties, solids mass and associated contaminants are buried or scoured as
needed to restore all surficial sediment segments to their initial volumes:
(1.40)
1=1
where: JB = burial flux, g [M]
AM = change in particle mass of the surficial sediment, g [M]
AV = change in volume of the surficial sediment segment, m3 [L3]
ti - ti_i = time interval over which settling, resuspension occur, days [T]
vs = settling velocity, m/day [L/T]
vr = resuspension velocity, m/day [L/T]
GI = solids concentration in the water column, mg/L (g/m3) [M/L3]
€2 = solids concentration in the surficial sediment, mg/L (g/m3) [M/L3]
A = interfacial area across which settling or resuspension occurs [L2]
Whenever burial is computed, sediment (particle) mass is moved (advected) through all
subsurface sediment segments that underlie a given surficial sediment segment (forming a
vertical stack) as needed to restore the surficial sediment segment volume to its initial value.
When burial is positive, sediment mass is moved from the surficial sediments, through the all
subsurface sediments, and across the bottom boundary of the model while pore water and
associated chemicals are "squeezed out" of each layer and transported to the layer above. When
27
-------
burial is negative (scour), sediment mass is moved from the bottom boundary of the model, into
'the bottom segment, through the all subsurface sediments, and into the surficial sediments.
1.5.3.2 Semi-Lagrangian Option
. With the semi-Lagrangian sediment bed option, the model frame of reference is fixed to
the bottom of the sediment column (including any elements in the ghost stack) and the elevation
of the sediment-water interface is allowed to float (relocate) in response to sediment transport.
Particles and associated contaminants "move" between surficial and subsurface sediment
segments and the ghost stack by a segment re-indexing procedure analogous to burial and scour.
Segment re-indexing may be downward or upward in direction. The downward re-indexing of
segments is referred to as "push" (analogous to burial) and. the upward re-indexing as "pop"
(analogous to scour). Segment re-indexing is automatically computed in the framework based on
the surficial sediment segment volume (depth) conditions .(which is itself determined by the
difference between the particle settling and resuspension fluxes); the user does not specify
burial/scour velocities.
The solids concentration in all sediment segments is held constant. The volume (depth)
of all surficial sediment segments is allowed to vary from their initial values in response to
erosion and deposition. The physical characteristics (including volume) of all subsurface
sediment segments is held constant until a re-indexing event occurs. Re-indexing occurs when
the surficial sediment segment reaches the minimum or maximum volume (depth) condition
specified by the user. When a re-indexing event occurs, model sediment segments and ghost
stack elements exchange identity to account for the net increase or decrease in sediment bed
elevation. When segments exchange identity, the properties of each segment/ghost are assigned
to a neighboring segment (or ghost stack element) below' (push) or above (pop) throughout the
sediment stack. In response to reaching the maximum volume condition, the surficial sediment
segment is split to form two segments (one surficial and one subsurface) and all subsurface
sediment segments (and ghost stack elements) are "moved" one element down the stack. When
this occurs, the sediment stack at that location will contain one additional member and the total
sediment bed elevation at that location increases. In response to reaching the minimum volume
condition, the residual volume (if any) of the surficial sediment segment is joined with the
subsurface segment (or ghost stack element) immediately below to form one surficial segment
and all subsurface sediment segments (and ghost stack elements) are "moved" one element up
the stack. When this occurs, the sediment stack at that location will contain one less member and
the total sediment bed elevation at that location decreases. Note that after upward re-indexing,
28
-------
the formerly occupied bottom-most segment or ghost will be "empty" (unoccupied). In this
manner, mass is never transported out of or into the model network across the bottom barrier.
Sediment bed elevations at any location can change over time. In cases where the total
supply of sediments in the sediment column at any location is exhausted, the volume of all
segments (and ghost stack elements) can be zero. Zero sediment volume represents reaching the
barrier at the bottom of the sediment column at that location. When this occurs, no further loss
of sediments (i.e. resuspension) can occur, regardless of the resuspension velocity specified,
because no sediments are available to resuspend. If sediments deposit to the barrier at the bottom
of the sediment column, the sediment column is reconstructed (from the top down) once stack
element at a time. In cases where sediment bed elevation increases beyond the initial elevation
of the sediment column, elements are added to the ghost stack until all "empty" (unoccupied)
elements in the stack are filled.
1.5 A The Sediment Bted
The sediment bed plays an important role in the transport and fate of hydrophobic
contaminants. Particles and associated pollutants in the surficial sediments may enter deeper
sediment layers by burial or be returned to the surficial sediments by scour. The movement of
particles and contaminants is automatically computed in IPX; the user does not specify burial
velocities. Interactions between segments in the sediment column depend on specification of the
Eulerian or semi-Lagrangian sediment bed options as presented in Section 1.5.3.
1.5.4.1 Surficial Sediments
As described in Section 1.5.2. the surficial sediments are comprised of age layers. IPX
allows the user to define an arbitrary number of sediment age layers. Sediments settling to the
bed over any defined interval (generally one day) enter the surficial age layer of surficial
sediments. Sediments in each layer undergo an aging process characterized by changes in the
physical characteristics of the layer. These changes alter the erosion potential of the sediment
bed. Freshly deposited sediments are more resuspendable than those that have undergone some
degree of aging. Sediments in the final age layer of surficial sediment have the lowest
resuspension potential. Not all age layers within surficial sediment segments will necessarily
hold sediment at all times during a simulation. Frequently, the oldest layer will hold most of the
sediment mass.
As particles settle to surficial sediments, the particles enter the first (freshest) age layer of
the surficial sediments and the volume (depth) of those segments increases. However, the
particle concentration (density/porosity) in each surficial sediment segment is constant. At a user
29
-------
specified sediment aging interval, particles are "moved" to the next (older) age layer. A typical
sediment aging interval is 1 day. A typical number of sediment age layers is 7. If no new (fresh)
sediment enters the segment, all particles within the segment will eventually enter the final
(oldest) age layer (and all other layers will be empty). All particles in the final age layer have the
same resuspension potential regardless of how long they have been in the layer. Users should be
aware that particles are never transported to the subsurface sediments as a result of aging. Unless
subsequently resuspended, particles will remain in the surficial sediments until burial occurs.
As particles are resuspended from surficial sediments, particles are removed from one or
more age layers and the volume (depth) of those segments decreases. (Recall that the total
particle concentration of the surficial sediments will be constant.) The particle mass removed
from the surficial sediments as a result of resuspension will be equal to the amount e (as
computed from Equations 1.37 - 1.39), modified by the sediment age factor (Z) of the age layer
from which resuspension is occurring. If the mass of sediment to be resuspended exceeds the
mass present in an age layer, all particles in that layer are resuspended and the process proceeds
with the next age layer, sequentially from freshest to oldest, until the proper particle mass is
resuspended. During this process, several age layers may be temporarily unoccupied (empty)
until particles enter the empty layers by settling and subsequent aging. However, IPX is designed
so that if the mass of particles to be resuspended during any time-step exceeds the total mass
present at that time, resuspension will stop (i.e. resuspension stops if the available supply of
resuspendable sediments is exhausted). If at the end of any simulation time-step, the volume of
any surficial sediment segment is less than the minimum volume (depth) condition, burial is
automatically triggered.
The user should be aware that sediment age layers serve only to describe how the
resuspension potential of sediments change with time. Age layers are conceptual divisions
within a surficial sediment segment; the number of model segments is not affected by the number
of age layers simulated in the sediments. For the purpose of computing a mass balance, the
sediments in all age layers are considered. Alternatively stated, the sum of the sediment masses
in the age layer equals the total sediment mass in the segment; the sum of the depths of the age
layer equals the total depth, and therefore volume, of the segment. Particle concentrations in the
surficial sediments are constant.
1.5.4.2 Subsurface Sediments
Eulerian Sediment Bed Option: For Eulerian representations of the sediment bed, burial
computations are performed at a fixed time interval (regardless of sediment volume conditions)
30
-------
as specified by the input variable TDINTS in Data Group C. The interval defined should be
representative of the time required for physical, biological, and chemical processes, such as
bioturbation and compaction, to change sediment properties. In addition, a minimum volume
condition also exists for the Eulerian sediment bed option such that burial computations are
performed whenever any surficial sediment segment is reduced to 25% of its initial volume. Any
time burial computations are performed, the volumes of all surficial sediment segments are
globally restored to their initial volumes. A description of burial computations for the Eulerian
sediment bed option is presented in Section 1.5.3.
Whenever burial/scour occurs, particles and associated chemicals are transported through
each subsurface sediment segment within a vertical stack under a given surficial sediment
segment. If the volume of the surficial sediment segment at the top of the stack is greater than its
initial volume, burial occurs and particles and chemicals are transported downward to the first
subsurface sediment segment. However, the particle concentration and volume of each
subsurface sediment segment are constant. An equal mass of particles, and the chemical mass
associated with these displaced particles, is transported from the uppermost subsurface sediment
segment, through each subsequent subsurface layer, and across the bottom boundary of the
deepest subsurface sediment segment in the stack.
If the volume of the surficial sediment segment at the top of the stack is less than its
initial volume, scour occurs and particles and chemicals are transported upward across the
bottom boundary of the model into the deepest subsurface sediment segment of the stack.
(Recall that the particle concentration and volume of each subsurface sediment segment are
constant.) Particle and chemical concentrations at this and all other boundaries are specified by
the user. For the sediment bed, particle boundary conditions must be greater than zero.
Chemical boundary conditions for the sediments can have any positive value. However, the
model network should be developed so that chemical concentrations at the sediment
boundary of the system are zero. Users should be aware that if the chemical sediment
boundary condition is greater than zero, the boundary can act as an infinite source of the
chemical to the system. An equal mass of particles, and associated chemical mass, is then
transported from the deepest subsurface sediment segment, through each overlying subsurface
layer, and into the oldest age layer of the surficial sediment segment at the top of the stack.
Pore water and mobile-phase (dissolved and DOC-bound) chemicals are squeezed
upward through the sediments and into the water column whenever "burial or scour occur.
Although the particle concentration and volume of each subsurface sediment segment .remain
31
-------
constant, chemical concentrations will vary during a simulation. If burial exceeds scour in the
long term, chemical concentrations in the sediments will approach the chemical concentration on
particles in the water column. Conversely, if scour exceeds burial, chemical concentrations in
the sediments will approach the chemical concentration specified at the bottom boundaries of the
model network.
Semi-Lagrangian Sediment Bed Option: For semi-Lagrangian representations of the
sediment bed, burial computations (segment re-indexing) are performed when a surficial
sediment segment reaches either the minimum or maximum volume (depth) conditions as
specified by the input variables MINDEP and MAXDEP in Data Group L. The minimum
volume condition for a stack is computed as the product of MINDEP and the initial volume of
the surface sediment segment. Similarly, the maximum volume condition is computed as the
product of MAXDEP and the initial volume of the surface sediment segment. The minimum and
maximum volume conditions defined should be representative of the physical (and biological)
processes that alter sediment properties and contaminant distributions over time. Burial
t
computations are performed on an individual basis as each stack reaches a minimum or
maximum volume conditions. A description of burial computations (push and pop) for the semi-
Lagrangian sediment bed option is presented in Section 1.5.3.
Whenever burial (segment re-indexing) occurs, the identities and properties (volume,
area, thickness (depth), state variable concentrations) of each segments in the sediment stack are
exchanged. When the maximum volume condition is reached, the surficial sediment segment is
split to become two segments. The upper portion of the split surface sediment segment becomes
the new (smaller) surficial sediment segment, the lower portion of the split segment becomes the
next segment in the stack, and all other segments (and ghosts) are pushed down the stack and re-
indexed. The downward re-indexing of sediments is managed by Subroutine PUSH. When the
minimum volume condition is reached, the residual volume of the surficial sediment segment (if
any) is joined with the volume of the next element in the sediment stack (if that element is
occupied) to become one segment. The joined segment becomes the surficial sediment segment
and all other segments (and ghosts) are popped up the stack and re-indexed. The upward re-
indexing of sediments is managed by Subroutine POP.
During the course of a simulation, the entire sediment stack at any location, including any
ghost stack elements, can be empty (unoccupied). Conversely, the entire sediment stack and all
ghost elements can be occupied. If downward re-indexing occurs at a locations where all
elements of the stack are occupied, the simulation will proceed according to the ghost collapse
32
-------
option. The ghost collapse option is set by the input variable IGOPT specified for that stack in
Data Group L. A separate IGOPT value is specified for eaeh sediment stack that uses the semi-
Lagrangian sediment bed option. The ghost collapse option allows the bottom two ghost
elements in a stack to merge to create an empty element in the stack prior to downward re-
indexing. If IGOPT = 0, ghost stack elements do not collapse. When all stack elements are full,
the next downward re-indexing event for that stack causes the simulation to abort and display a
message indicating the location (surficial sediment segment number) and time that the condition
occurred. If IGOPT = 1, the bottom two ghost elements of the sediment stack merge to form a
single (larger) ghost. When all stack elements are full, the next downward re-indexing event for
that stack causes the bottom two ghosts to merge and the simulation continues. Subsequent
downward re-indexing events result in further ghost collapse and the (potentially infinite) growth
of the collapsing ghost element at the bottom of the sediment stack. Regardless of the value of
the ghost collapse option, mass is never transported out of or into the model network across
the bottom barrier of a sediment stack using the semi-Lagrangian sediment bed option.
1.5.5 Screening-Level IPX Sediment Transport Simulation Guidelines
1. Resuspension events are defined by the hydrograph. These events can occur at random but
typically occur with greater frequency in the Spring and Fall. The duration of resuspension
events is a function of the hydrograph, starting when flows begin to increase, through the
flow peak, and continuing until flows decrease to base levels. Depending on the system,
events can last from a few days (or less) to several weeks or even months.
2. High flow events disturb the sediments and alter their resuspension properties. Following a
resuspension event, physical, chemical, and biological processes restore sediment
resuspension properties to their typical values in 30 days. Additional high flow events (peak
flow must exceed the previous peak) that occur within this 30 day window should be treated
as continuations of the first event of the series that disturbed the sediments. The sediments
return to their typical resuspension properties in 30 days following the last peak (highest
overall flow). ,
3. Settling velocities of medium-grained, cohesive particles may vary by season. During
Winter, when particles are largely inorganic, settling velocities are greatest. During Summer,
when particles are most organic (biotic), the settling velocities are smallest. During Spring
and Fall, settling velocities are somewhere between the Winter and Summer values.
33
-------
4. Background resuspension of fine-grained sediments occurs because these materials are kept
in a quasi-disturbed state by water flowing over or through the surficial sediments. In fast-
moving river reaches where coarser-grained materials predominate, background resuspension
is minimal (less than in slower reaches) because fine-grained materials do not accumulate
there (i.e. are not available for resuspension) and under non-event conditions coarse-grained
materials do not easily resuspend.
5. The extent of apparent background resuspension may vary by river reach and by season as a
function of local water velocity/shear stress differences. Background resuspension may be
greater in high velocity reaches and during the Spring and Fall when flows tend to be larger.
6. High velocity river reaches have less cohesive sediments available for resuspension so
sediments in these reaches are more difficult to resuspend (greater Z and/or lower m values).
7. Easily resuspendable surficial sediments may build up over the Winter under ice-cover. The
first resuspension event in the Spring after ice-off may be strongly influenced by these easily
resuspendable sediments (smaller Z values in the early Spring for river reaches were
sediments may accumulate).
8. In areas where non-cohesive sediments predominate (often narrow, high velocity river
reaches) resuspension reaches its maximum rapidly but returns to background levels more
slowly than do areas with predominantly cohesive sediments (i.e. considerable erosion
continues during the declining limb of an event because non-cohesive sediments do not
armor to the extent that cohesive sediments do).
9. Substantial deposition can occur following the peak of a resuspension event as particles are
removed from the water column by settling. In regions where sediment accumulate, the
sediments may be relatively uncompacted and more resuspendable for a short period until
compaction or other processes occur that restore sediment properties to their typical values
(smaller Z and/or TC values for up to 30 days).
1.6 SIMULATING THE TRANSPORT AND FATE OF TOXIC CHEMICALS
USING IPX
In addition to physical transport due to the movement of water, physical, chemical, and
biological mass transfer and transformation processes can affect the transport and fate of toxic
chemicals in the aquatic environment. As presented in Figure 1.5. the transport and fate
34
-------
Kinetic Processes:
Biodegredation
Hydrolysis
Oxidation
Photolysis
... . Total C
Water
Advection ^
Dispersion
Atmospheric
Deposition
1 Volatilization '
ontaminant
DOC-bound
Contaminant
t
Surficial '
Sediments Porewatei
i
Subsurface
Sediments
>
•s
(sum of all phases)
binding
Dissolv<
Contamir
s"*^ ^
(^ Kinetic C
^ Processes *
r
•Transport
'
Porewate
•
k External Loads and
Boundary Conditions
> '
.
, ^ D rf" ! f ^ /\QvGCIIOn
^0 ^ r"^3riiCUi3IG ,^ '
iant partitioning & Contsminsnt p ^ >
Dispersion
\^'^/^^t^'t^^^^^^f^^^^
•
Besses
r
' 1 r
r Transport Settling Resuspension
\ t 1
1
Scour Burial
> Partitioning between dissolved, particulate and Ddb phases occurs in the sediments as conceptualized «:
«: in the water coiumn. Biodegradation, hydrolysis and oxidation can also occur in the sediments as
conceptualized in the water column. ^ >«'*-'<,.- "°
-------
Table 1.5. Concentration related symbols used in mathematical equations.
Symbol
q
C\yj
C" •
^* wj
Csj
Csj
CBI
CBJ
msj ,
Msj
Msj
B!
Bi
»j
Kps
KB
Definition
Concentration of total chemical in segment j.
Concentration of dissolved chemical in segment j.
Concentration of dissolved chemical, in water in segment j
on a water volume basis.
Concentration of chemical sorbed to sediment type s in
segment j.
Concentration of chemical sorbed to sediment type s in
segment j on a mass basis (C'Sj = Csj/Msj).
Concentration of DOC-bound chemical in segment j.
Concentration of DOC-bound chemical in segment j on a
mass basis (C'BJ = CBJ/BJ).
Concentration of sediment type s in segment j.
Concentration of sediment type s in segment j
(Mj = mj/106).
Concentration of sediment type s in water in segment j
(M's^Msj/n,-).
Concentration of DOC in segment j.
Concentration of DOC in water in segment j .
Porosity of segment j (Vwater/VWater + solids).
Partition coefficient of chemical for sediment type s.
DOC-binding coefficient of chemical.
Units
mgc/L
mgc/L
mgc/Lw
mgc/L
mgc/kgs
mgc/L
mgc/kgb
mgs/L
kgs/L
kgs/Lw
kgb/L
kgb/Lw
Lw/L
Lw/kgs
Lw/kgb
chemicals may be independent, such as different pesticides or coupled by reaction pathways and
yields, such as a parent compound-daughter product sequence. Users should be aware that unlike
WASP4 (TOXI4), IPX does not simulate ionic speciation. Chemicals can exist in five phases:
dissolved, bound to dissolved organic carbon (DOC), and sorbed to each of the three possible
solids types that can be simulated. Local equilibrium is assumed and the distribution of the
chemical between each phase is defined by partition coefficients. The concentration of a
chemical in any phase can be computed from the total concentration of the chemical.
In the aquatic environment, chemicals may be transferred between media (air, water, and
solids) and may be transformed by a number of physicochemical and biological processes.
Defined mass transfer processes include sorption (equilibrium partitioning) and volatilization.
Transformation processes include biodegradation, hydrolysis, photolysis, and oxidation.
Sorption is treated as an equilibrium partitioning process. All mass transfer and transformation
36
-------
processes are described by rate equations. Rate equations may be quantified by first-order
reaction constants or by second-order, chemical-specific and environment-specific reaction
constants that may vary in space and time.
IPX uses a finite difference form of Equation 1.2 (presented in more detail as Equations
1.124 and 1.125) to calculate sediment and chemical mass and concentrations for every segment
in a specialized network that may include surface water, underlying water, surface bed, and
underlying bed. In a simulation, sediment is treated as a conservative constituent that is advected
and dispersed through water segments, settles to and resuspends from surficial sediment
segments, and moves between surficial and subsurface sediment segments burial, unburial, or
bed load.
During a simulation, chemicals can undergo several physical, chemical, or biological
transformations. It is convenient to group these into fast and slow reactions. Fast reactions have
characteristic reaction times on the same order as the model time-step and are handled with the
assumption of local equilibrium. Slow reactions have characteristic reaction times much longer
than the model time-step. These are handled with the assumption of local first-order kinetics
using a lumped rate constant specified by the user, or calculated internally, based on summation
of several process rates, some of which are second-order. Thus, the effective first-order decay
rate can vary with time, and space, and is recalculated as often as necessary throughout a
simulation. The chemical is advected and dispersed among water segments, and exchanged with
surficial benthic segments by dispersive mixing. Sorbed chemical settles through water column
segments and deposits to or erodes from surficial benthic segments. Within the bed, dissolved
chemical migrates downward or upward through percolation and pore water diffusion. Sorbed
chemical migrates downward or upward through net sedimentation or erosion. Both rate
constants and equilibrium coefficients must be estimated in most toxic chemical studies.
Although these can be calculated internally from chemical properties and local environmental
characteristics, site-specific calibration or testing is desirable.
One important limitation should be kept in mind when developing simulations: chemical
concentrations should be near trace levels. A good guideline is that concentrations should be less
than half the chemical solubility limit. At higher concentrations, the assumptions regarding
linear partitioning and other transport and fate processes may not apply. Also, at very high
concentrations, chemical density may become important, particularly near the chemical source;
such might be the case when simulating a chemical spill. Extremely high concentrations can
37
-------
affect key environmental characteristics such as pH or bacterial populations, thereby altering
transformation rates. IPX does not allow simulation of such feedback phenomena.
1.6,1 Equilibrium Sorption
Dissolved chemicals in the water column and sediments interact with particles and non-
filterable organic materials, conceptualized as dissolved organic carbon or DOC, to form five
phases: dissolved, DOC-bound, and sediment-sorbed (particulate) to each of the three possible
sediment types ("s"). These interactions can be written with respect to a unit volume of water:
n
(1.41)
n
(1.42)
For particles, the forward reaction is sorption and the reverse reaction is desorption. For
DOC, the forward reaction is binding and the reverse reaction is unbinding. These reactions are
usually fast in comparison with the model time-step, and can be considered in local equilibrium.
The phase concentrations Cw, Cs, and CB are governed by the equilibrium partition and binding
coefficients Kpso and KB (L/kg):
(1.43)
_CB/n_CB
jtX n — —'»—--"—'• '
B C C
(1.44)
These equations give the linear form of the Freundlich isotherm, applicable when
sorption/binding sites on sediment and DOC are plentiful:
C =
^ •
(1-45)
(1.46)
The partition and binding coefficients depend upon characteristics of the chemical and the
sediments or DOC to which partitioning occurs. Many pollutants of public health interest are
non-polar, hydrophobic, organic compounds. The partition and binding coefficients of these
compounds correlate well with the organic carbon fraction (foc) of the sediment. Rao and
Davidson (1980) and Karickhoff et al. (1979) developed empirical expressions relating
38
-------
equilibrium coefficients to laboratory measurements, leading to fairly reliable means of
estimating appropriate values. Conceptually, dissolved organic materials are composed entirely
of organic carbon (foc = !)• However, measured DOC-binding constants are generally one to two
orders of magnitude smaller than predicted (Eadie et al. 1990; Landrum et al. 1985; Landrum et
al. 1987; Capel and Eisenreich, 1990), suggesting that not all non-filterable material measurable
as organic carbon is available for binding. The partitioning and binding expressions used in IPX
are:
f
J o
(1.47)
(1.48)
where: KQC = organic carbon partition coefficient, Lw/kgoc [L3/M]
focs = organic carbon fraction of sediment or DOC
DE = DOC-binding effectiveness constant = 0.01 - 0.1
If no KQC values are available, a value can be generated internally from the octanol-water
partition coefficient (Kow) (Lw/Loct) using the following correlation:
logKoc = a0+allogK0
(1.49)
If ao and ai are not specified, default values of log 0.6 (-0.2218) for ao and 1.0 for ai are used.
The value of the partition coefficient depends on numerous factors in addition to the
fraction organic carbon of the sorbing/binding matrix. Of these, perhaps the most potentially
significant and the most controversial is the effect of particle concentration, which was first
presented by O'Connor and Connolly (1980). Based on empirical evidence, O'Connor and
Connolly concluded that the partition coefficient was inversely related to the particle
concentration. Much research has been conducted to prove or disprove this finding. At present,
the issue remains contentious. A particle interaction model has been proposed (Di Toro, 1985)
which describes the effects of particle concentration. This model was shown to be in conformity
with observations for a large set of adsorption-desorption data. The expression defining the
partition coefficient is:
ps
[ + M,K _n/vv
(1.50)
39
-------
where: Kpso = limiting partition coefficient with no particle interaction = focs Koc, L/kg
DL3/M]
Ms = solids concentration, kg/L [M/L3]
vx = particle interaction parameter = the ratio of adsorption to particle-induced
desorption rate
Di Toro found that vx was of order 1 over a broad range of chemical and solids types.
This formulation is included in IPX. The user may include the effect of particle concentration on
adsorption by using a vx value of order 1 (see Di Toro, 1985 for more detail); the effect may be
eliminated by specifying a large value for vx (the default value is 1012). If vx is specified to be
1.0, the model will predict a maximum paniculate fraction in the water column of 0.5 for all
hydrophobic chemicals (Kpso Ms > 10).
For each chemical modeled, up to 5 partition coefficients can be defined representing the
dissolved and four sorbents (DOC and three types of solids) states. Sorption/binding of the
chemical to solids and DOC is defined by the foc of the sorbent (assumed to be 1 for DOC but
modifiable through selection of other model parameters), the organic carbon partition coefficient
of the chemical (KoC), or the user defined relationship between KQW and Koc, and the particle
interaction parameter vx value.
The total chemical concentration is the sum of the five phase concentrations:
(1.51)
Substituting Equations 1.47 and 1.48 into Equation 1.51, factoring, and rearranging terms, yields
the dissolved fraction fa:
fd -'
1
(1-52)
Similarly, the particulate (sediment-sorbed) and DOC-bound fractions are:
fb =•
KBB
(1.53)
(1-54)
40
-------
These fractions are determined in time and space throughout a simulation from the partition
coefficients, internally calculated porosities, simulated sediment concentrations, and -specified
DOC concentrations.
Given the total concentration and the five phase fractions, the dissolved, sorbed, and
bound concentrations are uniquely determined:
(for each possible solids type "s")
cb =.fbc
(1.55)
(1.56)
(1.57)
These five concentrations have units of mg/L, and can be expressed as concentrations within
each phase:
n
C' -
CB~
B
(1.58)
(1.59)
(1.60)
These concentration expressions have units of mg/Lw, mg/kgs, and mg/kgB, respectively.
In some cases, such as near discharges, the user may have to alter input partition
coefficients to describe the effect of incomplete (non-equilibrium) partitioning. As guidance,
Karickhoff and Morris (1985) found that typical sorption reaction times are related to the
partition coefficient:
(1.61)
where:
=• desorption rate constant, hr [T]
41
-------
Table 1.6. IPX partitioning/binding data.
Description
Suspended solids concentration
Sediment bed solids concentration
Dissolved organic carbon
Partition coefficient, phase i
Lumped metal distribution coefficient
Octanol-water partition coefficient
Organic carbon fraction, phase i
Particle interaction parameter
Notation
ms
MB
DOC,B
KDi
KD
•Kow
loci
Vx
Common
Range
10 - 100
0.5-2
0-10
10'1 - 105
lO0-^
10° - 106
0.005 - 0.5
1 - 101*
S.I.
Units
mg/L
kg/L
mg/L
L/kg
L/kg
-
_
-
Thus, compounds with high, medium, and low Kow values of 105, 103, and 10 sorbing
onto 2% organic sediment should have reaction times of 1 day, 1/2 hour, and seconds. Given
that the time to equilibrium is roughly three times the reaction time, the three compounds should
reach equilibrium within 3 days, 1 hour, and 30 minutes. IPX data specifications for
partitioning/binding are summarized in Table 1.6.
1.6.2 Kinetic Transformations
The various phases of a chemical in water column and sediment segments are subject to
several transformation processes. Several variables may influence each process, leading to a
multi-term, and often non-linear, lumped transformation rate. If a single process is dominant in a
homogeneous aquatic system, then a single rate constant may be sufficient to describe the kinetic
reaction:
o __ y /~* s-\ ^r\\
•Jfe — — AteC (1.62)
where: Sfcc = total kinetic transformation rate for chemical c g/m3-day [M/L3/T]
Kkc = first-order rate constant for process k, day-1 [1/T]
C = total concentration of chemical, mg/L (g/m3) [M/L3]
Reaction can also be specified by a half-life. If a half-life is entered, then it will be converted to
a rate constant:
0.693
(1.63)
42
-------
If multiple rate constants are entered, they are summed:
(1.64)
For nonhomogeneous aquatic systems where rates vary in space, the user may specify a spatially
variable, lumped first-order rate constant KTC(X), so that:
Skc=-KTc(x}C
(1.65)
For nonhomogeneous aquatic systems where rates may vary in space and time, or for
cases where rate constants are unknown or cannot be calibrated, IPX uses the strategy
implemented in the original Exposure Analysis Modeling System (Burns et al. 1982). Each
process is considered separately using mixed second-order kinetics:
c+teL
(1.66)
where: [E]k = the intensity of environmental property "E" that affects process "k", such as
light intensity or bacterial population
= transformation product for process k acting on chemical c
The corresponding reaction rate Skc in mg/L-day for process k acting on chemical c is:
S!K=kkc[E]kYkcC (1.67)
where: kkc
= second-order rate constant for process k on chemical c
= yield coefficient for production of chemical from process k acting on chemical
c; assumed to be -1 for the production of chemical c by itself
Given a local value for [E]k, a pseudo-first-order rate coefficient Kkc in day1 can be specified:
K -k \E\ Y • ' (1.68)
*- ~ K£!'J[ \i.«o^
For a compound undergoing several reactions, the lumped transformation reaction is:
(1.69)
k c
43
-------
Table 1.7. IPX general kinetic data.
Description
First-order rate constant for process "k"
Half-life for process "k"
Lumped first-order rate constant
Chemical solubility
Water temperature
Notation
Kk
tHk
KT
Sol
T
Range
0-10
0.07 - oo
0-10
10"6 - 106
4-30
Units
day1
days
day'1
mg/L
°C
The local first-order assumption is generally accepted to be accurate for most chemicals
at environmental concentrations. The assumption is invalid at concentrations near the solubility
limit, however. If the user does not specify a maximum concentration CMAX(l), IPX sets this
limit at half the solubility or 10"5 molar, whichever is less, and aborts the simulation if
concentrations exceed this value.
The individual transformation processes considered by IPX are hydrolysis, photolysis,
oxidation, and biodegradation. In addition, volatilization is calculated and added to the
transformation rate. Good discussions of these processes are presented in Smith et al. (1977),
Burns et al. (1982), Mill et al. (1982), Mabey et al. (1982), Mills et al. (1985), and others. The
following sections summarize how IPX calculates the local rate constant for each of these
processes. Input data requirements are given for each process. The general kinetic data required
by IPX are summarized in Table 1.7.
1.6.3 Hydrolysis
Hydrolysis, the reaction of a chemical with water, is a major degradation pathway for
many toxic organic contaminants. An example reaction is shown in Figure 1.6. The reaction can
be catalyzed by hydrogen ions or proceed by consuming hydroxide ions. Figure 1.7. illustrates
the effects of pH on the base hydrolysis of carbaryl, the neutral hydrolysis of chloromethane, and
the acid and base hydrolysis of 2,4-D.
Li IPX, hydrolysis by specific acid, neutral, or base-catalyzed pathways is considered for
the various phases of each chemical:
(1.70)
(1.71)
44
-------
Neutral
Acid-
Catalysis
Base-
Catalysis
Hydrolysis
C = H2O
C = H2O
C = H2O —
OH'
-*• P + P'
-*• P + P'
-*> P + P1
Example
carbary!
H2NCH3 + CO2
naphthanol + methylamine + carbon
dioxide
Figure 1.6. Hydrolysis reactions.
(1.72)
where: KHN
KHH
KHOH
kr,
neutral hydrolysis rate constant, day1 [1/T]
acid catalyzed hydrolysis rate constant, day1 [1/T]
base catalyzed hydrolysis rate constant, day1 [1/T]
specific acid catalyzed (a) and base (b) rate constants for chemical in phase j,
molar1 day1 [L3/moles/T]
neutral rate constant for chemical in phase j, day1 [1/T]
fraction of chemical in phase j
IPX hydrolysis data specifications are summarized in Table 1.8. The reaction coefficients
can be specified as constants, with activation energy constants left as 0. If non-zero activation
energies are specified, IPX can compute the rates based on the temperature-based Arrhenius
function for each rate constant kn as follows:
(1.73)
45
-------
1 —
o —
*J -i-ir.iiL
CD
T5
Carbaryl
D Chloromethane
A 2,4-D (2-butoxyethyl ester)
6
pH
Figure 1.7. pH dependence of example hydrolysis reactions.
where: Tk = water temperature, K
TR = reference temperature for which reaction rate is reported, K
EaH = Arrhenius activation energy for hydrolysis reaction, kcal/mole
R = 1.99cal/moleK
1000 = cal/kcal
Activation energies may be specified for each hydrolysis reaction (acid, neutral, base)
simulated. If activation energies are input, reaction rates are adjusted to ambient water
temperatures. If no activation energies are specified, then rate constants are not adjusted for
temperature.
46
-------
Table 1.8. IPX Hydrolysis Data
Description
Negative log of hydrogen ion activity [H+J
Acid hydrolysis rate constant chemical in
phase j
Neutral hydrolysis rate constant for
chemical in phase j
Base hydrolysis rate constant for chemical
in phase j
Water temperature
Activation energy for hydrolysis reaction
for chemical
Notation
pH
kfiAj
kHNi
kHBj
T
EaH
Range
5-9
0-107
0-102
0-107
4-30
15-25
Units
kcal
mole[#+]°C
day1
kcal
mole[0#-]°C
°C
kcal
mole°C
1.6.4 Photolysis
Photolysis is the transformation of a chemical due to absorption of light energy. An
example of several photochemical pathways is given in Figure 1.8. The first-order photolysis
rate at the water surface can be calculated as:
(1.74)
where:
4>s
Lk =
first-order photolysis rate coefficient at the water surface, day1 [1/T]
reaction quantum yield for chemical, mole/E
(E is Einstein, a molar unit for light)
extinction coefficient of wavelength k by chemical, 1/M/cm
rj
average light intensity of wavelength k, mE/cm^/day
Usually, kps (or, preferrably £k and (J>) is treated as a reference rate in the literature and chemical
databases. Adjustment to the water body and conditions of interest is made in the computation of
the depth-integrated photolysis rate constant, kp.
, The rate of photolysis over the depth of the water column depends upon the attenuation of
light in water. For many organic chemicals, light absorption occurs primarily in the UV-B
wavelengths (280-320 nm), which are strongly attenuated by water as well as dissolved organic
matter. As a result, essentially all photolysis takes place in the upper meter of the water column.
47
-------
Internal
conversion
hv
A0 + heat
Internal
conversion
Absorption
Intersystem crossing
A0 + hv*
Phosphorescence
>
Quenching
"On
Chemical reaction
Chemical
reaction
A0 - ground state of reactant molecule
A* - excited state
Q0 - ground state of quenching molecule
Q* - excited state
Figure 1.8. Photolysis reactions.
It may then be assumed that all of the light is absorbed in the surficial water column segments,
for which the depth-integrated photolysis rate constant can be calculated as:
a,.
(1.75)
where: kp
= depth-integrated photolysis rate, m/day [L/T]
= light attenuation coefficient at wavelength k, nr1 [1/L]
For UV-B wavelengths, representative measurements of cc^ are 1/m for large water bodies such
as the Great Lakes (Calbins, 1975; personal communication, T. Mill, SRI), and 10/m for rivers in
the southeastern U.S. (Zepp and Cline, 1977).
48
-------
Adjustment of kp for the difference in latitude between the water body of interest and the
reference location is made by a latitude correction factor:
.63c
-------
Table 1.9. Average light intensity Lk of wavelength k at 40°N, mE/cm2/day.
Wavelength (nm)
297.5
300
302.5
305
307.5
310
312.5
315
317.5
320
323.1
330
340
350
360
370
380
390
400
410
420
Spring
1.85E-05
1.06E-04
3.99E-04
1.09E-03
2.34E-03
4.17E-03
6.51E-03
9.18E-03
1.20E-02
1.48E-02
2.71E-02
9.59E-02
1.23E-01
1.37E-01 .
1.52E-01
1.63E-01
1.74E-01
1.64E-01
2.36E-01
3.10E-01
3.19E-01
Summer
6.17E-05
2.69E-04
8.30E-04
1.95E-03
3.74E-03
6.17E-03
9.07E-03
1.22E-02
1.55E-02
1.87E-02
3.35E-02
1.16E-01
1.46E-01
1.62E-01
1.79E-01
1.91E-01
2.04E-01
1.93E-01
2.76E-01
3.64E-01
3.74E-01
Fall
7.83E-06
4.76E-05
1.89E-04
5.40E-04
1.19E-03
2.19E-03
3.47E-03
4.97E-03
6.57E-03
8.18E-03
1.51E-02
5.44E-02
7.09E-02
8.04E-02
9.02E-02
9.77E-02
1.05E-01
9.86E-02
1.42E-01
1.87E-01
1.93E-01
Winter
5.49E-07
5.13E-06
3.02E-05
1.19E-04
3.38E-04
7.53E-04
1.39E-03
2.22E-03
3.19E-03
4.23E-03
8.25E-03
3:16E-02
4.31E-02
4.98E-02
5.68E-02
6.22E-02
6.78E-02
6.33E-02
9.11E-02
1.20E-01
1.24E-01
Table 1.10. IPX photolysis data.
Description
Depth-integrated chemical photolysis rate
Reaction quantum yield fraction for
chemical in phase j
Notation
kp
Fj
Range
0-10
0-1
Units
m/day
50
-------
A much more sophisticated and elaborate algorithm, based upon EXAMS H (Burns and
Cline, 19.85), was used to calculate near-surface and depth-integrated photolysis rates in the
WASP4 (TOXI4) framework. Substitution of the algorithm described in this section, was
motivated by the desire for a simpler, more readily verified method to estimate photolysis rates
than the EXAMS-based code found in the early version of WASP4 that served as the foundation
for IPX. The authors hope to restore this algorithm in a later version of the IPX framework.
1.6.5 Oxidation
Chemical oxidation of chemicals in aquatic systems can be caused by interactions
between free radicals and the pollutant. Free radicals can be formed as a result of photochemical
reactions. Free radicals that have received some attention in the literature include alkylperoxy
radicals, RC>2, hydroxide radicals, OH", and singlet oxygen, O.
In IPX, oxidation is modeled as a general second-order process for the various phases of
each chemical:
(1.79)
where: Ko =
[R02] =
k0j =
net oxidation rate constant, day1 [1/T]
molar concentration of oxidant, moles/L [M/L3]
second-order oxidation rate constant for chemical in phase j, L/mole-day
[L3/M/T]
The reaction coefficients may be specified as constants, with activation energy constants
left (default value) as 0. If non-zero activation energies are specified, IPX can compute the rates
based on the temperature-based Arrhenius function for each rate constant k0 as follows:
k(Tt) = k(TR
where: Eao = Arrhenius activation energy for oxidation reaction, kcal/mole
(1.80)
Activation energies may be specified for each chemical simulated. If no activation
energies are given, then rate constants will not be adjusted to ambient water temperatures.
Because of the large number of alkylperoxy radicals that potentially exist in the
environment, it would be impossible to obtain estimates of k0 for each species. Mill et al. (1982)
propose estimation of a rate coefficient using t-butyl hydroperoxide as a typical oxidizing agent.
51
-------
Table 1.11. IPX oxidation data.
Description
Oxidation rate constant for chemical phase j
Activation energy for oxidation of chemical
Water temperature
Concentration of oxidants
Notation
Ko
Eao
T
[R02]
Range
15-25
4-30
io-lx -10-*
Units
L/mole-day
kcal/mole
°C
moles/L
Mill et al. argue that other alkylperoxides exhibit similar reactivities to within an order of
magnitude. The second-order rate coefficients are specified in IPX.as constants.
In addition to estimating a rate coefficient, estimates of free radical concentrations must
be made to completely define the expression for free radical oxidation. Mill et al. (1982)
reported RC>2 concentrations on the order of 10~9 M and OH' concentrations on the order of 10'17
M for a limited number of water bodies. Zepp and Cline (1977) reported an average value on the
order of 10'12 M for singlet oxygen in water bodies sampled. The source of free radicals in
natural waters is photolysis of naturally occurring organic molecules. If a water body is turbid or
very deep, free radicals are likely to be generated only near the air-water interface, and
consequently, chemical oxidation will be relatively less important. In such cases, the
concentrations cited above are appropriate in only the near-surface zones of water bodies. The
molar oxidant concentrations are input to IPX using the parameter OXRADG(ISEG). IPX
oxidation data specifications are summarized in Table 1.11.
1.6.6 Bacterial Degradation
Bacterial degradation, sometimes referred to as microbial transformation, biodegradation,
or biolysis, is the breakdown of a compound by the enzyme systems in bacteria. Examples are
given in Figure 1.9. Although these transformations can detoxify and mineralize toxins and
defuse potential toxins, they can also activate potential toxins.
Two general types of biodegradation are recognized: growth metabolism and co-
metabolism. Growth metabolism occurs when the organic compound serves as a food source for
the bacteria. Adaptation times from 2 to 20 days were suggested by Mills et al. (1985).
Adaptation may not be required for some chemicals or in chronically exposed environments.
Adaptation times may be lengthy in environments with a low initial density of degraders (Mills et
al. 1985). For situations where biodegradation is limited by the population size of the
52
-------
(Potential Toxin)
O(CH2)3COOH
C!
(Less Toxic Substance)
OH
O CH2 CH2 OSO3 H
Cl
Cl
(Potential Toxin)
O CH2 COOH
Cl
Figure 1.9. Example microbial transformations (from Alexander, 1980).
microorganism responsible for degradation, adaptation is faster for high initial microbial
populations and slower for low initial populations. Following adaptation, biodegradation
proceeds at fast, first-order rates. Co-metabolism occurs when the organic compound is not a
food source for the bacteria. Although adaptation is seldom necessary, transformation rates are
slow compared with growth metabolism.
In IPX, first-order biodegradation rate constants or half-life for the water column and the
benthos may be specified. If these rate constants have been measured under similar conditions,
this first-order approach is likely to be as accurate as more complicated approaches. If first-order
rates are unavailable, or if they must be extrapolated to different bacterial conditions, then the
second-order approach may be used. It is assumed that bacterial populations are unaffected by
the presence of the compound at low concentrations. Second-order kinetics for chemical in the
water column and the bed are considered:
H4/J (L81)
«/, (1-82)
53
-------
Table 1.12. IPX Biodegradation Data.
Description
Observed first-order degradation rate in
water column
Observed first-order degradation rate in
benthos
Bacterial activity or concentration of
bacterial agent
Observed second-order rate coefficients for
chemical in water and benthos
Biodegradation temperature coefficients
for chemical in phase j in water and
benthos
Water temperature
Notation
KBW
KBs
Pbac
kfiwj
kfisj
Qxwj
QTSJ
T
Common
Range
0-0.5
0-0.5
102 - 107
0 - 10'b
1.5 - 2.5
4-30
Units
day1
day'1
cells/mL
mL/cell-day
°C
where: KBW = net biodegradation rate constant in water segment, day1
KBS = net biodegradation rate constant in sediment (benthic) segment, day1
kBwj = second-order biodegradation rate constant for chemical in phase j in water
segments, ml/cell-day
kfisj = second-order biodegradation rate constant for chemical in phase j in benthic
segments, mL/cell-day
Pbac(t) = active bacterial population density in segment, cell/ml
fj = fraction of chemical in phase j
IPX biodegradation data specifications are summarized in Table 1.12. The second-order
rate constants for water and for sediment segments can be specified as constants. Temperature
correction factors can be zero. If non-zero temperature correction factors are specified, IPX can
modify each rate constant RB as follows:
•Bwj
Bsj
(T~20)/l°
(1.83)
(1.84)
where: Qxwj = "Q-10" temperature correction factor for biodegradation of chemical in phase j
in water
54
-------
Table 1.13. Size of Typical Bacterial Populations in Natural Waters.
Water Body Type
Oligotrophic Lake
Mesotrophic Lake
Eutrophic Lake
Eutrophic Reservoir
Dystrophic Lake
Lake Surficial Sediments
40 Surface Waters
Stream Sediments
Rur River (Winter)
Bacterial Numbers (cells/ml)
50 - 300
450 - 1,400
2000 - 12,000
1000 - 58,000
400 - 2,300
8 x 109 - 5 x 1010 cells/100 g dry wt
500 - 1 x 106
107 - 10s cells/100 g dry wt
3xl04
Reference
a
a
a
a
a
a
b
c
d
a) Wetzel (1975); enumeration techniques unclear.
b) Paris et al. (1981); bacterial enumeration using plate counts.
c) Herbes & Schwall (1978); bacterial enumeration using plate counts.
d) Larson et al. (1981): bacterial enumeration using plate counts.
QTSJ = "Q-10" temperature correction factor for biodegradation of chemical in phase j
in benthic segments
T = ambient temperature in segment, °C
The temperature correction factors represent the change in biodegradation rate constants resulting
from a 10°C temperature increase. Values in the range of 1.5 to 2 are common.
Total bacterial populations. for water and benthic segments are input using parameter
BACTOG(ISEG). Typical population size ranges are given in Table 1.13. Time functions that
multiply the water and benthic segment populations are input using functions BACNW and
BACNS. The product of parameter BACTOG and these time functions gives a description of the
time- and space-variation of bacterial populations. If the time functions are omitted, populations
remain constant in time.
Environmental factors other than temperature and population size can limit bacterial
rates. Potential reduction factors must be considered externally by the user. Nutrient limitation
can be important in oligotrophic environments. The following reduction factor was used by
Ward and Brock (1976) to describe phosphate limitation of hydrocarbon degradation:
55
-------
fpot —
0.0277C
PO.
1+0.0277C,
(1.85)
where: Cpo4 = dissolved inorganic phosphorus concentration, ug/L
This adjustment must be made to the input rate constants by the user for situations of
nutrient limitation. Low concentrations of dissolved oxygen can also cause reductions in
biodegradation rates. Below DO concentrations of about 1 mg/L, rates decrease. When anoxic
conditions prevail, most organic substances biodegrade more slowly. Because biodegradation
reactions are generally more difficult to predict than physical and chemical reactions, site-
specific calibration may be important.
1.6.7 Volatilization
Volatilization is the gradient-driven movement of a chemical across the air-water
interface. Other mass transfer processes between the air and water (dry and wet particle
deposition and rain dissolution) are treated as loads in IPX. The dissolved chemical attempts to
equilibrate with the gas phase partial pressure, as illustrated in Figure 1.10. The equation in this
figure shows that equilibrium occurs when the dissolved concentration equals the partial pressure
divided by Henry's Law Constant. Volatile organic chemical concentrations in the atmosphere
are often much lower than partial pressures equilibrated with water concentrations. In such
cases, volatilization reduces to a first-order process with a rate proportional to the conductivity
and surface area divided by volume:
— s —
-
(1.86)
where: Kv = net volatilization rate constant, day1
kv = conductivity of the chemical through water, m/day
As = surface area of water segment, m2
V = volume of the water segment, m3
D = average depth of the segment, m
fd = dissolved fraction of the chemical
The value of kv, the conductivity, depends on the intensity of turbulence in the water
body and the overlying atmosphere. Mackay and Leinonen (1975) discussed conditions under
which the value of kv is primarily determined by the intensity of turbulence in the water. As the
Henry's Law constant increases, the conductivity tends to be increasingly influenced by the
56
-------
air
dt
°a \
H/RT/
Cw -
H =
R =
I* =
D =
dissolved concentration in water, ug/L
concentration in air, ug/L
Henry's law constant, atm/M
gas constant (8.206 x 10-5), atm/M-K
water temperature, K
depth, m
= rate constant, m/day (conductivity)
Figure 1.10. Volatilization.
intensity of water turbulence. As the Henry's Law constant decreases, the conductivity tends to
be increasingly influenced by the intensity of atmospheric turbulence.
Henry's Law constants generally increase with increasing vapor pressure and decrease
with increasing solubility of a compound. Highly volatile, low solubility compounds are most
likely to exhibit mass transfer limitations in the water (liquid phase resistance). Relatively
nonvolatile, high solubility compounds are more likely to exhibit mass transfer limitations in the
air (gas phase resistance). Volatilization rates are usually smaller in comparatively quiescent
lakes and reservoirs than in relatively turbulent rivers and streams.
57
-------
In cases where the volatilization rate is regulated by turbulence in the water phase,
estimates of volatilization can be obtained from results of laboratory experiments. As discussed
by Mill et al. (1982), small flasks containing a solution of a pesticide dissolved in water that have
been stripped of oxygen can be shaken for specified periods of time. The amount of pollutant
lost and oxygen gained through volatilization can be measured and the ratio of conductivities
(KVOG) for pollutants and oxygen can be calculated. As shown by Tsivoglou and Wallace
(1972), this ratio should be constant and independent of the turbulence in the water body.
Therefore, if the reaeration coefficient for a receiving water body is known or estimated, and the
ratio of the conductivity for the pollutant to the reaeration coefficient has been determined, the
'pollutant conductivity can be directly estimated.
In IPX, the dissolved concentration of a compound in a surface water column segment
can volatilize at a rate determined by the two-layer resistance model (Whitman, 1923), where the
conductivity is the reciprocal of the total resistance:
(1.87)
where: RL
RG
KL
KG
= liquid phase resistance, day/m
= gas phase resistance, day/m
= liquid phase transfer coefficient, m/day
= gas phase transfer coefficient, m/day
The two-resistance model assumes that two "stagnant films" are bounded on either side
by well mixed compartments. Concentration differences are the driving force for the water layer
diffusion. Differences in partial pressures drive the diffusion for the air layer. At equilibrium,
the same mass must pass through both films, and the resistances combine in series. Another
resistance, the transport resistance between the two interfaces, is assumed to be negligible.
However, this may not be true in two cases: very turbulent conditions and in the presence of
surface active contaminants. Although this two-resistance method, the Whitman model, is rather
simplified in its assumption of uniform layers, it has been shown to be as accurate as more
complex models (for example, surface renewal and penetration theories). Laboratory studies of
volatilization of organic chemicals confirm the validity of the method as an accurate predictive
tool (Burns et al. 1982).
IPX allows the user considerable flexibility in specification of the volatilization rate. The
volatilization rate may be input directly or calculated from a combination of semi-theoretical
58
-------
relationships for the liquid and gas phase mass transfer rates, KL and KG, that follow. These
options are selected by specification of the constants XKL and XKG.
Several formulations are available for computing liquid phase mass transfer rates. For
segments with depths less than 0.61 m the Owens formula can be used to calculate an oxygen
reaeration rate, which is then converted to a chemical-specific liquid phase mass transfer rate:
u
0.67
KL= 5.349 .„
*- T-\O.S5
(MW02
\0.5
(1.88)
where: u =
D
MW02 =
MWC =
velocity of the water, m/s
segment depth, m
molecular weight of diatomic oxygen, g/mole = 32.
molecular weight of chemical, g/mole
For segments where the velocity is less than 0.52 m/s or the depth is greater than 13.584
U2.9135 (m)5 a modified form of the O'Connor-Dobbins formulation can be used to compute a
chemical-specific liquid phase mass transfer rate:
D
N0.5
w02
D
0.5
, ^2. 8.64xl04
(MWC }
(MWO
(1-89)
where: Dwo = diffusivity of oxygen in water, cm2/s
For other cases, the Churchill formulation can be used to calculate a liquid phase mass
transfer rate:
^,=5.049 .673
M°-969 f MW02
D(
\0.5
(1.90)
For cases where mass transfer is controlled by the wind velocity, typical of lakes and
reservoirs, formulations presented by O'Connor (1983b) and Mackay and Paterson (1986) can be
used. Two forms of the O'Connor liquid phase mass transfer formulation are available. The
O'Connor "short" form, which is applicable for hydrodynamically-smooth (viscous) conditions
is: . ' •
59
-------
\
limiting roughness length, cm = 0.35
2.5 '
transitional shear velocity, cm/s =10
(1.92)
(1.94)
60
-------
The shear velocity for this formulation is computed iteratively by successive substitution, using the
relationship which expresses the drag coefficient (CdoL) in terms of the roughness height (Zo):
10
In
/•
1000
1 , \&1U*OL
z* L va J
-Mil
e\**)
(1.95)
The values of the semi-empirical coefficients u*oLC, TO, Ze, A,i, and u#oLt are those suggested by
O'Connor to be applicable to large water bodies (oceans or large lakes).
The Mackay formulation is:
KL =
KL = [lO'6
>0.3m/s
<0.3m/s
where:
= shear velocity, m/s = 0.01 Ww (6.1+0.63 Ww)'
0.5
(1.96)
If the rate of oxygen reaeration for the water body is known (for instance, if reaeration
studies have been conducted for dissolved oxygen modeling), then a chemical-specific liquid
phase mass transfer rate may be calculated as:
MW,
N0.5
02
(1.97)
where: KLQ, = oxygen reaeration rate, m/day
kvo = ratio of chemical volatilization to reaeration conductivities
This rate can be modified in space and time through proper specification of model segment
parameters and time functions for reaeration.
Several formulations are also available for computing gas phase mass transfer rates. The
O'Connor-Rathbun formulation can be used to compute a chemical-specific gas transfer rate
(O'Connor, 1983a; Rathbun, 1990):
KC =
aoQi
X<>.667
'w
8.64xl04
(1.98)
61
-------
where: Sca = Schmidt Number for air =
Da = diffusivity of chemical in air, cm2/s
va = kinematic viscosity of air, cm2/s
Mill et al. (1982) found that the gas phase mass transfer rate for a chemical can be
estimated from the evaporation rate of water:
-leg
\0.25
MWC
w
10
(1.99)
where: MWW = molecular weight of water, g/mole =18
The O'Connor gas phase mass transfer formulation (O'Connor, 1983b) is:
_0l33
K,
[,-0.33
„ E.
"*- ^
8.64x10"
(1.100)
The Mackay gas phase mass transfer formulation (Mackay and Paterson, 1985)js:
#c=[w.m0.0462Sc^66+Kr3]8.64xl04 (1.101)
For flowing systems, Ambrose et al. (1988) assume that the gas phase mass transfer rate
is constant at 100 m/day. As wind speeds decrease, mass transfer rates in the gas phase also
decrease. Under still-air conditions, Mackay and Leinonen (1975) suggested that KG is 86.4
m/day (0.001 m/s), representing a minimum rate of gas phase mass transfer. The user has the
option of specifying gas phase mass transfer rates below this value, however.
Chemical diffusivities in air (Da) and water (Dw) are calculated internally by IPX. At a
minimum, the chemical molecular weight must be input to allow calculation of diffusivities. If
the LeBas molar volume (Vb) is also input, then diffusivities will be calculated as functions of
temperature. The calculation of Vb based upon molecular structure, is presented in Lyman et al.
(1982).
The Henry's Law constant, HLC or H, of a chemical is defined as the ratio of the vapor
pressure, Pv to solubility, S:
62
-------
(1.102)
where: H = Henry's Law constant (HLC), arm m3/mol
Pv = vapor pressure, atm
S = solubility, mg/L
The Henry's Law constant for a chemical may be input directly, or calculated from vapor
pressure and solubility data. Mackay and Shiu (1981) and Schnoor et al. (1987) provide
extensive reviews of available HLC data.
Henry's Law constants vary with water temperature. IPX provides several options to
handle this temperature dependence. For many chemicals, the HLC can be directly expressed as
a function of temperature in the following form:
InH = an —
a.
(1.103)
T + 273.15
where: ao, ai = user input HLC temperature function factors
T -= water temperature, °C
For example, the Henry's Law constants for polychlorinated biphenyls (PCBs) (Tateya et al.
1988) and mirex (Yin and Hassett, 1986) expressed as functions of temperature are:
PCBs:
mirex:
InH =18.53 —
7868
InH = 29.27-
T + 273.15
10849
T + 273.15
Alternatively, Henry's Law constant may be input as a chemical constant, or calculated by
IPX according to Equation 1.102. For chemicals which are solids at ambient temperature,
calculations of Henry's Law constant using Equation 1.102 require subcooled liquid vapor
pressure and solubility. IPX can convert solid-phase input of Pv and S to subcooled liquid
values, using an algorithm described by McLachlan and Mackay (1987). This conversion is only
performed if the melting temperature is below ambient temperature. The conversion from solid
(Ps) to subcooled liquid vapor pressure (PL) is:
63
-------
(1.104)
where: As = Entropy of fusion, 56 J/mol K
R = Gas constant, 8.3 14 J/mol K
TM = Melting temperature, K
T = Temperature at which Ps is valid, K
PL is then corrected for water temperature using the Kistiakowski Linear Enthalpy equation:
(1.105)
where: PL
PB
TB
T
= Subcooled liquid vapor pressure at temperature T
= Vapor pressure at boiling point ( 1 atm)
= Boiling temperature (K) at 1 atm
= Water temperature (K)
Iterative solution of Equation 1.105 is performed to calculate the chemical's boiling temperature
using the subcooled liquid vapor pressure and its reference temperature (assumed to be 293.15
K). Then, PL is adjusted for temperature using Equation 1.105, and HLC is calculated using
Equation 1.102. Solubility is not corrected for temperature, because S is nearly constant over the
range of 0 to 25°C (at least in comparison to the variation in vapor pressure).
Alternatively, if the overall volatilization rate is computed for a temperature of 20°C, this
rate can be adjusted for temperature as follows:
r-20
(1.106)
where: 0 = user-specified temperature correction factor
T = water temperature, °C
This would be an appropriate option if HLC was input as a chemical constant.
Although many calculations may be necessary to compute chemical volatilization rates,
most are performed internally and require users to specify relatively few parameters. The number
of constants required depends on the options selected. For example, if the Henry's Law constant
64
-------
Table 1.14. IPX volatilization data.
Description
Measured or calibrated conductance
Henry's Law Constant
Vapor pressure
Solubility
LeBas molar volume
Melting temperature
Concentration of chemical in the
atmosphere
Molecular weight of chemical
Oxygen Reaeration coefficient
Experimentally measured ratio of
volatilization to reaeration
Current velocity
Water depth
Water temperature
Wind speed 10 m above water surface
Notation
kv
H
Pv
S
vb
TM
Ca
MWC
KLo2
kyo
U
D
T
Wio
Range
0-25
10'7 - 10'1
0 - 1000
10 - 103
0.6 - 25
0-1
0-2
4-30
0-20
Units
m/day
atmrn3
mol
atm
mg/L
cm3/mol
°C
vgfl
g/mol
m/day
m/sec
m
°C
m/sec
for a solid-phase chemical is calculated internally from vapor pressure and solubility, then
melting temperature becomes a required input. If kvo is not measured, it will be calculated
internally from molecular weight. IPX volatilization data specifications are summarized in Table
1.14.
1.6.8 Extra Reaction
IPX allows the user to' specify an additional second-order reaction for the various phases
of each chemical:
,/,. (1-107)
where: KE =
kej =
net extra reaction rate constant, day1
intensity of environmental property driving the reaction
second-order rate constant for chemical in phase j, [E]'1 day1
fraction of chemical in phase j
65
-------
An example of a kinetic process that may be modeled by this extra reaction is reduction.
For a reduction reaction, [E] may be interpreted as the concentration of
environmental reducing agents RH2, so that:
C+RH2=*P (1.108)
where: C = concentration of chemical C, mole/L
RH2 = concentration of RH2, (= [E]) moles/L
P = reduced product, moles/L
The reducing agent and the second-order rate constant must be identified and quantified
by laboratory kinetics studies. If both the environmental oxidizing and reducing agents are in
excess, then two chemicals may be simulated as a redox pair:
Cl+RO2 <=>C2+RH2
(1.109)
where: Cj = reduced chemical
C2 = oxidized chemical
RO2 = oxidizing agent
RH2 = reducing agent
Laboratory kinetics studies can control the concentrations of RO2 and RH2 to determine
rate constants for both oxidation and reduction. These may be specified as constants k0 (for
oxidation) and kn (for the "extra" reaction). In this example, since there is interconversion
between the two chemicals being simulated, yield coefficients must also be specified. The
spatially variable concentrations [RO2] and [RH2] must be specified as parameters.
1.7 SIMULATING THE TRANSPORT AND FATE OF HEAVY METALS USING
IPX
Although designed explicitly for organic chemicals, IPX can be used to simulate metals
with judicious specification of model parameters. Heavy metals in the aquatic environment can
form soluble complexes with organic and inorganic ligands, sorb/bind with organic and inorganic
particulates, and precipitate or dissolve. Because of the inherent complexity of metals behavior,
site-specific calibration is required.
Although IPX does not compute ionic speciation, geochemical models such as
MINTEQA1 (Brown et al. 1987) can be used to predict metal speciation for a given set of
66
-------
environmental conditions. These calculations can be used to parameterize an IPX simulation.
All soluble complexes with the free ion are summed to give the dissolved metal concentration.
Precipitated metal phased can be lumped with all sorbed phases to give particulate "sorbed"
metal concentration. A spatially variable lumped partition coefficient Kp describes the two
phases. There is no general consistency in reported Kp values for particular metals in the natural
environment, so site-specific values will likely be necessary. Partition coefficients may depend
upon the sorbent character, including mineralogy, chemical structure, composition and electrical-
chemical properties, presence of coatings, and the age and origin of humic substances. Table
1.15 summarizes Kp values reported in Delos et al. (1984) for eight metals. These values are
generally high, and are provided as a starting point for simulation. Spatially-variable Kp values
can be input to IPX using the parameter FOC(ISEG,I) and omitting all other partitioning
parameters. '
1.8
SIMULATION COMPLEXITY LEVELS
The chemical transport and fate processes in IPX can be implemented at various levels of
complexity. These different levels involve increasing sophistication in solids behavior,
equilibrium partitioning, and kinetic reactions. Solids behavior may be modeled at four levels of
complexity: 1. descriptive solids concentration field, 2. descriptive solids concentration field with
specified solids transport rates, 3. simulated total solids, and 4. three simulated solids types.
Equilibrium partitioning may be modeled at five levels of complexity: 1. constant partition
coefficient, 2. spatially variable partition coefficients, 3. hydrophobic sorption, and 4. solids-
dependent partitioning. Kinetic reactions may be modeled at four levels of complexity: 1.
constant half-life or rate constants, 2. spatially variable rate constants, 3. second-order rates, and
4. second-order rates with transformation products. Each complexity level is discussed in detail
in the following section.
The increasing complexity levels permit more realistic description and extrapolation of
the solids, equilibrium, and kinetic processes. However, more complex representations of
processes generally require that users specify more model parameters. Although quicker and
easier to construct simulations for lower complexity, the corresponding decrease of model
capabilities requires more judgment on the part of the user. For example, consider analyzing a
problem where the dissolved phase of a transformation product is extremely toxic. Useful
simulations of this problem at a low complexity level would be difficult. For best results, the
user must match the complexity of the simulation with the requirements of the problem.
67
-------
Table 1.15. Speciation of priority metals between dissolved and particulate
phases as a function of suspended solids concentrations in streams.
Metal
Arsenic
Cadmium
Chromium
Copper
Lead
Mercury
Nickel
Zinc
SS (mg/L)
1
10
100
1000
1
10
100
1000
1
10
100
1000
1
10
100
1000
1
10
100
1000
1
10
100
1000
1
10
100
1000
1
10
100
1000 ,
Kp(L/kg)
5xl05
9xl04
2xl04
3 x 103
4xl06
3x10^
2 x 104
2xlOj
3 x 10°
4xl05
5xl04
5x10"
IxlO6
2xlOs
3 x 104
6xl03
3x10"
2x10'
IxlO5
9xl04
3xl06
2x10"
2xl04
IxlO3
5x10'
IxlO'
4 x 104
9xl03
IxlO6
2x10*
5xl04
IxlO4
% Dissolved
70
50
30
24
20
25
30
40
25
20
17
15
50
30
25
14
75
30
10
1
25
30
30
45
70
50
20
10
40
30
17
10
% Particulate
30
50
70
76
80
75
30
60
75
80
83
85
50
70
75
86
25
70
90
99
75
70
70
55
30
50
80
90
60
70
83
90
68
-------
As each simulation complexity level is described, specific IPX input variable names and
options are mentioned. Users should refer to Chapter 2 for more detailed explanations.
1.8.1 Simulating Solids
Level l~The simplest description of solids is to specify an average solids concentration
field. This is accomplished by setting the initial conditions for system 1 to observed
concentrations and setting SYSBY(l), SYSBY(2), and SYSBY(3) to 1. It is not necessary to
specify settling and resuspension velocities. Concentrations, of the specified solids will remain
constant but will influence chemical partitioning and, indirectly, transport and transformation.
Particulate phase chemicals will not settle or resuspend.
Level 2-The next level for solids is to specify an average concentration field for a single
solids type along with settling, resuspension velocities. As before, initial conditions for system 1
are specified as average solids concentrations and SYSBY(l), SYSBY(2), and SYSBY(3) set to
1. Solids transport velocities are specified in transport fields 3 and 6. Solids concentrations
remain constant but will directly influence chemical partitioning; solids transport velocities will
allow particulate chemical transport between the water column and sediments.
Level 3-The third level for solids is to simulate total solids. Loads, boundary
concentrations, and initial conditions are specified for system 1. Settling and resuspension
velocities are specified in transport fields 3 and 6. Solids concentrations can be calibrated to
observed data, leading to more accurate calculations of particulate chemical transport.
Particulate chemicals are transported between water and sediments.
Level 4~The fourth level for solids is to simulate multiple sediment types. Loads,
boundary concentrations, and initial conditions are specified for systems 2, 3, and 4. Solids
settling and resuspension velocities are specified in transport fields 3-6. Concentrations of each
solids type can be independently calibrated to observed data, leading to more accurate
calculations of solids transport, chemical partitioning, and chemical transport.
1.8.2 Simulating Equilibrium Partitioning
Level l--The simplest equilibrium partitioning is described by a single, constant partition
coefficient. This is done by specifying a value for PIXC(1,1) and omitting all other partitioning
, information such as LKOW, LKOC, and FOC(ISEG,1). Although the partition coefficient is
constant, the dissolved and sorbed chemical fractions vary with solids concentrations:
69
-------
(1.110)
fv =
l + K
(1.111)
where: K,
•P
MI
n
= partition coefficient, L/kg
= solids concentration, kg/L
= porosity
Porosity is calculated internally in IPX from input sediment density and solids concentration:
(1.112)
DSED(l)
where: DSED(l) = sediment density of solids type 1, specified under initial conditions, kg/L
Level 2~The next level for equilibrium partitioning is to specify spatially-variable
partition coefficients. This is done by specifying Kp values with the parameter FOC(ISEG,1) and
omitting all other partitioning information such as LKOW, LKOC, and PIXC(1,1). The
equations used by IPX are as described for Equilibrium Level 1. This option allows improved
site-specific calibration of partition coefficients to observed ratios between dissolved and
particulate chemical concentrations.
Level 3~The third level of complexity for equilibrium partitioning is spatially-variable
Kp values for hydrophobic sorption. This is done by specifying values for FOC(1,ISEG) and
LKOC or LKOW. If values for DOC(ISEG) are specified, then IPX will also calculate binding to
dissolved organic carbon. If concentrations for sediment types 2 and 3 are specified or simulated
along with corresponding FOG(2,ISEG) and FOC(3,ISEG) values, then IPX will calculate the
fraction of chemical in five phases:
j = 1, 2, 3
(1.113)
/
-------
J v ~
1 + KM fn + KDOCI n
KocDOC/n
.Mj/n + KocDOC/n
= l,2,3
= l,2,3
(1.115)
(1.116)
This level allows improved description of chemical partitioning, binding, and transport as
well as improved description of hydrophobic sorption for different sediment regimes.
Level 4~The next level of complexity adds solids-dependent partitioning to hydrophobic
sorption. This is done by specifying values for NUX(l) (Vx), the solids-dependent partitioning
constant. IPX calculates the partition coefficient as:
K
po
K "V f
OC / i J OCJ
where: Kp =
K,
-po
. Iv
the solids dependent partition coefficient, L/kg
= the solids independent partition constant, L/kg
(1.117)
If no values are specified for NUX(l) (vx), then TOXI4 assigns a large default value (i.e. no
particle effect). A value of between 1 to 10 is representative of many surface waters. However,
solids-dependent partitioning is not applied to sediment bed segments.
1.8.3 Simulating Kinetic Reactions
Level 1—The simplest kinetic reaction is described by a constant half-life or rate constant.
If the user supplies first-order decay constants KV, KBW, KBS, KHN, KHH, KHOH, KO, KF,
or KE for the transformation reactions, then they will be used directly:
:;c (i.iis)
where: C
Ki
KHN
KHH
KHOH
chemical concentration, mg/L
first-order decay constants, day1, including:
neutral hydrolysis constant, day1
acid-catalyzed hydrolysis constant, day1
base-catalyzed hydrolysis constant, day1
biodegradation constants for the water column and
sediments, day1
71
-------
KF = photolysis constant, day1
KO • = oxidation constant, day1
Ky = volatilization constant, day1
KE = extra reaction constant, day i
As an alternative, half-life may be are specified for any transformation process. If half-
life is specified for the transformation reactions, it is internally converted to first-order rate
constants:
0.693
THi
(1-119)
where: THI
THHN
THHH
THOH
THBW,
THF
THO
T1
IHV
THE
half-life, days, including:
neutral hydrolysis half-life, days
acid-catalyzed hydrolysis half-life, days
base-catalyzed hydrolysis half-life, days
biodegradation half-life for the water column and sediments, day1
photolysis half-life, days
oxidation half-life, days
volatilization half-life, days
extra reaction half-life, days
The converted half-life rate constants are then used as described by Equation 1.118.
Level 2—Because environmental conditions may change throughout a water body, decay
rates are expected to vary. The second kinetic level allows spatially variable decay rate constants
TOTKG(ICHM, ISEG) to be specified by the user so that:
Skc = -KTc(x)C
(1.120)
where:
= spatially variable lumped first-order decay rate constant for chemical c,
days"1
For segments where a non-zero KT is supplied for a chemical, IPX uses this value and
bypasses further kinetic calculations. For those segments where KT is zero or not specified, IPX
will apply any specified process rate constants.
72
-------
Level 3—Because environmental conditions may vary temporally as well as spatially,
empirically-determined lumped decay-rate constants may be inappropriate when extrapolated.
The third level of kinetic complexity calculates decay rates based on second-order kinetics so
that:
(1.121)
where:
kjkc
[E]k
= overall first-order rate constant for chemical c, day1
= second-order rate coefficient for phase j, process k of chemical c
= intensity of environmental property affecting process k
= fraction of chemical c in phase j
Users may implement any given reaction by specifying values for the rate constants (by phase),
relevant environmental parameters, and time functions.
Level 4-The fourth level of kinetic complexity allows simulation of transformation
products. This level is implemented by simulating two or more chemicals (NOSYS > 4) that
react and interconvert by specifying appropriate yield coefficients for each process:
Skcacm - ~
for n = 4,NOSYS, n # m (1.122)
cm k
for m = 4, NOSYS, m * n (1.123)
where: SkcnCm
Skcmcn
Kkcn
Kkcm
Ykcncm
production rate of chemical cn from chemical cm undergoing reaction k,
mg/L-day
production rate of chemical cm from chemical cn undergoing reaction k,
mg/L-day
effective rate coefficient for chemical cn, process k, day1
effective rate coefficient for chemical cm, process k, day1
yield coefficients for production of chemicals cn from chemical cm
undergoing reaction k, gm cn/gm cm
Ykcmcn = yield coefficients for production of chemicals cm from chemical cn
. undergoing reaction k, gm cm/gm cn
73
-------
Yield constants can be specified for every possible combination of reactant chemicals,
product chemicals, and reaction pathways.
1.9
SUMMARY OF DATA REQUIREMENTS
Chemical partitioning and transformation pathways add detail to the general mass
equations presented in Section 1.2. However, these additional processes require specification of
the environmental parameters, chemical constants, and environmental time functions as
discussed in the preceding sections. This section provides a summary of data requirements.
The environmental data required for a chemical simulation depend upon which
transformation processes are important. Table 1.16 lists the environmental properties influencing
each process in IPX, and a range of expected values. For a series of simulations involving many
compounds, approximate values for all environmental properties should be specified. For those
processes found to be most important, better estimates of the relevant environmental .properties
can be provided in a second round of simulations.
The chemical properties of each compound control which transformation processes are
important in a particular environment. Table 1.17 summarizes chemical properties influencing
each process in IPX. Although the model allows specification of different rates for the dissolved,
particulate, and DOC-bound chemical phases, such data are not generally available. Measured
rate constants are often assigned to the dissolved chemical phase. The model also allows
specification of temperature correction parameters for each process. Such data are often difficult
to find without special studies, and need not be input except for very hot or cold conditions, or
where seasonal variability is being studied.
Time variable functions can be used to study diurnal or seasonal effects on pollutant
behavior. The 21 environmental time functions are summarized in Table 1.18. As shown, some
of these time functions are multiplied by spatially variable parameters within TOXI4 to simulate
environmental conditions that vary in time and space. If no temporal variability is required, time
functions may be omitted. Default values for time functions is 1.0. However, unused functions
can affect a simulation. For example, if left unspecified the default water temperature is 0 and all
temperature-dependent reactions (such as volatilization) will be computed for 0°C (273.15 K).
Although the amount and variety of data potentially required for a simulation is large,
actual data requirements for any particular simulation are generally small. Usually only sorption
74
-------
Table 1.16. Environmental properties affecting interphase
transport and transformation processes.
Environmental Property
Solids Concentrations:
Suspended, mg/L
Sediment Bed, kg/L
Organic Carbon Fraction:
Suspended Sediment
Benthic Sediment
Dissolved Organic Carbon, mg/L
Water Column Depth, m
Water Temperature, °C
Average Water Velocity, m/s
Wind Speed at 10 m, in m/sec
pH, Standard Units
Concentration of Oxidants, moles/L
Surface Light Intensity, Langleys/day
Cloud Cover, tenths of sky
Light Extinction Coefficient, 1/m
Active Bacterial Population:
Suspended, cells/ml
Benthic, cells/lOOg
Input
Value
5-500
0.4 - 1.7
0.01-0.10
0.01-0.10
0- 10
0.5 - 500
4-30
0-2
0-20
5-9
10-9 . io-12
300 - 700
0-10
0.1-5
103 - 106
103 - 106
Environmental Process
Kp
(1)
X
X
X
X
X
KV
(2)
X
X
-x
X
KH
(3)
X
X
Ko
(4)
X
X
KF
(5)
X
X
X
X
KB
(6)
X
X
X
(1) Sorption; (2) Volatilization; (3) Hydrolysis; (4) Oxidation; (5) Photolysis; (6) Bacterial
Degradation
75
-------
Table 1.17. Chemical properties affecting interphase
transport and transformation processes.
Chemical Property
Molecular Weight
LeBas molar volume
Solubility
Vapor Pressure
Melting temperature
Octanol-Water Partition Coefficient
Organic Carbon Partition Coefficient
Henry's Law Constant
Liquid Phase Volatilization/
Reaeration Ratio
Alkaline Hydrolysis Rate Constant
Neutral Hydrolysis Rate Constant
Acid-Hydrolysis Rate Constant
Activation Energy for Alkaline,
Neutral, and Acid Hydrolysis
Oxidation Rate Constant
Activation Energy for Oxidation
Depth-Integrated Photolysis Rate
Measured Biodeg. Rate Constant
Water Column Rate Constant
Benthic Rate Constant
Temperature Dependence Multiplier
(for 10°C change)
Input Units
g/mole
cm^/mol
mg/L
atmospheres
°C
LW/LO
L
-------
Table 1.18.
(combinations
Time variable environmental forcing functions
of IPX time functions and constants/parameters).
Time Function
TEMPN(TMPFN)
VELN(VELFN)
WINDN
PHNW
PHNS
REARN
AIRTMPN
CHLN
PHTON
BACNW
BACNS
X
X
X
X
•X
X
X
X
X
X
Constant or
Parameter
TEMPOSEG)
VELOCG(ISEG)
WVEL(ISEG)
PH(ISEG)
PH(ISEG)
REAER(ISEG)
AIRTMP
CHPHL(ISEG)
BAC(ISEG)
BAC(ISEG)
=
=
=
=
=
**•"*
»
=
^:
™
=
Environmental Property
Water temperature (x,t), °C
Water velocity (x,t), m/sec
Wind speed at 10 m above the water
surface (x,t), m/sec
Water column pH (x,t), log activity
Benthic pH (x,t), log activity
Reaeration or volatilization rate (x,t), .
m/day
Air temperature (t), °C
Chlorophyll concentration, mg/L
Seasonal adjustment to light intensity
(dimensionless)
Water column bacteria population
(x,t), cells/mL
Sediment bacteria population (x,t),
cells/mL
and one or two transformation processes will significantly affect a particular chemical. To
simulate the transport of many soluble compounds in the water column, even sorption can often
be disregarded. For empirical screening-level studies most environmental parameters, chemical
constants, and time functions, can be ignored except the user-specified transformation rate
constant TOTKG(ICHM,ISEG), the partition coefficient LKOC, and organic carbon fraction
FOC(ISEGJ). Thus, IPX can be used as a first-order water pollutant modeling framework to
conduct standard simulations of dye tracers, salinity intrusion, or coliform die-off. What is
gained by the second-order process functions and resulting input data burden is the ability to
extrapolate more confidently to future conditions. The user must determine the optimum amount
of empirical calibration and process specification for each application.
1.10 SUMMARY OF MODEL EQUATIONS
The equations implemented in IPX account for the transport of dissolved and particulate
materials in the water column and sediments as summarized below:
77
-------
(1.124)
flow, precipitation, and evaporation
water column and pore water advection
n
n
water column and pore water dispersion
solids transport
L N B
point, non-point, and boundary loads
k c
kinetic transformations
where: j = segment index
= adjacent segment index
= solids transport field index
= point source index
= non-point source index
= boundary source index
= kinetic transformation index
= chemical index
Vj = volume of segment j,m3
Cj = concentration of the water quality constituent in segment j, g/m3
t = time, day
Ej = evaporation rate from segment j, m/day
Pj = precipitation rate into segment j, m/day
Aj = surface area of segment j, m2
l
s
L
N
B
k
c
78
-------
0)
w.
psj
advective flow between segments i and j, defined as positive when leaving
segment j, and negative when entering, m3/day
pore water flow between segments i and j, defined as positive when leaving
segment j, and negative when entering j, m3/day
constituent concentration advected between i and j, g/m3
to Cj + (1 -n) Q when entering j
CO Q + (1 -n) Cj when leaving j
advection factor (numerical weighing factor), 0 - 0.5
solids transport velocity between segments i and j, defined as positive when
leaving segment j, and negative when entering, m/day
dissolved fraction of chemical in segment j
fraction of chemical sorbed to solid type s in segment j
dispersive flow between segments i and j, m3/day
R,
•py
WLj
WNj
Qjo
= dispersion coefficient between segments i and j, m2/day
= cross-sectional area between segments i and j, m2
= characteristic mixing length between segments i and j, m
= pore water diffusive exchange flow, m3/day
Skcj =
Lcij Tij
average tortuosity of segments i and j, m water/m
average porosity of segments i and j, m3 water/m3
point source loads into segment j, g/day
non-point loads into segment j, g/day
boundary inflows to segment j, m3/day
boundary concentrations for segment j, g/m3
kinetic transformation k for chemical c within segment j, g/m3-day
Water column advection and dispersion parameters are specified for transport between
water column segments. Pore water advection and diffusion parameters are specified for
transport between sediment segments as well as between adjacent water column and sediment
segments. Sediment transport parameters representing settling and resuspension also link
adjoining water column and sediment segments.
79
-------
1.11 OTHER MODEL INPUTS
This section summarizes the model input parameters that must be specified in order to
uniquely define and solve the mass balance equation and conduct a simulation. Detailed
information regarding proper specification of all model inputs is presented in Chapter 2.
1.11.1 Model Identification and Control Parameters
These parameters identify the model input file. They include the number of water quality
constituents simulated and the number of segments in the network. This group of parameters
controls the simulation and checks the numerical stability of the solution. Simulation parameters
include the initial and final times, integration time-steps, the advection factor, maximum
concentrations, and a negative solution option. Also provided are text lines used to describe the
water body simulated.
Initial simulation time. days-The time at the beginning of the simulation must be
specified in order to synchronize all the time functions. The day, hour, and minute can be input.
A simulation typically begins on day 1 (time-break zero).
Final simulation time. days-The desired time at the end of the simulation must be
specified. The final time is input in fractional form (i.e. including decimal fraction, if any). For
example, a one year simulation beginning on day one (time zero) would end at 365.0.
Integration time-step. days-A sequence of integration time-steps (At) must be specified,
along with the time interval over which they apply. Given specific network and transport
parameters, time-steps are constrained within a specific range to maintain stability and minimize
numerical dispersion (solution inaccuracies). To maintain numerical stability, the advected,
dispersed, and transformed mass of a segment at any time-step must be less than the resident
mass of that segment at the beginning of the time-step. Mathematically, this can be expressed as:
(L&CJ +5>AC; + ^SV})At < VjCj (1.126)
Rearranging Equation 1.118 to solve for At and applying the criterion over the entire
network gives the maximum stable time-step size:
=Min
V,
(1.127)
80
-------
. If reactions are linear, the last term in the denominator reduces to I^KVj. However, most
tributaries are advectively dominated and At will generally be controlled by advective flows.
Advection factor, dimensionless-IPX uses a finite difference approximation of the
general mass balance differential equation (see Appendix A). The advection factor co (described
in Appendix A) is used to modify the finite difference approximation used by IPX. The value of
this factor influences the stability of the solution at the expense of numerical accuracy (numerical
dispersion). For €0 = 0, a backward difference approximation is used. This approximation is the
most stable and is recommended for most applications. However, the backward difference
introduces the most numerical dispersion into the solution. For CO = 0.5, a central difference
approximation is used. For co = 1.0, a forward difference is used. However, as implemented in
IPX, central and forward differences are numerically unstable and are not recommended.
The use of finite difference approximations in place of continuous spatial derivatives
introduces numerical dispersion in the model solution. Numerical dispersion has the effect of an
additional dispersive transport term on the numerical solution computed for the mass balance
equation. If the advection factor co = 0, a backward difference approximation of spatial
derivative (dC/dx) is used for the advection term. The degree of numerical dispersion resulting
from the backward difference approximation of the spatial derivative is:
—
—
(1.128)
where: u = water (advective) velocity
L = length of the segment (Ax)
For the Euler numerical integration method (the solution method used in IPX), the
forward difference approximation of the temporal derivative (9C/3r) is used. The degree of
numerical dispersion resulting from the forward difference approximation of the temporal
derivative is:
u At
(1.129)
Therefore, the overall numerical dispersion resulting from the approximation of both the spatial
and temporal derivatives is:
81
-------
Enum =(L-
(1.130)
Note that increasing the time-step to L/u (= Ax/u or V/Q) decreases numerical dispersion
to zero. However, the conditions for numerical stability require a time-step somewhat less than
V/Q for most segments. So, to maintain stability and minimize numerical dispersion in a water
body subject to unsteady flow, the sequence of time-steps should be as large as possible, but
always less than At max given in Equation 1.127.
A non-zero advection factor is helpful in situations where the network size and time-step
produce large numerical dispersion. A non-zero advection factor reduces the numerical
dispersion produced by a particular velocity, length, and time-step combination. According to
Bella and Grenney (1970):
(1.131)
Note that a oo of 0 reduces this to Equation 1.130. Values of Enum for a segment length of
2000 meters and various combinations of velocities and time-steps are presented in Table 1.19.
For a particular velocity, such as 0.4 m/sec, numerical dispersion can be reduced by increasing
the time-step. For CO = 0, increasing the time-step from 1000 to 4000 seconds decreases Enum
from 320 to 80 m2/s. If the time-step must be 1000 seconds, however, numerical dispersion can
still be reduced by increasing oo. In this case, increasing o> from 0 to 0.4 decreases Enum from
320 to 0 m2/s.
Negative solution option-Normally, concentrations are not allowed to become negative.
Unless the negative solution option is selected, if the concentration predicted at time t + At is
negative, IPX maintains a positive value by halving the concentration that was present at time t
instead of allowing the value to become negative. The negative solution option lets the user
bypass this procedure, allowing negative concentrations. This may be desirable for simulating
dissolved oxygen deficit in the benthos, for example.
1.11.2 Boundary Parameters
Boundary concentrations. mg/L-Steady or time-variable concentrations must be specified
for each water quality constituent, at each boundary. A boundary is either a tributary inflow, a
downstream outflow, or an open water end of the model network across which dispersive mixing
82
-------
Table 1.19. Numerical Dispersion (m2/s) values resulting from
u and At for a Segment Length of L = 2000 m
Q)
0.0
0.1
0.2
0.3
0.4
0.0
0.1
0.2
0.3
0.4
0.0
0.1
0.2
0.3
0.4
0.0
0.1
0.2
0.3
0.4
u (m/s)
0.1
0.2
0.4
0.6
0.8
1.0
At = 1000 sec
95
75
55
35
15
180
140
100
60
20
320
240
160
80
0
420
300
180
60
—
480
320
160
0
—
500
300 '
100
—
—
At- - 2000 sec
90
70
50
30
10
160
120
80 .
40
0
.240
160
80
0
—
240
120
0
—
—
160
0
—
—
—
0
—
—
—
—
At' = 4000 sec
80
60
40
20
0
120
80
40
0
_..
80
0
—
—
•
—
—
—
—
—
—
—
—
—
—
—
.
—
—
—
At = 8000 sec
60
40
20
0
—
40
0
—
_..
__
—
—
—
—
—
—
—
—
—
~
—
—
—
__
--
.
— -
—
—
—
can occur. Advective and dispersive flows across boundaries are specified by the transport
parameters.
1.11.3 Point Source and Non-Point Source Loads
Point Source Loads. kg/day-Loads (external inputs) may be specified for each water
quality constituent in any segment. Specified loads can steady or time-variable. These loads
represent municipal and industrial wastewater discharges, urban and agricultural runoff,
precipitation, and atmospheric deposition of pollutants, or any other external contaminant source.
83
-------
1.11.4 Initial Conditions
Initial concentrations. mg/L-Concentrations of each constituent in each segment must be
specified for the time at which the simulation begins. (Note: the concentrations specified will be
used as the starting conditions for a simulation regardless of the start time specified.) For those
water bodies with low transport rates, the initial concentrations of conservative substances may
persist for a long period of time. Accurate simulation requires accurate specification of initial
contaminant concentrations.
Dissolved fractions. dimensionless~The initial fraction of chemical dissolved in the water
portion of a segment input as a fraction of total chemical concentration (0-1). The dissolved
fraction is important in. determining the mass of chemical transported by pore water flow and
dispersion, and by solids transport. For simulation where equilibrium partitioning and binding
are specified, dissolved fractions are computed from sorption parameters.
Solid densities, g/cm3 —The density of each solid type. This is needed to compute the
porosity of sediment bed segments. Porosity is a function of sediment concentration and the
density of each solid type.
Maximum concentrations. mg/L—Maximum concentrations must be specified for each
water quality constituent. The simulation is automatically aborted if a calculated concentration
falls outside these limits. This usually arises from computational instability, and indicates the
time-step must be reduced.
1.12 APPLYING THE MODEL
The first step when applying the model is to analyze the problem to be solved. What
questions are being asked? How can a simulation model be used to address these questions? A
water quality model can do three basic tasks: simulate present water quality conditions, provide
generic predictions, and provide site-specific predictions. The first, descriptive task is to extend
in some way a limited site-specific data base. Because monitoring is expensive, data seldom give
the spatial and temporal resolution needed to fully characterize a water body. A simulation
model can be used to interpolate between observed data; for example, locating the dissolved
oxygen sag point in a river or the maximum salinity intrusion in an estuary. Of course such a
model can be used to guide future monitoring efforts. Descriptive models also can be used to
infer the important processes controlling present water quality. This information can be used to
guide not only monitoring efforts, but also model development efforts.
84
-------
Providing generic predictions is a second type of modeling task. Site-specific data may
not be needed if the goal is to predict the types of water bodies at risk from a new chemical, or
the critical conditions or risk (potential occurrence of event) for a particular contaminant
exposure. A crude set of data may be adequate to screen a list of chemicals for potential risk to a
\
particular water body. Generic predictions may sufficiently address the management problem to
be solved (provided that the uncertainty of such predictions are explored and communicated with
the results), or they may be a preliminary step in detailed site-specific analyses.
Providing site-specific predictions is the most stringent modeling task. Calibration to a
good set of monitoring data is definitely needed to provide credible predictions. Because
predictions often attempt to extrapolate beyond the present data base, however, the model also
must have sufficient process integrity. Examples of this type of application include waste load
allocation to protect water quality standards and feasibility analysis for remedial actions, such as
tertiary treatment, phosphate bans, or agricultural best-management practices.
Analysis of the problem should dictate the spatial and temporal scales for the modeling
analysis. Division of the water body into appropriately sized segments was discussed in Section
1.3. Users should extend the network upstream and downstream beyond the influence of the
waste loads being studied to minimize the influence of uncertain boundary conditions. Users
should also design the model network so that sampling stations and points of interest (such as
water withdrawals) fall near the center (centroid) of a segment. Point source loads in streams and
rivers with unidirectional flow should be located near the upper end of a segment. In estuaries
and other water bodies with reversing flow, waste loads are best located near the centers of
segments.
Once the model network is defined, a modeling study will proceed generally through four
steps that involve, in some manner, hydrodynamics, mass transport, water quality
transformations, and environmental toxicology. The first step addresses the issue of where water
goes. This can be answered by a combination of gaging, special studies, and hydrodynamic
modeling. Flow data or outputs from hydrodynamic models can be interpolated or extrapolated
throughout the model network using the principle of continuity. The second step answers the
question of where the material in the water is transported. This can be answered by a
combination of tracer studies and model calibration. Dyes, salinity (conductivity), and heat are
often used as tracers. The third step answers the question of how materials in the water and
sediments are transformed and where they are ultimately transported. This is the main focus of
many studies. Answers will generally depend on a combination of laboratory studies, field
85
-------
monitoring, parameter estimation, calibration, and testing. The net result is sometimes called
model validation or verification, which are elusive concepts. The success of this step depends on
the skill of the user, who must combine specialized knowledge with common sense and
skepticism into a methodical process. The final step answers the question of how this material is
likely to affect target organisms, such as people or fish, or the ecological balance. Often,
predicted concentrations are simply compared with water quality criteria adopted to protect the
general aquatic community. Care must be taken to insure that the temporal and spatial scales
assumed in developing the criteria are compatible with those predicted by the model. Sometimes
principles of physical chemistry or pharmacokinetics are used to predict chemical body burdens
and resulting biological effects.
86
-------
CHAPTER 2
IPX INPUT DATA FILE STRUCTURE
2.1 IPX INPUT DATA FILE STRUCTURE OVERVIEW
The input required to run the IPX water quality modeling framework is divided into a
series of data groups. The purpose of these data groups is to arrange all required model inputs
into a consistent, logical structure. Each data group is used to enter model inputs related to a
specific aspect of a simulation such as reaction kinetics and transport processes, etc. Data within
a group are related to other data within the same group. An overview of each input file data
group, A through K, is presented. Each data group is further divided into data records.
Following this overview, the structure of each data group and its records are described in detail.
A - Model Identification and Simulation Control
Data Group A is used for model identification and specifying simulation control options.
This data group requires specification of the number of segments and the number of systems.
Numerical integration time-steps and print intervals for reporting simulation output are also
specified.
B - Exchange (Dispersion) Coefficients
Data Group B is used to input dispersive exchanges for the water column and sediment
pore water.
C - Volumes
Data Group C is used to input the volumes and hydraulic characteristics for each model
segment. •
D - Flows, Sediment Transport, and Precipitation/Evaporation .
Data Group D is used to input water column flows, pore water flows, settling velocities,
resuspension velocities, and precipitation/evaporation rates.
87
-------
E - Boundary Concentrations
Data Group E is used to input the concentration of each state variable at the model
boundaries. Boundary conditions must be specified for all state variables at each boundary.
F - Loads (Forcing Functions)
Data Group F is used to input loads (forcing functions) of each state variable to any
model segment. Loads can be positive (source) or negative (sink) in magnitude and represent
both point and non-point source inputs.
G - Environmental Parameters
Data Group G is used to input site-specific environmental characteristics of the modeled
system that are constant temporally but can vary spatially (from segment to segment). A list of
all parameters that may be included in a simulation are presented in Table 2.2.
H - Chemical Constants (Physicochemical Properties)
Data Group H is used to input the physicochemical properties of each state variable
(solids type and chemical) simulated. Reaction yields for chemicals that transform from one
modeled chemical to another modeled chemical are also specified. A list of all constants that
may be included in a simulation are presented in Table 2.3-2.12.
I - Time Functions
Data Group I is used to input environmental characteristics of the modeled system that are
constant spatially but can vary temporally. A list of all time functions that may be included in a
simulation are presented in Table 2.13.
J - Initial Conditions - >
Data Group J is used to input initial concentrations (concentrations at the start time of the
simulation) of each state variable in each model segment.
K - Surficial Sediment Age Layer Conditions
Data Group K is used to define surficial sediment aging characteristics and the initial
distribution of sediment resuspension properties in each age layer for each surficial sediment
segment in the model. The aging properties of all sediments in each age layer and the
resuspension properties of freshly deposited sediments are also specified.
88
-------
L - Initial Conditions for Deep Sediment Layers (Ghost Stack)
Data Group L is used to define initial conditions for deep sediment layers when the semi-
Lagrangian sediment bed option is used. The physical configuration and physicochemical
properties of each deep sediment layer (if any) in each sediment stack using the semi-Lagrangian
bed option are specified.
2.2 IPX INPUT FILE DATA GROUP DESCRIPTIONS
2.2.1 Data Group A; Model Identification and Simulation Control
All parameters that define the number of state variables simulated, the number of model
segments, time-step for numerical integration.
Record 1-Title of Simulation (A80)
HEADS = descriptive title of simulation (A80).
Record I—Description of Simulation (A80)
TYPE = description of simulation (A80).
Record 3--Names of Variables Listed in Record 4 (A80)
HEADER = names of Record 4 variables, positioned properly; for user
convenience only (A80).
Record 4-Simulation Control Parameters (815, 2F5.0, F3.0, F2.0, 315, F10.0)
KSIM = simulation type; 0 = dynamic. (15)
Previous versions of the WASP4 framework from which IPX was
derived allowed the user to chose either a dynamic or steady-state
solution. However, IPX will only perform a dynamic mass balance
so the value of KSIM must always be 0. This variable was left in
the IPX code so users could use previously developed WASP4
input data files with a minimum of effort.
NOSEG
NOSYS
= number of segments in model network. (15)
= number of model systems (state variables). (15)
Systems 1-3 are solids and Systems 4-n are chemicals. If only one
solids type is simulated with several chemicals, Systems 1, 2, or 3
89
-------
ICFL
MFLAG
JMASS
NEGSLN
could be used to simulate the single solids type while the two
unused Systems would be "placeholders". Most IPX simulations
will require that at least four systems be specified. For example, if
a user wishes to simulate one chemical and one solid, the number
of systems that must be specified is four: System 1 for all solids,
Systems 2-3 would be "dummy" systems for which all
computations should bypassed, and System 4 would be for the
chemical.
= flag controlling use of restart file; 0 = neither read from nor write
to restart file (initial conditions located in input'file); 1 = write
final simulation results to restart file (initial conditions located in
input file); 2 = read initial conditions from restart file created by
earlier simulation, and write final simulation results to new restart
file. (15) (Option Not Available. Set ICFL = 0)
= flag controlling messages printed on screen during simulation; 0 =
all messages printed; 1 = simulation time only printed; 2 = no
messages printed. (15)
; system number for which export computations will be performed.
(15)
IPX computes a detailed, segment by segment mass balance for all
state variables (solids type or chemical) simulated. Export is
defined as the sum of the mass transported by advection and
dispersion across any model interface (internal or at a model
boundary). In addition to computing export for the system (state
variable) specified by JMASS, IPX automatically computes export
for the sum all solids types simulated.
negative solution option; 0 = prevents negative solutions; 1 =
allows negative solutions. (15)
The authors recommend that NEGSLN = 0 for all simulations.
This option was originally developed to allow WASP4 users to
simulate dissolved oxygen (DO) etc. where negative concentrations
90
-------
INTYP
ADFAC
ZDAY
ZHR
MIN
IDSY
IDSG1
IDSG2
represent DO deficit conditions. Code operation for NEGSLN = 1
has not been verified.
= time-step option; 0 = user inputs time-steps; 1 = model calculates
time-step. (15)
The authors recommend that INTYP = 0 and model time-steps be
input by the user. Code operation for INTYP = 1 has not been
verified.
= advection factor; 0 = backward difference; 0.5 = central difference;
[1 = forward, difference but is unconditionally unstable;] 0
recommended. (F5.0)
= day at beginning of simulation; 1 = first day. (F5.0)
= hour at the beginning of simulation. (F3.0)
= minute at the beginning of simulation. (F2.0)
= system number for which results (concentrations) will be displayed
on screen during simulation. (15)
= "upstream" segment for which the mass export of JMASS will be
computed and for which results (concentrations) will be displayed
on screen. (15)
= "downstream" segment for which the mass export of JMASS will
be computed and for which results (concentrations) will be
displaj'ed on screen. (15)
IDSG1 and IDSG2 define the model segment interface for which
export computations will be performed. Export is computed as the
sum of the mass transported across a model interface by advection
and dispersion from IDSG1 to IDSG2. If the, interface for which
export is calculated is at a model boundary, one of the segments
specified will be 0. Care must be taken to choose IDSG1 and
- IDSG2 such that the pair specifies a meaningful model interface. If
91
-------
IDSG1 and IDSG2 do not define a model interface, all export
computation will be spurious.
TADJ = factor by which input kinetic rates are adjusted; 0 or 1.0 will cause
no adjustment; 24.0 will adjust rates input as hours"1 to days"1;
86400.0 will adjust rates input as seconds"1 to days"1. (F10.0)
Record 5-Time-Steps (15, F5.0)
NOBRK = number of different model time-steps (15)
TRST = frequency (time interval) in days at which output is written to the
restart file; default value (if TRST is zero or not specified) = 5
days. (F5.0)
Record 6~Time-Steps (4(F10.0, F10.0)
SIMFVALCO = simulation time-step in days that will be used until the simulation
time is greater than SMFTIM(I). (F10.0)
SIMPriMOQ = simulation time in days up to when time-step SIMFVAL(I) will be
used. (F10.0)
where I = 1, NOBRK
Record 7~Number of Print Intervals for Reporting Output to *.dmp File and Option for
Temporal Averaging (215)
NPRINT1 = number of print intervals for reporting output to *.dmp file. (15)
VFLAG = option for reporting time averaged output to an additional file; 0 =
instantaneous (point value) results to *.dmp file; 1 = time-averaged
results reported to *.dma file (which is generated in addition to the
*.dmpfile). (15)
Record 8-Print Intervals for Reporting Ouptput to *.dmp File (4(F10.0, F10.0))
PRINT1(T) = simulation results print interval in days that will be used until the
simulation time is greater than TPRNTl(I). (F10.0)
92
-------
TPRNTIQ) = simulation time in days up to when print interval PRINT(I) will be
used. (F10.0)
where I = l.NFRINTr
Record 9-Number of Print Intervals for Reporting Output to *.exp File (15)
NPRINT2 = number of print intervals for reporting to *.exp file. (15)
Record 10-Print Intervals for Reporting Output to *.exp File (4(F10.0, F10.0))
PRINT2(I) = simulation results print interval in days that will be used until the
simulation time is greater than TPRNT2(T). (F10.0)
TPRNT2(I) , = simulation time in days up to when print interval PRINT2(I) will
be used. (F10.0)
where I = 1.NPRINT2
Record 1 l~Number of Print Intervals for Reporting Output to *.dma File (15)
NPRINT3 = number of print intervals for reporting to*, dma file. (15)
Record 12-Print Intervals for Reporting Output to *.dma File (4(F10.0, F10.0))
PRINTS® = simulation results print interval in days that will be used until the
simulation time is greater than TPRNT3(I). (F10.0)
TPRNT3(I) = simulation time in days up to when print interval PRINT® will be
used. (F10.0)
where I = 1.NPRINT3
Record 13-System Bypass Options (1615)
SYSBYQSYS) = bypass option for system ISYS; 0 = system will be simulated; 1 =
system will be bypassed. (15)
where ISYS
= l.NOSYS
Organization of Records. Records 1 through 10 and 13 are entered once for Data Group
A. Records 11 and 12 are entered once only if the value of VFLAG entered in Record 7 is equal
93
-------
to 1. If the value of VFLAG is 0, Records 11 and 12 are skipped. Record 8 uses as many lines as
needed to input NPRBSfTl pairs of. PRINT1(I) and TPRNTl(I), 4 pairs occupying each line.
Record 10 uses as many lines as needed to input NPRINT2 pairs of PRJNT2(I) and TPRNT2(I),
4 pairs occupying each line. Record 12 (if needed) uses as many lines as needed to input
NPRINT3 pairs of PRINTS® and TPRNT3Q), 4 pairs occupying each line. Record 13 will have
NOSYS entries (including any "dummy" solids state variables).
2.2.2 Data Group B: Exchange (Dispersion) Coefficients
Exchange coefficients are computed from input dispersion coefficients, cross-sectional
areas, and characteristic lengths. Dispersion coefficients may vary in time according to
piecewise linear time functions, with groups of segment pairs having the same dispersion time
function. Exchange data are read for each exchange field. Field one contains dispersion
coefficients for water column exchanges. Field two contains exchange data for pore water
exchange. Fields three, four and five contain sediment exchange data, with a separate field
available for each solid type.
Record 1-Number of Exchange Fields (15, 75X)
NRFLD = number of exchange fields; 1 = total constituent (dissolved,;bound,
and paniculate) exchanges (usually water column), 2 = dissolved
constituent exchanges (usually water column and/or pore water);
NRFLD will generally equal 1 or 2. (K)
It is also possible to specify exchanges to represent sediment
mixing for each solids type simulated (Field 3 for Solids 1, Field 4
for Solids 2, and Field-5 for Solids 3). However, this is not
generally recommended.
Chemical mixing in sediment layers (exclusive of ghost stack
elements) can be specified in exchange file 1.
TITLE
= name of data group (for user convenience). (75X)
If no exchange rates are to be read, set NRFLD to zero and continue with Data Group C. If
NRFLD is greater than zero, Records 2 through 6 are read in and repeated as a group NRFLD
times.
94
-------
Record 2-Exchange Time Functions for Each Exchange Field (15, 2F10.0)
NTEX(NF) = number of exchange time functions for field NF. If no exchange
time functions are input for field.NF, set NTEX to zero and
continue with the next exchange field. (15)
SCALR
CONVR
where NF
= scale factor for exchange coefficients. All exchange coefficients
for field NF will be multiplied by this factor. (F10.0)
= conversion factor for exchanges in field NF. (F10.0)
= 1,NRFLD
Record 3-Exchange Data (15)
NORS(NT) = number of exchanges for field NF, time function NT. (15)
where NT
= 1,NTEX(NF)
Record 4-Areas, Characteristic Lengths (2F10.0, 215)
A(K) = area in m2 for exchange pair K defined by JR and IR. (F10.0)
EL(K) = characteristic length in m for exchange pair K. (F10.0)
JR(K),IR(K) = segment pair between which exchange occurs. The order of the
segments is unimportant. (215)
where K . = 1,NORS(NF,NT)
RecordS-Number of Time-Breaks in the Exchange Time Function (K)
NBRKR(NF,NT) = number of time-break points in the exchange time function. (15)
Record 6-Piece Linear Dispersion Time Function (4(F10.0, F10.0))
RT(K) = value of dispersion coefficient in m2/sec at time TR(K). (F10.0)
TR(K)
where K
= time in days at which the dispersion coefficient is RT(K). (F10.0)
= 1,NBRKR(NF,NT)
95
-------
Record 7—Exchange Bypass Options (1615)
RBY(ISYS) = exchange bypass option; 0 = exchange occurs for system ISYS; 1
bypass exchange for system K. (15)
where ISYS
= 1,NOSYS
Organization of Records. Record 1 is entered once for Data Group B. Records 2 through
6 are repeated for each exchange field, and Records 3, 4, 5, and 6 are repeated for each time
function in a given exchange field. Record 4 uses as many lines as necessary to input NORS sets
of A(K), EL(K), IR(K), and JR(K) (1 set per line). Record 6 uses as many lines as needed to
input NBRKR pairs of RT(K) and TR(K), (4 pairs occupying each line. After data for all
exchange fields are entered, Record 7 is input on the following line with NOSYS entries.
2.2.3 Data Group C: Volumes
Segment volumes, types, and hydraulic characteristics are input by the user. The
segments above and below each model segment are specified to compute sediment transport
between the water and sediments and burial through each sediment layer. The user also specifies
options for sediment segment volumes in response to sediment transport.
Record 1-Preliminary Data (215, F10.0, 60X)
IVOPT = water column volume option; 1 = constant water column volumes
(water column inflows and outflows specified in Data Group D
must balance); 2,3 = water column volumes adjusted (change with
time) to maintain flow continuity (water column flows specified in
Data Group D need not balance). (15)
IBEDV
TDINTS
= sediment bed volume option; 0 = constant sediment bed volumes; 1
= sediment bed volumes change in response to sediment transport.
05)
Surficial sediment segment volumes in IPX vary in response to
sediment transport so IBEDV must be 1. The porosity of all
sediment segments is always constant.
= sediment bed control options; time-step in days for porosity
computations when IBEDV = 0; time-step in days for burial and
96
-------
TITLE
unburial to and from deeper sediment layers when IBEDV
(F10.0)
= name of data group. (60X)
= 1.
Record 2-Scale Factors (2F10.0)
SCALV = scale factor for segment volumes. All volumes are multiplied by
this factor. (F10.0)
CONVV
= conversion factor for volumes. (F10.0)
Record 3-Segment Types and Volumes (110, 215,110, 5F10.0)
ISEG = segment number. (110)
ITOPSG(ISEG) = segment immediately above ISEG; 0 = that there is no segment
above (i.e. at a model boundary). (15)
IBOTSG(ISEG) = segment immediately below ISEG; 0 = that there is no segment
below (i.e. at a model boundary). (15)
HYPE(ISEG)
BVOL(ISEG)
VMULT(ISEG)
VEXP(ISEG)
= segment types;
1 = surface water segment,
2 = subsurface water segment,
3 = upper bed segment (Eulerian bed option),
4 = lower bed segment (Eulerian bed option)),
5 = upper bed segment (Semi-Lagrangian bed option),
6 = lower bed segment (Semi-Lagrangian bed option). (110)
= volume of segment ISEG in cubic meters. (F10.0)
= hydraulic coefficient "a" for velocity in ISEG as a function of flow:
v = aQb
If b = 0, VMULT is a constant velocity in m/sec. (F10.0)
= hydraulic exponent "b" for velocity in ISEG as a function of flow
(0-1). A value of 0.4 represents rectangular channels. (F10.0)
DMULT(ISEG) = hydraulic coefficient "c" for depth of ISEG as a function of flow:
97
-------
d = cQ
If d = 0, DMULT is a constant depth in m. (F10.0)
DXP(ISEG) = hydraulic exponent "d" for depth of ISEG as a function of flow (0-
1). A value of 0.6 represents rectangular channels. (F10.0). A
value of d = 0 should be used for sediment segments.
where ISEG = l.NOSEG
Organization of Records. Records 1 and 2 are entered once for Data Group C. Record 3
is repeated NOSEG times. If ICFL = 2 in Data Group A, volumes are read from the restart file (
*.RST, where * is the input data file name), and Records 2 and 3 should not be included in the
input data set.
2.2.4
Data Group D; Water Column Flows; Pore Water Flows; Settling Velocities;
Resuspension Velocities; and Precipitation/Evaporation
Data Group D is used to specify the flows, sediment transport velocities, and
precipitation/evaporation inputs used in the model and may be specified for any or all of several
fields. Field one is for advective flows in the water column. Field two is for pore water flows.
Fields three, four, and five are for settling velocities for solids types one, two, and three,
respectively. Field six is for resuspension velocities and is applied to all solids types simulated.
Field seven is for precipitation and evaporation. All time functions specified are treated as
piecewise linear time functions.
Overall Organization of Records for Data Group D. Record 1 is read first. Data Sub-
Group Dl is read next. Data Sub-Groups D2, D3, D4, D5, D6, and D7 follow in order for
NFIELD = 2, 3, 4, 5, 6, and 7, respectively. Following all specified Data Groups, Record 7 is
read.
Record 1—Data Input Options: Number of Flow Fields (215)
IQOPT = Surface Water flow (Field 1) input option:
1 = surface water flows are specified directly by user (qsurfl: mass
transported by the net outflow).
2 = water column flows are specified directly by user (qsurfZ: mass
transported by the difference between the gross inflow minus the
gross outflow).
98
-------
3 = water column flows are read from a formatted file created by
DYNHYD. (15) (Option Disabled)
The user is required to directly input the surface water flows in IPX
so IQOPT must be 1 or 2. All linkages to the DYNHYD
hydrodynamic model have been disabled.
NFffiLD = number of flow fields. Field 1 is for surface water (advective)
flows. Field 2 is for pore water flows. Field 3 is for Solids 1
settling velocities. Field 4 is for Solids 2 settling velocities. Field
5 is for Solids 3 settling velocities. Field 6 is for resuspension
velocities and is applied to all solids types. Field 7 is for
precipitation and evaporation. If no transport fields are to be input,
set NFIELD to zero and continue with Data Group E. (15)
TITLE = name of data group. (VOX)
DATA SUB-GROUP Dl: FIELD ONE (WATER COLUMN) FLOWS (IQOPT = 1,2)
Surface water flows are input as a separate time function for each inflow to the model.
Each inflow is routed from a boundary (segment 0) through each water column segment water
from a given source flows until it leaves the system through a boundary. All phases (dissolved,
DOC-bound, and paniculate) of each system simulated are transported by this field.
Record 2—Water Column Flow Time Function Specification (15, 2F10.0)
NINQ(l) = number of time functions for Field One. If water column flows are
not used, set NINQ to zero and skip to Data Sub-Group D2. (15)
SCALQ
CONVQ
= scaling factor. All flows in Field one are multiplied by SCALQ.
(F10.0)
= units conversion factor. (F10.0)
Record 3~Number of Flows (15)
NOQS(1,NI) = number of segment pair interfaces (including those at boundaries)
across which surface water from source NI flows . (15)
where NI
99
-------
Record 4-Flow Routing (4(F10.0, 2115))
BQ(1,NI,K) = decimal fraction of flow time function NI that flows between
segment pair K defined by JQ and IQ. (F10.0)
JQ(14SU,K) = segment water is flowing from. (K)
IQ(1,NI,K) = segment water is flowing to. (15)
where K = 1,NOQS(1,NI)
Record 5-Number of Breaks in Advective Time Functions (15)
NBRKQ(1,NI) = the number of flow and time pairs used to describe piecewise linear
time function NI. (15)
Record 6-Piecewise Linear Advective Time Function (4(2F10.0))
QT(1,NI,K) = water column (advective) flow in m3/s. (F10.0)
TQ(1,NI,K)
where K
= time in days.. (F10.0)
= 1,NBRKQ(1,NI)
Organization of Records for Data Sub-Group PI. Record 2 is input once. Records 3, 4,
5, and 6 are input once for each flow time function. Record 4 uses as many lines as needed to
input NOQS sets of BQ, JQ, and IQ, with four sets per line. Record 6 uses as many lines as
necessary to input NBRKQ sets of QT and TQ, with four sets on each line.
DATA SUB-GROUP D2: FIELD TWO (PORE WATER) FLOWS
Pore water flows are input as separate time functions for each inflow to the model. Each
inflow is routed from a boundary (segment 0) and through the sediments. If pore water are
routed to the water column, a surface water flow function must also be specified for field 1 Only
the dissolved and DOC-bound phases of each systems simulated are transported by this field.
Record 2--Number of Pore Water Time Functions (15, 2F10.0)
NINQ(2) = number of pore water time functions. If pore water flows are not
used, set NINQ(2) to zero and skip to Data Sub-Group D3. (15)
100
-------
SCALQ
= scaling factor for pore water flows. (F10.0)
CONVQ
= units conversion factor. (F10.0)
Records—Number of Hows (15)
NOQS(2,NT) = number of segment pair interfaces (including those at boundaries)
across which pore water from source NI flows.
where NI
= 1,NINQ(2)
Record 4-Flow Routing for Pore Water Hows (4(F10.0, 215))
BQ(2,NI,K) = decimal fraction of pore water flow for time function NI that flows
between segment pair K defined by JQ and IQ. (F10.0)
JQ(2,NI,K) = segment pore water is flowing from. (15)
IQ(2,NI,K) = segment pore water is flowing to. (15)
Record 5—Number of Breaks in Pore Water Time Function (K)
NBRKQ(2,NI) = number of pore water flow and time pairs used to describe
piecewise linear time function NI. (15)
Record 6~Piecewise Linear Pore Water Time Function (4(2F10.0))
QT(2,NI,K) = pore water flow in m3/s. (F10.0)
TQ(2,NI,K)
= time in days. (F10.0)
where K
= 1,NBRKQ(2,NI)
Organization of Records for Data Sub-Group D2. Record 2 is input once. Records 3, 4,
5 and 6 are input once for each pore water time function. Record 4 uses as many lines as
necessary to input NOQS sets of BQ, JQ, and IQ, with four sets on each line. Record 6 uses as
many lines as necessary to inputNBRKQ sets of QT and TQ, with four sets on each line.
101
-------
DATA SUB-GROUPS D3-D5: SEDIMENT TRANSPORT - SETTLING VELOCITIES
Settling velocities are input separately for each solids type simulated. A separate field is
specified for each solid type: Data Sub-Group D3 is used to input settling velocities for solids
type one; D4 for solids type two; and D5 for solids type three. Velocities may vary with time.
Only the specified solids type and the chemical sorbed to that solids type are transported by each
field.
Record 2~Number of Settling Velocity Time Functions (15, 2F10.0)
NINQ(NF) = number of settling velocity time functions for Field NF. If settling
velocities for a given solids type are not used, or if a solids type is
not simulated, set NINQ(NF) to zero and skip to the next Data
Sub-Group. (15)
SCALQ
CONVQ
where NF
= scaling factor for settling velocities. (F10.0)
= units conversion factor. (F10.0)
= 3,5
Record 3—Number of Segment Pairs (15)
NOQS(NF,NT) = number of segment surface water-surficial sediment pair interfaces
across which sediments settle at the velocities specified by time
function NL (15)
where NI
= 1,NINQ(NF)
Record 4~Areas for Settling (4(F10.0, 215))
BQ(NF,NI,K) = area in m2 between segment pair K defined by JQ and IQ. (F10.0)
JQ(NF,NI,K) = water column segment that solids type (NF-2) settles from. (15)
IQ(NF,NI,K) = surficial sediment segment that solids type (NF-2) settles to. (15)
where K
= 1, NOQS(NF,NI)
102
-------
Record 5-Number of Breaks in Velocity Time Function (15)
NBRKQ(NF,NT) = number of settling velocity and time pairs used to describe
piecewise linear time function NI. (15)
Record 6-Piecewise Linear Settling Velocity Time Function (4(2F10.0))
QT(NF,NI,K) = settling velocity in m/s. (F10.0)
TQ(NF,NI,K) . = time in days. (F10.0)
where K
= 1,NBRKQ(NF,NT)
Organization of Records for Data Sub-Groups D3-5. Records 2 through 6 are input for
each solid type. Records 3, 4, 5 and 6 are input for each time function in each field. Record 4
uses as many lines as needed to input NOQS sets of BQ, JQ, and IQ, with four sets on one line.
Record 6 uses as many lines as needed to input NBRKQ sets of QT and TQ, with four sets per
line.
DATA SUB-GROUPS D6: SEDIMENT TRANSPORT - RESUSPENSION VELOCITIES
Resuspension velocities are input once and applied to each solids type simulated. One
field is specified for all solid types; Data Sub-Group D6 is used to input all resuspension
velocities. Velocities may vary with time. All solids simulated and the chemical sorbed to each
solid type are transported by this field.
Record 2—Number of Resuspension Velocity Time Functions (15, 2F10.0)
NINQ(6) = number of resuspension velocity time functions for Field NF. If
resuspension velocities are not used, or if solids are not simulated,
set N1NQ(6) to zero and skip to Data Sub-Group D7. (15)
SCALQ
CONVQ
= scaling factor for resuspension velocities. (F10.0)
= units conversion factor. (F10.0)
Record 3-Number of Segment Pairs (15)
NOQS(6,NI) = number of surficial sediment-surface water segment pair interfaces
across which sediments resuspend at the velocities specified by
time function NI. (15)
103
-------
where NI
= 1,N1NQ(6)
Record 4-Areas for Resuspension (4(F10.0,215))
BQ(6,NI,K) = area in m2 between segment pair K defined by JQ and IQ. (F10.0)
JQ(6,NI,K) = surficial sediment segment that solids resuspend from. (15)
IQ(6,NI,K) = water column segment that solids resuspend to. (15)
where K = 1,NOQS(6,NT)
Record 5~Number of Breaks in Velocity Time Function (15)
NBRKQ(6,NT) = number of resuspension velocity and time pairs used to describe
piecewise linear time function NI. (15)
Record 6~Piecewise Linear Resuspension Velocity Time Function (4(2F10.0))
QT(6JMI,K) = resuspension velocity in m/s. (F10.0)
TQ(6,NI,K)
where K
= time in days. (F10.0)
= 1,NBRKQ(6,NI)
Organization of Records for Data Sub-Group D6. Records 2 is input once. Records 3, 4,
5 and 6 are input for each time function in each field. Record 4 uses as many lines as needed to
input NOQS sets of BQ, JQ, and IQ, with four sets on one line. Record 6 uses as many lines as
needed to input NBRKQ sets of QT and TQ, with four sets per line.
DATA SUB-GROUP D.7: PRECIPITATION AND EVAPORATION
Precipitation and evaporation data are input to represent rainfall, seasonal evaporation, or
other events that change surface water levels. In response to the specified time functions water
volumes are adjusted to maintain continuity. In addition to modifying water volumes, chemicals
and solids are transported into the model by precipitation. However, since only one boundary
condition can be specified per model segment, it is recommended that all contaminant
inputs resulting from precipitation be treated as loads. No chemical or solids mass is
transported out of the model by evaporation.
104
-------
Record 2--Number of Precipitation/Evaporation Time Functions (15,2F10.0))
N1NQ(7) = number of precipitation and evaporation time functions. (15)
SCALQ
CONVQ
= scaling factor for precipitation/evaporation. (F10.0)
= units conversion factor. (F10.0)
Record 3--Number of Segment Pairs (15)
NOQS(7,NI) = number of surface water segment-air interfaces across which water
precipitates/evaporates at the rate specified by time function ML
(15)
Record 4-Areas for Precipitation/Evaporation (4(F10.0, 215))
BQ(7,NI,K) = area in m2 between segment pair K defined by JQ and IQ. (F10.0)
JQ(7,NI,K) = segment water is transported from; use JQ = 0 to represent
precipitation. (15)
IQ(7,NI,K) = segment water is transported to; use IQ = 0 to represent
evaporation. (15)
where K - = 1,NOQS(7,NI)
Record 5~Number of Breaks in Precipitation/Evaporation Time Function (15)
NBRKQ(7,NT) = number of precipitation/evaporation rate and time pairs used to
describe piecewise linear time function NI. (15)
Record 6-Piecewise Linear Precipitation/Evaporation Time Functions (4(2F10.0))
QT(7,NI,K) = precipitation/evaporation rate in m/s; if more traditional units of
cm/day or cm/year are desired, then specify CONVQ = 1.1574E-7
or 3.169E-10, respectively. (F10.0)
TQ(7,NI,K)
where K
= time in days. (F10.0)
= 1,NBRKQ(7,NI)
-End of Data Sub-Groups For Data Group D-
105
-------
Record 7-Transport (Data Group D) Bypass Options (1615)
QBY(ISYS) = transport bypass option for system ISYS; 0 = system will be
transported; 1 = system will be not be transported. The bypass
option applies to transport by all fields. (15)
where ISYS
= 1,NOSYS
2.2.5. Data Group E; Boundary Concentrations
Data Group E is repeated, in its entirety, NOSYS times. No boundary conditions are
specified for the bottom segment in any stack of sediments using the semi-Lagrangian sediment
bed option.
Record 1—Number of Boundary Conditions (110, 70X)
NOBC(ISYS) = number of boundary conditions used for system ISYS. If no
boundary conditions are to be input for System ISYS, set
NOBC(ISYS) equal to zero and continue entering boundary
conditions for the next system. When boundary conditions for all
systems have been specified (or skipped by setting NOBC(ISYS)
to zero) skip to Data Group F. (110)
TITLE = name of data group. (70X)
where ISYS = 1, NOSYS
Record 2-Scale Factor for Boundary Conditions (2F10.0)
SCALE = scale factor for boundary conditions for System ISYS. All
boundary conditions for System ISYS will be multiplied by this
factor. (F10.0)
CONVB = unit conversion factor for boundary conditions for System ISYS.
Boundary conditions are expected to be in milligrams per liter
(mg/L). (F10.0)
Record 3-Boundary Conditions (215)
D3C(ISYS,NB) = boundary segment number. (15)
NOBRK(ISYS,NB)= number of boundary condition value and time pairs used to
describe the boundary condition. (15)
106
-------
where NB
= 1,NOBC
Record 4-Boundary Concentrations Time Function(4(2F10.0))
BCT(ISYS,NB,K) = boundary concentration in mg/L at time T(K) for segment IBC.
(FIO.O)
T(ISYS,NB,K)
where K
= time in days. If the length of the simulation exceeds
T(ISYS,NB,NOBRK), the function is assumed to be periodic and
repeated, starting at T(l) with period equation to
T(ISYS,NB,NOBRK). (FIO.O)
= l.NOBRK
Organization of Records for Data Group E. Records 1 and 2 are entered once. Records 3
and 4 are a set and are repeated NOBC times. Within each NOBC set, Record 3 is entered once
and Record 4 is repeated until NOBRK entries are input. Four entries (four BCT(K)-T(K) pairs)
will fit on each 80-space line. The whole group is repeated NOSYS times, once for each model
system. Boundary conditions are never specified for ITYPE = 5 or TTYPE = 6 sediment
segments (see description of semi-Lagrangian bed option).
2.2.6 Data Group F: Point and Non-Point Source Loads (Forcing Functions)
Data Group F is used to specify point and non-point source loads (forcing functions) used
in the model. Point source forcing functions are specified in Data Sub-Group Fl. Non-point
forcing functions are specified in Data Sub-Group F2.
Overall Organization of Records for Data Group F. Data Sub-Group Fl is read first and
is repeated once for each system simulated. Data Sub-Group F2 is read second.
DATA SUB-GROUP Fl: POINT SOURCE LOADS
Data Sub-Group Fl is used to input all point source waste loads in the model. Data Sub-
Group Fl is repeated NOSYS times.
Record I—Number of Forcing Functions (110, 70X)
NOWK(ISYS) = number of forcing functions for System ISYS. Forcing functions
specified as sources (+) or sinks (-) of a water quality constituent
107
-------
TITLE
(system). If forcing functions are not input, set NOWK(ISYS) to
zero, and continue with next system or go to next data group.(I10)
= name of data group. (70X)
Record 2-Forcing Function Scale and Conversion Factors (2F10.0)
SCALW = scale factor for forcing functions. All forcing functions will be
multiplied by this factor. (F10.0)
CONVW = unit conversion factor for forcing functions. Forcing functions are
expected to be in kilograms per day. If forcing functions are given
in English units (pounds per day), this .factor will be 0.4535.
(F10.0)
Record 3~Number of Point Sources (215)
IWK(NW) = segment number to which forcing function WKT(K) is applied.
(15)
NOBRK(NW)
where NW
= number of forcing function and time value pairs used to describe
the forcing function WKT(K). (15)
= 1,NOWK
Record 4~Point Source Time Function (4(2F10.0))
WKT(K) = value of the forcing function at time T(NW,K), in kg/day. (F10.0)
T(K)
where K
= time in days. If the length of the simulation exceeds T(NOBRK),
the approximation is assumed to be periodic and is repeated
starting at T(l) with a period equal to T(NOBRK). (F10.0)
= 1,NOBRK
Organization of Records for Data Sub-Group Fl. Records 1 and 2 are input once.
Records 3 and 4 are a set and repeated as a set NOWK times. Within each set, Record 3 is input
once and Record 4 is repeated until all NOBRK entries are input (four WKT(K), T(K) pair
entries per line). Data Sub-Group Fl is repeated NOSYS times, once for each system.
108
-------
DATA SUB-GROUP F2: NON-POINT SOURCE LOADS (FORCING FUNCTIONS)
Data Sub-Group F2 allows the user non-point source loads that are read in from an
external file named NPS.DAT that must be created prior to running the model.
The authors recommend that all non-point source loads be treated as point source loads
and input through Data Sub-Group Fl. This non-point load option was originally developed to
allow WASP4 users to link water quality simulation to storm water/runoff models such as
SWMM. However, code operation for computing non-point source forcing functions has not
been verified.
Record 1-Number of Runoff Loads, Initial Day (215)
NOWKS = number of segments receiving runoff loads. If non-point source
loads are not specified, set NOWKS equal to zero and skip to Data
Group G. If NOWKS >0, continue with records 2, 3, and 4. (15)
NPSDAY
= the time in the runoff file corresponding to the initial simulation
time, in days. (15)
Record 2~Scale Factor for Runoff Loads (2F10.0)
SCALN = scale factor for runoff loads. All runoff loads will be multiplied by
. this factor. (F10.0)
CONVN
= unit conversion factor for runoff loads. Runoff loads are expected
in kilograms per day. If runoff loads are given in English units
(pounds per day), this factor will be 0.4535. (F10.0)
Record 3-Runoff Segments (1615)
INPS(J) = segment number to which runoff load J is applied. (15)
where J
= 1,NOWKS
Record 4--Print Specifications (1615)
KT1 = initial day for which nonzero runoff loads from file NPS.DAT will
be printed. (15)
109
-------
KT2
KPRT(D
where I
= final day for which nonzero runoff loads from file NPS.DAT will
be printed. (15)
= indicator specifying whether nonzero runoff loads will be printed
for each system. If KPRT(I) is greater than zero, then runoff loads
will be printed for system I. (15)
= 1,NOSYS
Organization of Records for Data Sub-Group F2. Records 1 and 2 are input once.
Record 3 has NOWKS entries and uses as many 80-space lines as needed to enter all NOWKS
segment numbers (16 entries will fit on one line). Record 4 is input once.
2.2.7 Data Group G: Parameters
Data Group G is used to enter all environmental parameters required for a simulation.
The number of parameters needed for any simulation will vary as a function of the type of water
body modeled and the simulation complexity. The input format for all parameters, however, is
constant.
Record I—Number of Parameters (110, 70X)
NOPAM = number of parameters used in the simulation. If no parameters are
specified, set NOPAM to zero and skip to Data Group H. (110)
TITLE
name of data group. (70X)
Record 2-Scale Factors for Parameters (4(A5,15, F10.0)
ANAME(ISC) = descriptive name for parameter ISC. (A5)
ISC
PSCAL(ISC)
where ISC
= parameter number identifying parameter; all parameter numbers
are defined in Table 2.2. (15)
= scale factor for parameter ISC that will be applied to parameter ISC
in all segments. (F10.0)
= 1, NOPAM
110
-------
Record 3--Segment Number (110)
ISG = model segment number to which the following parameter values
will be applied. (110)
Record 4-Segment Parameters (4(A5,15, F10.0))
PNAME(ISC) = an optional one to five alphanumeric character descriptive name
used to identify PARAM(ISG,ISC). (A5)
ISC
= parameter number identifying the parameter; all parameter
numbers are defined in Table 2.2. (15)
PARAM(ISEG,K) = the value of parameter ISC in segment ISG. (F10.0)
where K
ISEG
= 1, NOPAM
= 1,NOSEG
Organization of Records for Data Group G. Record 1 is input once and occupies one line.
Record 2 has NOPAM entries. Four entries will fit on one line; thus, Record 2 uses as many 80-
space lines as needed to enter all NOPAM entries. Records 3 and 4 are entered NOSEG times,
once for each segment. For each segment, Record 4 uses as many lines as needed to enter all
NOPAM entries.
2.2.8 Data Group H; Physicochemical Constants
Data Group H is used to enter all chemical specific physicochemical constants required
for a simulation. The number of constants needed for any simulation will vary as a function of
the type of water body modeled and the kinetic complexity of the simulation.
The structure of Data Group H varies. This data group is subdivided into constants that
are applied to individual chemicals or solids (systems) and global constants that are applied to all
systems (thus NOSYS+1 sub-groups of constants are read). Each sub-group can be further
subdivided into a user specified arbitrary number of sets that contain similar kinds of data that
allows the user to logically organize all Data Group H inputs.
Record 1-Header (SOX)
TITLE = name of data group. (SOX)
111
-------
Record 2-Data Fields in Group K (A10,110)
CHNAME(ISYS) = ten-character descriptive name for System ISYS. (A10)
NFLD
where ISYS
= number of groups of constants for System ISYS; 0 = no constants
for this System (skip to next System); the user may subdivide the
constants into any number of arbitrary groups. (110)
= l.NOSYS+1
Groups for IS YS=1-NOS YS are for input of constants that apply to
an individual system (chemical or solids. The group for
IS YS=NOS YS+1 is for input of global constants.
Record 3-Number of Constants in Group (A10,110)
FLDNAME = ten-character name identifying field of constants. (A10)
NCONS
= number of constants to be entered in this group; 0 = no constants
for this group (skip to next field). (110)
Record 4~Constants (2(A10,110, F10.0))
TNAME(ISC) = name identifying constant ISC. (A10)
ISC
= number identifying constant; all constant numbers are defined in
Tables 2.3-2.12. (110)
CONST(1SC) = value of constant ISC in units specified in Tables 2.3-2.12. (F10.0)
Record 5~Reaction Yields (110)
NYLDS = total number of reaction yields (110)
Reaction yields must be specified anytime a simulated chemical
undergoes a reaction through which it is transformed into a
chemical that is also simulated. If a chemical reacts and is
transformed into a state that is not being simulated, no reaction
yield will be specified.
112
-------
Record 6--Yield Constants (2(315, F10.0))
FROM = reacting chemical number; Chemicals 1 to N are Systems 4 to
(N+3). (15)
TO = product chemical number. (15)
The chemical number into which reacting chemicals are
transformed is never zero (to represent transformation to a
chemical state not simulated). If chemical reactions transform a
model state variable into a product that is not also being simulated,
no reaction yield is specified.
PROCESS = reaction process pathway number. (15)
1 = water column biodegradation
2 = benthic biodegradation
3 = alkaline hydrolysis
4 = neutral hydrolysis
6 = oxidation
7 = photolysis
8 = extra (user defined) reaction
Y = reaction yield constant in gram product chemical produced/gram
reacting chemical reacted. (F10.0)
Organization of Records for Data Group H. Record 1 is input once. Records 2 through 4
are input for NOSYS +1 groups, one for each modeled system (1 to NOSYS) and one for global
constants (NOSYS+1). For each group, Records 3 and 4 are entered NFLD times. For each
field, Record 4 uses as many lines as needed for NCONS entries (2 entries per line). Record 5 is
input once. Record 6 uses as many lines as needed for NYLDS entries (2 entries per line).
2.2.9 Data Group I; Kinetic Time Functions
Data Group I is used to enter all kinetic time functions required for a simulation. The
number of kinetic time functions needed for any simulation will vary as a function of the
simulation complexity. The input format for all time functions, however, is constant.
113
-------
Record 1 -Number of Time Functions (110, 70X)
= number of time functions used in the simulation. If no time
functions are specified, set NFUNC to zero and skip to Data Group
J- (110)
TITLE
= name of data group. (70X)
Record 2~Time Function Descriptions (A5,215)
ANAME(ISC) = an optional one to five alphanumeric character descriptive name
for the time function ISC(I). (A5)
ISC(I)
NOBRKOBC)
where I
= number identifying the time function; all time function numbers
are defined in Table 2.13. (15)
= number of breaks used to describe the time function ISC(I). (15)
= 1, NFUNC
Record 3-Time Functions (4(2F10.0))
VALT(K) = value of time function ISC(1) at time T(K). (F10.0)
T(K)
where K
= time in days. If the length of the simulation exceeds T(NOBRK),
the time function is assumed to be periodic and is repeated starting
at T(l) with a period equal to T(NOBRK). (F10.0)
= 1,NOBRK
Organization of Records for Data Group T. Record 1 is input once. Records 2 and 3 are a
set and repeated as a set NFUNC times. Within each set, Record 2 is input once and Record 3 is
repeated until all NOBRK entries are input (four VALK(K)-T(K) pair entries per line).
2.2.10 Data Group J: Initial Concentrations
Data Group J is used to specify the concentrations of each system in each model segment
at the start of the simulation. Conditions specified are the segment concentrations and densities
for the state variables at the simulation start time (ZDAY+ZHOUR/24+ZMIN/1440) specified in
Data Group A. If ICFL = 2 in Data Group A, initial concentrations are read from the restart file.
114
-------
Even if initial conditions are to be read from a restart file, it is still necessary to include Data
Group J in the input data set to allow proper read in of Data Group K.
Record 1-System Information (A40,15, F5.0, F10.0, 20X)
CHEML = system name (chemical or solid). (A40)
IFIELD
DSED
CMAX
TITLE
= solids field (3, 4, or 5) that transports this system in its pure or
sorbed form (15).
Previous versions of the WASP4 framework from which IPX was
derived were developed to simulate eutrophication. For the
eutrophication version of the framework, UKLbLD was used to
specify which the settling velocity field in Data Group D was used
to transport planktonic solids. For IPX simulations, set IFIELD to
3 for solids type 1, 4 for solid 2, and 5 for solid 3; set IFIELD to 1
for all chemical systems simulated.
= density of system in kg/L; choose 0.5-2.5 for solids; choose 0.0 for
chemicals. (F5.0).
= maximum concentration of system permitted during the simulation
in mg/L; if the computed concentration ever exceeds CMAX, the
simulated is aborted. (F10.0)
= name of data group. (20X)
Record 2-Mtial Conditions (3(A5, 2F10.0)
ANAME(ISEG) = an optional one to five alphanumeric character descriptive name or
number identifying segment ISEG. (A5)
C(ISYS,ISEG)
DISSF
= initial concentration of system ISYS in segment ISEG in mg/L.
(F10.0)
= dissolved fraction of chemical in segment K. (F10.0)
The value specified for DISSF is unimportant. The fractions of
chemical in each phase (dissolved, DOC bound, or paniculate) are
115
-------
where ISEG
ISYS
computed as a function of partition coefficients specified in Data
Group H. (If no partitioning constants are specified, all chemical
mass will be in the dissolved phase.)
= 1,NOSEG
= l.NOSYS
Organization of Records for Data Group J. Records 1 and 2 are a set and repeated as a set
NOSYS times. Within each set Record 2 is repeated until NOSEG entries are input (three
ANAME-C-DISSF entries per line).
2.2.11 Data Group K; Surficial Sediment Age Laver Conditions
Data Group K is used to specify the surficial sediment aging control parameters used
throughout the simulation and the initial distribution of sediments in each surficial sediment age
layer at the start of the simulation.
Record 1 - Sediment Aging Information (A40,2I5,F 10.0,1 OX)
SOLNAME = Data Group name for user reference. (A40)
NSURFSEDS
= number of surficial sediment segments in the model (segments
with ITYPE = 3 or 5). (K)
NSEDAGES = Number of age layers in the sediments. (15)
TSEDAGES = time for sediment aging computations, in days. (F10.0)
Record 2 - Sediment Age Layer Initial Conditions (A5,(NSEDAGES(F10.0))
BNAME(I) = identifier of surficial sediment segment I for user reference. (A5)
LAYERDEP(I,J)
depth of sediment in age layer J of segment I, in cm.
(F10.0). The sum of the nSedAges layerDep values should
equal the segment depth (d) defined for each surficial
sediment segment in Data Group C (volumes):
nSedAges
116
-------
where I
where J
= 1,NSURFSEDS
= l.NSEDAGES
Record 3 — Sediment Resuspendability Control Parameters (A5,F10.0)
CNAME(J) = Name of age layer J for user reference. (A5)
TDEPOSITION(J) = sediment resuspendability factor for age layer J (the value of
parameter Z in the epsilon equation, Equation 1.25). (F10.0)
where J
= l.NSEDAGES
Record 4 - Fresh Sediment Resusperision Control Parameters (3F10.0)
AO = fresh sediment resuspension coefficient (fresh sediments only).
(F10.0)
M
TAUCRIT
fresh sediment resuspension exponent (fresh sediments only).
(F10.0)
critical shear for fresh sediment resuspension, in dynes/cm^ (fresh
sediments only). (F10.0)
Organization of Records for Data Group K. Record 1 is input once. Record 2 is repeated
for each surficial sediment segment (ITYPE = 3 or 5) in the model with nSedAges initial age
layer depths appearing on each line. Record 3 is repeated for each sediment age layer. Record 4
is input once.
Example: Data Group K. Data Group K, from an IPX input data file used to model
sediment transport in the Fox River, is listed below to illustrate the specification of sediment
resuspension parameters:
Surficial Sediment Age Layer Conditions 48
Seg25 11 1 1
1.0
(lines for Segments 26 thru 71 deleted)
Seg72
1.
1
117
-------
Layrl
Layr2
Layr3
Layr4
LayrS
LayrS
Layr7
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0.008
3.
2.2.12—Data Group L; Initial Conditions for Deep Sediment Layers (Ghost Stack) rSemi-
Lagrangian Bed Option)
Data Group L is used to specify control parameters and initial conditions for sediment
compartments comprising the ghost stack using the semi-Lagrangian sediment bed option. If
there is at least one surficial sediment segment specified as ITYPE = 5, the semi-Lagrangian
sediment bed option is enabled and Data Group L must be entered. Data Group L is skipped only
if there are no segments specified as ITYPE = 5. "Note that if no segments are specified as
ITYPE = 5, then there should not be any segment specified as ETYPE = 6.
Record 1 ~ Minimum and Maximum Thickness (Volume) Control Parameters (2F10.0)
MINDEP = Minimum depth factor to trigger upward sediment re-indexing
MAXDEP
= Maximum depth factor to trigger downward sediment re-indexing
Record 2 — Ghost Stack Configuration Information (315)
SEGNUM = Surface sediment segment number (dummy variable just for user
reference). (15)
NGHOSTO(ISEG) = Number of elements in ghost stack under surface sediment segment
ISEG at the start of the simulation (do not include subsurface
segments specified in Data Group C). (15)
IGOPT(ISEG) = Ghost collapse option flag. (15)
0 = Ghost collapse option off. Ghosts do not collapse.
1 = Ghost collapse option on. Ghost collapse enabled.
118
-------
Record 3 - Ghost Stack Element Geometry (2F10.0)
AGO(ISEG,IG) = Surface area of ghost stack element IG under segment ISEG in m2.
(F10.0)
DGO(ISEGJG) = Thickness of ghost stack element IG under segment ISEG in
meters. (F10.0)
Record 4 — Ghost Stack Element State Variable Concentrations (8F1Q.O)
CG(ISYS,ISEG,IG)= Concentration of state variable ISYS in ghost stack element IG
under surface segment ISEG in mg/L. (F10.0)
where ISYS = l.NOSYS
IG = 1,NGHOSTO
ISEG = 1, NOSEG for all ISEG where ITYPE = 5
Organization of Records for Data Group L. Record 1 is input once. Record 2 is input
once for each semi-Lagrangian sediment stack (surficial sediment segment specified as ITYPE =
5) in the model with nSedAges initial age layer depths appearing on each line. Records 3 and 4
are repeated once for each ghost element in the stack. Record 4 takes as many lines as are
needed to input values for all state variables with a maximum of 8 values per line. If the initial
number of occupied ghost elements in any stack is zero, only Record 2 is input for that stack.
2.3 IPX MODEL PARAMETERS, CONSTANTS, AND TIME-FUNCTIONS
OVERVIEW
The environmental parameters, physicochemical constants, and time functions that may
be specified in Data Groups G, H, and I are described in the following sections. The number of
parameters, constants, and time functions that the user must specify depends on the complexity
with which chemical reactions and the interactions between chemicals and solids are to be
simulated. For convenience, the complexity necessary to represent environmental dynamics has
been conceptually divided into the five equilibrium and four kinetic complexity levels presented
in Table 2.1. These complexity levels are best conceptualized as points along a spectrum of
possible simulation complexities; not all chemicals simulated will, or need be, simulated at the
same complexity. However, more complex representations of environmental dynamics generally
require more constants than simpler representations.
119
-------
Table 2.1. Conceptual levels of simulation complexity.
Complexity Level
Equilibrium 1
Equilibrium 2
Equilibrium 3
Equilibrium 4
Kinetic 1
Kinetic 2
Kinetic 3
Kinetic 4
Explanation
Constant partition coefficient (Kp = constant)
Variable partition coefficients (Kp = f(x, y, z, t)
Hydrophobic partitioning (Kp = KOC foc)
Hydrophobic, solids-dependent partitioning
(Kp=.f(Koc,foc,mi))
Constant process half-life or rate constants
Variable process rate constants
Second-order process rates (rates depend on specified
characteristics of the modeled system)
Chemical reactions that yield transformation products
simulated chemical state variables
physical
that are
2,3.1 Environmental Parameters; Data Group G
The 16+ environmental parameters the user can specify are listed in Table 2.2. These
parameters permit the user to spatially vary the environmental properties that control the rates at
which simulated chemicals react. The user need input only those parameters required to
represent the particular reactions simulated. It is not necessary to specify any parameters unless
required to represent a specific mass transport pathway. The parameters specified in Data Group
G are coupled with the time functions in Data Group I to permit simultaneous spatial and
temporal variation of simulation conditions.
No parameters are necessary for Equilibrium Level 1. At this complexity level, all
partitioning coefficients for each chemical are constant and best specified in Data Group H. For
equilibrium level 2 FOC(ISEG,1), FOC(ISEG,2), and FOC(ISEG,3) can be used as
dimensionless multipliers of the partition coefficients input in Data Group H. For equilibrium
levels 3 and 4, FOC(ISEG,1), FOC(ISEG,2), and FOC(ISEG,3) are the fraction organic carbon of
each solids class. DOC(ISEG) may be entered to represent three-phase partitioning.
No parameters are necessary for Kinetic Level 1. At this complexity level, all reaction
rates for each chemical are constant and best specified in Data Group H. At kinetics level 2,
TOTKG(ISEG,1), TOTKG(ISEG,2), and TOTKG(ISEG,3) can be used as dimensionless
120
-------
Table 2.2. Data Group G environmental parameters.
lisc
i
2
3
4
5
6
7
8
9
10
11
12
13
PARAM (ISEGJSC)
VELFN (ISEG)
TMPFN (ISEG)
TEMP (ISEG)
WVEL (ISEG)
REAR (ISEG)
DOC (ISEG)
FOC(ISEG,1)
FOC (ISEG,2)
FOC (ISEG,3)
CHPHL (ISEG)
PH (ISEG)
XKE2 (ISEG)
OXRAD(ISEG)
Definition, Units, Reactions Affected
Pointer to water velocity time function (1-4) that will be
specified in Data Group I; V.
Pointer to water temperature time function (1-4) that will be
specified in Data Group I; ALL.
Water temperature in segment (°C) and/or multiplier for the
selected water temperature time function (TEMPN); ALL.
Wind velocity over segment (m/s) and/or multiplier for
WIND time function; V.
Oxygen reaeration in segment (m/d) and/or multiplier of
REARN time function; usage depends on volatilization
option selected; V.
Dissolved organic carbon concentration in segment (mg/L);
S,P.
Fraction organic carbon of solids type 1 in segment and/or
multiplier of MFOCd) time function; S.
Fraction organic carbon of solids type 2 in segment and/or
multiplier of MFOC(2) time function; S.
Fraction organic carbon of solids type 3 in segment and/or
multiplier of MFOC(3) time function; S .
Chlorophyll concentration in segment and/or multiplier for
CHLN time function (mg/L); P.
pH in segment and/or multiplier for pH time functions
(PHNW and PHNS); H.
Light extinction coefficient of photochemically active light
in segment (1 /meter); P.
Concentration of oxidants in segment, such as O3 for H2O2
(moles/L); O.
S = sorption; V = volatilization; B = biodegradation; H = hydrolysis; O = oxidation; P
photolysis; E = extra (user defined) reaction
121
-------
Table 2.2 (continued). Data Group G environmental parameters.
ISC
PARAM (ISEGJSC)
Definition, Units, Reactions Affected
14
BAG (ISEG)
Density of active bacteria (cells/100 cc) in segment and/or
multiplier of bacteria time functions (BACNW and
BACNS). Units for BAG must be consistent with those of
KBIO20 (constants 146-160); the product of BAG and
KBIO20 must yield units of day-1; B.
15
EXENV (ISEG)
Environmental property that affects the user-defined "extra
reaction." Units for EXENV must be consistent with those
of KE20 (constants 576-590); the product of EXENV and
KE20 must yield units of day1; E.
16-N
TOTKG (ISEG,N)
Total lumped first-order decay rate constant for chemical (N-
15) (day-1).
where N= 15 + number of chemicals simulated
S = sorption; V = volatilization; B = biodegradation; H = hydrolysis; O = oxidation; P
photolysis; E = extra (user defined) reaction
multipliers of the overall reaction rates input in Data Group H. Kinetics level 3 may. require the
other parameters, depending on the kinetic processes of importance. For example, if water
temperatures differ significantly from 20°C, then it may be necessary to input TEMP(ISEG) to
properly represent reaction processes. If volatilizing chemicals are simulated, it may be
necessary, or required, to input REAR(ISEG), DEPTH(ISEG), VELOC (ISEG) or WVEL
(ISEG), etc. 'In general, if a chemical reaction pathway is a function of a physical property of the
water body, it may be necessary to specify one or more parameters.
2.3.2 Physicochemical Constants: Data Group H
The physicochemical constants the user can specify are listed in Tables 2.3-2.12. These
constants permit the user to specify chemical-specific properties that control the rates at which
simulated chemicals react. The user need input only those constants required to represent the
particular equilibria and reactions simulated. It is not necessary to specify any parameters for any
solids class simulated. For each chemical it is almost always necessary to specify the molecular
weight (constant 81). Beyond the molecular weight, it is not necessary to specify any additional
constants unless required to represent a specific mass transport or reaction pathway. The
122
-------
Table 2.3. Partition coefficients and rate constants for simple reactions.
Constant
=====
111
116
121
140
141
142
181
182
1183
256
287
571
143
144
252
253
J254
257
*"* —
1 289
1 572
Variable
=
PIXC(1,1)
PKC(2,1)
PKC(3,1)
Ki
KV
KBW
KBS
KHOH
KHN
KHH
KO
KF
KE
THi
THBW
THBS
THHOH
THHN
THHH
THO
THF
THE
Constant partition coefficient (Kpl) for sorption to
class l,Lw/kgs
Constant partition coefficient (Kp2) for sorption to
class 2, Lw/kgs
Constant partition coefficient (Kp3) for sorption to
class 3, Lw/kgs
solids
solids
solids
General first-order loss rate constants, day1
Volatilization
Water column biodegradation
Benthic (sediment) biodegradation
Alkaline hydrolysis
Neutral hydrolysis
Acid hydrolysis
Oxidation
Photolysis
Extra (user defined) reaction
dpne.ra.1 reaction half-life, day
Water column biodegradation J
Benthic biodegradation _|
Alkaline hydrolysis
Neutral hydrolysis
Acid hydrolysis
Oxidation J
Photolysis J
Extra (user defined) reaction 1
123
-------
Table 2.4. General chemical constants.
Constant
3
4
9
81
82
83
84
Variable
VB
TMELT
TDINT
MOLWT
SOLG
VAPRG
LKOW
Definition
LeBas molar volume, cm^/mol
Melting temperature, "C
Time interval at which reaction rates are recomputed,
days
Molecular weight, g/mole
Solubility, mg/L
Vapor pressure, atm
Log octanol-water partition coefficient, Loct/Lwater
Table 2.5. Sorption constants.
Constant
84
101
102
103
106
111
116
121
Variable
LKOW
LKOC
AO
Al
NUX(l)
PKC(1,1)
PKC(2,1)
PKC(3,1)
Definition
Log octanol-water partition coefficient, Lwater/L0ct
Log of the organic carbon partition coefficient,
Lwater/kgoc; for hydrophobic solids dependent
partitioning, this is log Koc at zero solids concentration.
Intercept in the KQW - KQC correlation: log Koc = AQ; log
KOW + AI; default AQ = log 0.6
Slope in the Kow - Koc correlation; default =1.0
Solids-dependent partitioning parameter (vx) of the
chemical onto solids; default = 1012 makes Kp
independent of solids concentration
Solids-independent (limiting) partition coefficient Kpo for
sorption to solid 1, (Lwater/kgsolids); if PKC is not
specified, Kpo for the chemical will be calculated from
LKOC and parameter FOC
Solids-independent (limiting) partition coefficient Kpo for
sorption to solid 2, (Lwater/kgsolids)
Solids-independent (limiting) partition coefficient Kpo for
sorption to solid 3, (Lwater/kgsoiids)
124
-------
Table 2.6. Hydrolysis constants.
Constant
603
Variable
XHYDRO
Definition
Hydrolysis Option:
0 = No hydrolysis
1 = Hydrolysis occurs
Neutral pH(pH~ 7.0)
184
201
206
211
236
TREFH
KH2O(2,1,1)
KH2O(2,2,1)
KH2O(2,3,1)
EHN(l)
Reference temperature at which hydrolysis rates were
measured, °C
20°C neutral hydrolysis rate constant for dissolved
chemical, day1
20°C neutral hydrolysis rate constant for DOC-bound
chemical, day1
20°C neutral hydrolysis rate constant for particulate
chemical, day1
Activation energy for neutral hydrolysis, kcal/mole
AcidpH(pH<7)
216
221
226
241
KH20(3,1,1)
KH20(3,2,1)
KH20(3,3,1)
EHH(l)
Second-order, 20°C acid hydrolysis rate constant for
dissolved chemical, L/mole-day
Second-order, 20°C acid hydrolysis rate constant for
DOC-bound chemical, L/mole-day
Second-order, 20°C acid hydrolysis rate constant for
particulate chemical, L/mole-day
Activation energy for acid hydrolysis, kcal/mole
Alkaline pH (pH > 7)
186
191
196
231
KH20(1,1,1)
KH2O(1,2,1)
KH20(1,3,1)
EHOH(l)
Second-order, 20°C alkaline hydrolysis rate constant
for dissolved chemical, L/mole-day
Second-order, 20°C alkaline hydrolysis rate constant
for DOC-bound chemical, L/mole-day
Second-order, 20°C alkaline hydrolysis rate constant
for particulate chemical, L/mole-day
Activation energy for alkaline hydrolysis, kcal/mole
125
-------
Table 2.7. Photolysis constants.
Constant
286
291
551
556
561
Variable
XPHOTO
KDPG(l)
QUANTG(l)
QUANTG(2)
QUANTG(3)
Definition
Photolysis option:
0 = no photolysis
1 = depth-integrated photolysis rate incut bv user
Depth-integrated photolysis rate (m/day)
Quantum yield fraction of dissolved chemical
Quantum yield fraction of paniculate chemical
Quantum yield fraction of DOC-bound chemical
Table 2.8. Global constants (generally not used at this time).
Constant
1
3
4
6
7
11-23
24-36
37-49
50-62
63-75
Variable
TO
ELEVG
LATG
XLITE
DFACG
CLOUDG(l)
AIRTYG(l)
RHUMG(l)
ATURBG(l)
OZONEG(l)
Definition
Julian date at beginning of simulation
Average ground surface elevation, m
Latitude of water body, degrees
Water surface light intensity option; 0 = do not
compute light; 1 = annual average; 2 = average for
month indicated by TO; 3 = monthly step function
Ratio of optical path length to vertical depth: 1 17
Mean monthly cloudiness, in tenths of full sky
coverage (0-10)
Mean monthly air mass type; 1 = rural, 2 = urban, 3 =
maritime, 4 = tropospheric
Mean monthly daylight relative humiditv. percent
Mean monthly atmospheric turbidity, in equivalent
aerosol layer thickness km
Mean monthly ozone content of atmosphere, in cm
NTP (0.2 - 0.3)
126
-------
Table 2.9. Oxidation constants.
Constant
604
258
261
266
271
276
Variable
XOXID
TREFO
KOX20(1,1)
KOX2O(2,1)
KOX2O(3,1)
EOX(l)
Definition
Oxidation Option:
0 = No oxidation
1 = Oxidation occurs
Reference temperature at which oxidation rates were
measured, °C
Second-order, 20°C oxidation rate constant for
dissolved chemical, L/mole-day
Second-order, 20°C oxidation rate constant for DOC-
bound chemical, L/mole-day
Second-order, 20°C oxidation rate constant for
particulate -chemical, L/mole-day
Activation energy for oxidation, kcal/mole
Table 2.10. Biodegradation constants.
Constant
602
146
151
156
161
166
171
Variable
XBIO
KBI020(1,1)
KBI020(2,1)
KBI020(3,1)
QlODIS(l)
QlODOC(l)
QlOPAR(l)
Definition
Biodegradation Option
0 = No biodegradation
1 = Biodegradation occurs
Second-order 20°C biodegradation rate constant for
dissolved chemical, mL/cells-day
Second-order 20°C biodegradation rate constant for
DOC-bound chemical, mL/cells-day
Second-order 20°C biodegradation rate constant for
particulate chemical, mL/cells-day
Temperature correction factor for dissolved chemical
biodegradation; multiplication factor for 10°C
temperature increase
Temperature correction factor for DOC-bound
chemical biodegradation; multiplication factor for
10°C temperature increase
Temperature correction factor for particulate chemical
biodegradation; multiplication factor for 10°C
temperature increase
127
-------
Table 2.11. Volatilization constants.
Constant
Variable
Definition
136
XVL
Volatilization option for liquid phase mass transfer:
no volatilization
Owens
Modified O'Connor-Dobbins
Churchill
O'Connor "short" form
Mackay
Reaeration rate input by user and converted to a
liquid phase mass transfer rate
Volatilization rate directly input by user
O'Connor "lone" form
0
1 =
2 =
3 =
4 =
5 =
6 =
7 =
8 =
XVG
Volatilization option for gas phase mass transfer:
0 =
1 =
2 =
3 =
4 =
5 =
6 =
No gas phase resistance/no volatilization
(default if XVL=0)
O' Connor-Rathbun
Mills et al.
O'Connor
Mackay
Ambrose et al. for flowing systems
Mackay still-air
7 = Liss
XVH
Henry's Law constant option:
0 = HLC is a chemical constant (default if XVL=0)
1 = Direct expression of temperature
2 = Calculated from temperature-dependent vapor
pressure and solubility
137
HLC25
Henry's Law constant at 25°C, atm-m3/mole
600
HLCAO
Henry's Law constant temperature factor ao,
dimensionless. Used in temperature-based HLC
estimation expression:
InH = an —
a,
T + 273.15
128
-------
Table 2.11 (continued). Volatilization constants.
Constant
601
138
139
5
8
Variable
HLCA1
KLT
KVOG
AIRTMP
ATMOS
Definition
Henry's Law constant temperature factor ai,
dimensionless
Volatilization temperature correction factor,
dimensionless
Measured ratio of volatilization to reaeration rates
Multiplier for air temperature time function, °C or
normalized (dimensionless)
Atmospheric concentration of chemical, j-ig/L
Table 2.12. "Extra reaction" constants.
605
573
576
581
586
591
Variable
XERXN
TREFE
KE2O(1,1)
KE20(2,1)
KE2O(3,1)
EEX(l)
Definition
Extra Reaction Option:
0 = No extra reaction
1 = Extra reaction occurs
Reference temperature at which extra reaction rates were
measured, °C
Second-order, 20°C extra reaction rate constant for
dissolved chemical, l/[E]-day
Second-order, 20°C extra reaction rate constant for DOC-
bound chemical, l/[E]-day
Second-order, 20°C extra reaction rate constant for
particulate chemical, l/[E]-day
Activation energy for extra reaction, kcal/mole
129
-------
constants specified in Data Group H are coupled with the parameters in Data Group G and time
functions in Data Group I to permit simultaneous spatial and temporal variation of chemical-
specific equilibria and reactions.
There are many possible combinations of constants, parameters, and time functions that
can be specified to represent simulated equilibria or kinetics. As a result it is not possible to
enumerate every possible option for using each constant. However, the definitions and use of the
range of constants in IPX is consistent with typical environmental engineering usage as described
in references texts such as those by Thomann and Mueller (1987) and Chapra and Reckhow
(1983), and Chapter 1 of this manual. However, users are strongly advised to review the IPX
source code as well as Tables 2.3-2.12 to ensure that constants for desired equilibria and
reactions are properly specified. With this caveat in mind, a brief overview of commonly used
constants follows.
At equilibrium level 1, it is only necessary to also input a partition coefficient to represent
partitioning to each sorbing solids class using PIXC (constants 111, 116, and 121). For
equilibrium level 2, in addition to specifying PIXC, the user may use Data Group G parameters
numbers 6-9 and Data Group I time functions 19-21 to permit spatial and temporal variation of
partition coefficients; at this level, the parameters and time function should be treated as spatial
and temporal, dimensionless scale factors of PIXC. At equilibrium level 3, the user must specify
LKOC (constant 101), or other constants required to approximate LKOC, and enter the FOC-
related parameters and time functions in Data Groups G and I. At equilibrium level 4, the user
must specify NUX (constant 106) in addition to those inputs required for equilibrium level 3.
As previously noted, it is difficult to enumerate all possible options for specifying
constants. Although this is particularly true for representing chemical reaction kinetics, some
general guidelines are provided. At kinetic level 1, reaction rates are constant and can be entered
as either reaction half-life (t1/2) or rate constants (Ki). Reaction half-life or rates for each
chemical may be specified for each relevant process (hydrolysis, oxidation, biodegradation,
volatilization, photolysis, or user-defined) or lumped into a single "decay" constant. In general,
it is preferable to specify reaction rates for each relevant reaction pathway individually. At
kinetic level 2, reaction rates can be spatially and temporally varied through the use of Data
Group G parameters and Data Group I time functions. At kinetics level 3, reaction kinetics are
"second-order" and can depend on some characteristic of the water body (such as ozone or
bacterial concentrations) in addition to the specified reaction rate constant (which itself can be
spatially and temporally variable). Reactions need not be second-order to be simulated. For
130
-------
example, consider a simple oxidation reaction; if the concentration of the oxidant is constant, it
could be "lumped" with the rate constant and represented as an overall "first-order" reaction. At
kinetic level 4, reaction transformation products are simulated and yield coefficients for each
relevant transformation process must be specified.
To specify volatilization, the user must generally choose both a liquid and gas phase
resistance option, and a Henry's Law Constant option. If XVL is specified but XVG is not, the
gas phase resistance defaults to zero when computing the overall volatilization rate. If XVG is
specified but XVL is not, XVL has a default value of zero and no volatilization is computed for
the chemical. If XVH is not specified, then HLC will be treated as a constant, either specified by
input of HLC25 or calculated from VAPRG and SOLG. The user is reminded to verify that the
volatilization options specified are appropriate for the range of environmental conditions
encountered. Liquid and gas phase mass transfer options are described in Table 2.11. Some
caveats regarding the selection of volatilization options in conjunction with parameter and time
function requirements follow.
To use XVG options 1-5 and 8, the user must specify both WINDN, time function 9 in
Data Group I, and WVEL, parameter 4 in Data Group G of the input data set. The variable
SWIND is the product of WVEL and WINDN and is used for XVG options 1-4 and 7 (SWIND is
used to compute USTAR for options 3, 4, and 7). The default value for WVEL is zero and the
default for WINDN is 1. Unless both values are specified, SWIND will have a value of zero.
The use of REARN, tune function 12 in Data Group I, and REAR, parameter 5 in Data
Group G of the input data set, depend on the XVL option chosen. For XVL options 1-5 and 8,
REARN does not directly enter into the volatilization computations but is used to control
whether volatilization occurs or not. For XVL options 6-7, the internally computed variable
SREAER, defined as the product of REARN and REAR, enters directly into the computations.
For XVL options 1-5 and 8, REARN is used to toggle the volatilization computation on
or off to represent conditions such as ice cover. To use REARN as a toggle, the user should
specify a value of 1 (any non-zero value will work) for ice free periods and zero for periods of ice
cover. No volatilization is computed for any period when the value of REARN is zero. The
default value of REARN is 1. So, if REARN is not specified for XVL options 1-5 and 8,
volatilization computations will always be computed (ice free conditions). The parameter REAR
does not need to be specified for XVL options 1-5 and 8.
131
-------
For XVL option 6, the use of REARN depends on the value specified for KVOG. If
KVOG is specified as a chemical-specific reaeration rate or liquid phase mass transfer
coefficient, REARN is a dimensionless time function through which the reaeration rate is varied
with time. If KVOG is given a value of 1, REARN is a time series of reaeration/mass transfer
coefficients (units in m/day). The user should note that if REARN is used to specify a time series
of mass transfer coefficients, the same coefficient value will be applied to all chemicals
volatilizing. For either case, REAR is a segment-specific normalization constant through which
the liquid phase mass transfer coefficient can be varied segment by segment. The default value
of REARN is 1. So, if REARN is not specified the liquid mass transfer coefficient will not vary
with time. The parameter REAR must be specified for this option. The default value of REAR
is zero. So, if REAR is not specified the value of SREAER and the resultant liquid phase mass
transfer coefficient will be zero and a divide by zero error will occur.
For XVL option 7, the use of REARN depends on the value specified for KVOG. If
KVOG is specified as a chemical-specific volatilization rate, REARN is a dimensionless time
function through which the volatilization is varied with time. If KVOG is given a value of 1,
REARN is a time series of volatilization rates (units in m/sec). The user should note that if
REARN is used to specify a time series of volatilization rates, the same rate will be applied to all
chemicals volatilizing. For either case, REAR is a segment-specific normalization constant
through which the volatilization rate can be varied segment by segment. The default value of
REARN is 1. So, if REARN is not specified the volatilization rate will not vary with time. The
parameter REAR must be specified for this option. The default value of REAR is zero. So, if
REAR is not specified the value of SREAER and the resultant volatilization rate will be zero and
no volatilization will occur.
Potentially all volatilization options and rates depend on temperature. Specifically, the
contaminant diffusivity in water and air, air-water partition coefficient (including the Henry's
Law constant), and the density and viscosity of water and air are functions of temperature. To
account for temperature variation in volatilization, the water and air temperature must be
specified.
The water temperature, TMP, is internally computed as the product of a water
temperature time function TEMPN(l-4), time functions 1-4 in Data Group I, TMPFN, and
TEMP, parameters 2 and 3 in Data Group I. Up to four different water temperature time
functions can be specified. Unless environmental conditions suggest otherwise, it is generally
not necessary to specify more than one water temperature time function. TMPFN is a pointer
132
-------
that indicates which time function (1-4) is used to compute water temperature for a given
segment. TEMP is a segment-specific normalization constant through which the water
temperature for a given time function can be varied segment by segment. The default value of
TEMP is zero and the defaults of TMPFN and TEMPN are 1. So, unless all three values are
specified, TMP will have a value of 0°C (273.15 K).
The air temperature, ATMP, is internally computed as the product of AIRTMPN, time
function 13 in Data Group I, and AIRTMP, constant 5 in Data Group H of the input data set. The
default value of AIRTMP is zero and the default of AIRTMPN is 1. So, unless both values are
specified, ATMP will have a value of 0°C (273.15 K).
Due to the large number of possible XVL and XVG combinations, it is not feasible to
specifically list every possible combination or use of volatilization options that may require
specification of the water and/or air temperature. The user is advised to examine the IPX source
code to confirm whether specification of the water and/or air temperature is required for the
volatilization options chosen.
2.3.3 Kinetic Time Functions: Data Group I
The 21 kinetic time functions the user can specify are listed in Table 2.13. These time
functions permit the user to temporally vary the environmental properties that control the rates at
which simulated chemicals react. The user need only specify those functions required to
represent the particular reactions simulated. It is not necessary to specify any time functions
unless required to represent a specific mass transport pathway. The parameters specified in Data
Group G are coupled with the time functions in Data Group I to permit simultaneous spatial and
temporal variation of simulation conditions.
For simple simulations, no time functions need be specified. For simulations in which
contaminants are transferred between media (air-water, water-sediment) or react at variable rates,
time functions for each relevant process may be specified. TEMPN can affect all reactions.
Volatilization may require VELN(l-4), REARN, WINDN and AIRTMPN. Photolysis may
require PHTON. Hydrolysis may require PHNW and PHNS. Biodegradation may require
BACNW and BACNS. Values for functions not specified default to 1.0.
133
-------
Table 2.13. Data Group I kinetic time functions.
ISC
1
2
3
4
5
6
7
8
9
10
11
12
ANAME(ISC)
TEMPN(l)
TEMPN(2)
TEMPN(3)
TEMPN(4)
VELN(l)
VELN(2)
VELN(3)
VELN(4)
WINDN
«
PHNW
PHNS
REARN
VALT(ISC)
Water temperature function 1. Values for
functions TMPN(l-4) can be either normalized
(dimensionless) or actual temperatures in °C,
depending upon the use of parameter
TEMP(ISEG) specified in Data Group G.
Water temperature function 2, normalized or °C.
Water temperature function 3, normalized or °C.
Water temperature function 4, normalized or °C.
Water velocity function 1, m/sec. The velocities
input through functions VELN(l-4) are added to
the net water velocities VELOCG(ISEG)
computed from flow and the hydraulic parameters
specified in Data Group C.
Water velocity function 2, m/sec.
Water velocity function 3, m/sec.
Water velocity function 4, m/sec.
Wind velocity function, m/sec. This time function
is multiplied by the segment specific wind
parameter WVEL(ISEG) specified in Data Group
G.
Water column pH function; pH units or
dimensionless (normalized). This function is
multiplied by the segment specific pH parameter,
PH(ISEG), entered in Data Group G.
Benthic pH function; pH units or normalized
(dimensionless). This function is multiplied by
the segment pH parameter, PH(ISEG), for benthic
segments.
Reaeration coefficient, day1 or normalized. This
function is multiplied by the segment specific
reaeration parameter REAR(ISEG) entered in the
Data Group G.
134
-------
ISC
13
Table 2.13 (continued). Data Group I kinetic time functions.
ANAME(ISC)
.-s=^=^^=^^^=s
AIRTMPN
VALT(ISC)
Air temperature, °C or normalized. Used for
computing reaeration/volatilization.
14
CHLN
Phytoplankton chlorophyll concentration, mg/L.
This variable is multiplied by the segment specific
chlorophyll parameter CHPHL(ISEG) entered in
Data Group G.
15
PHTON
~
Normalized light intensity, dimensionless. This is
used for photolysis option 1 to adjust the
measured rate constant under controlled light
intensity to a predicted rate constant under
ambient light intensity.
16
BACNW
Time variable bacteria concentration in the water
column, mg/L. This is multiplied by the segment
specific parameter BAG entered in the parameter
section. • ___
17
BACNS
Normalized benthic bacteria function,
18
19
20
I21
ATCHEM
MFOC(l)
MFOC(2)
MFOC(3)
dimensionless. This is multiplied by the segment
bacteria parameter B AC(ISEG) for benthic
segments.
Normalized atmospheric chemical concentration
function, dimensionless. This is multiplied by the
atmospheric chemical concentration constant
ATMOS for each chemical.
Time variable fraction organic carbon content for
solids type 1 in the water column. This is
multiplied by the segment foc parameter
FOC(ISEG,1). Applies to the water column only.
Time variable fraction organic carbon content for
solids type 2 in the water column. This is
multiplied by the segment foc parameter
FOC(ISEG,2). Applies to the water column only.
Time variable fraction organic carbon content for
solids type 3 in the water column. This is
multiplied by the segment foc parameter
FOC(ISEG,3). Applies to the water column only. |
135
-------
CHAPTER 3
IPX SIMULATION OUTPUT AND POST-PROCESSING
3.1 IPX OUTPUT OVERVIEW
IPX produces four output files that hold a variety of simulation results. These files use
the root file name of the input data file with a unique identifying extension. The most important
of these is the *.DMP file (where * is the root name of the input data file), a direct access file that
contains all kinetic display variables for each state variable for each model segment at each print
interval throughout the simulation. Other files created by an IPX simulation include *.DMA,
*.EXP, *.MSB, and *.RST, *.OUT.
3.2 THE *.DMP FILE
The *.DMP file is a direct-access file in which instantaneous (point-in-time) simulation
results are stored. Results are stored for all water column and sediment segments as well as all
elements in the ghost stack (if any). The results stored in this direct access file can be retrieved
with the W4DIS274 post-processing program. The W4DIS274 program prompts the user to
input desired display options and creates a sequential-access (ASCII text) file that contains the
specified simulation results. Display options for solids include the total solids concentration and
the concentration of each solids type. Display options for chemicals include the total (sum of all
phases) chemical concentration, the dissolved, DOC-bound, and paniculate concentrations of
each chemical on a mass per volume basis, and the particulate chemical concentration on a mass
per mass basis. Descriptions of the W4DIS274 post-processing program and the available
display options are presented in Section 3.8.
3.3 THE *.DMA FILE
The *.DMA file is a direct-access file in which time averaged simulation results are
stored. Results are stored for all water column and sediment segments. No results for the
sediment ghost stack are stored in the *.DMA file. The results stored in this direct access file can
be retrieved with the W4DIS274 post-processing program. The process for retrieving simulation
results from the *.DMA file using the W4DIS274 post-processor are identical to the process used
for the *.DMP file.
136
-------
3.4 THE *.EXP FILE
The *.EXP file is a sequential-access file that contains the times series of simulated
export for total solids and any one additional state variable (solid or chemical) defined in Data
Group A. Export is the mass transported by advection and dispersion across the system interface
defined in Data Group A. The values displayed on each line of the *.EXP file are the cumulative
masses exported across the defined interface since the previous print interval.
3.5 THE *.MSB FILE
The *.MSB file contains a mass balance record for all state variables (solid or chemical)
simulated. The mass balance is computed for all segments in the model network over the entire
simulation period. In this file, the accumulated masses transported into and out of each model
segment by each fate pathway are reported in kilograms. The fate pathways include loading,
advection, dispersion, pore water transport, settling (in), resuspension (in), burial (out), scour
(out), biodegradation, hydrolysis, oxidation, photolysis, and volatilization. Also reported are the
total resident mass in each segment and the estimated residual (unaccounted for) mass, In
addition, a report of final conditions in the sediment ghost stack is reported element by element
for each sediment stack. Reported conditions include: ghost volume, surface area, thickness,
elevation within the sediment bed, and the concentration of each state variable.
3.6 THE *.OUT FILE
The *.OUT file contains an "echo" of all information in the model input data file and a
record of any input data structure or simulation errors. As an input file is read at the start of an
IPX simulation, all model input data are written to the *.OUT file. If any input data file structure
errors are found at the time the input file is read, model execution stops and a message is written
to the *.OUT file at the point where the error was discovered. IPX input data files can be very
long and difficult to read; the *.OUT file provides a simple means for the user to confirm that all
input data are correctly read by IPX. If conditions become unstable during a simulation, model
execution stops and an error message is written to the *.OUT file. However, the user should be
aware that IPX does "trap" arithmetic faults such as division by zero errors; arithmetic faults that
cause simulation execution to abort may not be reported.
3.7 THE *.RST FILE
The *.RST file contains a "snapshot" of the volume and concentrations of each system in
every model segment. Writing data to and/or reading data from the *.RST file is controlled
137
-------
through the variable ICFL that is input in Data Group A of the input file. Values are written to
the *.RST file every TRST days and at the end of the simulation, replacing (overwriting) any
previous values (if any) stored in the file. This file can be read by IPX to conduct a continuous
series of simulations or resume a long simulation simulation that did not complete.
3.8
W4DIS274: THE *.DMP AND *.DMA FILE POST-PROCESSOR
The W4DIS274 post-processing program is used to retrieve simulation output from the
direct-access *.DMP and *.DMA files. Display variables for water column and sediment
segment retrievals are listed in Table 3.1. Display variables for sediment archive stack retrievals
are listed in Table 3.2. The sediment archive stack is defined below. The results of W4DIS274
retrievals are stored in a tab-delimited, ASCH text format in a file with the default name *.TBL.
To retrieve output for water column and sediment segments, the user must differentiate
between chemical-independent and chemical-dependent display variables. Display variables 1-7
(include solids concentrations, segment volume, segment depth, and segment type) are chemical-
independent Display variables 8-19 are chemical-dependent.
When prompted for a display option list, the user must precede each chemical-dependent
display variable number with the chemical number and a colon (:). For convenience, ranges of
chemicals, display variables, and segment numbers can be specified using square brackets ([, ])
and colons (:). For clarity, four examples follow:
Example 1 - Specifying chemical-independent and chemical-dependent display variables.
1. Total Solid
2. Solid 1
3. Solid 2
4. SolidS
5. Segment Volu
6. Segment Dpth
7. Segment Type
8. Total Chm
9. Dis. Chm
10. Doc Chm
11. Part Chm
12. PartChms
13. Ionic Chm
14. BiodegChm
15. HydrolChm
16. PhotlChm
17. VolatChm
18. OxidChm
19. XtraChm
Enter display variable list -> 1, 5, 6, 1:8, 2:8
138
-------
Table 3.1. W4DIS274 display variables.
Number
===^===
o
3
4
f-
7
9
10
11
12
13
14
15
16
17
18
19
Display Variable
Total Solid
Solid 1
Solid 2
Solid 3
Segment Volume
Segment Depth
Segment Type
Total Chm
Dis. Chm
Doc Chm
Part Chm
Part ChmS
Ionic Chem
Biodeg Chm
Hydrol Chm
Photl Chm
Volat Chm
Oxid Chm
Xtra Chm
Total solids concentration, mg/L
Solids type 1 concentration, mg/L
Solids type 2 concentration, mg/L
Solids type 3 concentration, mg/L
Segment volume, m3
Segment depth (thickness), m
Segment type (1, 2, 3, 4, 5, or 6)
Total chemical concentration, ug/L
Dissolved chemical concentration, ng/L
DOC-bound chemical concentration, fig/L
Particulate chemical concentration, ng/L
Particulate chemical concentration on a mass basis,
Ug/kgtotal solids '• —
Total ionic chemical concentration, |ig/L
Biodegradation rate constant, day"1
Total hydrolysis rate constant, day
Photolysis rate constant, day"
Volatilisation rate constant, day l
OvirlfltioTi rate constant, day"1
Extra rate constant, day"1
This example input will select total solids, segment volume, segment depth, total
concentration of chemical 1 and total concentration of chemical 2.
Example 2 - Specifying ranges of chemical-dependent display variables with the [ ] and :
separators^ _
Enter display variable list -> 1,[1:3]:8,4:[9:12]
This example will select total solids, total chemical concentration for chemicals 1, 2, and
3, and dissolved, DOC-bound, paniculate and sorbed particulate concentrations for chemical 4.
139
-------
Example 3 - Specifying ranges of chemical-dependent display variables with the [ ] and :
separators.
Enter display variable list-> [1:10]:[8:19]
This example will select all chemical-dependent display variables [8:19] for chemicals 1
through 10 [1:10].
Example 4 - Specifying ranges of segments with the [ ] and: separators.
Enter in segment numbers -> 1,2,3,[10:17]
This example will select segments 1, 2, 3,10, 11, 12, 13,14, 15, 16 and 17.
The sediment stack archive is a report of conditions for each sediment stack including the
surficial sediment segment, all specified subsurface sediment segments, and all ghost elements in
a stack. As data are retrieved, each element in the archive stack is assigned a layer number
according to its vertical position above the stack bottom. The bottom-most element of the stack
is assigned the largest possible layer number value. Each layer above that is assigned the next
largest layer number. This process continues until all ghost stack and subsurface sediment layers
are assigned a layer number. Unoccupied elements of the sediment stack will appear as zeros
(nulls) over all occupied stack elements. The surficial sediment segment is always assigned a
layer number of zero (0). Users should be aware that unoccupied stack elements will appear as
empty layers between the last occupied ghost/subsurface sediment element and the surficial
sediment element. Empty layers contribute zero to the total elevation of the sediment stack. If
the surface sediment element is unoccupied, this indicates that no sediment is present. If any
sediment is present at a location, at least the surficial sediment element of the stack (Layer 0) will
be occupied.
Retrieving output for the sediment stack archive follows procedures similar to those used
for the water and sediment column. To reach the menu for archive stack retrievals, the user must
select Option 5 from the main menu. To retrieve output from the archive stack, the user must
again differentiate between chemical-independent and chemical-dependent display variables.
Display variables 1-8 (element volume, surface area, thickness, and solids concentrations) are
chemical-independent. Display variable 9 is chemical-dependent. When prompted for a display
option list, the user must precede each chemical-dependent display variable number with the
140
-------
Table 3.2. W4DIS274 sediment archive stack display variables.
1
2
3
4
e
6
7
Q
9
Volume
Surface Area
Thickness
Elevation
Total Solid
Solid 1
Solid 2
SolidS
Total Chm
Definition
Volume of stack element, m3
Surface area of stack element, m
Thickness (depth) of stack element, m
Distance from the bottom of. the sediment stack to the
top of the layer, m
Total solids concentration, mg/L
Solids type 1 concentration, mg/L
Solids type 2 concentration, mg/L
Solids type 3 concentration, mg/L
Total chemical concentration, ug/L
chemical number and a colon (:). As previously described, ranges of chemicals, display
variables, and segment numbers can be specified using square brackets ([, ]) and colons (:).
141
-------
CHAPTER 4
IPX INPUT DATA FILE PRE-PROCESSING
4.1 IPX Input Data File Pre-Processing Program Overview
A difficult but critical aspect of developing a water quality model for any river is to
estimate settling and resuspension velocities for each solids type simulated. Four input data pre-
processing programs are provided to help simplify this task: REDUCE, SETTLE, RESUSPND,
and TIMESTEP. These pre-processing programs are designed to help the user define model
inputs for the simulation hydrograph and also settling and resuspension velocities and numerical
integration time-steps as a function of the hydrograph. A description of each pre-processing
program follows. However, the user is advised to examine the FORTRAN source code for these
programs to fully understand their proper operation and capabilities.
4.2 REDUCE: The Simulation Hydrograph Pre-Processor
Flow data are available through a number of resources. Of particular note is the USGS
flow file portion of the USEPA's STORET database. Daily flows for many streams are available
from this data source. The REDUCE pre-processing program reads a user-provided input file
(from STORET for example) that may contain years of daily flow observations for a stream and
produces a file that can be used as input for Data Group D of the IPX framework. The maximum
number of data points that can be used to.define time series inputs in IPX is controlled by the
value of the FORTRAN parameter MB in the file WASPJPX.CMN and is ultimately
constrained by available computer memory. It is often desirable, and sometimes necessary due to
the memory limitations, to reduce the number of data points that define the simulation
hydrograph. For example, using daily observations, the hydrograph for a 25 year simulation
would be defined by 9125 points. It is sometimes undesirable to develop model inputs based on
such a large number of data points. This pre-processing program reduces the number of points
used to define the simulation hydrograph by eliminating those values that do not differ
significantly from the previous value. This procedure permits the user to manage the size of the
simulation hydrograph as well as settling and resuspension velocities and numerical integration
time-steps.
142 - .
-------
When using REDUCE, the user is interactively prompted for information that controls
how the program develops the simulation hydrograph. The first prompt requests the name of the
data file that contains the flow data. Input files for the REDUCE pre-processor should typically
be of the form *.PRE_HYDRO. It is important that the user consistently organize the input file
such that all flow observations for each year are in a column; each column (representing one year
of flow data) must have 365 rows (one row for each day of the year). For example if Day 1 of the
simulation is January 1, the value representing January 1 of the second year of simulation will be
in the first row of the second column. The second prompt requests the number of years for which
flow data are input. The remaining prompts request information used to reduce the number of
flow values that define the hydrograph. Although REDUCE can be used independent of the units
(ft3/s, m3/s) of the flow values, units of m3/s are recommended due to the input requirements of
the SETTLE and RESUSPND pre-processing programs.
The next two prompts request that the user specify thresholds that divide the input flows
into three classifications: high, medium, and low flows. The remaining prompts request that the
user specify tolerance levels for each flow classification. The tolerance levels are fractions (0 <
tolerance < 1) that are used to assess whether a flow value is sufficiently similar to neighboring
flow values that it need not be included as a time-break in the simulation hydrograph. After
reading a flow value from the *.PRE_HYDRO file, the program determines whether to include
the value as a time-point in the simulation hydrograph by comparing it to the preceding flow
value selected as a time-break. If the current flow value is within the tolerance (previous
value*(l - tolerance) < current value < previous value*(l+tolerance)), then it is not included as a
time-break in the simulation hydrograph. Jf the flow value is outside this range, then it is
included in the simulation hydrograph as a time-break. This process in repeated for each input
flow value. However, the last value input is always included as a time-break to insure proper
operation of the IPX framework.
Output from the REDUCE pre-processor is contained in a file named *.POST_HYDRO.
The number of time-breaks, k, in the simulation hydrograph is given on line 1 of the output file
(i5 format), followed by k lines of paired values that are the simulation time (in days) and flow
defining the hydrograph (i5,al,ell.4 format). The simulation hydrograph defined in the
*.POST_HYDRO file can be used as input to estimate settling and resuspension velocities,
numerical integration time-steps, as well as input for IPX Data Group D after simple formatting.
143
-------
4.3 SETTLE: The Settling Velocity Pre-Processor
The SETTLE pre-processor is used to estimate TSS settling velocities as a function of the
simulation hydrograph when simulating all solids as a single state variable (TSS). SETTLE
estimates the distribution of particle into three categories: fine, medium, and coarse. Particles
within each category are assigned a settling velocity. An overall TSS settling for each time-break
in the hydrograph is computed as the sum of the products of the fraction of particles in each
category and settling velocity of the particle category.
When using SETTLE, the user is interactively prompted for the name of the data file that
contains the hydrograph (from REDUCE) and settling velocities of the three particle categories at
each simulation time-break. All other information necessary to control how settling velocities
are estimated is read directly from the input file. Input files for the SETTLE pre-processor
should typically be of the form *.PRE_VS. The first line of the input file must specify whether
flows in the river system are controlled or uncontrolled (0 represents a controlled system, 1
represents an uncontrolled system). The user must then provide information describing the
frequency distribution of river flows. For controlled rivers, the user must specify the river flow
at the SO*, 70*, 80*, and 98* percentile in m3/s, each on a separate line. For uncontrolled
systems, the user must specify the river flow at the 80*, 98* and 99.9* percentile in nP/s, each
on a separate line. The remainder of the file must be organized such that the simulation time,
flow, fine particle and settling velocities for the fine, medium, and coarse particle categories
appear on a single line for each time-break in the hydrograph. All flows must be in m3/s. All
settling velocities must be in m/s.. All times must be in days (recall that all time output from
REDUCE is in days). SETTLE then calculates the fraction of fine, medium and coarse particles
in the water column at that time based on the flow-rate. These fractions are then used to
calculate the average settling speed for all particles based on the three input settling speeds.
Output from the SETTLE pre-processor is contained in two files. The first file is named
*.POST_VS. This file has k lines, where k is the number of time-breaks in the simulation
hydrograph. Each line contains the fraction of particles that are fine, medium, and coarse (the
sum of these fractions will be 1), the computed overall TSS settling velocity, and the simulation
time. The second file is named *.VS. This file includes the number of entries, k, on the first
line, followed by k pairs of TSS settling velocity and simulation time values. After the first line,
each line of the *.VS file will have up to four pairs of entries per line and will have as many lines
as are needed for k pairs of settling velocities and simulation times (4(el0.2,f 10.0) format). This
file can be used as input for IPX Data Group D without additional formatting.
144
-------
4.4 RESUSPND: The Resuspension Velocity Pre-Processor
The RESUSPND pre-processor is used to estimate resuspension velocities as a function
of the simulation hydrograph, the time history of resuspension events, and sediment armoring.
RESUSPND estimates the shear stress exerted by water flowing over the sediment bed from
which the mass of sediment resuspended during an event is computed. A resuspension velocity
for each time-break in the hydrograph is then computed from the mass resuspended as a function
of the sediment bulk density.
When using RESUPEND, the user is interactively prompted for the name of the data file
that contains the hydrograph (from REDUCE) and settling velocities of the three particle
categories at each simulation time-break. All other information necessary to control how
resuspension velocities are estimated is read directly from the input file. Input files for the
RESUSPND pre-processor should typically be of the form *.PRE_VR. The first line of the input
file must specify the sediment properties that control resuspension (as described in Section
1.4.5.2.): taucrit, AO, m, Ax, density, te, and NCFLG. The first three constants (taucrit, AO, m)
are parameters in the epsilon (£) equation (Equation 1.25.). The next constant (Ax) is a
parameter in the shear stress equation (Equation 1.26.). The shear stress equation relates stream
flow to the shear exerted at the sediment-water interface. The next two constants (density and te)
are parameters used to express the mass of sediment resuspended as a resuspension velocity
(Equation 1.27.). The last constant (NCFLG) indicates whether the sediment bed experiences
non-cohesive resuspension (0 indicates that non-cohesive resuspension does not occur, 1
indicates that non-cohesive resuspension does occur). The epsilon equation defines the total
mass of cohesive sediment that can be resuspended for a given stress applied at the sediment-
water interface. In general, once the mass of sediment resuspended during an event equals
epsilon, resuspension decreases to background levels. However, non-cohesive sediments (such
as coarse sands) do not armor and will continue to resuspend as long as shear stress exerted at the
sediment water interface exceeds the critical shear of the non-cohesive sediments. If NCFLG =
1, the second line of the input file must specify taucrit_nc, the critical shear stress for
resuspending non-cohesive sediment, in dynes/cm2. The remainder of the input file must be
organized such that the simulation time, flow, the sediment age factor of the surficial age layer
(Z), and the background resuspension velocity appear on a single line for each time-break in the
hydrograph. All flows must be in m3/s. All background resuspension velocities must be in m/s.
All times must be in days (recall that all time output from REDUCE is in days). RESUSPND
then calculates the resuspension velocity at each time, based on the flow-rate and the computed
time history of sediment erosion potentials.
145
-------
Output from the RESUSPND pre-processor is contained in two files. The first file is
named *.POST_VR. This file has k lines, where k is the number of time-breaks in the simulation
hydrograph. Each line contains the value of epsilon for the given flow-rate, the depth of
sediments resuspended, the event resuspension velocity (does not include background
resuspension), the overall resuspension velocity (includes background resuspension), and the
simulation time. The second file is named *.VR. This file includes the number of entries, k, on
the first line of the file, followed by k pairs of resuspension velocity and simulation time values.
After the first line, each line of the *.VR file will have up to four pairs of entries per line and will
have as many lines as are needed for k pairs of resuspension velocities and simulation times
(4(el0.2,flO.O) format). This file can be used as input for IPX Data Group D without additional
formatting.
4.5 TIMESTEP: The Simulation Time-Step Pre-Processor
The TIMESTEP pre-processor is used to estimate the optimum model time-step as a
function of the simulation hydrograph and volume of the limiting (critical) segment in the model
network. TIMESTEP estimates the hydraulic residence time (HRT) from the critical volume and
flow. A suggested time-step history for the given hydrograph is then computed.
When using TIMESTEP, the user is first interactively prompted for the volume of the
critical segment in the model network. Next, the user is prompted for a hydraulic residence time
divisor factor (HRTDF). Finally, the user is prompted for the name of the data file that contains
the hydrograph (from REDUCE). For advectively dominated tributary systems, the critical
model segment will typically be the smallest water column segment in the network. The HRTDF
is used as a "safety" factor, to ensure that the computed model time-steps do not exceed the
maximum stability criterion for the numerical integration approach IPX uses to solve the mass
balance equation. HRTDFs will typically have a value of 3.0. However, users should be aware
that IPX treats model time-steps as a step, rather than a linear, function so it may be necessary to
use a larger factor during rapidly changing portions of the hydrograph. Input files for the
TIMESTEP pre-processor should typically be of the form *.PRE_DT. The input file must be
organized such that the simulation time and flow appear on a single line for each time-break in
the hydrograph. TIMESTEP will read all lines from the input file until the end of the file is
reached.
146
-------
Output from the TDVLESTEP pre-processor is contained in two files. The first file is
named *.POST_DT. This file has k lines, where k is the number of time-breaks in the simulation
hydrograph. Each line contains the minimum hydraulic residence time, the maximum model
time-step, the "suggested", model time-step, and the simulation time. The second file is named
*.DT. This file includes the number of entries, k, on the first line of the file, followed by k pairs
of time-step and simulation time values. After the first line, each line of the *.DT file will have
up to four pairs of entries per line and will have as many lines as are needed for k pairs of time-
steps and simulation times (4(el0.2,flO.O) format). This file can be used as input for IPX Data
Group A without additional formatting.
147
-------
CHAPTERS
IPX PROGRAMMER'S GUIDE
5.1 IPX PROGRAMMING OVERVIEW
This chapter provides a brief overview of the important programming aspects of the IPX
framework. In addition to its inherent flexibility, one demonstrable strength of IPX and its
WASP4 predecessor is that, with a basic degree of familiarity, users can customize most
functions of the framework to create highly specialized applications to address a broad range of
water quality simulation issues. However, a word of caution is in order: many users assume that
there are no errors in "official" software distributed or supported by USEPA. IPX, like WASP4
and its antecedents, is a complex program. No guarantees are made regarding either the
suitability of IPX for constructing a water quality model or the computational performance of the
code. Although substantial efforts have been taken to examine all aspects of the code, the user is
strongly advised to carefully examine both the IPX source code and all model outputs to ensure
proper operation.
The authors would appreciate receiving notification of any problems encountered with the
IPX program. This may be done by mail, phone, or email to:
Mark Velleux, Wisconsin Department of Natural Resources WT/2, P.O. Box 7921,
Madison, Wisconsin 53707-7921; Phone (608) 267-5262; email: vellem@dnr.state.wi.us.
AVAILABILITY OF IPX PROGRAMS AND PORTING TO OTHER
PLATFORMS
5.2
The IPX programs, and this User's Guide, will be provided upon request to the authors.
IPX is not officially supported by either the USEPA or WDNR. This means that users should not
expect support beyond receiving the programs and User's Guide, unless special arrangements are
made with USEPA or WDNR.
148
-------
IPX and its supporting suite of pre- and post-programs were initially developed in
FORTRAN?? in a VAX/VMS environment. When porting the IPX source code to other
environments, users should be aware that many VAX FORTRAN programming conventions
(such as the treatment of uninitialized or undeclared variables, etc.) were used. Many compilers
have options to accommodate VAX conventions and VAX FORTRAN extensions. These
options should be used to help ensure that IPX operates consistently on all platforms to which it
is ported. IPX has been successfully ported to Sun (Solaris) and Compaq/DEC Alpha (Digital
UNIX, OSF/1) UNIX workstations. IPX has also been successfully ported to Cray YMP and
C90 supercomputers as well as LINUX-based and Windows NT-based personal computers.
5.3 UNIX ENVIRONMENT DEVELOPMENT TOOLS FOR IPX
Following its first release, the IPX program source code was ported to a UNIX computing
environment to facilitate development efforts. Three basic modifications were made at the time
the source code was ported. The first modification was splitting the program source code into a
separate file for each subprogram unit (subroutine). The second modification was
implementation of the UNIX C preprocessor utility cpp to allow flexible inclusion of the
platform-dependent codes used by different computer operating systems (i.e. VMS, OSF/1,
Solaris, etc.). The third modification was implementation of a UNIX Makefile to provide swift
and organized compilation and loading of the IPX source code.
The IPX source code is now contained in over seventy FORTRAN source files containing
the IPX subprograms and include files. Files that contain source code for the IPX subprograms
are identified by the file extension ".F'. Files that contain included common blocks are
identified by the file extension ".inc". Other files include the message file messfile_rl.dat
(needed during' program execution), the platform.h file (used to define platform specific
conditions for source code compilation), and the UNIX Makefile (used to manage source code
compilation).
UNIX cpp is a program run prior to compilation of program source code that reads the
source, interprets specially prepared directives imbedded in it, and outputs another version of the
source whose modifications are determined by the cpp directives and predefined variables named
cpp macros that are used both as flags to steer execution of the preprocessor and as symbols
representing some syntax of the program. These symbols can be defined either within the source
code file of the program itself, or in an independent file made known to cpp at the time that the
preprocessing occurs. (Such symbols can also be defined at compile time, but instances of this
149
-------
practice are more the exception than the rule and are not currently being used as a part of IPX
development.) The current implementation of IPX uses a file named platform* to define cap
macro flags. A copy of this file is provided in the Appendix, together with a description of the
individual flags currently being used. A simple example demonstrating the use of cpp is shown
below:
#define FLAG 1
PROGRAM MESSAGE
*
#if(FLAG)
PRINT *, DThis message is only printed if FLAG is setD
#endif
*
END
IPX source code development is assisted by the use of the make utility which is
responsible for organizing compilation and loading of applications in the UNIX environment
Use of the make command is very straightforward. Simply go into the UNIX directory
containing both the Makefile and the UNIX source and type make at the UNIX prompt. This will
force compilation of any source files that have not already been compiled since they were last
updated, and link the resulting object files into the application executable program; i.e., ipx-
«•
The instructions for compiling and loading or linking the application are contained in a
file named Makefile and are unfortunately complex. Most alterations to this file should be made
by the developers, however, some simple modifications to the make procedure of IPX are
possible Changing the compilation flags is accomplished by editing Makefile, locating the string
FFLAGS, and changing its arguments as desired. Similarly, new subprograms can be added to
the application by placing their source and object file names in the respective SOURCE and
OBJECTS macro lists of Makefiles.
5.4 IPX PROGRAM DESCRIPTION
Figure 5.1 is a flowchart of IPX and illustrates the functional relationships among the
subroutines. IPX is a modular program. Its many subroutines can be grouped into the functional
categories of input, process, output, and utility. The main program calls input subroutines calls
subroutine EULER, and closes files at the end of the simulation. The input subroutines are called
150
-------
IPX
WASP1
| WASP2 |
WASPS |
WASP4
j WASPS |
WASP6
READ.KINETIC
| WASP9 | | EU1
.ER
WAS6A
TINIT
WAS13 [.["pirny"] ISEDIMENT.AGJI |BURIAL| |PUSH[ |POP| 1WAS14 |
[ WAS8B | [aCALC [ pWASPB | | WAS12
WASPS
All
Kinetic Subroutines
QSURF1/QSURF2
QPORE
QSED
QEVAP
WA12A
BEDSED
TOPT
Figure 5.1. Hierarchy of major subprogram units in IPX.
sequentially, as shown. Subroutine EULER controls the actual simulation and calls DERIV each
time-step to recalculate mass derivatives. The output subroutines are called sequentially as
shown after the simulation is completed. Utility subroutines are called by the other subroutines
as needed. Data are shared among the subroutines primarily through common blocks.
5.4.1 Main Program
IPX_R1:
The IPX_R1 main program is the control module. It operates the calling sequence for the
input, calls subroutine EULER, and closes files at the end of the simulation.
5.4.2 Input Subroutines
WASP1
WASP1 opens the input and output units, then reads Data Group A for model
identification and system bypass options. Information is printed and values and arrays are
initialized. The subroutine then returns control to the main program.
151
-------
WASP2 ,
WASP2 reads Data Group B for sets of dispersion coefficients, cross-sectional areas, and
characteristic lengths. These are converted to bulk exchanges, and information is stored in
memory and printed. The subroutine then returns control to the main program.
WASPS
WASPS reads Data Group C for volumes. If indicated, volumes are read from the *.RST
restart file. Information is stored in memory and printed. The subroutine then returns control to
the main program.
WASP4
WASP4 reads Data Group D for advective flows, which are converted to internal units of
cubic meters per day. Information is stored in memory and printed. The subroutine then returns
control to the main program.
WASPS
WASPS reads Data Group E for boundary concentrations for each model system.
Information is stored in memory and printed. The subroutine then returns control to the main
program.
WASP6
WASP6 reads Data Group F for waste loads for each model system. Information is stored
in memory and printed. If indicated, WASP6 calls subroutine WAS6A to read non-point source
loads from an external file. The subroutine then returns control to the main program.
WAS6A
WAS6A opens an unformatted non-point source loading file NPS.DAT. The structure of
the NPS.DAT file is presented in Table 5.1. The runoff day corresponding with the initial IPX
simulation day is read. Input segment numbers corresponding to each runoff load are read.
Actual runoff loads from the file are printed as specified. Finally, the file pointer is positioned
properly to begin the simulation. The subroutine then returns control to subroutine WASP6.
e
READJKINETIC
READ_KINETIC reads Data Group G for parameters for each segment. It then reads
Data Group H for chemical constants. Finally, it reads Data Group I for kinetic time functions.
152
-------
Table 5.1. Structure of the non-point source load file NPS.DAT.
Record Number
Contents of Record
1
NWKS
((NPSWK(U), I=1,NOSYS), J=1,NWKS)
((NPSWKOJ), I=1,NOSYS), J=1,NWKS)
N+l
((NPSWK(I,J), I=1,NOSYS), J=1,NWKS)
Variable
Type
Definition
NWKS
1*4
The number of runoff loads
NPSWK
R*4
Runoff loads averaged over day, kg/day
NOSYS
1*4
Number of water quality variables (or systems)
1*4
Water quality variable counter
1*4
Runoff load counter
N
1*4
Number of days for which loads are available
Information is stored in memory and printed. The subroutine then returns control to the main
program.
WASP9
WASP9 reads Data Group J for bulk densities, maximum concentration, initial
concentrations, and dissolved fractions in all segments for each model system. If indicated,
initial concentrations are read from the *.RST restart file. Information is stored in memory and
printed. Then WASP9 reads Data Group K for surficial sediment age layer characteristics. This
information is also stored in memory and printed. The subroutine then returns control to the
main program.
5.4.3 Output Subroutines
AVERAGING
AVERAGING is called by subroutine TOXIDU to perform time averaging of simulation
results for output to the *.DMA file.
153
-------
TOXIDU
TOXDDU is called by subroutine WASPB at specified intervals to write simulation output
to the *.DMP and *.DMA files.
WRTTEEXPDATA
WRTTEEXPDATA is called in WAS 13 at specified intervals to write export results to the
*.EXPfile.
WRITEMSBDATA
WRTTEMSBDATA is called in EULER at the end of a simulation to write mass balance
results to the *.MSB file.
WRITERSTDATA
WRITERSTDATA is called in EULER at specified intervals to write simulation restart
information to the *.RST file.
^A Simulation Control and Process Subroutines to Compute
EULER
EULER is the computational engine of IPX. EULER performs the first-order numerical
integration of the mass balance equations. First, counters and time functions are initialized to
TZERO as needed by subroutine TINIT. Subroutine WAS13 is called to set up initial printouts.
Subroutine DERIV is called to compute initial mass derivatives. Input data are initially screened
for fatal error conditions and then the integration (simulation) proceeds, time-step by time-step.
For each time-step, EULER loops through each system and segment, calls subroutines
DERIV, SEDIMENT_AGE, and BURIAL, and computes the new mass as follows:
new mass = old mass + mass derivative x time-step
Each new concentration is set to the new mass divided by the new volume, and the mass
derivative is reset to zero. Unless the negative solution option is active, any negative
concentrations are replaced by one-half of the old mass divided by the new volume. Next,
EULER increments the time and adjusts the new day counter if necessary. When intermediate
print times are reached, EULER calls subroutine WAS 13 to produce intermediate printouts and
triggers storage of all display variables (by returning IDISK = 1). Volumes are stored if IDISK =
1. The final task for each time-step is to check for a new time-step and for the end of the
simulation. Subroutine WAS 14 is called if it is time to use a new time-step.
154
-------
When the final time for the simulation is reached, EULER triggers a final storage of
display variables, then stores final volumes and concentrations in the *.RST file. At the end of
the simulation, the subroutine returns control to the main program.
DERIV .
DERTV is called by EULER to calculate mass derivatives. Subroutines WAS8B and
QCALC is called to update flow and exchange time functions; Subroutine WASPB is called to
obtain kinetic derivatives. Subroutine WAS 12 is called to obtain transport and loading
derivatives.
WASPB
WASPB is the subroutine mat controls the computation of kinetic mass derivatives.
WASPB calls subroutine WASPS and all other subroutines necessary to compute the kinetic
portion of the mass derivative. An overview subroutine WASPB and all kinetic subroutines is
presented in Section 5.3.5.
WAS 12
WAS 12 is called by DERIV to obtain the transport and loading derivatives. Upon entry
to WAS 12, only the kinetic portion of the mass balance derivative has been evaluated by
WASPB. WAS12 calls subroutines WA12A, WAS8A, and BEDSED to determine the proper
upstream and downstream concentrations for advective flow and dispersive exchange and to
update time functions and sediment conditions. Then WAS 12 calls subroutines QSURF1 or
QSURF2, QPORE, QSED, and QEVAP and calculates the mass derivatives due to advective
flow, dispersive exchange, point source waste loading, and runoff loading, and adds them to the
kinetic derivative. WAS 12 goes through the following steps:
a. Using the IQ and JQ vectors as drivers, WAS 12 computes advective transport.
For each system, variable boundary concentrations are updated by calling WAS8A
if necessary. For each flow, Q, proper upstream and downstream concentrations
are assigned by calling WA12A. The advected concentration CSTAR is
determined, and mass derivatives for the downstream and upstream segments are
adjusted by ±Q CSTAR.
b. Using the IR and JR vectors as drivers, WAS 12 computes dispersive transport.
For each system and each exchange flow, R, proper upstream and downstream
155
-------
c.
d.
concentrations C2 and d are assigned by calling WA12A. Mass derivatives for
the downstream and upstream segments are adjusted by ± R (C2-Ci).
Using the IWK vector as a driver, WAS 12 computes point source loading. For
each system, variable loadings are updated by calling WAS 8A if necessary. For
each load L (in kg/day), the mass derivative for the affected segment is adjusted
by + L.
Using the INPS vector as a driver, WAS 12 computes diffuse source loading if
appropriate. New loads are read from file NPS.DAT at the beginning of each new
day. For each load L' (in kg/day), the mass derivative for the affected segment is
adjusted by + L'.
WA12A
WA12A is called by. WAS 12 to determine the proper upstream and downstream
concentrations C2 and d for advective flow from segment JQ to segment IQ or dispersive
exchange between segments JR and IR. For flows or exchanges with a downstream boundary,
the proper boundary concentration is located for d.' For flows or exchanges with an upstream
boundary, the proper boundary concentration is located for €2-
WASPS
WASPS is called by subroutine WASPB and updates piecewise linear time functions, if
any, for kinetic time functions. This means computing new slopes and intercepts, and setting a
variable to indicate the next simulation time that the functions are to be updated. The following
convention is used for the ife update:
slope
- /(4-.-/C4
intercept = f(t)M
next update time = tM
WAS8A
WAS8A updates piecewise linear time functions, if any, for boundary conditions and
forcing functions. This means computing new slopes and intercepts for any system or state
variable that requires an update, and setting a variable to indicate the next simulation time that
156
-------
the piecewise linear functions are to be updated. The same conventions used in WASPS are used
in WAS8A for computing slopes and intercepts.
WAS8B
WAS8B updates piecewise linear time functions for dispersion coefficients and flows.
Updated dispersion coefficients for each exchange field and time function are stored in the array
BRINT(NF,NT). Updated flows for each field and time function are stored in the array
QINT(NF,NT).
QCALC
QCALC, called in DERTV, sums the flows from all time functions for each segment.
QSURF1/QSURF2
QSURF1 or QSURF2, called in WAS 12, compute the mass derivative due to surface
water (advective) flow. QSURF1 is called only if IQOPT =1. QSURF2 is called only if IQOPT
= 2. Surface water flow options (IQOPT) are described in Section 1.4.1.
QPORE
QPORE, called in WAS 12, computes the mass derivative due to pore water flow.
QSED
QSED, called in WAS 12, computes the mass derivative due to settling and resuspension.
QEVAP
QEVAP, called in WAS 12, computes the mass derivative due to precipitation and
evaporation.
BEDSED
BEDSED, called in WAS 12, computes changes in volumes and porosities for sediment
bed segments, depending upon the sediment bed option used. For surficial sediment bed
segments, porosity is constant but volumes change in response to net difference between the
settling and resuspension fluxes. When the sedimentation (compaction) time interval is reached,
mass is transported between surficial and subsurface sediment segments. The porosity and
volume of all subsurface sediment segments are constant.
157
-------
SEDIMENT.AGE
SEDIMENT_AGE, called in EULER, manages the aging process for solids and
associated chemicals in each age layer of the surficial sediments.
BURIAL
BURIAL, called in EULER, computes the mass derivative due to upward (scour) and
downward (burial) movement of solids and associated chemicals through the surficial and
subsurface sediments for the Eulerian sediment bed option.
PUSH
PUSH, called in EULER, performs ghost collapse as specified by the ghost collapse
option and computes the downward re-indexing of sediment segments and ghost stack elements
for the semi-Lagrangian sediment bed option.
POP
POP, called in EULER, computes the upward re-indexing of sediment segments and
ghost stack elements for the semi-Lagrangian sediment bed option.
WAS13
WAS 13 is called every print interval by EULER to print intermediate concentrations or
mass checks on a designated constituent. At this time, the solution stability is checked by
comparing the maximum concentrations specified by the user with calculated concentrations. If
any concentrations exceed the maximum, the simulation is stopped (aborted).
WAS14
WAS 14 is called by EULER to adjust the integration step size (time-step) as specified by
the user in Data Group A.
TINT!
TINT! is called by EULER at the beginning of the simulation to adjust time functions to
the initial time TZERO. TINTT checks and adjusts time functions for exchanges, flows, kinetic
time functions, boundary concentrations, and loads.
AVEnSfTT
AVEINIT is called by TOXIDU to initialize program variables used to perform time
averaging of simulation results.
158
-------
TOPT
TOPT can be called by the user WASPB subroutine to maximize the time-step subject to
the flow and dispersion stability constraints. This should reduce numerical dispersion, but is not
unconditionally stable. The time-step is calculated by TOPT every 0.5 days will fall between
0 01 and 0.5 days. Note: the operation of subroutine TOPT has not been verified. The
authors recommend that users input the time-step in Data Group A.
5.4.5 Utility Subroutines
BRKERR
BRKERR prints an error message to output file and screen concerning the number of data
points in a time function; the simulation is aborted.
FMTER
FMTER prints an error message to output file and screen concerning input data formats;
the simulation is aborted.
GETMES
GETMES retrieves text messages from the messfile_rl,dat file for subsequent output to
screen.
PRTSCR
PRTSCR prints text messages to screen.
SETCA
SETCA sets a character array to a specified character value.
SETIA
SETIA sets an integer array to a specified integer value.
SETRA
SETRA sets a real array to a specified real value.
SETRB
SETRB sets a real three-dimensional array to a specified real value.
159
-------
SETXA
SETXA sets a double precision array to a specified double precision value.
TSTAMP
TSTAMP write a time and date at the state of a simulation to the *.OUT, *.DMP, *.EXP,
and *.DMA files.
WERR
WERR writes error messages for improper segment designations and missing boundary
conditions; the simulation is aborted.
WEXTT
WEXTT aborts the simulation when error conditions are encountered.
WMESS
WMESS prints a message when stability criteria are violated; the simulation is aborted.
5.4.6 Toxic Chemical Kinetic Subroutines
Source code for IPX toxic chemical kinetics are written in a modular structure that
includes numerous subroutines as presented in Figure 5.2. Each transformation and transfer
process is separated into one or several subroutines. This structure allows for convenient
addition or modification of the kinetic descriptions.
WASPB
WASPB is the central water quality subroutine that calculates the kinetic mass derivative
and stores the proper display variables for later printout. WASPB calls other subroutines as
presented in Figure 5.2. WASPB is called every time-step by subroutine EULER and perfc
the following tasks:
rorms
At the start of a simulation, subroutine TOXINT is called to initialize parameters
for kinetic processes. Chemical diffusivities at 25°C, basic partitioning
parameters, and transformation rates and are initialized.
The current values for the piecewise-linear time functions of time are calculated.
WASPS is called if a time-break has been reached for any of the functions.
160
-------
WASPB
ENV SOLID | CHEM | |RXPROD| | TOXIDU
PRTITON] |FRCTION| | BIODEG | HYDROL
VOLAT
OXID
PHOTO
EREACT
Figure 5.2. Hierarchy of major kinetic process subroutines in IPX.
3. For each segment of the model the following tasks are performed:
a. Subroutine ENV is called to evaluate current environmental characteristics.
b. Subroutine CHEM is called to evaluate the kinetic portion of the derivative
describing a chemical. Effective rate constants for kinetic processes are
stored in memory.
c. Subroutine SOLE) is called to evaluate the kinetic derivative for each solid
type. (There are no kinetic for solids in the present code.) Locations in the F
array that the transport of solids in the transport fields are set.
d. If a print interval has been reached, subroutine TOXIDU is called to store
simulation results in the *.DMP and *.DMA files.
CHEM(I)
CHEM determines the kinetic portion of the derivative describing a chemical in segment
I. Tasks executed are:
1. Subroutine PRTION is called to determine the partition coefficient of the
chemical.
161
-------
2.
3.
4.
5.
6.
The Henry's Law constant and air-water partition coefficient for the chemical are
computed.
Subroutine FRCTION is called to determine the fraction of total chemical in each
of the 3 phases (dissolved, DOC-bound, and sorbed to solids).
Subroutines BIODEG, HYDROL, PHOTO, VOLAT, OXID and EREACT are
called to determine the transfer and transformation rates that make up the kinetic
portion of the derivative.
The individual rates returned by the subroutines called are summed to yield the
total kinetic portion of the derivative.
The derivative is multiplied by the segment volume to be consistent with the
transport derivative calculation.
ENV
ENV computes segment-specific physical and chemical characteristics of each segment.
PRTTTONCICHEM)
PRTITON computes the equilibrium partition coefficients defining sorption of chemical
ICHEM to each of the three solids types in segment I. For the neutral chemical, the classical
partition coefficient (PKC) is computed from the organic carbon partition coefficient and the
fraction organic carbon of the solids. The partition coefficients (PART) are calculated from
PIXC using the particle interaction model for sorption (see Section 1.5).
FRCTION(ICHEM, DISTRIB, PIDOC)
FRCTION computes the fraction of chemical ICHEM in segment I present in the
dissolved, DOC-bound and particulate (solids-sorbed) phases. This computation consists of 3
parts:
1.
Distribution factors are computed for each phase. A distribution factor is the ratio
of component concentration to dissolved neutral chemical concentration. It is
calculated for sorbed chemical. The distribution factors for dissolved ionic
chemical are passed to FRCTION as subroutine argument DISTRIB.
162
-------
2. The distribution factors are summed. As part of the summation the factor for
dissolved components are multiplied by the porosity of the segment (PORE).
3. The fractions are computed as the distribution factor divided by the sum of
distribution factors.
4. Chemical concentrations in each phase are computed from the total concentration
and the distribution coefficient.
5. Locations in the F array that define the transport of the chemical in the transport
fields are set.
BIODEG(ICHEM)
BIODEG computes the loss rate of chemical ICHM in segment I due to biodegradation.
Separate loss rates are computed for solids sorbed chemical (BIOS) and dissolved and DOC-
bound chemical (BIOW) using 20°C rates passed to the routine in array KBIO20 and temperature
correction factors passed in arrays Q10DIS (dissolved chemical), Q10DOC (DOC-bound
chemical) and Q10PAR (solids sorbed chemical).
HYDROL(ICHEM)
HYDROL computes the loss rate of chemical ICHM in segment I due to chemical
hydrolysis. Separate loss rates (ug/L-d) are computed for alkaline hydrolysis (ALKH), neutral
hydrolysis (NEUTH) and acid hydrolysis (ACIDH) using 20°C rates passed to the routine in
array KH20 and Arrhenius temperature correction factors passed in arrays EHOH (alkaline
hydrolysis), EHN (neutral hydrolysis) and EHH (acid hydrolysis).
PHOTO(ICHEM)
PHOTO computes the loss rate of chemical ICHEM in segment I due to photolysis. If the
user has selected photolysis option 1 (IPHOTO = 1), the specified depth-integrated photolysis
rate constant (KDPG) and quantum yield (QUANTG) for each photolyzing phase are used to
compute as the overall photolysis loss as the product of the rate constant, the quantum yield, and
chemical concentration in each phase.
163
-------
VOLAT(ICHEM, VLT)
VOLAT computes the rate of loss by volatilization of a chemical ICHEM in segment I.
The calculation approach used depends on the volatilization options (XVL and XVG) selected by
the user.
OXID(ICHEM)
OXID computes the loss rate of chemical ICHM in segment I due to chemical oxidation.
Separate oxidation rates are computed for dissolved, DOC sorbed and solids sorbed chemical
fraction using 20°C rates passed to the routine in array KOX20 and Arrhenius temperature
correction factors passed in array BOX. These rates are summed and an overall oxidation rate is
computed.
EREACT(ICHEM)
EREACT computes the loss rate of chemical ICHEM in segment I due to the user-defined
extra reaction.
SOLID
SOLID computes transport fractions and kinetic derivatives for solids state variables.
TOXINT
TOX3NT is called by WASPB at the start of a simulation to initialize parameters for
kinetic processes. Chemical diffusivities in water and air at 25°C, basic partitioning parameters,
and transformation rates and are initialized. This subroutine is only called one time, at the start
of a simulation.
5.5 INPUT AND OUTPUT UNITS
All the input/output (I/O) units can be reassigned an integer value in the WASP MAIN
subroutine. However, new users should first develop a firm understanding of the structure and
function of the program before reassigning I/O units. A brief description of each I/O-related
integer and its default value follows:
AUX: Default value is 4. AUX refers to the use of an auxiliary flow file. This file has
been created outside the WASP programs and is used to input flows and volumes Example-
READ(AUX).
164
-------
IN: Default value is 2. The value 2 refers to the input data set. Input data set is a
sequential formatted file. May also be used to represent the integer 2. Example: READ(IN).
OUT: Default value is 5. OUT refers to the output file *.OUT, a sequential formatted
file. Example: WRTTE(OUT). .
RESRT: Default value is 9. RESTRT refers to the *.RST file that contains-a snapshot of
model conditions which can be used as initial conditions for a sequence of simulations or
restoring a simulation that did not complete. Example: WRTTE(RESTRT).
MESS: Default value is 6. MESS is used to write inquiry messages to the screen and
display run time status messages. For more information see MFLAG variable in Data Group A.
IEXP: Default value is 16. IEXP is used to write the mass exported across a model
interface during the'simulation to the *.EXP file for any one state variable and total solids as
specified in the input data file.
MASS: Default value is 20. MASS is used to write the mass balance file for all
systems.
IDMP: Default value is 15. IDMP is used to write all instantaneous (point-in-time)
simulation results that can be latter recalled and printed out using the W4DIS274 post-processing
program.
IDMA: Default value is 25. IDMA is used to write all time averaged simulation results
that can be latter recalled and printed out using the W4DIS274 post-processing program.
5.6 COMMON BLOCKS
Common blocks permit program subroutines to share data. The common blocks also
contain the parameter statement assignments used to control the dimension of variables that
constrain the maximum possible size of a simulation. The user can change the maximum
dimension of a simulation by changing the values assigned to each parameter, compiling, and
linking the source code. A brief description of each common block follows:
165
-------
AVERAGE.CMN: Common block used by output subroutines to compute time averaged results
from simulation calculations.
CONC.CMN: Common block used by kinetic process subroutines to store chemical
concentrations from simulation calculations.
ENVIRON.CMN: Common block used by kinetic process subroutines to store environmental
information.
IPX.CMN: Common block used to declare general variables used throughout the code. This
common block also contains the PARAMETER statements that control the dimension of
simulation constraint parameters, such as number of segments and number of systems.
KNETTC.CMN: Common block used by kinetic process subroutines to store degradation rates
and constants.
PARAM.CMN: Common block used by kinetic process subroutines to store information related
to the environmental parameters and time functions
PHYSCHM.CMN: Common block used by the kinetic process subroutines to store information
related to the physical-chemical calculations.
166
-------
REFERENCES
Abramowitz, M. and I. A. Stegun. 1972. Handbook of Mathematical Functions. National
Bureau of Standards, Applied Mathematics Series 55, Washington, D.C.
Alexander, M. 1980. Biodegradation of toxic chemicals in water and soil. In: Dynamics,
Exposure,'and Hazard Assessment of Toxic Chemicals; R. Haque, (ed.), Ann Arbor Science
Publishers, Ann Arbor, Michigan.
Ambrose, R.B. et al. 1983. User's manual for the chemical transport and fate model
(TOXIWASP), Version 1. EPA-600/3-83-005. U.S. Environmental Protection Agency, Office
of Research and Development, Environmental Research Laboratory, Athens, Georgia.
Ambrose, R.B., T.A. Wool, J.P. Connolly, and R.W. Shanz. 1988. WASP4, A hydrodynamic
and water quality model - Model theory, user's manual, and programmer's guide. EPA-600/3-
87-039. U.S. Environmental Protection Agency, Office of Research and Development,
Environmental Research Laboratory, Athens, Georgia.
Ambrose, R.B., J.L. Martin, and T.A. Wool. 1993. WASPS, A hydrodynamic and water quality
model - Model theory, user's manual, and programmer's guide. U.S. Environmental Protection
Agency, Office of Research and Development, Environmental Research Laboratory, Athens,
Georgia.
Bella, D.A. and WJ. Grenney. 1970. Finite-difference convection errors. J. Sanitary Engin.
Div., ASCE, 96(SA6): 1361-1375.
Bowie, G.L., W.B. Mills, D.B. Porcella, C.L. Campbell, J.R. Pagenkopf, G.L. Rupp, K.M.
Johnson, P.W.H. Chan, S.A. Gherini and C.E. Chamberlin. 1985. Rates, constants, and kinetics
formulations in surface water quality modeling. Second Edition. EPA-600/3-85-040. U.S.
Environmental Protection Agency, Office of Research and Development, Environmental
Research Laboratory, Athens, Georgia.
167
-------
Brown, L.C. and T.O. Barnwell. 1987. The enhanced stream water quality models QUAL2E
and QUAL2E-UNCAS: Documentation and user manual. EPA-600/3-87-007 US
Environmental Protection Agency, Office of Research and Development, Environmental
Research Laboratory, Athens, Georgia.
Brown, D.S. and J.D. Allison. 1987. MINTEQA1, An equilibrium metal speciation model-
Users manual. EPA-600/3-87-012. U.S. Environmental Protection Agency, Office of Research
and Development, Environmental Research Laboratory, Athens, Georgia.
Burban, P.Y., Y. Xu, and W. Lick. 1990. Settling speeds of floes in fresh and sea waters J
Geophysical Res., 95(C10): 18213-18220.
Bums, L.A., D.M. Cline, and R.R. Lassiter. 1982. Exposure analysis modeling system
(EXAMS): User manual and system documentation, EPA-600/3-82-023. U.S. Environmental
Protection Agency, Office of Research and Development, Environmental Research Laboratory
Athens, Georgia.
Burns, L.A. and DM. Cline. 1985. Exposure analysis modeling system, reference manual for
EXAMS H. EPA-600/3-85-038. U.S. Environmental Protection Agency, Office of Research and
Development, Environmental Research Laboratory, Athens, Georgia.
Capel, P.D. and SJ. Eisenreich. 1990. Relationships between chlorinated hydrocarbons and
organic carbon in sediment and pore water. J. Great Lakes Res., 16(2):245-257.
Chapra, S.C. 1997. Surface Water-Quality Modeling. McGraw-Hill Companies Inc New
York, New York. '
Chapra, S.C. and K.H. Reckhow. 1983. Engineering Approaches For Lake Management
Volume 2: Mechanistic Modeling. Butterworth Publishers. Woburn, Masschusetts.
Connolly, J.P. and R. Winfield. 1984. A user's guide for WASTOX, A framework for modeling
the fate of toxic chemicals in aquatic environments. Part 1: Exposure concentration. EPA-
600/3-84-077. U.S. Environmental Protection Agency, Office of Research and Development,
Environmental Research Laboratory, Gulf Breeze, Florida.
168
-------
Delos, C.G., W.L. Richardson, J.V. DePinto, R.B.. Ambrose, P.W. Rodgers, K. Rygwelski, J.P.
St. John, W.L. Shaughnessy, T.A. Faha, and W.N. Christie. 1984. Technical guidance manual
for performing waste load allocations. Book H. Streams and rivers, Chapter 3, Toxic substances.
EPA-440/4-84-022. U.S. Environmental Protection Agency, Washington, D.C.
Di Toro, D.M., J.J. Fitzpatrick, and R.V. Thomann. 1981, rev. 1983. Water quality analysis
simulation program (WASP) and model verification program (MVP) - Documentation.
Hydroscience, Inc., Westwood, New York, for the U.S. Environmental Protection Agency, Office
of Research and Development, Environmental Research Laboratory, Duluth, Minnesota, Contract
No. 68-01-3872.
Di Toro, D.M. 1985. A particle interaction model of reversible organic chemical sorption.
Chemosphere, 14(10): 1503-1538.
Eadie, B.J., N.R. Moorehead, and P.F. Landrum. 1990. Three-phase partitioning of hydrophobic
organic compounds in Great Lakes waters. Chemosphere, 20(1-2): 161-178.
Filkins, J.C., Smith, V.E., Rathbun, J.E., and Rood, S.J. 1993. ARCS toxicity/chemistry
workgroup sediment assessment guidance document: Chapters 3-5. U.S. Environmental
Protection Agency, Office of Research and Development, ERL-Duluth, Large Lakes Research
Station, Grosse He, Michigan.
Gailani, J., C.K. Ziegler, and W. Lick. 1991. The transport of suspended solids in the lower Fox
River. J. Great Lakes Res., 17(4):479-494.
Green, A.E.S., K.R. Cross, and L.A. Smith. 1980. Improved analytical characterization of
ultraviolet skylight. Photochem. Photobio., 31:59-65.
Herbes, S.E. and L.R. Schwall. 1978. Microbial transformation of polycyclic aromatic
hydrocarbons in pristine and petroleum-contaminated sediments. Appl. Environ. Microbiol.,
35(2):306-316.
Karickhoff, S.W., D.S. Brown, and T.A. Scott. 1979. Sorption of hydrophobic pollutants on
natural sediments. Water Res., 13:241-248. .
169
-------
Karickhoff, S.W. and K.R. Morris. 1985. Sorption dynamics of hydrophobic pollutants in
sediment suspensions. Environ. Toxicol. Chem., 4:469-479.
Landrum, P.F., M.D. Reinhold, S.R. Nihart, and B.J. Eadie. 1985. Predicting the bioavailability
of organic xenobiotics to Pontoporeia hoyi in the presence of humic and fulvic acids and natural
dissolved organic matter. Environ. Toxicol. Chem., 4:459-467.
Larson, R.J., G.G. Clinckemaillie, and L. VanBelle. 1981. Effect of temperature and dissolved
oxygen on biodegradation of nitrilotriacetate. Water Res., 15:615-620.
Lee, D. Y., W. Lick, and S.W. Kang. 1981. The entrainment and deposition of fine-grained
sediments in Lake Erie. J. Great Lakes Res., 7(3):264-275.
Leopold, L.B. and T. Maddox. 1953. The hydraulic geometry of stream channels and some
physiographic implications. Professional Paper 252, U.S. Geological Survey, Washington, D.C.
Leopold, L.B., M.B. Wolman and J.P. Miller. 1964. Fluvial Processes in Geomorphology. W.
H. Freeman and Co., San Francisco, California.
Lick, W., McNeil, J., Xu, Y., and Taylor, C. 1995. Measurement of the Resuspension and.
Erosion of Sediments in Rivers. Department of Mechanical and Environmental Engineering,
University of California-Santa Barbara, Santa Barbara, California, Dated June 2, 1995.
Liss, P.S. 1973. Deep-sea research. Volume 20. pp. 221-239.
Lyman, W.J., W.F. Reehl and D.H. Rosenblatt. 1982. Handbook of Chemical Property
Estimation Methods. McGraw-Hill Book Co., New York, New York.
Mabey, W.R., J.H. Smith, R.T. Podell, H.L. Johnson, T. Mill, T.W. Chou, J. Gates, I.W.
Partridge, H. Jaber, and D. Vandenberg. 1982. Aquatic fate process data for organic priority
pollutants. EPA-440/4-81-014. U.S. Environmental Protection Agency, Washington, D.C.
Mackay, D. and PJ. Leinonen. 1975. Rate of evaporation of low-solubility contaminants from
water bodies to atmospheres. Environ. Sci. Technol., 7:611-614.
170
-------
Mackay, D. and S. Paterson. 1986. Model describing the rates of transfer processes of organic
chemicals between atmosphere and water. Environ. Sci. Technol., 20 (8):810-816.
Mackay, D. and W.Y. Shiu. 1981. A critical review of Henry's Law constants for chemicals of
environmental interest. J. Phys. Chem. Ref. Data., 10 (4): 1175-1199.
McLachlan, M. and D. Mackay. 1987. A model of organic chemicals in the Niagara River.
Contract Report KA401-06-0195, Environment Canada, Toronto, Canada.
Mill, T., W.R. Mabey, P.C. Bomberger, T.W. Chou, D.G. Herdry, and J.H. Smith. 1982.
Laboratory protocols for evaluating the fate of organic chemicals in air and water. EPA-600/3-
82-022. U.S. Environmental Protection Agency, Office of Research and Development,
Environmental Research Laboratory, Athens, Georgia.
Miller, G.C. and R.G. Zepp. 1979. Effects of suspended sediments on photolysis rates of
dissolved pollutants. Water Res., 13:453-459.
Mills, W.B., D.B. Porcella, M.J. Ungs, S.A. Gherini, K.V. Summers, L. Mok, G.L. Rupp, G.L.
Bowie, and D.A. Haith. 1985. Water quality assessment: A screening procedure for toxic and
conventional pollutants, Parts 1 and 2. EPA-600/6-85-002a,b. U.S. Environmental Protection
Agency, Office of Research and Development, Environmental Research Laboratory, Athens,
Georgia.
O'Connor, D.J. and J.P. Connolly. 1980, The effect of concentration of absorbing solids on the
partition coefficient. Water Res., 14:1517-1523.
O'Connor, DJ. 1983. Wind effects on gas-liquid transfer coefficients. J. Environ. Engin.,
109(9):731-752.
O'Connor, DJ. and R.V. Thomann. 1972. Water quality models: Chemical, physical and
biological'constituents. In: Estuarine Modeling: An Assessment. EPA Water Pollution Control
Research Series 16070 DZV, Section 702/71. pp. 102-169.
O'Connor, DJ. 1983a. The significance of gas exchange in water quality assessment and
planning. In: Gas Transfer at Water Surfaces. D. Reidel Publishing Co., Boston, Massachusetts.
171
-------
O'Connor, DJ. 1983b. Wind effects on gas-liquid transfer coefficients. J. Environ Engin
109(9):731-752. • S •»•
O'Connor, D.J., J.A. Mueller, and KJ. Barley. 1983. Distribution of kepone in the James River
estuary. J. Environ. Engin. Div., ASCE. 109(EE2):396-413.
Paris, D.R, W.C. Steen, G.L. Baughman and J.T. Barnett, Jr. 1981. Second-order model to
predict microbial degradation of organic compounds in natural waters. Appl. Environ.
Microbiol., 4(3):603-609.
Partheniades, E. 1992. Estuarine sediment dynamics and shoaling processes. In: Handbook of
Coastal and Ocean Engineering, Volume 3: Harbours, Navigation Channels, Estuaries, and
Environmental Effects, pp. 985-1071. Herbich, J. B., Ed. Gulf Publishing Company, Houston
Texas.
QEA (Quantitative Environmental Analysis). 1999. PCBs in the Upper Hudson River, Volume
2: A Model of PCB Fate, Transport, and Bioaccumulation. Montvale, New Jersey. Report
prepared for General Electric
Rao, P.S.C. and J.M. Davidson. 1980. Estimation of pesticide retention and transformation
parameters required in non-point source pollution models. In: Environmental Impact of Non-
point Source Pollution, pp. 23-67. Ann Arbor Science Publishers, Ann Arbor, Michigan.
Rathbun,R.E. 1990. Prediction of stream volatilization coefficients. J. Environ Engin ASCE
116:615-631. "
Robinson, N., Ed. 1966. Solar Radiation. Elsevier Publishing, Amsterdam, London, and New
York. 347pp.
Schnoor, J.L., C. Sato, D. McKechnie, and D. Sahoo. 1987. Processes, coefficients, and models
for simulating toxic organics and heavy metals in surface waters. EPA-600/3-87-015. EPA-
600/3-87-015. U.S. Environmental Protection Agency, Office of Research and Development,
Environmental Research Laboratory, Athens, Georgia.
Smith, J.H., W.R. Mabey, N. Bohonos, B.R. Hoh, S.S. Lee, T.W. Chou, D.C. Bomberger and T.
Mill. 1977. Environmental pathways of selected chemicals in freshwater systems. Part I:
172
-------
Background and experimental procedures. EPA-600/7-77-113. U.S. Environmental Protection
Agency, Office of Research and Development, Environmental Research Laboratory, Athens,
Georgia.
Tateya, S., S. Tanabe, and R. Tatsukawa. 1998. PCBs on the globe: Possible trend of future
levels in the open ocean environment. In: Toxic ^Contamination in Large Lakes, Volume m, pp.
237-282, N.W. Schmidtke (ed.), Lewis Publishers, Chelsea, Michigan.
Thomann, R.V., R.P. Wmfield, and J.J. Segna. 1979. Verification analysis of Lake Ontario and
Rochester Embayment three dimensional eutrophication models. EPA-600/3-79-094. U.S.
Environmental Protection Agency, Office of Research and Development, ERL-Duluth, Large
Lakes Research Station, Grosse He, Michigan.
Thomann, R.V. and J.A. Meuller. 1987. Principles of Surface Water Quality Modeling and
Control. Harper and Row Publishers, Inc., New York, New York.
Tsai, C.H. and W. Lick. 1987. Resuspension of sediments from Long Island Sound. Water Sci.
Technol., 21(6/7): 155-184.
Tsivoglou, E.E. and J.R. Wallace. 1972. Characterization of stream reaeration capacity. EPA-
R3-72-012. U.S. Environmental Protection Agency, Washington, D.C.
Velleux, M. and D. Endicott. 1994a. Development of a Mass Balance Model for Estimating
PCB Export from the Lower Fox River to Green Bay. Journal of Great Lakes Research,
20(2):416-434.
Velleux, M., J. Gailani and D. Endicott. 1994b. A User's Guide to IPX, the In-Place Pollutant
Export Water Quality Modeling Framework. U.S. Environmental Protection Agency, Office of
Research and Development, ERL-Duluth, Large Lakes Research Station, Grosse He, Michigan.
Velleux, M. and D. Endicott, J. Steuer, S. Jaeger, and D. Patterson. 1995. Long-term Simulation
of PCB Export from the Fox River to Green Bay. Journal of Great Lakes Research, 21(3):359-
372.
173
-------
Velleux, M., J. Gailani and D. Endicott. 1996. Screening-level approach for estimating
contaminant export from tributaries. ASCE Journal of Environmental Engineering, 122(6):503-
514.
Ward, D.M. and T.D. Brock. 1976. Environmental factors influencing the rate of hydrocarbon
oxidation in temperate lakes. Appl. Environ. Microbiol., 31 (5):764-772.
Wetzel,R.G. 1983. Limnology (Second Edition). Saunders College Publishing. Philadelphia,
Pennsylvania. 767pp.
Wetzel,R.G. 1975. Limnology. W.B. Saunders Co. Philadelphia, Pennsylvania. 743pp.
Whitman, R.G. 1923. A preliminary experimental confirmation of the two-film theory of gas
absorption. Chem. Metallurg. Eng., 29:146-148.
Xu, Y. 1991. Transport properties of fine-grained sediments. Ph.D. dissertation, University of
California-Santa Barbara, Santa Barbara, California.
Yang, C. S. 1996. Sediment Transport: Practice and Theory. McGraw-Hill, Inc. New York,
New York.
Yin, C. and J.P. Hassett. 1986. Gas-partitioning approach for laboratory and field studies of
mirex fugacity in water. Environ. Sci. Technol., 20(12): 1213-1217.
Zepp, R.G. and D.M. Cline. 1977. Rates of direct photolysis in aquatic environment. Environ.
Sci. Technol., 11:359-366.
Ziegler, C. K., W. Lick, and J. Lick. 1988. The transport of fine-grained sediments in the
Trenton Channel of the Detroit River. University of California-Santa Barbara, Santa Barbara,
California. Report to the U.S. Environmental Protection Agency, Office of Research and
Development, ERL-Duluth, Large Lakes Research Station, Grosse He, Michigan.
174
-------
APPENDIX A
IMPLEMENTATION OF THE MASS BALANCE EQUATION
IPX solves a finite difference approximation of Equation 1.1 for a model network that
represents the important characteristics of the real water body. This section explains the
derivation of the finite difference mass balance equation using the one-dimensional form for
convenience. Regrouping the terms in Equation 1.2 for mathematical convenience gives:
AST (A.D
where: ST = total source/sink rate = SL + SB + SK, g/m3-day
Q = volumetric flow = A Ux, m3/day
Assuming that derivatives of C are single-valued, finite, continuous functions of x, as in Figure
Al, then the Taylor's series expansion gives:
3C ^A 2
IS. +AX &'
1*0
J_
2'
--Ax3
2
(A.2)
(A.3)
Assuming that terms containing third and^ higher order powers of Ax are negligible in
comparison with the lower powers of Ax, Equations A.2 and A.3 can be subtracted to give:
2Ax
+Errorterm
(A.4)
The error term is of order Ax2. Referring to Figure Al, Equation A.4 states that the slope of the
line AB is equal to the slope of the tangent centered at P. This is known as the central-difference
175
-------
X0-2AX X0-AX X0 X0+AX Xn+2AX
X0-2AX X0-AX X0 X0+AX Xn+2AX
j-1
Figure Al. Definition sketch for finite difference equation.
176
-------
approximation. The slope at P may also be approximated by the slope of the line PB, giving the
forward-difference formula:
(A.5)
Ax
Similarly, the slope at P may be approximated by the slope of the line AP, giving the backward-
difference formula:
(A.6)
1*0
Ax
Equations A.4 and A.6 can be obtained from A.2 and A.4, respectively, by assuming that the
second and higher order powers of Ax are negligible. The error term for both the forward and
backward difference approximations is of order Ax.
Substituting the central difference approximation into the advection term of Equation A.1
gives:
2Ax
Similarly, the dispersion term becomes:
.2 Ax
(A.7)
(A.8)
'
Substituting the central difference approximation for — in Equation A.8 gives:
2Ax
2Ax
2Ax
(A.9)
. 177
-------
When applying the difference approximations to segment j of a network as in Figure Al
x0 corresponds to the center of j, XQ + A* to the interface between j and j+1, *0 - £ to the
interface between j-1 andj, *0+2A* to the center of j+l, and *0-2A* to the center of j-1. The
mass balance equation for segment j can be written:
_
+
(A. 10)
Multiplying through by Lj gives:
- c , )- /?._,. (c . - c,,,
(A. 11)
where: Vj = volume of segment j = Aj Lj, m3
R = dispersive flow = E A/LC, m3/day
LC = characteristic length, m
The interfacial concentrations CjJ+1 and Cj.y must be expressed in terms of the
concentrations:
segment
(A12)
(A13)
where: o> = numerical weighting factor (advection factor) between 0 and 1
Specifying co = 0 gives a backward difference approximation for the advective term, oo = 0.5
gives a central difference approximation, and oo = 1 a forward difference
Equation A. 11 can be extended to the multi-dimensional form used in IPX. Consider i
segments adjoining segment j. Interfaces are denoted ij. The general equation becomes:
178
-------
(A.14)
>Kj
where: Qij
= flow, defined as positive when leaving segment j, and negative when entering
j, nrVday
Equation A.14 is the general expression used in IPX to evaluate the mass derivatives for
every segment j during each time-step t between the initial time to and the final time tF. Given
concentrations and volumes at time t, IPX calculates new masses at t + At using the one-step
Euler scheme:
u (A-15)
where: At = time-step (typical values between 0.001 and 0.5 days), days
Given new masses at time t + A/.IPX computes the new concentrations by dividing the new
masses by the new volumes:
'j,t+&t
(A.16)
The new volumes are calculated internally from the specified (or computed) flow fields using the
principle of continuity.
During normal simulations, IPX prohibits segment concentrations from going negative
and causing numerical instability of the solution. A negative concentration might be calculated
for constituents with low concentrations in the vicinity of significant spatial gradients. If a
calculated mass derivative would drive a segment concentration below zero at t + Al, IPX
maintains a positive segment concentration by halving the mass that is present at time t.
Experience shows that this procedure is generally acceptable. Users can activate or disable this
negative solution procedure through the use of negative solution option in Data Group A. If
negative concentrations and instability occur when this option is disabled, the simulation can be
repeated with a smaller time-step.
179
-------
-------
-------
United States
Environmental Protection Agency/ORD
National Health and Environmental
Effects Research Laboratory
Research Triangle Park, NC 27711
Official Business
Penalty for Private Use
$300
Please make all necessary changes on the below label,
detach or copy, and return to the address in the upper
left-hand corner.
If you do not wish to receive these reports CHECK HEREQ ;
detach, or copy this cover, and return to the address in the
upper left-hand corner.
EPA/600/R-01/074
------- |