&EPA
United States
Environmental Protection
Agency
Environmental Research
Laboratory
Gulf Breeze FL 32561
EPA-600/3-84-077
August 1984
Research and Development
A User's Guide for
WASTOX, a
Framework for
Modeling the Fate of
Toxic Chemicals in
Aquatic
Environments
Pa rt 1 :
Exposure
Concentration
-------
EPA-600/3-84-077
Auqust 1984
A USER'S GUIDE FOR
WASTOX, A FRAMEWORK FOR MODELING THE FATE OF TOXIC CHEMICALS
IN AQUATIC ENVIRONMENTS
PART 1: EXPOSURE CONCENTRATION
BY
John P. Connolly
Richard P. Winfield
Environmental Engineering and Science
Manhattan College
Bronx, New York 10471
Cooperative Agreement No. R807827/R807853
Project Officer
Parmely H. Pritchard
Environmental Research Laboratory
Gulf Breeze, Florida 32561
William L. Richardson
Large Lakes Research Station
Environmental Research Laboratory
Grosse lie, Michigan 48138
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
GULF BREEZE, FLORIDA
AND
GROSSE ILE, MICHIGAN 48138
-------
DISCLAIMER
The information in this document has been funded wholly or in part
by the U.S. Environmental Protection Agency (EPA) under contract R807827
with J. P. Connolly of Manhattan College, the Bronx, New York. It has been
subject to the Agency's peer and administrative review, and approved for
publication as an EPA document. Mention of trade names or commercial
products does not constitute endorsement or recommendation for use.
The WASTOX computer code has been tested against other computer
programs to verify its computational accuracy. Nevertheless, errors are
possible. The EPA assumes no liability for either misuse of the model or
errors in the code. The user should verify the code before using it.
ii
-------
FOREWORD
The protection of estuarine and freshwater ecosystems from damage
caused by toxic organic pollutants requires that regulations restricting
the introduction of these compounds into the environment be formulated on
a sound scientific basis. Accurate information describing the potential
exposure of indigenous organisms and their communities to these toxic
chemicals under varying conditions is required. The Environmental Research
Laboratory, Gulf Breeze, contributes to this information through research
programs aimed at determining:
. the effects of toxic organic pollutants on individual species,
communities of organisms, and ecosystem processes.
. the fate and transport of toxic organics in the ecosystem.
. the application of methodologies which integrate fate and effects
information to predict enviromental hazard.
The magnitude and significance of chemical contamination of aquatic
environments are increasingly evident. The potential persistence and
possible accumulation of those chemicals in aquatic food chains means
that the impact on the health and activities of man is more direct.
Therefore, the ability to predict exposure concentration, bioaccumulation,
and chronic exposure is critical to our efforts in hazard assessment.
Mathematical models provide a basis for quantifying the inter-relationships
among the various physical, chemical, and biological variables that
affect the fate, transport, and bioaccumulation of toxic chemicals. They
also provide a mechanism for extrapolating laboratory information to the
environment and a rationale and conceptually relevant basis for decisionmaking.
The modeling framework WASTOX described in this manual permits the
user to examine the transport of a toxic chemical dissolved in water or
sorbed to sediments, its transfer from one environmental medium (air,
water, sediment, biological tissue) to another, and its reactivity relative
to the processes of biodegradation, hydrolysis, photolysis, and oxidation.
(Part I deals with the exposure concentration component of the WASTOX
model. Part II will describe a food chain bioaccumulation model, and Part
III will relate exposure and bioaccumulation to toxic effects.)
F. Enos
Director
Environmental Research Laboratory
Gulf Breeze, Florida
iii
-------
ABSTRACT
A computer program was developed for modeling the fate of toxic
chemicals that are discharged to natural water systems. The program
permits the user to model the water and sediment transport in a natural
water system and the movement and decay of chemicals discharged to that
system. Either the equilibrium distribution of chemicals continually
discharged to the system or the concentrations in the system as a function
of time may be computed. From one to three types of solids may be
considered.
The reaction of the chemical and its transfer among phases are
computed from specified characteristics of the chemical and environmental
parameters of the system. The processes considered include photolysis,
hydrolysis, biodegradation, volatilization and adsorption. Adsorption to
the solids types included in the model is described as a local equilibrium
process defined by a partition coefficient and the local solids
concentration. All other processes are defined in terms of reaction rates.
WASTOX does not explicitly specify each of the transport processes
that may affect the chemical or solids. Transport is considered either
as an advective process defined by a flow or a mixing process defined by
a dispersion or exhange. Specification of separate transport processes
is made by the user by defining up to nine sets of flows and dispersions
(termed fields). Each field is applied to either dissolved chemical or
adsorbed chemical and solids, or both. For example, dispersion within the
stationary sediment is limited to dissolved chemical; therefore a field
of dispersions would be specified by the user and applied only to this
component. The user of such a non-specific transport structure permits
construction of models consistent with the understanding of the particular
natural water system and the question being addressed.
WASTOX is sufficiently general to be applied to all types of natural
water systems. It has been successfully applied to the James River
Estuary, the Great Lakes, and the U.S. Environmental Protection Agency
(EPA) experimental stream channels at Monticello, Minnesota.
iv
-------
CONTENTS
Foreword iii
Abstract iv
Figures vii
Tables viii
Acknowledgements ix
1. Introduction 1
2. Fundamental Equations 4
2.1 Water Column 4
2.2 Sediment Bed 7
2.3 Applications to a Specific Water Body 10
3. Modeling Framework 13
3.1 Finite Difference Approximation 13
3.2 Kinetics 16
4. Structure of Computer Code 27
4.1 Overview 27
4.2 Subroutines 28
5. Preparation of Data Input 33
5.1 Introduction 33
5.2 Summary of Card Groups 34
5.3 Card Group A - Model Identification and System Execution
Options 36
5.4 Card Group B - Exchange Coefficients 38
5.5 Card Group C - Volumes 52
5.6 Card Group D - Flows 55
5.7 Card Group E - Boundary Conditions 61
5.8 Card Group F - Forcing Functions 66
5.9 Card Group G - Parameters 71
5.10 Card Group H - Constants 72
5.11 Card Group I - Miscellaneous Time Functions 75
5.12 Card Group J - Photolysis 77
5.13 Card Group K - Initial Conditions 79
5.14 Card Group L - Stability Criteria 80
5.15 Card Group M - Print Control 80
5.16 Card Group N - Integration Information 80
5.17 Card Group 0 - Display Parameters 81
v
-------
6. Example Applications 82
6.1 Kepone - James River, Va 82
6.2 Plutonium, PCB - Great Lakes 88
6.3 Pentachlorophenol - EPA MERS Channels 90
7. Operational Considerations 96
7.1 Acquisition Procedures 96
7.2 Installation Procedures 96
7.3 Testing Procedures 96
7.4 Machine Limitations 97
References 98
Appendix 1, Glossary 100
Appendix 2, Listing of Kinetic Subroutine WASPB 108
Appendix 3, Test Program Input and Output. ............ .117
VI
-------
FIGURES
Number Page
1. Fluxes of toxic chemical associated with the bed 11
2. Typical characteristics of natural water systems:
a) Typical transport regimes, b) Typical bed conditions 12
3. Segmentation of the James River used in the Kepone analysis 84
4. Sequence of steps in James River Kepone analysis (23) 86
5. Comparison of observed and computed suspended solids and Kepone
for the 1000 cfs freshwater flow 88
6. Segmentation of the Great Lakes and Saginaw Bay 90
7. Comparison of 1971-1977 observed and calculated water column
239 240
' Pu concentration for three conditions of the particulate
settling velocity 92
8. Calculated surface sediment PCB concentration (ng/g) for condi-
tions on external bed and volatilization rate and comparison to
observed data 93
9. Model Comparison with September 16 and 17 Rhodamine WT Dye Data.. 95
10. Time Variable Model Comparison with MERS Channel J Total PGP
data - Stations 1,2,3,4,5,6,7,8,13,14. June 8 to 10, 1982 96
vii
-------
TABLES
Number Page
1 Subroutines comprising WASTOX 29
viii
-------
ACKNOWLEDGEMENTS
This work is a fruit of the toxic chemical modeling research conducted
at Manhattan College. As such, it reflects the expertise and direct assis-
tance of our colleagues in the water quality modeling group: Dominic Di
Toro, Donald O'Connor and Robert Thomann. The assistance of Paul Rontanini
in computer programming is also gratefully acknowledged.
The efforts of the project officers involved in this work were signifi-
cant. Parmely Pritchard provided constant support and review for both the
development and application of the model. William Richardson's long term
commitment to Great Lakes water quality modeling was of great importance to
the undertaking and completion of this work.
WASTOX was developed from a version of WASP used at Manhattan College
and we acknowledge the efforts of those involved in WASP's development,
particularly James Fitzpatrick of HydroQual, Inc. and Gregory Clemens of
Manhattan College.
The dedicated and friendly assistance of Eileen Lutomski and Margaret
Cafarella who typed this report is greatly appreciated.
IX
-------
SECTION 1
INTRODUCTION
In problems concerning toxic chemicals in aquatic environments, one of
two regulatory questions is generally posed: (1) what is the maximum allow-
able discharge of a toxic chemical to a water system for a specified use of
the system, and (2) how long will it take a contaminated natural water system
to recover to some specified level? The former question acknowledges that
chemical discharges are an unavoidable component of our technological society
and addresses the issue of allocation of the wasteload capacity of a natural
water system relative to a given use and associated criteria. For this case,
the relevant chemical concentration for regulatory evaluation is the equilib-
rium or steady-state value. The latter question deals with an unregulated or
accidental discharge whose adverse effects must be assessed in terms of future
use of the contaminated system. For this case, the time dependent behavior
of the chemical in the system must be estimated.
To address these questions the physical, chemical and biological charac-
teristics of both the system and the chemical must be specified. These char-
acteristics determine the fate of the chemical by defining the rates of trans-
port, transfer, and reaction.
Transport is the physical movement of the chemical caused by the net ad-
vective movement of water, mixing, and the scouring and deposition of solids
to which the chemical may be adsorbed. It is specified by the flow and dis-
pension characteristics of the natural water system and the settling velocity
and resuspension rate of the solids in the system.
Transfer is the movement of the chemical between the air, water, and
solid phases of the system. It includes the processes of volatilization, ad-
sorption, and bioconcentration-bioaccumulation.
Reaction is the transformation or degradation of the chemical. It in-
cludes biodegradation and the chemical reaction processes of photolysis -
hydrolysis, and oxidation.
The characteristics of the chemical are determined by a laboratory eval-
uation of the transfer and reaction processes. These controlled experiments
provide data from which the kinetics of reaction and transfer may be deter-
-------
mined. This information may then be combined with a theoretical development
for transport to form an analysis framework or mathematical model. The model
is the means by which the questions posed previously may be answered for a
specific water body.
As part of projects to determine the fate of toxic chemicals in estuaries
and the Great Lakes, a general model has been developed. This model calcu-
lates the time-variable or steady-state concentrations of dissolved and ad-
sorbed chemical in any natural water system. It is named WASTOX, an acronym
for Water quality Analysis Simulation of TOXics, and was developed from the
Hydroscience, Inc. water quality model WASP (20).
The equations that comprise the model are developed using the principle
of mass conservation, including the inputs with the transport, transfer and
reactions components. The general expression for the mass balance equation
about a specified volume, V, is:
V ^- = J + £R + £T + £W (1-1)
dt
in which
c = concentration of the chemical
J = transport through the system
R = reactions within the system
T = transfer from one phase to another
W = inputs
Equation (1-1) describes the mass rate of change of the chemical due to the
net effect of the various fluxes and transformations. The purpose of ex-
pressing the transfer rate (T), distinct from the transport (J) and reaction
(R), is to provide a basis for the development of the equations, which des-
cribe more fully the relevant phenomena. The transport, reaction and trans-
fer terms may be positive or negative depending on the direction of kinetic
routes.
To obtain a time variable solution of Equation (1-1), it is necessary to
employ numerical methods. The approach used is to approximate the spatial
derivatives that define the transport component of the equation and the temp-
oral derivative with difference expressions. This is equivalent to separating
-------
the natural water system being modeled into a series of completely-mixed sec-
tions. The size of the sections is controlled by the time step constraints
imposed by the stability limitations of the numerical scheme and the extent
of spatial resolution desired.
To obtain the steady state solution of Equation (1-1) difference expres-
sions are used to approximate the spatial derivatives and the time
derivative, dc/dt is set equal to zero. The resulting set of simultaneous
algebraic equations are solved using the Gaussian Elimination technique.
The model is general in nature. The natural water system under study
may be sectioned or segmented as the user chooses. Transport is defined by
dispersion coefficients and flows that are specified in a series of fields
that are individually applied to both the particulate and dissolved chemical
(e.g., hydrodynamic flow), the particulate chemical only (e.g.,
resuspension) , or the dissolved chemical only (e.g., interstitial water
diffusion).
The reaction and transfer components are computed in a single subroutine
separate from the transport and integration computations. This kinetic sub-
routine has been developed to provide the user with the capability to alter
the transfer and reaction formulations or to add additional kinetic processes.
The kinetic subroutine returns to the main program that portion of the
dc
time derivative Orr) contributed by the transfer and reaction components.
Changes may be made by adding terms to this derivative computation. The
variables that define a specific kinetic process are transmitted to the
kinetic subroutine through a segment dependent array of parameters and an
array of constants. These may be altered by the user through respecification
of the equivalence statements in the kinetic subroutine. A copy of the
present kinetic structure is included in Appendix A.
The simple Euler-Cauchy method is used for the integration of the time
derivative. Results are filed at a user specified print interval and dis-
played at the conclusion of the program execution.
-------
SECTION 2
FUNDAMENTAL EQUATIONS
One of the most distinguishing characteristics of toxic substances is
the partitioning between dissolved and particulate components. Thus equa-
tions are developed for each of these components and in addition for those
solids which provide sites for the adsorption of the substance. The analysis
involves therefore, the solution of at least three simultaneous equations,
describing the concentration of the various components in the water column.
Furthermore, for interaction with the bed, additional sets of equations are
developed to account for distribution in the benthal layer and its effect on
water column concentrations. Given these concentrations, the dissolved and
the particulate in the water column and in the bed, the distribution through
the food chain is then analyzed.
2.1 WATER COLUMN
Consider the concentration, c, to be the dissolved component of the chem-
ical in water. It interacts with the particulate concentration., p, through
an adsorption-desorption reaction with the solids. The particulate concentra-
tion is defined as:
p = rm (1)
3
p = particulate concentration M/L
r = mass of chemical/unit of solids mass M/M
3
m = concentration of the solids M/L
The equations governing the distribution of the dissolved and particu-
late components in a natural water system may be written as follows:
9c 8 ,„ 9c, 9 ,„ 9c. ^ 9 ,_, 9c. 9
3t = 31 (Ex 9^ + I? (Ey 3F> + 3F (Ez 9^ ' fcF V
~ "97 Uyc " fc V ' Kom° + K2P * VX'y'Z't}
Uyp - V + Wsp + Komp - K2P ±
-------
in which
3
c = concentration of the dissolved component [M/L ]
2
E = dispersion coefficient [L /T]
u = velocity [L/T]
w = settling velocity of the particulates [L/T]
3
/
-1
s 3
K = adsorption coefficient [L /(M«T)]
-
K,, = desorption coefficient [T ]
x,y,z = coordinate directions
t = time [T]
S = sources and sinks of the component due to reactions and
phase transfers
The first three terms in each equation represent dispersion or mixing
due to temporal and spatial velocity gradients and density differences within
the natural water system. The next three terms represent the longitudinal,
lateral and vertical advection, respectively. The seventh term in equation
(3) accounts for the vertical advection of the particulate component due to
settling. The following two terms define the rates of adsorption and desorp-
tion, respectively. The last term accounts for the chemical and biological
reactions and volatilization that may produce or degrade the component.
The vertical boundary conditions are assigned at the air-water and
water-bed interfaces, at which concentration and flux conditions are speci-
fied. The flux may be negative or positive, representing net deposition or
scour at the bed, respectively, or gas-liquid exchange at the air-water sur-
face. An additional source due to precipitation may be effective under cer-
tain conditions.
Note that since the distributions of the dissolved and particulate com-
ponents depend on the concentration of the solids, m, an expression equivalent
to Eqs. 2 and 3 must be written for the distribution of the solids:
3m _ 9 frf 3mN 9 (v 9iiK 9 ,„ 9m.. 9
9t ~ 9x <• x 3x; 9y ^ y 9y' 9z l z 9z' 9x
3 3 9
- — u m - — u m + — w m
9y y 3z z 9z s
(4)
-------
If more than one solid is considered, i.e., several solid classes, then two
additional equations must be included for each additional component to
describe the particulate and solids concentrations.
The particulate and dissolved concentrations can be summed to yield the
total concentration of chemical, c :
CT = 6c + p (5)
or if several solids types are considered:
n
c - 6c + Zp (6)
i=lx
where 6 is porosity (water volume/total volume) .
The particulate concentration is related to the concentration of solids, m,
as shown by equation (1):
p = rm
Equation (6) may therefore be written as:
n
c = 6c + Z r.m. (7)
1=1 1 X
If adsorption equilibrium is assumed then the relationship between the dis-
solved chemical concentration and the mass concentration of adsorbed chemi-
cal, r (M chemical/M solid), is specified by the sorption isotherm. For many
chemicals in natural water systems the isotherm is linear and is expressed in
terms of a partition coefficient 1F (L /M) :
r = He (8)
Combining equations (7) and (8), the total concentration may be ex-
pressed as:
CT = c[6 + ZILnu] (9)
The fraction of total chemical that is dissolved, f ,, may then be written as:
d
P CT e + ZILIIK
do)
-------
Similarly, the fraction of total chemical adsorbed to a particular
solids type j, f . , is given as:
p f m
f = -JL = - U - (11)
rPj cx 0 + ZILnu <• ;
It is evident from equations (10) and (11) that the distribution of the dis-
solved and particulate components may be determined directly from the distri-
butions of total chemical concentration and solids concentration. A reduc-
tion in the number of equations included in the modeling framework is realized
since only one equation for the chemical and one for each solids type is re-
quired whereas previously an equation for dissolved chemical and two equations
for each solids type, one for particulate chemical and one for the solids,
were required.
An equation describing the distribution of total chemical is obtained by
summing the equations for dissolved and particulate (equations (2) and (3)):
3c_
(12)
2.2 SEDIMENT BED
The transport components of the bed may include a longitudinal advection
and a vertical dispersion or mixing. The advective component of transport in
the bed is significant in rivers and estuaries where it is induced by shearing
stress of the overlying flowing water, with which it is in contact. The hori-
zontal velocity of the bed is a maximum at the water interface and diminishes
in depth to a point at which the bed may be regarded as stationary. The ver-
tical dispersion follows a similar pattern. A mass balance about an elemental
volume within the bed results in the following equation defining the solids
distribution:
3m 3 f f , 3 . u(z) 3Bm
= (e(x) '
-------
in which
2
e = vertical mixing coefficient (L /T)
u = horizontal velocity (L/T)
B = width of the bed (L)
The boundary condition in the vertical planes are specified at the water-
bed interface and at that depth of zero horizontal velocity,, which is the
interface between the transport layer and the stationary bed. The former is
the net flux, J, due to the settling of suspended solids from the overlying
water and the resuspension from the transport or surface layer of the bed:
b
= w m - w m, (14)
s w u b
in which m, is the concentration of particles in the sediment layer, m is
b w
the concentration of particles in the water column and w is the entrainment
u
velocity of these sediment particles across the interface. A net flux from
the water to the sediment layer implies that the mean settling is greater
than the resuspension of particles and conversely a net flux from the sedi-
ment layer to the water stipulates that uplift is greater than settling. The
net effect of this flux is a change in the elevation of the water-bed inter-
face, which rises or falls depending on the relative magnitude of the two
terms in the above equation. A volumetric change results in the bed that
must be taken into account in the mass balance.
Assuming vertical uniformity of solids concentration in the surface bed
layer, a mass balance results in the following equation defining the solids
distribution:
3Vm
-r-5- = w m A - w m, A + -£i_ -2- (Bm. ) (15)
8t s w u b B(x) 3x b
in which V is the bed volume and A is the surface area of the bed. The ver-
tical dispersion term of equation (14) does not appear because; by definition
there is no vertical solids gradient. It is replaced by the flux condition
at the boundary.
Replacing the volume term by the product of the length and cross-sec-
al area, A , of the layer, expanding the
dividing through by the cross-sectional area:
tional area, A , of the layer, expanding the left side of the equation, and
-------
3m, m, 9A w w .
b + b c_ _ s m s_m +u.f_ /Bm \
c
in which h = depth of the layer. Assuming the width of the bed is constant
in time:
, 3A ... w,
i c _ 1 on _ d
A~~ '3t h" dT = h~
c
•7— = w, = the rate of change of the bed elevation
dt d
The rate of change of the bed elevation is referred to as a sedimenta-
tion velocity. It may be envisioned as the velocity at which the water-bed
interface approaches the stationary bed if erosion is predominant or at which
the interface moves away from that datum if settling is the significant term.
Assuming the concentration of the bed, m, , is approximately constant, the
above reduces to, under steady-state:
w m m,
0 = -~-2- - T-5- (w + w.) + urn, b (18)
n h u d b
in which
1 dB
For a stationary bed u is zero and the last term in equation (18) drops out.
Intrepreting the bed elevation change as described, the flux about a lower
bed layer is:
Vb
in which
w, = sedimentation velocity in the lower bed layer
m = average concentration of the solids in the lower bed layer.
5
A mass balance on particulate chemical in the bed follows that for solids
directly. For dissolved chemical a flux at the sediment-water interface, J,,
specified in terms of a dispersion coefficient, E, is incorporated:
-------
J = -E — 1 (20)
d 9z 'z=sediment-water interface
Assuming a constant concentration of solids within a layer of depth H adjacent
to the water column, the equation describing the distribution of total chemi-
cal in that layer is:
9CTK " f
-IT • r fpcT - / (wu + VCT, + Ub<\ - IT <21)
w b b bbww
where £ is the average depth of the surface bed layer and the adjacent water
column layer. Dispersion of dissolved chemical between the surface and lower
bed layers adds a term to equation (21) that is analogous to that for disper-
sion with the water column. Figure 1 illustrates the components of equation
(21).
2.3 APPLICATION TO A SPECIFIC WATER BODY
The water column and bed equations are general in nature and are usually
simplified when applied to a specific water body. Typical transport regimes
for the various types of water bodies are shown in Figure 2. In rivers advec-
tion is dominant and the water column is fairly well-mixed vertically and lat-
erally. Therefore, the lateral and vertical spatial derivatives in equation
(12) may be dropped reducing the equation to a single dimension. Settling is
included as a boundary condition at the water-bed interface.
In estuaries vertical concentration gradients may be important and only
the lateral spatial derivative is dropped. Equation (12) is therefore reduced
to a two-dimensional equation.
In lakes, the bed equation is generally simplified because the low velo-
cities do not cause significant horizontal movement of the bed. The horizon-
tal velocity term in equation (21) is dropped reducing the equation to a
single dimension.
Further simplifications are sometimes made depending on the space and
time scales of interest. Several of these are illustrated in the applications
of WASTOX presented in Section 6.
10
-------
t
WATER COLUMN
SEDIMENT
PROCESS
DEPOSITION OF SOLIDS
MASS FLUX TERM
SCOURING OF SOLIDS
Q) NET SEDIMENTATION
HORIZONTAL BED MOVEMENT
DIFFUSION OF DISSOLVED CHEMICAL
(fdbCTb~
FIGURE 1. FLUXES OF TOXIC CHEMICAL ASSOCIATED WITH THE BED
11
-------
LAKES AND
COASTAL WATERS
FREE-FLOWING
STREAMS t RIVERS
WIND-DRIVEN
CURRENT
ESTUARIES AND
LAKE TRIBUTARY
TIPAU.Y
AVERAGE!?
VELOCITY
T7>
rrrrrrr
U - HORIZONTAL VELOCITY
U - VERTICAL VELOCITY
6 - VERTICAL MIXING
TYPICAL TRANSPORT REGIMES
WATER
BED
STATIONARY
SEDIMENTING BED
T>
MIXED
• LAYER
MIXED LAYER
SEDIMENTING BED
U. \ \ Uu
ffff- Ub TRANSPORT LAYER
MIXED LAYER
STATIONARY IED
LIST OF SYMBOLS
I -0. - SETTLING
t ^u - RESUSPENSION
< ^d - SEDIMENTATION
1C1 Ub- MIXING
TYPICAL BED CONDITIONS
BED TRANSPORT
MIXED LAYER
SEDIMENTING BED
FIGURE 2. TYPICAL CHARACTERISTICS OF NATURAL WATER SYSTEMS:
(a) Typical transport regimes, b) Typical bed
conditions
12
-------
SECTION 3
MODELING FRAMEWORK
3.1 FINITE DIFFERENCE APPROXIMATION
WASTOX is in essence a computer program that solves the equations pre-
sented in Section 2. The solution of these equations is accomplished using
finite difference approximations to the spatial derivatives in identical
fashion to the Water Quality Analysis Simulation Program (WASP) (20). The
water body is divided into completely-mixed segments in which the concentra-
tions of total chemical and solids are calculated by the finite-difference
approximations of their mass balance equations. For total chemical the equa-
tion for segment i is as follows:
dc E A
'i 3T 'I[Vic3 - f i.x'i + ? -iJ-T1 «: - V
J K J 1,J
(22)
E. A
- w A. f c. + w A f c + -~ '— (f c - f, c.) ± S
"i P ' ~f P ~i ' P "i P P "f
where
V = volume of segment i
c. = total chemical concentration in segment i (subscript t dropped
for convenience)
Q. . = hydrodynamic flow from segment j to segment i
3 »!
E. . = dispersion coefficient between segment i and segment j
1 > J
A. . = cross-sectional area between segments i and j
x»3
L. . = characteristic length or mixing length of segments i and j
i»j 2
w = settling velocity of particulates from segment i to segment H
w = resuspension velocity of particulates from segment Jl to segment i
H,i (included only if the i,£ interface is a water-bed interface)
(Note: E. is fundamentally the product of porosity and the dispersion
coefficient and would be included only if i,SL is a water-bed or bed-bed inter-
face)
13
-------
WASTOX applies equation (22) to each of the segments the water body is
divided into. However, the explicit separation of hydrodynamic flow, set-
tling, resuspension, sedimentation, water column dispersion and interstitial
water diffusion indicated by equation (22) is not formalized in the computer
code. The program considers all transport as either a flow or a dispersion.
Including different transport processes is accomplished through the use of
multiple sets of flows and dispersions, termed fields. These fields are
applied to either dissolved chemical or adsorbed chemical and solids, or both.
For example, dispersion within the bed is limited to dissolved chemical and a
separate grouping or field of dispersions is inputted by the user and applied
only to this component. The type and number of transport processes included
is the decision of the user based on the fundamental processes; indicated in
Section 2, the characteristics of the water body being modeled, and the par-
ticular segment, i.e., water column or sediment, at or away from the water-bed
interface.
For time variable solutions equation (22) is solved using a forward dif-
ference approximation for the time derivative. The concentration at time n+1
is equal to the concentration at time n plus the derivative evaluated at time
n multiplied by the time step, At;
(23)
To maintain a stable solution for the numerical scheme the size of the
time step must be restricted. The stability criteria and numerical error
associated with the solution scheme are the same as for WASP (20).
For computational stability a necessary condition is that
At ' Mln S.j * 4,j + Kivi} (24>
where
(25)
R. . = exchange coefficient = T
1>J i,j
K. = reaction rate of the chemical
14
-------
However, since K is, for many applications, non-linear and
time-dependent itself, and therefore may be difficult to
evaluate, the following criteria may be used for choosing
the integration step-size
At
J J
providing that
Max (K.V.) « IQ.. + ZR...
The use of the completely mixed finite segment approximation
in conjunction with a backward difference spatial approxima-
tion, introduces a numerical error (sometimes referred to as
numerical or pseudo-dispersion) into the model. The extent
of this effect is given by
E
num 2
where E is the numerical error expressed as a pseudo-dis-
persion coefficient and Ax is the length of the segment in
the direction of the velocity u. For some applications, es-
pecially in well mixed estuaries, where the time scale of
importance is on the order of days to seasons rather than
hours, and where u is the net advective freshwater velocity
(usually small), the effect of the numerical error, E , is
generally not significant.
For steady state solutions the time derivative in Equation (22) is set
equal to zero yielding an algebraic equation for each segment. This set of
simultaneous equations may be written as:
[A](c) = (w) (26)
where [A] is an nxn matrix (n = number of segments) including the transport
and reaction terms of Equation (22) , (c) is the vector of segment concentra-
tions, and (w) is the vector of mass loadings to each segment. The steady
state concentrations can be obtained by inversion of the [A] matrix, i.e.,
(c) = [Al'w) (27)
15
-------
Direct inversion of the [A] matrix is valid only if the reaction terms in-
cluded in the model are linear, as they are in WASTOX. If non-linear reac-
tions were included an iterative solution scheme would have to be used.
3.2 KINETICS
Kinetic expressions that define the time rate of transfer of chemical
between phases and transformation to degradation products comprise the terms
S, and S« in equations (2) and (3), respectively. The processes of volatili-
zation, photolysis, chemical hydrolysis and microbial degradation are in-
cluded.
3.2.1 Volatilization
The transfer of a volatile substance between air and water occurs when
the concentration in one medium is not in equilibrium with the concentration
in the other. The rate of transfer per unit area in each phase ±s propor-
tional to the driving force, which is the partial pressure gradient in the
gas phase and concentration gradient in the liquid phase. Under steady state
conditions the unit transfer is:
J - VCU - C£] = KG[Gg - Cgi] (28>
The subscripts £ and g refer to the liquid and gas phases, respectively. c0 .
Jci
and c . are the concentrations at the air-water interface, c. and c are the
gi *• g
bulk concentrations, and K. and K are the mass transfer coefficients between
* g
the interface and the bulk phase. Assume the concentrations are at equilib-
rium at the surface of contact between the two phases, i.e., c . = He.., in
which "H" equals Henry's constant. Solving for c. and substituting this ex-
pression into the flux equation yields:
in which
j = k£[c /H - c£] = k [c - Hc£] (29)
& O &
k. and k are the overall transfer coefficients
16
-------
For exchange between the atmosphere and natural water systems, the ratio
of the liquid film coefficient, IL , to the gas film coefficient, KG> is in
the range of 0.0001-0.005. Therefore for substances with Henry's constant
greater than 10, the resistance to transfer is localized in the aqueous film
-4
and for Henry's constant less than 10 , in the gas film.
The gas phase concentration of most toxic chemicals is essentially zero
and the flux may be written as:
j - v, <31>
It is evident that to evaluate the transfer rate, k , values for the
x«
liquid phase transfer coefficient, Hi., the gas phase transfer coefficient,
K , and the Henry's constant, H, are required.
(j
In the classic two-film theory of transfer (1), from which the above
equations were derived, the transfer coefficients may be defined as:
K = | (32)
iv 6n.iv!'ri 6 is the thickness of the diffusional or viscous layer and D is the
molecular diffusivity of the chemical. Since the relationship applies to
both air and water, the subscripts are omitted.
An approach based on the concept of surface renewal has been presented
(2). The analysis, postulated on a random replacement of the elements at the
interface, leads to a definition of the transfer coefficient as follows:
K = vfi? (33)
in which r = rate of surface renewal. A definition of the renewal rate has
been proposed (3) as the ratio of the intensity to the scale of turbulence
17
-------
in the vicinity of the interface. The vertical velocity fluctuation, v, is
representative of the intensity and the mixing length, &, of the scale, and
the renewal rate is the ratio of these:
r = Y (34)
Surface renewal implies turbulent flow extending to the interface with
no diffusional or viscous layer. It is applicable to the liquid phase of
flowing waters such as streams and estuaries. For such systems the rate of
surface renewal, r, has been shown (3) to be approximated by the ratio of the
velocity, u, to the depth, h:
r = ^ (35)
The depth, h, is the depth of the well-mixed layer, which is the total depth
in streams and the depth of the surface layer in stratified estuaries. Sub-
stituting equation (35) into equation (33) yields:
n
<> (36)
For more quiescent systems such as lakes, where surface renewal does not
extend to the interface, a viscous sublayer exists. In these systems mass
transfer is controlled by diffusion through the liquid and gas sublayers as
indicated by equation (32) . Because the structure of the sublayers is af-
fected by the shearing action of winds on the water surface the transfer
coefficients are related to wind activity.
The liquid mass transfer coefficient due to winds may be expressed as
(4):
DL >a W10
w w
i
L
where
D = diffusion coefficient of the chemical in water
LJ
v = viscosity of water
18
-------
c, = drag coefficient
d
p = density of air
3
p = density of water
w
W1. = wind velocity at 10 meters above the water surface
A = coefficient related to the thickness of the liquid viscous
Li
sublayer
The gas term K is primarily affected by wind action and the structure
d
of the diffusional sublayer in the interfacial region. By analogy to the
liquid mass transfer coefficient relationship to winds, K is expressed as:
G
a w G
where
D., = diffusion coefficient of the chemical in the atmosphere
(j
v = viscosity of air
3
A = coefficient related to the thickness of the gas viscous sublayer
The gas mass transfer coefficient for a flowing system like an estuary
when there is no wind is undefined. It is estimated at 50-100 m/d by con-
trast to values in the order of 1000 m/d for large open water systems under
moderate wind speeds.
Winds are most significant for open water bodies such as lakes and bays.
Most flowing systems are sheltered and are only significantly affected by
winds parallel to the longitudinal axis of the system. The lower portion of
an estuary, i.e., the harbor area, is generally more open than the upstream
area and winds may be significant to air-water transfer in these regions.
If the model is directed to stream or estuarine environments (flowing
systems) , the liquid mass transfer coefficient is specified by equation (36).
Wind effects are not considered in specifying the gas phase mass transfer
coefficient. It is set at a constant value of 100 m/d, acknowledging that
this may not be valid for harbor regions or systems subjected to strong winds
along their main axis. For lake and ocean environments the mass transfer
coefficients are calculated by equations (37) and (38).
19
-------
3.2.2 Photolysis
Photodecomposition (photolysis) is the transformation or degradation
of a compound that results directly from the adsorption of light energy. It
is a function of the quantity and wavelength distribution of incident light,
the light absorption characteristics of the compound and the efficiency at
which adsorbed light produces a chemical reaction. Photolysis is classified
into two types that are defined by the mechanism of energy absorption.
Direct photolysis is the result of direct absorption of photons by the com-
pound. Indirect or sensitized photolysis is the result of energy transfer to
the compound from some other molecule that has adsorbed radiation.
Direct Photolysis
In natural environments, the rate of light absorption at any wavelength
X, R (photons/Jl-s-nm), is directly proportional to the concentration of the
absorbing compound, c(mole/£).
R = K ,c (39)
a aA
where
K , = specific adsorption rate at wavelength X (photons/mole-nm-s)
The specific absorption rate at any given wavelength is dependent on the rate
at which light energy of that wavelength is incident on the compound [the
2
sealer irradiance, I (X) (photons/cm -s-nm))], and the extent: to which the
compound absorbs light of that wavelength [the molar absorptivity, e,(!l/mole-
cm) ] :
K . = 1000 e.I (X) (40)
aA A o
In natural waters, UV-light of wavelength greater or equal to 295 nm and
visible light up to 700 to 800 nm may contribute to photolysis. However, for
the majority of compounds light absorption is most significant in the UV
range (295-400 nm). The total specific adsorption rate, K , is determined by
3
integrating K , over the wavelength spectrum of the incident light:
K = 1000 / e,I (X)dX (41)
a X o
20
-------
The total rate of light absorption is then the product of K and the
3
molar concentration of the compound. If all the light absorbed were con-
verted to photoreaction, the rate of decay due to photolysis would simply
equal the rate of absorption:
• -V
In fact, the absorption-reaction step has a low efficiency and equation (42)
must be modified by the fraction of absorbed light that causes chemical
change; the quantum yield, <|>:
ft • -*v (43)
The quantum yield is an inherent property of the compound that may be
readily computed from kinetic data obtained in the laboratory. It may be
lowered by interactions of the compound with molecules called quenchers that
remove the excess energy from the compound returning it from its electron-
ically excited state to its ground state. Oxygen is the predominant quencher
in natural water systems. The quantum yield may be increased by chemical
reaction of the electronically excited compound with another molecule. It
may also be altered by adsorption of the compound on suspended solids. In
general, because of the short excited state lifetime (usually < 10 ns) and
low concentrations of the pollutants of concern, natural photolysis and photo-
lysis in pure air-saturated water are assumed equivalent (5).
Sensitized Photolysis
In addition to absorbing energy directly from photons, photoreactive
compounds may take up energy from certain compounds called photosensitizers
that have absorbed light. Photosensitizers are compounds capable of achieving
a relatively long-lived electronically excited state called the triplet state.
Humic materials are the most common sensitizers in natural waters.
The quantitative expression of sensitized photolysis includes the rate
of light absorption by the sensitizer and quantum yields for the three steps
involved in the photodegradation: (1) formation of the triplet state photo-
sensitizer ( ), (2) transfer of energy to the photoreactive compound ( )
s t
21
-------
and, (3) degradation of the electronically excited compound (<)>). The re-
sulting rate equation is:
where
s = concentration of sensitizer
K = total specific absorption rate of the sensitizer
3
s
The quantum yield defining the fraction of triplet energy transferred to the
photoreactive compound, A, is directly proportional to the concentration of
the photoreactive compound:
4»t = Ktc (45)
If the quantum yields and the concentration of sensitizer are assumed
constant, Eq. (44) reduces to a form identical to the rate expression for
direct photolysis (Eq. 43):
£=-K'c (46)
where
K' = s
-------
Ultraviolet light intensity at the water surface is computed from re-
cently published analytical expressions relating light attenuation to solar
altitude, ozone and aerosol light absorption, and aerosol and air molecule
light scattering (6,7). Following Zepp and Cline (8) surface intensity of
visible light is computed using equations described by Leighton (9).
Light attenuation in water is a first-order process described by an equa-
tion of the form:
where
-Ke z
I (X,z) = I (A,o)e (47)
Ke = attenuation coefficient for light of wavelength X (m )
A
z = depth (m)
Calculation Procedure
The model currently considers only direct photolysis. The calculation
procedure follows Zepp and Cline (8). From a user specified location and
season, atmospheric light intensities are computed. Using user specific
light attenuation coefficients, the average light intensity in each water
column segment of the model is then computed. Molar absorptivities of the
chemical under study are combined with the segment light intensities to
compute the rate of light absorption by the compound in accordance with
equation (35). Given the quantum yield, photolysis rate constants are then
computed. The average daily photolysis rate is determined by averaging com-
puted rates for solar altitude increments up to the maximum solar altitude
for the specified location and season.
3.2.3 Hydrolysis
Hydrolysis is a reaction in which a cleavage of a molecular bond of
the compound and the formation of a new bond with the hydrogen and hydroxyl
components of the water molecule results. Hydrolytic reactions are usually
mediated by an acid or base and are identified in that fashion. To a more
limited degree, there may also be a neutral reaction with water. By its
nature, the rate of the reaction is a function of pH and, as with most chem-
ical reactions, temperature.
23
-------
It is a common phenomenon in many chemical and biochemical processes.
In the biochemical processes, the fundamental food substances, such as car-
bohydrates, lipids and proteins are decomposed in this manner by various
enzymes. These reactions are discussed under the Biodegradation Section
which follows. In chemical processes, certain salts undergo hydrolysis.
Salts formed in neutralization of a strong acid by a weak base or of a strong
base by a weak acid yield acidic and alkaline reactions respectively, because
of the higher concentrations of hydrogen or hydroxyl ion. The hydrolysis of
organic chemicals, in some cases, is brought about by both biochemical as
well as chemical factors and it is difficult to distinguish between them in
prototype and certain laboratory conditions. It is pertinent to discuss each
separately and to examine the factors which affect the individual processes.
The reaction equation, fundamentally second order, containing the
product of the concentration of the chemical and either the acid, or alkaline
concentration, in general, may be expressed as:
|| - -Kc (48)
in which
K = K [H]+ for acid hydrolysis
3
for alkaline hydrolysis
= K [H?0] for neutral hydrolysis
The effect of pH on the reaction may be quite pronounced. Due to the
acid and base catalytic action, greater reaction coefficients are observed at
the acid and alkaline pH extremes and minimum values in the neutral zone.
The effect of temperature is typical of chemical and biochemical re-
actions with the rate increasing with a rise in the temperature. The effect
is usually expressed in terms of the Arrhenius equation, which relates the
reaction coefficient to the activation energy (calories/mole), the absolute
temperature and the gas constant. Based on this relation, the customary
procedure of expressing the rate coefficient is as follows:
24
-------
T-?n
K = K 61 ZU (49)
Reported laboratory and field (in soil experiments) indicate a range of 6 of
about 1.05-1.15 with an upper value of 1.2, i.e., temperature differences of
between 5 and 15°C causing a doubling of the rate.
In addition to the effect of temperature and pH, there appear to be
other factors which may influence the reaction rate, as is borne out by com-
paring rates of degradation in soil and in water. In general, the experimen-
tal data indicate the rate is greater in soil in some cases, by an order or
more. This medium apparently catalyzes the degradation, perhaps due to the
combined activity of both biological and chemical hydolysis.
The model considers acid, alkaline, and neutral hydrolysis. The user is
required to input the rate constant for each and the pH in each segment of the
model.
3.2.4 Biodegradation
Biodegradation encompasses the broad and complex processes of enzymatic
attack by organisms on organic chemicals. Bacteria, and to a lesser extent
fungi, are the mediators of biological degradation in natural systems. Dehalo-
genation, dealkylation, hydrolysis, oxidation, reduction, ring cleavage, and
condensation reactions are all known to occur either metabolically or via co-
metabolism. Co-metabolism refers to degradation of pesticides by microorgan-
isms when the microbe is not capable of utilizing the pesticide as a substrate
for growth. Sometimes the microbe does not even derive carbon or energy from
the degradation, rather the pesticide is "caught up" in the overall metabolic
reactions (10).
A chemical that is metabolized is used a a source of carbon, some other
nutrient element or energy. When such compounds are introduced into a natural
ecosystem the heterotrophs able to grow on them proliferate, the active popu-
lation increasing in size with time.
In cometabolism the population density of the responsible species does
not increase because it has no selective advantage in the presence of the
compound and cannot use it as a nutrient. Because an increase in population
25
-------
density of species growing on a chemical is generally paralleled by an in-
crease in degradation rate, cometabolism is characterized by the lack of an
increase in disappearance rate with time after introduction of the chemical.
Metabolic reactions follow Monod or Michaelis-Mention enzyme kinetics.
The degradation rate and bacterial growth rate can be expressed as:
, V c B
dc _ max
IF - - 1+ c
.
dt dt
c = pesticide concentration
t = time
V = maximum rate of substrate utilization
max
B = bacteria concentration
k^ = Michaelis half saturation constant
a = yield coefficient for bacteria utilizing pesticide.
Values of k typically range from 0.1 to 10 mg/£.
In general, chemical concentration is much less than the half-saturation
constant and equation (50) reduces to:
. V
dc _ max _ /,-9>.
_ --- _ cB _ . ^cB (52)
This equation is used in the model to describe biodegradation and there-
fore the user is required to specify the rate constant, K, , and J:or each seg-
ment, the bacterial concentration.
Specification of the proper biomass is critical since the degradation
rate is directly proportional to it. Obviously, it is the concentration of
bacteria capable of degrading the chemical under study that must be specified.
However, for chemicals whose degradative pathway is a nonspecific enzymatic
reaction such as hydrolysis use of total biomass may be sufficient.
26
-------
3.2.5 Adsorption
As discussed in the development of the equations describing transport,
sorption kinetics are not considered in the modeling framework. Instead,
local equilibrium is assumed between the dissolved and particulate phases.
This assumption implies that in comparison to any of the other processes
occurring in the system, the sorption-desorption equilibrium is achieved so
quickly that in describing the other processes the sorption-equilibrated con-
centrations can be found.
A linear isotherm relationship (equation (3)) is used to describe the
equilibrium sorption. This type relationship has been observed in numerous
adsorption studies (e.g., 11, 12, 13, 14) and is the most frequently used
isotherm for describing organic chemical adsorption in soils and sediments
(15, 16).
The partition coefficient describing the linear isotherm at a fixed
solids concentration has been shown for numerous chemicals to depend on the
concentration of adsorbing solids in the following fashion (17, 18):
If = ~ + H. (53)
m
where
m = concentration of solids
a,8 = empirical constants
11^ = limiting partition coefficient at high solids concentration
For most chemicals H was found to be zero and is not considered in the
GO
model.
The model uses this solids-dependent partitioning to compute the par-
ticulate and dissolved components of the chemical under study. A constant
partitioning condition may be simulated by setting b equal to the partition
coefficient and a to zero.
27
-------
SECTION 4
STRUCTURE OF COMPUTER CODE
4.1 OVERVIEW
WASTOX was developed as a general purpose computer program for modeling
the fate of toxic chemicals in any water body. It is based on the Water Ana-
lysis Simulation Program (WASP)(20) which was developed for modeling conven-
tional pollutants in water bodies. The major differences between WASTOX and
WASP are in the way transport is handled and in the reaction kinetics. In
WASP flow and dispersion are assumed to be the same for all constituents being
modeled. To account for the differences in flow and dispersion between dis-
solved chemical and adsorbed chemical and the multiple transport processes
affecting each, WASTOX considers groups of flows and dispersions that are
applied to either dissolved chemical or adsorbed chemical and solids, or both.
In addition, the complete freedom allowed the user of WASP in specifying the
reaction kinetics does not exist in WASTOX. Because WASTOX is directed to
toxic chemicals only, the kinetics and associated input have been incorp-
orated. The user may add to the kinetic processes but may not completely
respecify them.
WASTOX consists of a mainline and twenty-eight subroutines. The mainline
controls the assignment of logical units to the disk files used to store cer-
tain inputted information and the computational results and calls the subrou-
tines. The subroutines may be separated by function into the categories
input, computation, output, and utility. A list of the subroutines by cate-
gory is given in Table 1.
28
-------
TABLE 1. Subroutines comprising WASTOX
Mainline
RTWASP
Input
WASP 2
WASPS
WAS 3 A
WAS PA
WASPS
WASP6
WASP 7
PHOTO
WASP9
WAS 10
WAS 11
Computation
GREEN
SPECIR
ISLAM
WAS 12
WAS12B
WAS 15
WASPS
WASPB
WTXSS1
WTXSS2
WTXSS3
WTXSS4
Output Utility
WAS 16 ABORT
WSTSS3 SCALF
BWRITE
BREAD
FILEOC
FILEOP
RESET
4.2 SUBROUTINES
WASP 2
WASP2 initializes several variables and reads model identification
information, the number of segments into which the water body is divided, the
number of constituents being modeled, and system by-pass options. (Card
Group A)
WASP 3
WASP3 reads the input related to dispersion between segments. Depending
on which of the six possible options is chosen, dispersions (E) and
associated areas (A) and characteristic lengths (L) or exchanges (EA/L) are
read in as constant values or variable in time as broken line approximations.
(Card Group B)
29
-------
WAS 3 A
WAS3A reads the volume of each segment either as a constant or as a
timevariable function represented by a broken line approximation. (Card
Group C)
WASP 4
WASP4 reads the input related to flow between segments. Depending on
which of three possible options is chosen, flows are read in as constants or
as time-variable functions represented by broken line approximations. Flows
are converted from their input units of cubic feet per second to million gal-
lons per day. (Card Group D)
WASPS
WASPS reads the concentrations of each constituent modeled at the boun-
daries of the model where a flow or dispersion (exchange) occurs. Depending
on which of three possible options is chosen, boundary concentrations are
read in as constants or as time-variable functions represented by broken line
approximations. (Card Group E)
WASP6
WASP6 reads the mass rate of discharge into the water body of the
constituents modeled (forcing functions). Depending on which of three
possible options is chosen, forcing functions are read in as constants or as
timevariable functions represented by broken line approximations. (Card
Group F)
WASP7
WASP7 reads the parameters for each segment, the kinetic constants, and
any miscellaneous time functions. The number of parameters and constants
read in is computed internally for the current kinetic structure of
subroutine WASPB. To add parameters and constants for any modifications of
the kinetics requires changing the calculation of NPAM and NCONS. Note that
for NPAM > 11 or NCONS > 30 the size of arrays PARAM and CONST in
COMMON/MAIN/ must be changed. (Card Groups G,H and I)
30
-------
PHOTO
PHOTO reads the coefficients necessary to compute the first-order pho-
tolysis rate constant and calls the subroutines GREEN, SPECIR, and TSLAM that
compute the rate constant. The photolysis input and rate constant calculation
may be bypassed using a switch. (Card Group J)
WASP9
WAS10 reads the print interval. (Card Group M)
WASH
WASH reads the integration time step and the total time of the run.
(Card Group N)
GREEN
GREEN computes the time of day for a given solar zenith angle, latitude,
longitude, and ephemeride data. This information is used in subroutine SPECIR
and TSLAM.
SPECIR
SPECIR computes the direct and diffuse UV spectral irradience at the
water surface for use in calculating the photolysis rate constants.
TSLAM
TSLAM computes the direct and diffuse visible irradiance at the water
surface for use in calculating the photolysis rate constants.
WAS12-WAS12B
WAS12 and WAS12B compute the transport component of the time derivative
and add in any discharges or forcing functions.
WAS 15
WAS15 performs the numerical integration using an Euler scheme. It also
is the time counter for the model.
31
-------
RESET
RESET recomputes the break times for time variable functions to repeat
the function for times greater than the time of the last break in the func-
tion.
WTXSS1
steady-state solution subprogram which writes flow and bulk dispersion
coefficients onto a disk file.
WTXSS2
accesses WASPB to determine reaction coefficients and then combines these
with the transport coefficients to form steady-state solution subprogram which
calculates each element of the coefficient matrix .
WTXSS3
sequentially computes the concentration vector for each system, and
prints out steady-state concentration for each system.
WTXSS4
solves the set of simultaneous equations defined by the coefficient
matrix and the concentration vector.
WASPS
WASP8 evaluates the slopes of the broken-line approximations used to
represent time variable flows, dispersions, boundary conditions, and miscel-
laneous kinetic functions at each change in slope.
WASPB
WASPB computes the transfer and reaction components of the derivative.
Rates of photolysis, hydrolysis, biodegradation and volatilization are
computed. The fraction of chemical dissolved and adsorbed is computed from
the sorption partition coefficient and the porosity and solids concentration
of each segment. Segment concentrations are written to disk files at
specified intervals.
32
-------
WAS 16
WAS 16 reads the names, variable numbers and segment numbers for the
variables to be printed. The times and concentrations are read from the disk
files generated in WASPB during the simulation and printed.
ABORT
ABORT prints the reason for the termination of the run due to a user
generated error.
SCALF
SCALF modifies the values in an array by a scale factor. The scaling
operation may be a multiplication, division, addition or subtraction. ,
BWRITE
BWRITE writes information for VAL and T onto a specified disk file
beginning with a specified record.
BREAD
BREAD reads information into arrays VAL and T from a specified disk file
beginning with a specified record.
DEC PDF Subroutines
Two special subroutines were written for the DEC computer system. These
subroutines FILEOP and FILEOC were necessitated due to the way the DEC opera-
ting system handles disk output, i.e., requiring separate core buffers for
each disk file. FILEOP and FILEOC permit the disk files to share a common
disk buffer, reducing the excessive core requirements required for separate
buffers, at little cost to execution time.
33
-------
SECTION 5
PREPARATION OF DATA INPUT
5.1 INTRODUCTION
The computer code requires that all advective transport be specified as
flows. The flows resulting from the settling and resuspension velocities in
equation (22) are the product of the velocity and the interfacial area of the
segments associated with the velocity. Flows are inputted in groups (fields)
that are applied to dissolved chemical, or particulate chemical and solids,
or both. Most models will have separate flow fields for the hydrodynamic
flows, settling flows, resuspension flows, and bed sedimentation flows.
Dispersion coefficients are also inputted in groups or fields that are
applied to dissolved chemical, or particulate chemical and solids, or both.
Most models will have separate dispersion fields for water column dispersion
and interstitial water diffusion.
The state variables in the model are termed systems. System 1 is total
chemical and systems 2 through N are solids types. The maximum number of
systems is 4 (i.e., a maximum of 3 solids types).
To properly account for certain kinetic processes the model must know
which segments are water column segments. This is accomplished by requiring
that in the segment numbering scheme chosen by the user the water column seg-
ments be numbered from 1 to the number of water column segments and the bed
segments be numbered from one plus the number of water column segments to the
total number of segments.
34
-------
5.2 SUMMARY OF CARD GROUPS
Card Group
A. Model Identification and System Execution Options
1. Program Option
2. Model Identification Numbers
3. Title Card
4. Solution Type
5. System Execution Options
B. Exchange Coefficients
1. Input Option Number: Number of Exchange Fields
2. Number of Exchanges; Number of Systems Affected; Title
3. System Scale Factors
4. Exchange Coefficients
Card groups 2-4 are read for each exchange field
C. Segment Volumes
1. Number of Volumes; Input Option Number
2. Volumes
D. Flows
1. Input Option Number; Number of Flow Fields
2. Number of Flows: Number of Systems Affected; Title
3. System Scale Factors
4. Flows
Card groups 2-4 are read for each flow field
E. Boundary Conditions
1. Number of Boundary Conditions; Input Option Number
2. Boundary Conditions
Cards 1-2 are input for each system in the model
F. Forcing Functions
1. Number of Forcing Functions; Input Option Number
2. Forcing Functions
Cards 1-2 are inputted for each system in the model
35
-------
SUMMARY OF CARD GROUPS (cont'd)
G. Parameters
1. Parameters
Card group is inputted for each segment in the model
H. Constants
1. Constants
I. Miscellaneous Time Functions
1. Number of Time Functions
2. Function Name; Number of Breaks in Function
Cards 2-3 are inputted for each time function in the model
J. Photolysis
1. Execution Switch
2. Number of Wavelengths
3. Toxicant Molar Extinction Coefficient
4. Toxicant Quantum Yield
5. Latitude; Longitude; Season
6. Extinction Coefficients
7. Cloud Cover Effect
K. Initial Conditions
1. Initial Conditions for each system in the model
L. Stability Criteria
1. Stability Criteria
M. Print Control
1. Print Interval
N. Integration Information
1. Integration Interval; Total Time
0. Display Parameters
1. Variable Names
2. Dump parameters
Cards 1 and 2 are inputted for each system in the model
36
-------
5.3 CARD GROUP A - MODEL IDENTIFICATION AND SYSTEM EXECUTION OPTIONS
The variables which appear on each card are as follows:
5.3.1 Exposure Concentration - Food Chain Option
PRGOPT
FORMAT(15)
PRGOPT = 1; execute exposure concentration component of WASTOX only
PRGOPT = 2; execute food chain component of WASTOX only
PRGOPT = 3; execute exposure concentration and food chain components
of WASTOX
5.3.2 Model Identification Numbers
5 10 15 20 25
MODEL ISER IRUN NOSEG NOSYS
FORMAT(515)
MODEL - model designation
ISER - series designation
IRUN - run number
NOSEG - number of model segments (maximum of 75)
NOSYS - number of systems (maximum of 4. Note that system 1 is toxi-
cant and system 2 through N are solids types. Minimum number
of systems is 2).
MODEL, ISER, IRUN, although not actually used by the program, can
assist the user in maintaining a log of computer simulations.
5.3.3 Title
80
REAL TIME WASP TOXICS MODEL
FORMAT(2OA4)
Card columns 1-80 contain any information the user feels would be
helpful in describing the run and identifying the output for later reference,
37
-------
5.3.4 Time Variable or Steady State Solution
Type
FORMAT(A4)
Data input options:
TYPE - TIME VARIABLE; time variable solution requested
TYPE - STEADY STATE; steady state solution requested
Card columns 1 to 4 only are used to establish run type (e.g. TIME or
STEA)
5.3.5 System Execution Options
SYSEXC(l) SYSEXC(2) ... SYSEXC(NOSYS)
FORMAT(411)
SYSEXC(I) - 0, perform the kinetic and transport phenomena associated
with system I (numerically integrate or iterate the dif-
ferential equations).
1, bypass all kinetic and transport phenomena associated
with system I (concentrations read as initial conditions
for system I apply throughout simulation period).
38
-------
5.4 CARD GROUP B - EXCHANGE COEFFICIENTS
Exchange coefficients may be inputted in one of two forms, actual
exchange coefficients or, they may be calculated from inputted dispersion
coefficients and accompanying crosssectional areas and characteristic lengths.
Exchange coefficients for bed segments are assumed to be the product of poro-
sity and the exchange coefficient.
5.4.1 Data Input Option Number; Number of Exchange Fields
10
IROPT NRFLD
FORMAT(2I5)
Data input options:
IROPT - 1, constant exchange coefficients.
- 2, all exchange coefficients proportional to one piecewise linear
approximation.
- 3, each exchange coefficient represented by its own piecewise
linear approximation.
- 4, constant exchange coefficients calculated from the dispersion
coefficient, cross-sectional area, and characteristic lengths
specified for each interface.
- 5, all exchange coefficients proportional to one piece linear
approximation, calculated from a piecewise linear dispersion
coefficient approximation and respective cross-sectional areas,
and characteristic lengths.
- 6, each exchange coefficient proportional to its own piecewise
linear approximation, calculated from a piecewise linear approx-
imation for the dispersion coefficients, cross-sectional area,
and characteristic length specified for each area.
NRFLD - number of exchange fields.
If no exchange coefficients are to be read, set IROPT equal to zero, and
continue with Card Group C.
39
-------
Card types 2 through 4 are repeated NRFLD times, once for each
separate exchange field incorporated in the model.
5.4.2 Number of exchanges, Number of systems affected, Exchange field title
5 10 21 80
NOR NOSE RTITL
FORMAT(2I5,10X,15A4)
NOR = number of exchange coefficients in this field
NOSE = number of systems that this exchange field is to be applied to
RTITL = any alphanumeric descriptor to describe this exchange field
examples:
WATER COLUMN EXCHANGES
INTERSTITIAL WATER DIFFUSION
5.4.3 System Scale Factors
10 15 25 75_
SCALES (1) RSYS(l) SCALES (2) RSYS(NOSE)
FORMAT(4(E10.3,15))
RSYS(I) = system affected by this exchange field
SCALES(I) = scale factor for system RSYS(I). All exchanges in this
field are multiplied by SCALES(I) when the field is
applied to system RSYS(I).
5.4.4 Exchange Coefficients
The data input format is determined by the option selected.
40
-------
Option 1
Each card in this package contains the exchange coefficient information
for four interfaces. The number of exchange coefficients read is equal to
NOR. The last card is a scale factor for the exchange coefficients. The
information on each card is described below:
10
BR(I)
50
BR(l+2)
15
IR(D
55
IP(I+2)
20
JRCD
60
JR(I+2)
30
BR(I+1) I
70
BR(I+3) I
35
R(I+D
75
R(I+3)
40
JR(I+1)
80
JR(I+3)
FORMAT(4(F10.0,215))
BR(I) ~ exchange coefficient between segments IR(I) and JR(I)
in million cubic feet per day.
IR(I),JR(I) - segments between which exchange takes place the
order of the segments is not important; if a segment
exchanges with a boundary, the boundary is specified
as zero.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value
UNITS - the units of the exchange coefficient after FACTOR
has been applied.
NOTE: To request no scale factor put the value 1. after column 5
but before column 16 on the scale factor card. Also leave
columns 1 to 5 blank.
41
-------
Option 2
The card package consists of two sub-packages. Subpackage I contains the
exchange coefficient data, while sub-package II contains a detailed specifica-
tion of the piecewise linear approximation to which all exchange coefficients
are proportional.
Sub-Package I - Exchange Coefficients
Each card in this sub-package contains the exchange coefficient informa-
tion for four interfaces. The number of exchange coefficients read is equal
to NOR. The last card of this subpackage contains the scale factor for the
exchange coefficients. The information on each card is described below:
10
BR(I)
50
BR(I+2) I
15
IR(I)
55
R(I+2)
20
JR(D
60
JR(I+2)
30
BR(I+1) IR
70
BR(I+3) I
35
(I+D
, 75
R(I+3)
40
JR(I+1)
80
JR(I+3)
FORMAT(4(F10.0,215))
BR(I) - ratio of the exchange coefficient between segments IR(I) and
JR(I) to the piecewise linear approximation.
IR(I),JR(I) - segments between which exchange takes place
NOTE: the order of the segments is not important; if a segment
exchanges with a boundary, the boundary is specified as zero.
15
35
OP
FACTOR
UNITS
FORMAT(A2.3X.F10.2,2A10)
OP
FACTOR
UNITS
the operation (* / + or **) of the scale factor.
the scaling value.
the units of the exchange coefficient after FACTOR has been
applied.
42
-------
Sub-Package II - Piecewise Linear Approximation
The number of breaks required to describe the broken line approximation
is followed by a time series describing the broken line approximation. Each
time series element consists of two parts; an exchange value and a time (nor-
mal time scale is days). The last card of this subpackage contains the
scale factor for exchange value. The input is as follows:
NOBRK
FORMAT(15)
NOBRK - number of values and times used to describe the piecewise linear
approximation.
10
20
30
40
50
60
70
80
RT(I)
RT(I+2)
(T+2) RT(I+3) T(I+3)
FORMAT(8F10.0)
RT(I)
•JP
FACTOR
UNITS
value of the approximation at time T(I), in million cubic feet
per day.
time in days; if the length of the simulation exceeds T(NOBRK),
the piecewise linear approximation will repeat itself, starting at
time T(l); i.e., the approximation is assumed to be periodic with
period equal to T(NOBRK), this holds true for all piecewise linear
functions of time.
15
35
OP
FACTOR
UNITS
FORMAT(A2.3X.F10.2,2A10)
the operation (* / + - or **) of the scale factor
the scaling value
the units of the exchange value after FACTOR has been applied.
43
-------
Option 3
Each exchange coefficient is defined by a package of cards consisting of
two sub-packages. The first sub-package identifies the two segments between
which the exchange will take place, and the number of values comprising the
piecewise linear approximation. The second subpackage defines the piecewise
linear approximation which describes the exchange coefficient and a scale
factor. The input is as follows:
SubPackage 1
5 10 15
JR(I)
FORMAT(213,14)
IR(I),JR(I) - segments between which exchange takes place;
NOTE: for exchange only, order of segments is not important.
If segment exchanges with a boundary, the boundary is speci-
fied as zero.
NOBRK - number of values and times used to describe the piecewise
linear approximation.
All exchanges must have the same number of breaks.
SubPackage II - Piecewise Linear Approximation
This consists of a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; an exchange value,
and a time (normal time units days). The last card is a scale factor for
the exchange value. The input is as follows:
RT(I) T(I) RT(I+1) T(I+1) RT(I+2) T(I+2) RT(H-3)
FORMAT(8F10.0)
44
-------
RT(I) - value of the piecewise linear approximation at time T(I) in million
cubic feet per day.
T(I) - time in days. All break times must agree for all segments, i.e.,
T(l) must be the same for all exchanges, T(2) must be the same for
all exchanges, etc.
15
35
40
OP
FACTOR
UNITS
IFLG
FORMAT(A2,3X,F10.2,2A10,I5)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the exchange value after FACTOR has been applied.
IFLG - flag set equal to 1, if this scale factor is to apply to all the
remaining exchange coefficients in this field.
If IFLG = 1, then scale factor cards are not included for the
remaining exchange coefficients in this field.
45
-------
Option 4
Each card in this package contains the information to calculate the ex-
change coefficient for two interfaces. The number of dispersion coefficients
is equal to NOR. The last card is a scale factor for the dispersion coeffi-
cient. The information on each card is described below:
10
E(I)
50
20
A(I)
60
25
IL(I)
65
30
JL(I)
70
35
IR(I)
75
40
JR(D
80
E(I+1) A(I+1) IL(I+1) JL(I+1) IR(I+1) JR(I+1)
FORMAT(2(2F10.0.2F5.0,215))
E(I) - dispersion coefficient for the interface between segment IR(I) and
JR(I) in square miles per day.
A(I) - the interfacial crosssectional area between segments IR(I) and
JR(I), in square feet.
IL(I) - the length of segment IR(I), with respect to the IR(I)-JR(I) inter-
face, in feet.
JL(I) - the length of segment JR(I) in relation to the IR(I)-JR(I) inter-
face, in feet. If a segment exchanges with a boundarj', the charac-
teristic length of the boundary should be set equal to the length
of the segment with which it is exchanging.
IR(I),JR(I) - segments between which exchange takes place;
NOTE: for exchange only, order is not important—if a segment ex-
change with a boundary, the boundary is specified as 2:ero.
2 15 35
OP FACTOR UNITS
FORMAT(A2.3X.F10.2,2A10)
46
-------
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the dispersion coefficient after FACTOR has been
applied.
47
-------
Option 5
The card package consists of two sub-packages. Sub-Package I contains
the information necessary to calculate the exchange coefficients, while sub-
package II contains a detailed specification of the piecewise linear approxi-
mation to which the dispersion coefficients contained in sub-package I are
proportional.
Sub-Package I
Each card in this sub-package contains the information necessary to cal-
culate the exchange coefficients for two interfaces. The number of dispersion
coefficients is equal to NOR. The last card is a scale factor for the disper-
sion coefficient. The information on each card is described below:
10 20
E(I) A(I)
50 60
25
IL(I)
65
30
JL(I)
70
35 40
IR(I) JR(I)
75 80
E(I+1) A(I+1) IL(I+1) JL(I+1) IR(I+1) JR(I+1)
FORMAT(2(2F10.0.2F5.0,215))
E(I) - the ratio of the dispersion coefficient between segment IR(1) and
JR(I) to the piecewise linear approximation.
A(I) - the interfacial crosssectional area between segments; IR(I) and
JR(I), in square feet.
IL(I) - the length of segment IR(I) in relation to the IR(I)-JR(I) interface,
in feet.
JL(I) - the length of segment JR(K) in relation to the IR(I)-JR(I) interface,
in feet. If a segment exchanges with a boundary, the characteristic
length of the boundary should be set equal to the length of the seg-
ment with which it is exchanging.
IR(I),JR(I) - segments between which exchange takes place;
NOTE: for exchange only, order is not important.
48
-------
_2 15 35^
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the dispersion coefficient after FACTOR has been
applied.
Sub-Package II - Piecewise Linear Approximation
The number of breaks required to describe the piecewise linear approxima-
tion is followed by a time series describing the piecewise linear approximation.
Each time series element consists of two partsj a dispersion coefficient and a
time (normal units are days). The last card is the scale factor for the disper-
sion coefficient. The input is as follows:
NOBRK
FORMAT (15)
NOBRK - number of values and times used to describe the piecewise linear
approximation .
10 20 30 40 50 60 70 80
_
RT(I) T(I) RT(I+1) T(I+1) RT(I+2) T(I+2) RT(I+3)
FORMAT (8F10.0)
RT(I) - value of the piecewise linear approximation at time T(I), in square
miles per day.
T(I) - time in days.
_2 _ 15 _ 35
OP _ FACTOR UNITS
FORMAT (A2,3X,F10. 2, 2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the dispersion coefficient after FACTOR has been
applied.
49
-------
Option 6
Each exchange coefficient is defined by a package of cards consisting of
three sub-packages. The first sub-package defines the interfacial cross-sec-
tional area and the characteristic length. The second sub-package identifies
the two segments between which the exchange will take place, and defines the
number of values comprising the piecewise linear approximation. The third
sub-package defines the piecewise linear approximation which describes the
dispersion coefficient.
Sub-Package I
The sub-package defines the interfacial cross-sectional area and the
characteristic length of the segments involved, also the scale factors.
10
20
A(I) A(I+1). . .A(NOR1) A(NOR)
FORMAT(8F10.0)
the interfacial cross-sectional area between segments in square
feet.
15
35
OP FACTOR UNITS
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the area after FACTOR has been applied.
10
20
XL(I)
XL(NOSEG-1)XL(NOSEG)
FORMAT(SFIO.O)
XL(I) - the length of the segment in the direction of the dispersion, in
feet. If a segment exchanges with a boundary, the characteristic
length of the boundary should be set equal to the length of the
segment with which it is exchanging.
50
-------
__2 _ 15 _ 35
OP FACTOR _ UNITS
FORMAT (A2,3X,F10. 2, 2A10)
OP - the operation (* / + - or **) of the scale factor)
FACTOR - the scaling value
UNITS - the units of the characteristic length after FACTOR has been
applied.
Sub-Package II
5 10 15
IR(I) JR(I) NOBRK
FORMAT(2I3,I4)
IR(I),JR(I) - segments between which exchange takes place;
NOTE: for exchange only, order is not important
NOBRK - number of values and times used to describe the piecewise
linear approximation. All exchanges must have the same number
of breaks.
Sub-Package III - Piecewise Linear Approximation
This consists of a time series describing the piecewise linear approxima-
tion. Each time series element consists of two parts; a dispersion coefficient,
and a time (consistent with the normal time scale of the model). The last card
is a scale factor for the dispersion coefficient. The input is as follows:
10 20 30 40 50 60 70 80
_ _ _ _ _ _
RT(I) T(I) RT(I+1) T(I+1) RT(I+2) T(I+2) RT(I+3)
FORMAT (8F10.0)
RT(I) - value of the piecewise linear approximation at time T(I), in square
miles per day.
T(I) - time in days; all break times must agree for all segments, i.e.,
T(l) must be the same for all exchanges, T(2) must be the same for
all exchanges, etc.
51
-------
OP
15
FACTOR
35_
UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the dispersion coefficient after FACTOR has been
applied.
IFLG - flag set equal to 1 if this scale factor is to apply to all the
remaining exchange coefficients in this field.
If IFLG = 1, then scale factor cards are not included for the
remaining exchange coefficients in this field.
52
-------
5.5 CARD GROUP C - VOLUMES
5.5.1 Data Input Option Number; Number of Volumes
10
IVOPT NOV
FORMAT(2I5)
Data input options:
IVOPT - 1, constant volumes
2, all volumes proportional to one piecewise linear approximation.
3, each volume is represented by its own piecewise linear approxima-
tion.
NOV - number of volumes; normally NOV is equal to NOSEG, the number of
segments;
5.5.2 Volumes
The data input format is determined by the option selected.
Option 1
Each card in this package contains the volume information for eight
segments. The number of volumes is equal to NOV. The last card is the scale
factor for the volumes. The information on each card is described below:
10 20 30 70 80
BVOL(I)BVOL(I+l)BVOL(I+2) .... BVOL(I+6)BVOL(I+7)
FORMAT(SFIO.O)
BVOL(I) - volumes of segment I, in million cubic feet. The volumes are
listed consecutively starting with segment 1, and ending with
segment NOV.
53
-------
OP
15
FACTOR
35
UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor
FACTOR - the scaling value
UNITS - the units of the volume after FACTOR has been applied.
Option 2
The card package consists of two sub-packages. Sub-package I contains
the volume data while sub-package II contains a detailed specification of the
piecewise linear approximation to which all volumes are proportional.
Sub-Package I - Volumes
Each card in this subpackage specifies the volume ratios of eight seg-
ments. The number of volumes is equal to NOV. The last card is the scale
factor for the volume ratios. The input is described below:
10
20
30
80
BVOL(I)BVOL(I+l)BVOL(I+2)
BVOL(I+7)
FORMAT(8F10.0)
BVOL(I) - ratio of the volume in segment I to the piecewise line approxima-
tion. The ratios are listed consecutively, starting with segment
1 and ending with segment NOV.
OP
FORMAT
UNITS
FORMAT(A2.3X.F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor,,
FACTOR - the scaling value.
UNITS - the units of the volume ratios after FACTOR has been applied.
54
-------
Option 3
This package reads a separate piecewise linear volume function for each
segment, from 1 through NOV. The last card is a scale factor for the volumes.
NOBRK
FORMAT(15)
NOBRK - number of values and times used to describe the piecewise linear
approximation. NOBRK must be the same for each broken line approxi-
mation.
10
20
30
40
50
60
70
80
VT(I)
VT(I+2)
VT(I+3)
FORMAT(8F10.0)
VT(1) - value of the piecewise linear approximation at time T(I), for
segment I in million cubic feet.
T(I) - time in days; if the length of the simulation exceeds T(NOBRK), the
approximation will repeat itself starting at time T(l), i.e., the
approximation is assumed to be periodic with period T(NOBRK). All
break times must agree for all segments, i.e., T(l) must be the same
for all volumes, T(2) must be the same for all volumes, etc.
15
35
40
OP
FACTOR
UNITS IFLG
FORMAT(A2,3X,F10.2,2A10,15)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value. UNITS the units of the volumes after FACTOR
has been applied.
IFLG - flag set equal to 1 if this scale factor is to apply to all the
remaining volumes. If IFLG = 1 then scale factor cards are not
included for the remaining volumes.
55
-------
5.6 CARD GROUP D - FLOWS
5.6.1 Data Input Option Number; Number of Flow Fields
5 10
IQOPT NOQF
FORMAT(2I5)
Data Input Options:
IQOPT - 1, constant flows.
2, all flows proportional to one piecewise linear approximation.
3, each flow is represented by its own piecewise linear approximation.
NOQF - number of flow fields.
If no flows are to be inputted, set IQOPT to zero and go to Card
Group E. The maximum number of flow fields is 9.
Card types 2 through 4 are repeated NOQF times, once for each
separate flow field incorporated in the model.
5.6.2 Number of flows, Number of systems affected, Flow direction switch,
Flow field title
5 10 15 21 80
NOQ NOSE IDIRSW QTITL
FORMAT(3I5,5X,15A4)
NOQ = number of flows in this field
NOSE = number of systems that this flow field is to be applied to
IDIRSW = switch that when set equal to 1 reverses the direction of all
the flows in this field
QTITL = any alphanumeric descriptor to describe this flow field
examples:
HYDRODYNAMIC FLOW
SETTLING OF SOLIDS
RESUSPENSION OF SOLIDS
56
-------
5.6.3 System Scale Factors
10 15 25 7_5
SCALES(I) QSYS(I) SCALES(2) . . . QSYS(NOSE)
FORMAT(5(E10.3,I5))
QSYS(I) = system affected by this flow field
SCALES(I) = scale factor for system QSYS(I). All flows in this field
are multiplied by SCALES(I) when the field is applied to
system QSYS(I).
5.6.4 Flows
The data input format is determined by the option selected.
Option 1
Each card in this package contains the flow information for four
interfaces, the number of flow specifications is equal to NOQ. The last card
is a scale factor for the flows. The information on each card is described
below:
10 15 20 30 35 40
BQ(I) IQ(I) JQ(I) BQ(I+1) IQ(I+1) JQ(I+1)
50
BQ(I+2)
55
IQCI+2)
60
JQ(I+2)
70
BQ(I+3)
75
IQCI+3)
80
JQ(I+3)
FORMAT(4(F10.0,215)
BQ(I) - flow between segment IQ(I) and JQ(I) in cfs. Convention: if the
flow value is positive, then flow is from segment JQ(I) to IQ(I).
IQ(1) - downstream segment
JQ(I) - upstream segment
If flow is from a segment to a boundary, then IQ(I) is set equal to zero;
if a flow is from a boundary to a segment, then JQ(I) is set equal to zero.
57
-------
_2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the flows after FACTOR has been applied.
Option 2
The card package consists of two sub-packages. Sub-package I contains
the flow routing while sub-package II contains a detailed specification of the
piecewise linear approximation to which all the flows are proportional.
Sub-Package I - Flows
Each card in this sub-package contains the flow information for four
interfaces. The number of flow specifications is equal to NOQ. The last card
is a scale factor for the flows. The information on each card is described
below:
10
BQ(I)
50
BQ(I+2)
15
IQ(D
55
IQ(I+2)
20
JQ(D
60
JQ(I+2)
30
BQ(I+1)
70
BQ(I+3)
35
IQ(I+1)
75
IQCI+3)
40
JQO+D
80
JQCI+3)
FORMAT(4(F10.0,215)
BQ(I) - ratio of the flow between segments IQ(I) and JQ(I) to the piecewise
linear flow approximation.
IQ(I) - downstream segment
JQ(I) - upstream segment
If flow is from a segment to a boundary, then IQ(I) is set equal to zero,
if a flow is from a boundary to a segment, then JQ(I) is set equal to zero.
58
-------
__2 15 35^
OP FACTOR UNITS
FORMAT(A2,3X,F102,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the flow ratios after FACTOR has been applied.
Sub-Package II - Piecewise Linear Flow
The number of breaks required to describe the piecewise linear approxima-
tion is followed by a time series describing the piecewise linear flow approx-
imation. Each time series element consists of two parts; a flow and a time.
The last card is a scale factor for the flows. The input is as follows:
N 0 B R K
FORMAT(I5)
NOBRK - number of values and times used to describe the piecewise linear
approximation .
10 20 30 40 50 60 70 80
_ _ _ _ _
QT(I) T(I) QT(I+1) T(I+1) QT(I+2) T(I+2) QT(I+3)
FORMAT(8F10.0)
QT(I) - value of the piecewise linear approximation at time T(I), in cubic
feet per second.
T(I) - time in days, if length of the simulation exceeds T(NOBRK), the broken
line function will repeat itself, starting at time T(I) , i.e., the
approximation is assumed to be periodic, with period equal to T (NOBRK) .
59
-------
_2 15 3_5
OP FACTOR UNITS
FORMAT(A2.3X.F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the flow after FACTOR has been applied.
Option 3
Each flow is defined by a package of cards consisting of two sub-packages.
The first subpackage identifies the two segments between which the flow
occurs, and the number of values comprising the piecewise linear flow approxi-
mation. The second subpackage defines the piecewise linear approximation
which describes the flow. The input is as follows:
SubPackage I
5 10 15
IQ(I) JQ(I) NOBRK
FORMAT(315)
IQ(I) - downstream segment, flow from segment JQ(I) to IQ(I), assuming
positive flow.
JQ(I) - upstream segment flow from segment IQ(I), assuming positive flow.
NOBRK - number of values and times used to describe the broken line
approximation. All flows must have the same number of breaks.
SubPackage II
Sub-package II is a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; a flow and a time.
The last card is a scale factor for the flows. The input is as follows:
60
-------
10 20 30 40 50 60 70 80
QT(I) T(I) QT(I+1) T(I-H) QT(I+2) T(I+2) QT(I+3) T(I+3)
FORMAT(8F10.0)
QT(I) - value of the piecewise linear flow approximation at time T(I) in
cfs.
T(I) - time in days, if the length of the simulation exceeds T(NOBRK)
the broken line function will repeat ifself, starting at time
T(l). All break times must agree for all segments; i.e., T(l)
must be the same for all flows, T(2) must be the same for flows,
etc.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the flow after FACTOR has been applied.
IFLG - flag set equal to 1 if this scale factor is to apply to all the
remaining flows in this field.
If IFLG = 1 then scale factor cards are not included for the
remaining flows in this field.
61
-------
5.7 CARD GROUP E - BOUNDARY CONDITIONS
All input is read NOSYS times; once for each system of the model.
5.7.1 Data Input Option Number; Number of Boundary Conditions!
5 10
IBCOP(I) NOBC(I)
FORMAT(2I5)
Data Input Options:
IBCOP(I) - 1, constant boundary conditions
2, all boundary conditions proportioned to one piecewise
linear approximation.
3, each boundary condition represented by its own piecewise
linear approximation.
NOBC(I) - number of boundary conditions used for system I.
If no boundary conditions are to be inputted, set NOBC(I) equal to zero,
and continue with the next system, or go to the next card group.
5.7.2 Boundary Conditions
The data input format is determined by the option selected.
Option 1
10 15 25 30 40
BBC(I) IBC(I) BBC(I+1) IBC(I+1) BBC(I+2)
45 55 60 70 75
IBC(I+2) BBC(I+3) IBC(I+3) BBC(I+4) IBC(I+4)
FORMAT(5(F10.0,15))
BBC(I) - boundary condition of segment IBC(I) in mg/£.
IBC(I) - segment number to which boundary condition BBC(I) is to be
applied.
62
-------
_2 15 35_
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the boundary condition after FACTOR has been
applied.
Option 2
This card package consists of two sub-packages. Sub-package I contains
boundary condition data, while sub-package II contains a detailed specifica-
tion of the piecewise linear approximation to which all boundary conditions
are to be proportional.
Sub-Package I
Each card in this sub-package contains the boundary condition information
for five segments. The number of boundary condition specifications is equal
to NOBC. The last card of the package is the scalefactor for the boundary con-
ditions. The information on each card is described below:
10 15 25 30 40
BBC(I) IBC(I) BBC(I+1) IBC(I+1) BBC(I+2)
45 55 60 70 75
IBC(I+2) BBC(I+3) IBC(I+3) BBC(I+4) IBC(l+4)
FORMAT(5(F10.0,15))
BBC(I) - ratio of the boundary condition for segment IBC(I) to the piece-
wise linear approximation.
TBC(I) - segment number.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
63
-------
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of boundary condition specifications after FACTOR has
been applied.
Sub-Package II - Piecewise Linear Boundary Condition Approximation
The number of breaks required to describe the piecewise linear boundary
condition approximation is followed by a time series describing the boundary
approximation. Each time series element consists of two parts; boundary con-
centration, and a time. The last card is a scale factor. The input is as
follows:
5
N 0 B R K
FORMAT(15)
NOBRK - number of values and times used to describe the piecewise linear
approximation.
10 20 30 40 50 60 70 8£
BCT(I) T(I) BCT(I+1) T(I+1) BCT(I+2) T(I+2) BCT(I+3) T(I+3)
FORMAT(8F10.0)
BCT(I) - value of the broken line approximation at time T(I) mg/£.
T(I) - time at breaks in broken line approximation in days.
If the length of the simulation exceeds T(NOBRK), the piecewise linear
approximation is repeated, starting at T(l), i.e., the approximation is assumed
to have period equal to T(NOBRK).
2 15 35
OP FACTOR UNITS
FORMAT(A2.3X.F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling values.
UNITS - the units of the boundary concentration after FACTOR has been
applied.
64
-------
Option 3
Each boundary condition is defined by a package of cards consisting of two
sub-packages. The first sub-package identifies the segment associated with the
boundary condition and the number of values comprising the piecewise linear
approximation. The second sub-package defines the piecewise linear approxima-
tion which describes the boundary conditions. All boundary condisions for one
system must have the same number of breaks. The input is as follows:
SubPackage I
5 10
IBC(I) NOBRK
FORMAT(2I5)
IBC(I) - boundary segment number.
NOBRK - number of values and times used to describe the broken line
approximation. The number of breaks must be equal for all boun-
dary conditions within a system.
Sub-Package H - Piecewise Linear Boundary Condition Approximation
The segment number and the number of breaks required to describe the
broken line approximation is followed by a time series describing the broken
line approximation. Each time series element consists of two parts: a boun-
dary concentration, and a time (consistent with the normal time scale of the
model). The number of breaks must be the same for all boundary approximations.
The last card is a scale factor for the boundary concentrations. The input is
as follows:
10 15 25 30 40 45 70 75
BCT(I) T(I) BCT(I+1) T(I+1) BCT(I+2) T(I+2) BCTd+4) T(I+4)
FORMAT 5(F10.0,F5.0)
BCT(I) - value of the boundary approximation at time T(I) in mg/£.
65
-------
T(I) - time in days; if the length of the simulation exceeds T(NOBRK),
the broken line approximation is repeated, starting at T(l),
i.e., the approximation is assumed to be periidic, with period
equal to T(NOBRK). All break times must agree for all segments,
i.e., T(l) must be the same for all exchanges, T(2) must be the
same for all exchanges, etc.
2 15 35 40
OP FACTOR UNITS IFLG
FORMAT(A2,3X,F10.2,2A10,I5)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the boundary concentration after FACTOR has been
applied.
IFLG - flag set equal to 1 if this scale factor is to apply to all the
remaining boundary conditions.
If IFLG = 1 then scale factor cards are not included for the
remaining boundary conditions.
66
-------
5.8 CARD GROUP F - FORCING FUNCTIONS
All input is read NOSYS times, once for each system of the model,
5.8.1 Data Input Option Number; Number of Forcing Functions
5 10
IWOP(I) NOWK(I)
FORMAT(2I5)
Data Input Options:
IWOP(I) - 1, constant forcing functions
- 2, all forcing functions are proportioned to one piecewise
linear approximation.
- 3, each forcing function represented by its own piecewise
linear approximation.
NOWK(I) - number of forcing functions used for system I. NOTE: forcing
functions may also be considered as sources (loads) or sinks
of a water quality constitutent. If no forcing functions are
to be inputted, set NOWK(I) to zero, and continue with next
system or go to next card group if this is the last system.
5.8.2 Forcing Functions
The data input format is determined by the option selected.
67
-------
Option 1
10 15 25 30 40
BWK(I) IWK(I) BWK(I+1) IWK(I+1) BWK(I+2)
45 55 60 70 75
IWK(I+2) BWK(I+3) IWK(I+3) BWK(I+4) IWK(I+4)
FORMAT(5(F10.0,I5))
BWK(I) - forcing function of segment IWK(I), in pounds per day.
IWK(I) - segment number to which forcing function BWK(I) is to be
applied.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the forcing function after FACTOR has been applied.
Option 2
The card package consists of two sub-packages. Sub-package I contains
the forcing function data, while sub-package II contains a detailed specifica-
tion of the piecewise linear approximation to which all forcing functions are
proportional.
Sub-Package I
Each card in this sub-package contains the forcing function information
for five segments. The number of specifications is equal to NOWK(I). The
last card is a scale factor for the forcing function information. The infor-
mation on each card is described below:
68
-------
10 15 25 30 40^
BWK(I) IWK(I) BWK(I+1) IWK(I+1)BWK(I+2)
45 55 60 70 75_
IWK(I+2) BWK(I+3) IWK(I+3) BWK(I+4) IWK(I+4)
FORMAT(5(F10.0,15))
BWK(I) - ratio of the forcing function for segment IWK(I) to the piecewise
linear approximation
IWK(I) - segment number to which forcing function BWK(I) is to be applied.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the forcing function ratio after FACTOR has been
applied.
Sub-Package II - Piecewise Linear Forcing Function Approximation
The number of breaks required to describe the piecewise linear forcing
function approximation is followed by a time series describing the forcing
function. Each time series element consists of two parts; a function value
and a time. The last card is the scale factor for the function value. The
input is as follows:
N 0 B R K
FORMAT(15)
NOBRK - number of values and times used to describe the piecewise linear
approximation.
69
-------
10 20 30 40 50 60 70 80
WKT(I) T(I) WKT(I+1) T(I+1) WKT(I+2) T(I+2) WKT(I+3) T(I+3)
FORMAT(8F10.0)
WKT(I) - value of the forcing function at time T(I), in pounds per day.
T(I) - time in days; if the length of the simulation exceeds T(NOBRK),
the forcing function approximation is repeated, starting at
T(l), i.e., the approximation is assumed to be periodic, with
period equal to T(NOBRK).
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor
FACTOR - the scaling value.
UNITS - the units of the forcing function after FACTOR has been applied.
Option 3
Each forcing function is defined by a package of cards consisting of two
sub-packages. The first sub-package identifies the segment associated with
the forcing function and the number of values comprising the piecewise linear
approximation. The second sub-package defines the approximation which des-
cribes the forcing function. The input is as follows:
SubPackage I
5 10
IWK(I) NOBRK
FORMAT(2I5)
IWK(I) - segment number which has forcing function BWK(I)
NOBRK - number of breaks used to describe the forcing function approxi-
mation. The number of breaks must be equal for all forcing
functions within a system.
70
-------
Sub-Package II - Piecewise Linear Forcing Function Approximation
The segment number and the number of breaks required to describe the
piecewise linear forcing function approximation is followed by a time series,
describing the forcing function. Each time series element consists of two
parts: a function value and a time. The last card is scale factor for the
forcing function. The input is as follows:
10
20
30
40
50
60
70
80
WKT(I) T(I) WKT(I+1) T(H-l) WKT(I+2) T(I+2) WKT(I+3) T(I+3)
FORMAT(5(F10.0.F5.0))
WKT(I) - value of the forcing function at time T(I), in pounds per day.
T(I) - time in days; if length of the simulation exceeds T(NOBRK), the
approximation is assumed to be periodic with period equal to
T(NOBRK). All break times must agree for all segments; i.e.,
T(l) must be the same for all forcing functions, T(2) must be
the same for all forcing functions, etc.
15
35
40
OP
FACTOR
UNITS IFLG
FORMAT(A2,3X,F10.2,2A10,I5)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value.
UNITS - the units of the forcing function after FACTOR has been applied,
IFLG - flag set equal to 1 if the scale factor is to apply to all the
remaining forcing functions.
If IFLG = 1 then scale factor cards are not included for the
remaining forcing functions.
71
-------
5.9 CARD GROUP G - PARAMETERS
All input is read NOSEG times, once for each segment.
The definition of the parameters is set in the kinetics subroutine
(WASPB). The number of parameters used in the model is dependent on the number
of solids types considered.
5.9.1 Segment Parameters
15
20
30
ANAME(I) PARAM(ISEG.I) ANAME(I+1) PARAM(ISEG,1+1)
35
45
ANAME(I+2) PARAM(ISEG,I+2)...
FORMAT(5(1X,A4,F10.0))
ANAME(I) - a one to four alphanumeric character descriptive name
for parameter PARAM(ISEG,I)
PARAM(ISEG,I) - the value of parameter ANAME(I) in segment ISEG.
Parameter
1
2
5
6
Fortran Variable
Name
DEPTH
TDEPTH
DWNSEG
AREA
PH
TEMP
WVEL
Required Parameters
Suggested Description Units
ANAME
DPTH Segment depth meters
TDPT Depth of segment sur- meters
face
DSEG Segment downstream of
current segment
2
AREA Area between current ft
segment and down-
stream segment
PH pH of segment
TEMP temperature of seg °C
ment
WVEL Wind velocity 10 m/s
meters above the
water surface at the
location of the cur-
rent segment
72
-------
Parameter Fortran Variable
Name
8 BACW
BAG SI
Required Parameters (cotit'd)
Description
Suggested
ANAME
BACW concentration of com-
pound degrading bac-
teria in water
BCS1 concentration of com-
pound degrading bac
teria on solids type
1
Units
num/A (or units
compatible with
degradation
rate constant)
nutn/mg (or units
compatible with
degradation rate
constant)
IF ADDITIONAL SOLIDS TYPES ARE INCLUDED IN THE MODEL:
for each additional solids type add a parameter analogous to parameter 9.
For example, if 3 solids types are considered (NOSYS=4) then a 10th and
llth parameter, BACS2 and BACS3, respectively, roust be included.
5.10
CARD GROUP H - CONSTANTS
The definition of the constants is set in the kinetics subroutine (WASPB),
The number of constants used in the model is dependent on the number of solids
types considered.
5.10.1 Constants
15
20
30
35
ANAME(I) CONST(I) ANAME(I+1) CONST(I+1) ANAME(I+2)
45
50
60
65
75
CONST(I+2) ANAME(I+3) CONST(I+3) ANAME(I+4) CONST(I+4)
FORMAT(5(1X,A4,F10.0))
ANAME(I) - a one to four alpha-numeric character descriptive name for
constant CONST(I).
CONST(I) - the value of constant ANAME(I).
73
-------
Required Constants
Constant
Fortran Variable
Name
KD20
Suggested
ANAME
KD20
8
9
10
11
12
KDT
KP20
KDT
KP20
KPT
KHOH20
KHN20
KHH20
KHT
XWSEG
HENRY
MOLWT
KLT
KPT
KHOH
KHN
KHH
KHT
WSEG
HNRY
MLWT
KLT
Description
2nd Order biodegradation
rate constant for dis-
solved toxicant at 20°C
Temperature correction
factor for biodegrada-
tion of dissolved toxi-
cant
2nd Order biodegrada-
tion rate constant for
sorbed toxicant at 20°C
Temperature correction
factor for biodegrada-
tion of sorbed toxicant
2nd Order alkaline
hydrolysis rate constant
at 20°C
1st Order neutral
hydrolysis rate constant
at 20°C
2nd Order acid hydroly-
sis rate constant at
20°C
Temperature correction
factor for hydrolysis
Number of water column
segments
Henry's Constant (set to
zero to skip volatiliza-
tion)
Molecular weight of tox-
icant
Temperature correction
factor for volatiliza-
tion
Units
&/d-num (or
units com-
patible with
bacterial
concentra-
tion BACW)
5,/dnum (or
units compat-
ible with bac-
terial concen-
tration BAGS)
£/mole-d
1/d
£/mole-d
Unitless
74
-------
Required Constants (cont'd)
Constant
Fortran Variable
Name
Suggested
ANAME
Description
Units
13
14
15
16
17
ATMOS
IFLOW
AIRTMP
Al
Bl
18
RH01
ATMS Concentration of chemical
in the air immediately
above the water
IFLW Switch to indicate a flow-
ing (0) or non-flowing (1)
natural water system
ATMP Air temperature
Al Solids dependent parti-
tioning exponent for
solids type 1
Bl Solids dependent parti- l/mg'(mg/£)
titioning coefficient for
for solids type 1
RH01 Density of solids type 1
Al
IF ADDITIONAL SOLIDS TYPES ARE INCLUDED IN THE MODEL:
for each additional solids type add 3 parameters analogous to parameters 16,
17 and 18.
For example, if 2 solids types are considered (NOSYS=3) then constants
19 (A2), 20 (B2), and 21 (RH02) must be included.
75
-------
5.11 CARD GROUP I - MISCELLANEOUS TIME FUNCTIONS
The definition of miscellaneous piecewise linear time functions will
vary depending upon the structure and the kinetics of the systems comprising
each model. No time functions are currently considered and NFUNC should be set
equal to zero. The user may wish to include time functions in a modified WASPB
(e.g., temperature) and the input format is detailed below:
5.11.1 Number of Time Functions
NFUNC
FORMAT(15)
NFUNC - number of time functions required for the model. If no time
functions are to be inputted set NFUNC equal to zero and go
to card group J.
5.11.2 Time Functions
The following package of cards is required for each time function.
The first subpackage defines the function name and number of breaks in the time
function. The second subpackage contains a detailed specification of the piece-
wise linear time function.
SubPackage I
10
ANAME(I) NOBRK
FORMAT(A4,16)
ANAME(I) - a one to four alpha-numeric character descriptive name for
time function I.
NOBRK - number of breaks used to describe the time function I.
76
-------
Sub-Package II
Each time series element consists of two parts; a function value and a
time. The last card is a scale factor for the function value.
10 20 30 40 50 60 70 8£
VALT(I) T(I) VALT(I+1) T(I+1) VALT(I+2) T(I+2) VALT(I+3) T(I+3)
FORMAT(8F10.0)
VALT(I) - value of the function at time T(I).
T(I) - time in days; if the length of the simulation exceeds T(NOBRK),
the time function will repeat itself, starting at T(l), i.e.,
the approximation is assumed to be periodic, with period equal
to T(NOBRK). All time functions must have the same number of
breaks and all break times must agree for all functions; i.e.,
T(l) must be the same for all functions, T(2) must be the same
for all functions, etc.
2 15 35
OP FACTOR UNITS
FORMAT(A2,3X,F10.2,2A10)
OP - the operation (* / + - or **) of the scale factor.
FACTOR - the scaling value
UNITS - the units of time function after FACTOR has been applied,
77
-------
5.12 CARD GROUP J - PHOTOLYSIS
1st order photolysis rate constants are computed for each water
column segment in the model from the input indicated below:
The calculation is adapted from Zepp and Cline (E,S & T, 11:359.
1977)
5.12.1. Execution switch
ISW
FORMAT(I5)
ISW =0 - Photolysis is not considered and input is bypassed
1 - Photolysis included in model
5.12.2 Wavelengths for which data will be inputted
NOWAV
FORMAT(15)
NOWAV = number of wavelengths considered.
Wavelengths are counted from a minimum wavelength of 297.5 nm in the
following sequence:
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, 430, 440, 450, 460, 470, 480, 490, 500,
525, 550, 575, 600, 625, 650, 675, 700, 750, 800
NOWAV is thus controlled by the maximum wavelength at which the toxicant
considered is photochemically active.
78
-------
5.12.3 Toxicant Molar Absorptivity
10 20
EPSILON(l) EPSILON(2) . . . EPSILON(NOWAV)
FORMAT(8F10)
EPSILON(I) = molar absorptivity of the toxicant at wavelength
number I £/mole-cm
5.12.A Toxicant quantum yield
10
OY
FORMAT(F10.)
QY = quantum yield of the toxicant
5.12.5 Location
10 20 30
XLAT XLONG SEASON
FORMAT(3F10.)
XLAT = latitude of the natural water system being modeled (degrees)
XLONG = longitude of the natural water system being modeled (degrees)
SEASON = time of year: 1 - spring
2 - summer
3 - fall
4 - winter
5 - annual average
Note: if the system being modeled encompasses a wide range of
latitude or longitude use average values in the model.
5.12.6 Water column light extinction
XKE(l) XKE(2) . . . XKE(NOWAV)
XKE(I) = light extinction coefficient (I/meter) for wavelength
number I
79
-------
5.12.7 Effect of Cloud Cover
FCLOUD BETA
FORMAT(2F10.)
FCLOUD = fraction of sky that is cloud covered
BETA = light reduction factor based on cloud type
Cloud type BETA (avg. values)*
Cirrus 0.8 - 0.89
Cirrostratus 0.65 - 0.84
Altocumulus 0.45 - 0.52
Altostratus 0.41
Stratocumulus 0.31 - 0.35
Stratus 0.24 - 0.25
Fog 0.17 - 0.18
Leighton, Photochemistry of Air Pollution, Academic Press, 1961.
5.13 CARD GROUP K - INITIAL CONDITIONS
5.13.1 Initial Conditions
5
ANAME(I)
35
ANAME(I+2)
15
20
30
C(ISYS,I) ANAME(I+1) C(ISYS,I+1)
45
C(ISYS,I+2)
50
ANAME(I+4)
60
C(ISYS,I+4)
FORMAT(4(A5,F10.0))
ANAME(I) - a one to five alpha-numeric character descriptive name the
initial condition in segment I of system ISYS.
C(ISYS,I) initial condition in segment I of system ISYS in the appro-
priate units (normally mg/£ or ppm)
The user will be required to input initial conditions for each system even
if the system is bypassed or if the initial conditions are zero. The initial
conditions are read in one system at a time (from system 1 through system
NOSYS), with concentrations being read from segment 1 through NOSEG within a
segment package.
80
-------
5.14 CARD GROUP L - STABILITY CRITERIA
10 20 30 80_
CMAX(I) CMAX(I+1) CMAX(I+2) CMAX(NOSYS)
FORMAT(8E10.0)
CMIN(I) - stability criteria for system I; i.e., the concentration
(normal units mg/Jl or ppm) above which program execution
will be terminated because of probable numerical instabil-
ity.
5.15 CARD GROUP M - PRINT CONTROL
5.15.1 Print Interval
10
P R N T
FORMAT(F10.0)
PRNT - print interval in days. (Set to zero for steady state option)
5.16 CARD GROUP N - INTEGRATION INFORMATION
5.16.1 Integration Interval, Total Time
10 20
DT NTCHA
FORMAT(2F10.0)
For time variable option
DT - integration step size (normal units - days)
NTCHA - total time of time variable simulation
For steady state option
DT - set equal to 0.0
81
-------
5.17 CARD GROUP 0 - DISPLAY PARAMETERS
All input is read NOSYS times, once for each system of the model,
For each variable-segment combination chosen a time history of the segment
will be displayed (dumped). For steady-state option no input is required;
skip this card group.
5.17.1 Variable Names
8 16 80_
ANAME(l) ANAME(2) ANAME(IO)
FORMAT(10A8)
ANAME(K) = a one to eight alphanumeric character descriptive
name for display variable K. The order of these
names is determined via the appropriate disk file
WRITE in the users kinetic subroutine 2.
5.17.2 Variable Number, Segment Numbers
27 31
VARNO SEG(K) SEG(K+1) ... SEG(K+7) ALLDMP
FORMAT (913, A4)
VARNO = the position of the desired variable, to be displayed,
in the WRITE file statement in the kinetic subroutine
(see previous note) .
SEG(K) = segment number to be displayed (note: order of display
unimportant, i.e., need not be sequential).
ALLDMP = flag set equal to ALL if all segments in the model are to
be displayed.
A blank card terminates display for system, ISYS. Then, another Variable Name
card, followed by Variable Number, Segment Number card(s) ±s read until system
NOSYS has been read.
NOTE; ^ thg current version of the WASPB subroutine the
variables written to
output are:
VARNO SYSTEM 1 SYSTEM 2
1 Total toxicant (yg/£) Solids concentration (mg/£)
2 Dissolved toxicant (yg/il)
3 Particulate toxicant yg/g
82
-------
SECTION 6
EXAMPLE APPLICATIONS
6.1 KEPONE-JAMES RIVER, VA. (19)
The discharge of Kepone to the James River estuary from 1966 through
1975 resulted in widespread contamination of the estuarine system, extend-
ing the 120 kilometers from Hopewell, Virginia to the mouth of the James
and into Chesapeake Bay. A model was developed to analyze the distribu-
tion of Kepone and project the time required to reduce concentrations to
various levels. WASTOX was not available for use in the original work but
was later applied to verify the original calculation.
To predict Kepone concentrations in the water column and sediment,
the actual geomorphology, hydrology and tidal phenomena must be ade-
quately reflected. To accomplish this, the water body must be divided
into a sufficient number of segments so that each segment represents
localized parameters of the system. In this way, variations of any para-
meter along the length of the river can be taken into account. Segment
lengths roust be short enough so that expected gradients in water and sedi-
ment concentrations can be accurately calculated.
With the above criteria in mind, the WASTOX version of the James
River model was constructed of 74 completely mixed segments. The channel
of the river consists of 56 segments and the four side bays and two tribu-
taries of the river are represented by a total of 18 segments. Four
layers of segments are used for the main channel, two in the water column
and two in the sediment. In the more shallow side bays and tributaries,
three vertical layers of segments are used, one for the water and two for
the bed.
A schematic of an elevation view of the segments in the main channel
is shown in the upper portion of Figure 3. Note that segments 1 through
28 define the two layers of the water column, segments 35 through 48 are
in the first sediment layer and segments 55 through 68 constitute the
second sediment layer.
The top sediment layer, subjected to horizontal movement by the lower
water layer velocities, is referred to as the "moving sediment" layer or
83
-------
MAIN CHANNEL
WATER
COLUMN
BED
1
2
35
55
3
4
36
56
5
6
37
57
7
8
38
58
9
10
39
59
11
12
40
60
13
14
41
61
15
16
42
62
17
18
43
63
19
20
44
64
21
22
45
65
23
24
46
66
25
26
47
67
27
28
48
68
c
H
E
S
A
P
E
A
K
E
B
A
Y
134 124 114 104 94 84 74 64 54 44 34 24
BAYS & TRIBUTARIES
14 4
WATER
COLUMN
BED
1
29
49
69
3 5
30
50
70
31
51
71
11
32
52
72
15
33
53
73
19 21
34
54
74
BAILEYS CHICKAHOMINY
APPOMATTOX BAY TAR RIVER COBHAM
RIVER BAY BAY
BURWELL'S
BAY
-1 : KM
POINT
FIGURE 3. SEGMENTATION OF THE JAMES RIVER USED IN THE KEPONE ANALYSIS
84
-------
the "moving bed" and the second sediment layer is the "stationary bed".
As shown in the lower portion of the figure, segments 29 through 34 are
the water seg ments, sequentially numbered, in the Appotnattox River,
Baileys Bay, Tar Bay, Chickahominy River, Cobham Bay and Burwell Bay. For
the same locations, segments 49 through 54 and 69 through 74 represent the
moving and stationary bed segments, respectively.
In the main channel, segment lengths are all approximately 10 kilo-
meters long. With 14 segments end-to-end, the model represents 135 kilo-
meters of river from its mouth to above Hopewell.
The water column was divided into two layers to account for the net
tidal flow profile in the estuary in which high salinity water moves up-
stream in the lower layer and less dense low salinity water moves down-
stream in the upper layers. The depths of the two layers were defined by
the plane of no net motion, i.e., the location of zero velocity where the
flow direction changes. Based on available sediment data, depths of 10
and 30 centimeters were selected for the moving and stationary beds,
respectively.
The analysis procedure, which is appropriate for those cases in
which the system has been previously subjected to an input of a toxic
substance, essentially consisted of four steps, shown diagrammatically in
Figure 4. First, equations of continuity, momentum and state were
employed in a steady-state tidally averaged mode to generate the
horizontal and vertical velocities of the estuarine circulation. A
steady state model of the water column was then used to determine the
vertical dispersion coefficient between the two water layers using ocean
salt as a tracer. Secondly, settling velocities of solids to the bed and
resuspension of solids from the bed were estimated to match observed
solids concentration in the water column. A settling velocity was
assigned and assuming a surface bed solids concentration of 50,000 mg/£
the resuspension velocities needed to match suspended solids
concentrations were determined. These velocities were checked by
calculating the dissolved and particulate organic Kepone concentrations
in the water column using the observed Kepone concentrations in the
surficial sediments as a boundary condition. From this procedure a
settling velocity of 4 ft/d (1.2 ro/d) was assigned as characteristic of
the James River solids.
85
-------
UPSTREAM
V
1-!
n
2-
3
3 g 7 - .
1 2 vv'
v
^ [5 7 ^=r
^ 2x^
v
^J^l j ^~
•^J" ^_
A-
*
€ ^(yy f
172 g
3,
\
wu-' '•'•'•'•'
SUSPENDED SOLIDS
Q 3 TRANSPORT LAYER ,,
i.. ^
1
— V r>
00 4 STATIONARY BED
i,.,
wd
(SEDIMENTATION)
•'
M/ V
Ul / ^™
r^»
f | 2 |
5 ••••••••••••
ki .'.'.'.».'.'.'.'.
OQ 4 •:•:•:•:•:•:•:•:•
!
N
^2-5
"°J-4
~l>xMx
SUSPENDED AND BED SOLIDS
DISSOL VED AND PARTICULA TE ORGANIC CHEMICAL
t
t
FIGURE 4. SEQUENCE OF STEPS IN JAMES RIVER USED IN THE KEPONE
ANALYSIS(23)
86
-------
In the third step, the solids transport in the bed was addressed.
Net horizontal velocities of the transport layer were assigned and the
flux of sed iment between the two sediment layers was determined by mass
balance. With selected values of solids concentrations for both the
transport and stationary beds of 50,000 and 500,000 mg/Jl respectively and
the settling and resuspension rates from step two, the net sedimentation
rate of the bed was calculated and compared to observed rates in the
estuary.
The last step consisted of calculating the dissolved and particulate
con centrations of Kepone in the water and bed. The diffusivity of the
dissolved component in the interstitial water was assigned as 5 x 10
2
cm /s. A solids dependent sorption partition coefficient was used.
The Kepone distribution in the water and the bed were first calcu-
lated for various freshwater flows held constant over the period of analy-
sis. The actual hydrograph for the time period beginning in 1965 was then
approximated by six constant flow conditions and the Kepone distribution
calculated. WASTOX has been applied only to the constant flow condition
although it could also be used for the variable flow simulation.
The results of the model for a freshwater flow of 1000 cfs are shown
in Figure 5.
87
-------
I
CO
g
co 40
a
iii
o
z
LU
Q.
W 0
CO
01
o»
UJ
O
Q.
LU
10.0
5.0
O)
O)
=5
tiJ
O
O.
UJ
0.4
0.2
LOWER WATER LAYER
SUSPENDED SOLIDS
UPPER •"•
WATER LAYER
LOWER WATER LAYER
ADSORBED KEPONE
STATIONARY SEDIMENT
ADSORBED KEPONE
' I
/ a
125
100
75
50
25
DISTANCE FROM MOUTH-km
FIGURE 5. COMPARISON OF OBSERVED AND COMPUTED SUSPENDED SOLIDS
AND KEPONE FOR THE 1000 CFS FRESHWATER FLOW
88
-------
6.2 PLUTONIUM-239, PCB - GREAT LAKES (21)
The Great Lakes represent the major fresh water resource of the
United States as well as one of the major resources for Canada. Because
of the significance of these bodies of water, there is a unique and long
terra concern over the fate of chemicals that may be potentially toxic to
the aquatic ecosystem as well as to the general public health. The dis-
charge of chemicals such as the heavy metals including mercury, and
organic chemicals such as polychlorinated biphenyls (PCB) and mirex have
resulted in the closing of commercial and sports fishing opportunities.
As a result of these concerns, a need exists to develop a framework for
the Great Lakes that would permit estimates to be made of the fate of
chemicals discharged to the Lakes.
WASTOX was used for the calculations associated with a physico-chemi-
cal model of the Great Lakes system. Full details of the model, calibra-
tion and application to chemical fate in the Great Lakes is given in
Thomann and Di Toro (21).
The time scale of the model is considered to be long term, i.e. year
to year. The physical segmentation of the model considers the Lakes to
be completely mixed with the exception of Lake Erie (Figure 6). This Lake
is divided into three basins; west, central, and east to reflect varying
regions of solids deposition and water column solids concentrations. In
addition, Saginaw Bay is included as a separate embayment from Lake Huron
to represent a more local region interacting with a large lake. Three
sediment segments of 2 cm each in depth are included under each of the
lakes or region of lake. This results in a model with eight water column
segments and twenty-four sediment segments for a total of thirty-two seg-
ments.
The calibration procedure was as follows: (a) From a review of data
on fine grain solids loading to the Lakes, net depositional flux of
solids, and water column suspended solids concentrations, the net loss
rate of solids from the water column was estimated. From assigned poros-
ity in the surface sediment layer, particle density and net flux of solids
to the sediment, the net sedimentation rate is computed. (b) With the
estimate of the net loss of solids, w , a range of particle settling velo-
89
-------
lil
W
&
8
M
ft.
90
-------
cities were assigned and the resuspension velocity necessary to maintain
the solids balance was computed. (c) Since there is an infinite number
of combinations of settling and resuspension velocities that will result
in the same solids balance, a time variable plutonium-239 model was used
to provide the tracer calibration. All decay mechanisms and sediment
diffusion were assumed to be zero and a sensitivity analysis using three
values of plutonium solids partitioning was conducted.
Following calibration of the model using plutonium-239 as a tracer,
the model was applied to PCB, benzo(a) pyrene and cadmium.
The comparisons of observed plutonium-239 data in the 1970's to model
calculations are shown in Figure 7 (w = particle settling velocity, w
3 DG t
= net loss of solids, w = w indicates zero resuspension). The cali-
bration indicated that in general the sediments are interactive with the
water column in the Great Lakes through resuspension and/or horizontal
transport.
The results of the application of the model to PCB using a 20 year
calculation indicate that a load level ranging from 640 to 1390 kg/yr
with volatilization (at an exchange rate of 0.1 m/d) appears to be repre-
sentative of observed surface sediment data for the open lake waters (Fig-
ure 8). Fifty percent response times for PCB following cessation of load
and including volatilization varied from less than 5 years to 10-20 years
for the other lakes without volatilization. Comparison of these response
times to decline of concentrations of PCB in Lake Michigan indicates that
at least for that lake volatilization is occurring at an exchange rate of
about 0.1 m/d.
6.3 PENTACHLOROPHENOL - EPA MERS CHANNELS (22)
Pentachlorophenol (PCP) and its salts have historically been used for
the preservation and treatment of wood. However, their antimicrobial,
antifungal, herbicidal and insecticidal properties have led to widespread
application of PCP. Increased knowledge of the behavior and fate of PCB
in the aqueous environment will enable its impact to be effectively man-
aged. During 1982 EPA executed a PCP fate and effects study utilizing the
semi-natural stream channels at the Monticello Environmental Research Lab-
oratory. The 1400 ft. channels consist of alternative pool and riffle
91
-------
* 1.0
o
3 0.5
of
«
CM
0
19
1.0
™
U
/ 0.5
t
CM
of
«
CM
o
\ SUPERIOR
\
\
\ S1-"
\ "
:*%. \ ~
X 105
0)
«
CM
1 I I .!_ 1... .J 0
Nx MICHIGAN
\
\
\
h*i^ X
i i i i i i
71 72 73 74 75 76 77 1971 72 73 74 75 76 77
Year Year
HURON
= 1'°\
5 \
wa-wnet
=2.5 m/d
-*- =5.0 m/d
~ ^^**«B=£— A— A TA
0*0.5. X^l ^
0)
n
M 0 i i ' i i >
1971 72 73 74 75 76 77
Year
ERIE 1-0
\
U
V. c/0-5
i "T.rrr ~i \ i o
ONTARIO
[ T s=Sept.
A= Aug.
^=^_e^
^r^^31"
1971 72 73 74 75 76 77 1971 72 73 74 75 76 77
Year Year
FIGURE 7. COMPARISON OF 1971-1977 OBSERVED AND CALCULATED WATER COLUMN
039 ?40
' Pu CONCENTRATION FOR THREE CONDITIONS OF THE PARTICU-
LATE SETTLING VELOCITY.
92
-------
800i
800
800
Extended Load
Condition
Upper Level
Lower Level
Volatilization Rate
(m/d)
0.0
1
B
0.1
n
*
- Mean & Range
MICHIGAN
SAGINAW BAY
ONTARIO
FIGURE 8. CALCULATED SUERFACE SEDIMENT PCB CONCENTRATION (ng/g) FOR CON-
DITIONS ON EXTERNAL BED AND VOLATILIZATION RATE AND COMPARISON
TO OBSERVED DATA
93
-------
3
sections and during the study a constant 0.45 cfs (0.013 m /s) flow was
maintained. Process experiments to quantify photolysis and biodegradation
decay rates, two comprehensive survey data sets and the results of a dye
study were utilized to construct and calibrate a mathematical model of the
transport and fate of PGP in the channels. Manhattan College's toxic sub-
stance computer program (WASTOX) was used as the modeling framework.
Since photolysis was a major degradation mechanism for PGP, a time
variable model was constructed which utilized observed hourly solar radia-
tion measurements. The results of a detailed travel time/dispersion study
were used to verify that the transport component of the time variable
model. Each of the pool and the riffle sections of the channel was
divided into 3 model segments (42 total model sections). Comparison of
the first model run with the dye data, as shown in Figure 9, indicates
good agreement. Once the transport component of the model was judged
adequate the kinetic calibration of the model was initiated. Using the
observed biodegradation and photolysis rates from the process experiments
and the synoptic survey data as shown in Figure 10, the model was success-
fully calibrated using the June 1982 survey data. The results of the ana-
lysis indicate that photolysis is a major degradation pathway for PCP in
a natural environment. As shown in Figure 10, significant diurnal varia-
tion of PCP was observed in the channel which was successfully character-
ized by the model. As the MERS study proceeded, the importance of biodeg-
radation increased as the channel became acclimated to the PCP dosing.
94
-------
FIGURE 9: Model Comparison with September
16 and 17 Rhodamine WT Dye Data.
95
-------
MERS CHANNEL 5 STATION MW 1
MERS CHANNEL 5 STATION MW 2
180
^ 160
a
O 120
a.
-J 100
0 °"
k-
GO
• •
ktf • o 0
kd-0.6/d
180
^ 160
•3 140
0.
O 120
Q.
< 10°
0 80
K
60
\ . .'"1 •':
w^^-^-
_ _ 1 1- • ' ' 1, I 1
-10 0 10 20 30 40 50 60 70 BO
TIME (hours)
MEFI3 CHANNEL 5 STATION MW 3
180
5 16° I
?
-^ 140
TIME (hours)
MERS CHANNEL 5 STATION MW 4
180
^ 160
3-140 |
a.
O 120
a.
-i 100
<
o 80
0 10 20 30 40 50 60 70 80
TIME (hours)
MERS CHANNEL 5 STATION MW 7
180
^ 160
•5 140
a.
o 120
0.
tiU II,,
-10 0 10 20 30 40 50 60 70 60
TIME (hours)
MERS CHANNEL 5 STATION MW 8
180
C 160
0.
O 120
Q.
_l 100
O 80
-10 0 10 2'l JO to SO G.j la 80
TIME (hours)
CO
-10 0 10 20 30 40 50 60 70 80
TIME (hours)
MERS CHANNEL 5 STATION MW 13 MERS CHANNEL 5 STATION MW 14
J,4C
180
5 '60
40
a.
o 120
o.
< lo°
o 80
H
60
180
^ 160
•JM40
O.
O 120
0-
-I 100
O 80
60
U I • , .
-10 0 10 20 00 40 50 CO 70 80
TIME (hours)
-10 0 10 20 30 40 50 60 70 80
TIME (hours)
FIGURE 10: Time Variable Model Comparison with MERS Channel J
Total PCP data - Stations 1, 2, 3. 4, 5. 6, 7, 8, 13, 14
June 8 to 10, 1982
96
-------
SECTION 7
OPERATIONAL CONSIDERATIONS
This chapter describes how to obtain the computer program WASTOX,
how to install it on a DEC PDF mini computer, how to test the program
with a sample dataset, and what machine limitations limit the program.
7.1 AQUISITION PROCEDURES
To obtain the program WASTOX along with a sample dataset and support
software, write to:
Center for Water Quality Modeling
Environmental Research Laboratory
U.S. Environmental Protection Agency
College Station Road
Athens, GA 30613
A nine-track magnetic tape will be mailed to you. Please copy the con-
tents and return the tape.
7.2 INSTALLATION PROCEDURES
The subroutines that comprise WASTOX must be compiled and linked into
a task image. This is accomplished on the PDP IAS operating system by
running the command file "WXTCMP.CTL." If the compilation succeeds, then
linkage is automatically attempted with the command file "WXTLNK.CTL."
7.3 TESTING PROCEDURES
Once WASTOX is installed, the sample input dataset should be run and
compared with the sample output dataset to verify that the program is cal-
culating correctly. To perform a simulation on the PDP, submit the batch
input sequence "WXTRUN.CTL."
97
-------
7.4 MACHINE LIMITATIONS
Currently, TOXIWASP is set up for the following configuration:
PDF 11/70 Hardware
RSTS/E operating system
FORTRAN IV
60 segments - steady-state
75 segments - time-variable
4 systems
The PDP 11/70 computer utilizing RSTS/E operating system allocates a
32k word (64k byte) user area for execution of programs. WASTOX occupies
at least 31k words of memory. Any enlargement of this program may result
in an overflow of the user area.
98
-------
REFERENCES
1. Lewis, W.K. and W.C. Whitman. Principles of Gas Absorption, Ind.
& Eng. Chem. 16. 1924.
2. Danckwert, P.V. Significance of Liquid Film Coefficients in Gas
Absorption. Ind. & Eng. Chem. 42. June 1951.
3. O'Connor, D.J. and W.E. Dobbins. Mechanism of Reaeration in Natural
Streams. Transactions, ASCE, 641:123. 1958.
4. O'Connor, D.J. The Effect of Winds on the Gas-Liquid Transfer
Coefficient. Manhattan College, Environ. Eng. & Sci. 1981.
5. Zepp, R.G. and G.L. Baughman. Prediction of Photochemical Tranformation
of Pollutants in the Aquatic Environment ±n_ Aquatic Pollutants: Trans-
formation and Biological Effects, edited by 0. Hutzinger, L.H. Van
Lelyveld and B.C.J. Zoeteman. Pergamon Press. 1978.
6. Green, A.E.S., T. Sawada, and E.P. Shettle. The Middle Ultraviolet
Reaching the Ground. Photochem. & Photobio. 19:251-259. 1974.
7. Green, A.E.S., K.R. Cross, and L.A. Smith. Improved Analytic
Characterization of Ultraviolet Skylight. Photochem. & Photobio.
31:59-65. 1980.
8. Zepp, R.G. and D.M. Cline. Rates of Direct Photolysis in the Aquatic
Environment. Environ. Sci. & Technol. 11:359. 1977.
9. Leighton, P.A. Photochemistry of Air Pollution. Academic Press, N. Y.
1961.
10. Horvath, R.S. Microbial Co-Metabolism and the Degradation of Organic
Compounds in Nature. Bacteriological Reviews. 36(2):145-155. 1972.
11. Lambert, S.M., P.E. Porter, and H. Schieferstein. Movement and Sorption
of Chemicals Applied to the Soil. Weeds 13:185-190. 1965.
12. Haque, R., F.T. Lindstrom, V.H. Freed and R. Sexton. Kinetic Study
of the Sorption of 2,4-D on Some Clays. Environ. Sci. Technol. 2(3):
207-211, 1968.
14
13. Garnas, R.L., A.W. Bourquin, and P.H. Pritchard. The fate of C-Kepone
in Estuarine Microcosms. In Kepone in the Marine Environment, Appendix
C to the EPA Kepone Mitigation Project Report. EPA 440/5-78-004C. 1978.
14. Karickhoff, S.W., D.S. Brown, and T.A. Scott. Sorption of Hydrophobic
Pollutants on Natural Sediments. Water Res. 13:241-248. 1979.
99
-------
REFERENCES (continued)
15. Hamaker, J.W. and J.M. Thompson. Adsorption. Ir\_ Goring, C.A.I, and
J.W. Hamaker (Editors), Organic Chemicals in the Soil Environment pp.
49-142. Marcel Dekker, Inc. New York. 1972.
16. Poinke, H.B. and G. Chesters. Pesticide-Sediment-Water Interactions.
J. Environ. Qual. 2(1):29-45. 1973.
17. Connolly, J.P. The Effect of Sediment Suspension on Adsorption and Fate
of Kepone. Ph.D. Thesis, The University of Texas at Austin. 1980.
18. O'Connor, D.J. and J.P. Connolly. The Effect of Concentration of
Adsorbing Solids on the Partition Coefficient. Water Res. 14:1517-1523.
1980.
19. O'Connor, D.J., K. Farley, and J. A. Mueller. Mathematical Models of
Toxic Substances in Estuaries with Application to Kepone in the James
River. Draft Final Report, EPA Grant No. R804563. 1981
20. Di Toro, D.M. , J.J. Fitzpatrick, and R.V. Thomann. Water Quality Analy-
sis Simulation Program (WASP) and Model Verification Program (MVP) - Doc-
umentation. Hydroscience, Inc., Westwood, N.J. for U.S. Environmental
Protection Agency, Duluth, MN, Contract No. 68-01-3872. 1982.
21. Thomann, R.V. and D.M. Di Toro. Physio-Chemical Model of Toxic Sub-
stances in the Great Lakes. Final Report to USEPA, Grosse lie, Michigan.
1983.
22. Winfield, R.P., M. Labiak, D.M. Di Toro, and W.L. Richardson. Mathemati-
cal Model of the Fate of Pentachlorophenol in the EPA MERS channels.
Final Report to USEPA, Grosse He, Michigan. 1983.
23. O'Connor, D.J., J.A. Mueller, and K.J. Farley. Distribution of Kepone
in the James River Estuary. ASCE, J. Environ. Eng. Div., 109(2):396-413.
1983.
100
-------
APPENDIX 1
GLOSSARY OF VARIABLES
I. Variables in Common Block MAIN
BBC - Matrix of order (NOSYS, Max. No. of boundaries per system)
Contains the boundary concentrations in mg/&.
BQ - Vector of length (Max. No. of Flows)
ft 3
Contains flows in 10 ft /d.
BR - Vector of length (Max. No. of exchanges)
ft O
Contains the exchanges in 10 ft /d.
BVOL - Vector of length (NOSEG)
/• n
Contains the segment volumes in 10 ft
BWK - Matrix of order (NOSYS, Max. No. of loads per system)
Contains the loadings or discharge rates in Ib/d.
C - Matrix of order (NOSYS, NOSEG)
Contains the current concentrations for each system and segment
in mg/£-d.
CD - Matrix of order (NOSYS, NOSEG)
Contains the current time derivative for each system and segment
in mg/£-d.
CMAX - Vector of length (NOSYS)
Contains the maximum allowable concentration for each system in
mg/2..
CONST - Vector of length (Max. No. of Constants)
Contains the kinetic coefficients used in WASPB
DT - Integration interval in days.
IBC - Matrix of order (NOSYS, Max. No. of boundaries per system)
Contains the segment numbers for the corresponding boundary
concentrations in BBC.
IBCOP - Vector of length (NOSYS)
Contains the option numbers for inputting boundary
concentrations.
101
-------
IQ - Vector of length (Max. No. of Flows)
Contains the numbers of the segments to which flow is directed
for the corresponding flows in BQ.
IQOPT - Option number for inputting flows.
IR - Vector of length (Max. No. of Exchanges)
Contains the segment numbers for the first segment: corresponding
to the exchanges in BR
IROPT - Option number for inputting exchanges.
ITIMB - Vector of length (NOSYS)
Counter for time variable boundary conditions.
ITIMD - Counter for time variable dispersions.
ITIMF - Counter for time variable straight line functions.
ITIMQ - Counter for time variable flows.
ITIMR - Counter for time variable exchanges.
ITIMV - Counter for time variable volumes.
ITIMW - Vector of length (NOSYS)
Counter for time variable forcing functions.
IVOPT - Option number for inputting volumes.
IWK - Matrix of order (NOSYS, Max. No. of loads per system)
Contains the numbers of the segments in which the corresponding
discharges in BWK are occurring.
IWKOP - Vector of length (NOSYS)
Contains the option number for inputting discharges or (forcing
functions).
JQ - Vector of length (Max. No. of Flows)
Contains the numbers of the segments from which flow is leaving
for the corresponding flows in BQ.
JR - Vector of length (Max. No. of Exchanges)
Contains the segment numbers for the second segment corresponding
to the exchanges in BR.
102
-------
MBC - Matrix of order (NOSYS, Max. No. of boundaries per system)
For a time-varying boundary concentration MBC contains the
current slope of the broken line approximation. For a time-
varying option BBC contains the intercept of the current line
segment of the broken line approximation.
MQ - Vector of length (Max. No. of flows)
Slope for time-varying flows analogous to MBC.
MR - Vector of length (Max. No. of Exchanges)
Slope for time-varying exchanges analogous to MBC.
MVOL - Vector of length (NOSEG)
Slope for time-varying segment volumes analogous to MBC.
MWK - Matrix of order (NOSYS, Max. No. of Loads per system)
Slope for time-varying loads or (forcing functions) analogous
to MBC.
NFUNC - No. of user defined broken line functions.
NFUNT - Time of next break in user defined broken line functions in days,
NOBC - Vector of length (NOSYS)
Contains the number of boundary conditions for each system.
NOQ - Contains the number of flows for a particular flow field.
NOR - Contains the number of exchanges for a particular exchange field.
NOWK - Vector of length (NOSYS)
Contains the number of loads (forcing functions) for each system.
NQT - Time of next break in flow broken-line function in days.
NRT - Time of next break in exchange broken-line function in days.
NTCHA - Total time of run in days.
NVOLT - Time of next break in segment volume* b"roken-line function in
days.
NWKT - Vector of length (NOSYS)
Contains the time of the next break in the forcing function
broken-line function for each system in days.
103
-------
PARAM - Matrix of order (NOSEG, Max. No. of parameters per segment)
Contains the values for each segment of the parameters used in
WASPB.
PRNT - Print interval in days.
TIME - Current time in days.
II. Variables in Common Block SYSTRN
FIELDQ - Vector of length (Max. No. of flow fields)
Contains the number of flows read in for the current and all
prior flow fields.
FIELDR - Vector of length (Max. No. of exchange fields)
Contains the number of exchanges read in for the current and all
prior exchange fields.
FPART - Matrix of order (NOSYS, NOSEG)
Contains the fraction of total chemical in each segment that is
dissolved or adsorbed to the solids of each system.
SCALEQ - Matrix of order (NOSYS, Max. No. of flow fields)
Contains the scale factors by which the flows in each flow field
are multiplied for each system.
SCALER - Matrix of order (NOSYS, Max. No. of exchange fields)
Contains the scale factors by which the exchanges in each flow
field are multiplied for each system.
III. Variables in Common Block MAINA
IN - Unit number designation for input.
IREC - Counter for number of times results are written to disk.
IRECL - Vector of length (NOSYS)
Contains the record number of the first record of the results on
disk for a specific time.
ISYS - Current system being processed.
104
-------
NOSEG - Number of segments in the model.
OUT - Unit number designation for output.
SYSEXC - Vector of length (NOSYS)
Contains the system bypass option for each system.
IV. Variables in Common Block PHOTO
DAYK - Vector of length (NOSEG)
Contains the photolysis rate constants for each segment.
V. Variables in Common Block COMP
BFUNC - Vector of length (Max. no. of user defined broken-line
functions)
IEV - Switch to indicate to WASPS which input broken-line function is
being evaluated.
ITIME - Counter for user defined broken-line functions.
MFUNC - Vector of length (Max. no. of user defined broken-line functions)
Contains the current slope of each broken-line function.
NBCT - Vector of length (Max. no. of boundary conditions)
Time of next break in each boundary condition broken-line
function.
NDEVL - Switch to indicate time for writing results on disk.
VI. Variables in Common Block DUMP
MXDMP - Number of variables written to disk for each system each time
results are dumped.
MXSEG - Maximum number of segments allowed in the model.
105
-------
VII. Local variables in the kinetic subroutine WASPB
A1,A2,A3 - Solids dependent adsorption experiment for solids type 1,2, or
3.
ATMOS - Concentration of chemical in the air above the water body in
AREA
- Vector of length (NOSEG)
2
Area of each segment at its downstream face in it .
B1,B2,B3 - Solids dependent adsorption coefficient for solids type 1,2,
or 3.
BACS1
BACS2
BACS3
BACW
BIOS
BIOW
DEPTH
DIFF
DISTOX
DTIME
DVMMY
DWNSEG
H
Vectors of length (NOSEG)
Concentration of compound degrading bacteria on solids type.
1,2, or 3 in num/mg.
Vector of length (NOSEG)
Concentration of compound degrading bacteria in water in
num/£.
Vector of length (NOSYS-1)
Computed rate of biodegradation on each solids type in ug/Jl-d.
- Computed rate of biodegradation in water in
- Vector of length (NOSEG)
Depth of each segment in meters
2
- Molecular diffusivity of chemical in m /s
- Concentration of dissolved chemical in yg/£
- Time between current time and next break in user defined
broken-line approximation.
- Extra variable used when results are written to disk to keep
the number of variables written for each system equal to MXDMP
when less than MXDMP variables are needed for a system.
- Vector of length (NOSEG)
Number of segment downstream of each segment.
- Hydrogen ion concentration in moles/Jl.
106
-------
HENRY - Henry's Constant - unitless.
HYDROL - Computed hydrolysis rate in vg/H for dissolved chemical
IFLOW - Switch to indicate a flowing(O) or non-flowing(l) water body.
KD - Biodegradation second order rate constant for dissolved
chemical in Jl/d-num.
KDI - Temperature correction factor for biodegradation of dissolved
chemical.
KD20 - Biodegradation second order rate constant for dissolved
chemical at 20°C in Jt/d-num.
KH - Total hydrolysis rate constant in 1/d.
KHH - Acid hydrolysis second order rate constant in Jt/d-mole.
KHH20 - Acid hydrolysis second order rate constant at 20°C in
£/d-mole.
KHN - Neutral hydrolysis first order rate constant in 1/d.
KHN20 - Acid hydrolysis second order rate constant at 20°C in
fc/d-mole.
KHOH - Alkaline hydrolysis second order rate constant in Jl/d-mole.
KHT - Temperature correction factor for hydrolysis.
KL - Volatilization rate constant in m/d.
KLA - Volatilization rate constant in 1/d.
KLT - Temperature correction factor for volatilization.
KP - Biodegration second order rate constant for adsorbed
chemical in £/d-num.
KPT - Temperature correction factor for biodegradation of adsorbed
chemical.
KP20 - Biodegradation second order rate constant for adsorbed
chemical at 20°C in Jl/d-num.
MOUNT - Molecular weight of chemical
OH - Hydroxide ion concentration in moles/a
107
-------
PART - Vector of length (NOSYS)
Partition coefficient for each system in £/mg
PARTOX - Vector of length (NOSYS-1)
Concentration of adsorbed chemical for each solids type
in yg/Jl
PH - Vector of length (NOSEG)
pH of each segment
PHOTOL - Photolysis rate in ug/£-d
PORE - Porosity of a segment
RHO - Vector of length (NOSYS-1)
density of each solids type in g/m£
TDEPTH - Vector of length (NOSEG)
Depth from the water surface of the top of the
segment in meters
TEMP - Vector of length (NOSEG)
Temperature of each segment in °C
VEL - Flow velocity in m/s
VOLAT - Volatilization rate in ug/&-d
XKGH - Gas film transfer coefficient
XKL - Liquid film transfer coefficient
108
-------
APPENDIX 2
Listing of Kinetic Subroutine WASPB
109
-------
FORTRAN IV
0001
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
V02.5-2 Tue 13-Sep-83 16:45:
SUBROUTINE WASPB
38 PAGE 001
THIS MODEL IS SET UP TO COMPUTE THE SPATIAL AND TEMPORAL
DISTRIBUTION OF TOXIC CHEMICALS IN
TOTAL TOXICANT AND UP TO 3 CLASSES
OF TOXICANT BETWEEN THE DISSOLVED
PHASES IS COMPUTED FROM THE TOTAL
ESTUARIES. IT MODELS
OF SOLIDS. THE DISTRIBUTION
AND VARIOUS PARTICULATE
TOXICANT ASSUMING
INSTANTANEOUS EQUILIBRIUM LINEAR PARTITIONING.
SYSTEM 1 IS TOTAL TOXICANT
SYSTEMS 2 THROUGH N ARE SOLIDS
>***********************************
TTTT^TT~~~^~^~~^~^~TTVTVT^T^TTT^~T~
THIS MODEL IS DIMENSIONED FOR
4 SYSTEMS
75 SEGMENTS
100 EXCHANGES
200 FLOWS
50 BOUNDARIES
17 FORCING FUNCTIONS
11 SYSTEM PARAMETERS
30 CONSTANTS
*****************************
^TTV~T^TT^TT^TTV^TT^TT^^^^^T^
THE FOLLOWING:
C * 4 STRAIGHT LINE FUNCS
C * 9 FLOW FIELDS
C * 9 EXCHANGE FIELDS
C *
C *
C *
C *
0002
COMMON /MAIN/
* BBC( 4,50), BCK200), BR(IOO),
* CC 4, 75), CDC 4, 75), CMAXC 4)
BVOLC 75), BWKC 4, 17),
, CMINC 4), CONSTC30),
* DT, IBCC 4,50), IBCOPC 4), 10(200),
* IQOPT, IRC100), IROPT,
* ITIMF, ITIMQ, ITIMR,
* IVOPT, IWKC 4, 17), IWKOP(
* MBCC 4,50), MQC200), MR(IOO),
* NFUNC, NFUNT, NOBCt 4)
* NOWKC 4), NOT, NRT,
ITIMB( 4), ITIMD,
ITIMV, ITIMWC 4),
4), JOC200), JRC100),
MVOLC 75), MWKC 4, 17),
, NOQ, NOR,
NTCHA, NVOLT,
* NWKTC 4), PARAMC 75,11),PRNT,
* TIME
C *
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
0003
REAL
110
-------
FORTRAN IV
¥02,5*2
Tue 13-Sep-83 16|45|38
PAGE 002
0004
OOOS
0006
0007
0008
0009
0010
c
c
c
c
c
c
c
c
* MBC,
* NFUNT,
* NWKT
*
*
COMMON /MAINA/
* IN,
* NOSEG,
*
INTEGER
* OUT,
*
*
COMMON /COMP/
* BFUNCC 4),
* NDEVL
*
REAL
* MFUNC,
*
*
MO, MR, MVOL, MWK,
NOT, NRT, NTCHA, NVOLT,
IREC, IRECLC 4), ISYS,
NOSYS, OUT, SYSEXC( 4)
SYSEXC
IEV, ITIME, MFUNCC 4), NBCT( 4),
NBCT, NDEVL
COMMON /SYSTRN/
C
C
C
C
C
C
* SCALEQC 4,
* ,FPARTC 4,
*
INTEGER
* FIELDOC 9),
*
*
*
9),SCALER( 4, 9 ) , FIELDQ , FIELDR
75)
FIELDRC 9)
********************************************************************
COMMON /DUMP/
MXDMP, MXSEG
C*************#*****n*** 4C****************************»*##****#****#***#*»
C
C DESCRIPTION OF PARAMETERS AND CONSTANTS
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
PARAM NAME
DEPTH
TDEPTH
DWNSEG
AREA
PH
TEMP
WVEL
BACW
BACSI
DESCRIPTION
DEPTH OF THE SEGMENT
DEPTH FROM THE WATER SURFACE OF THE
TOP OF THE SEGMENT
SEGMENT DOWNSTREAM OF THE CURRENT SEGMENT
AREA BETWEEN THE CURRENT SEGMENT AND THE
DOWNSTREAM SEGMENT
PH OF THE SEGMENT
TEMPERATURE OF THE SEGMENT
WIND VELOCITY ABOVE THE SEGMENT
CONCENTRATION OF BACTERIA IN THE WATER
CONCENTRATION OF BACTERIA ON SOLIDS
TYPE I
UNITS
METERS
METERS
FT**2
CENTIGRADE
M/S
NUM/L
NUM/MG
111
-------
FORTRAN IV
V02.5-2
Tue 13-S6P-83 16:45:38
PAGE 003
0011
0012
0013
0014
0015
OOlb
0017
0018
0019
0020
0021
0022
0023
C
C
C CONST NAME DESCRIPTION UNITS
C
C
C FOR DISSOLVED TOXICANT AT 20 DEGREES C L/D-NUM
C
C
C
C FOR SORBED TOXICANT L/D-NUM
C
C
C
C CONSTANT AT 20 DEGREES C L/MOLE-D
C
C CONSTANT AT 20 DEGREES C I/DAY
C
C AT 20 DEGREES C L/MOLE-D
C
C
C
C VOLATILIZATION) UNITLESS
C
C
C
C
C
C
C
C
C
C
C
C RHOI DENSITY OF SOLIDS TYPE I GM/ML
C
£***************************************************** *****************
C
C
CONST NAME
KD20
KDT
KP20
KPT
KHOH20
KHN20
KHH20
KHT
X-SEG
HENRY
MOLWT
KLT
ATMDS
IFLDW
AI
RI
RHOI
DESCRIPTION
2ND ORDER BIODEGRADATION RATE CONSTANT
FOR DISSOLVED TOXICANT AT 20 DEGREES C
TEMP. CORRECTION FACTOR FOR BIODEGRA-
DATION OF DISSOLVED TOXICANT
2ND ORDER BIODEGRADATION RATE CONSTANT
FOR SORBED TOXICANT
TEMP. CORRECTION FACTOR FOP BIODEGRA-
DATION OF SORBFD TOXICANT
2ND ORDER ALKALINE HYDROLYSIS RATE
CONSTANT AT 20 DEGREES C
1ST ORDER NEUTRAL HYDROLYSIS RATE
CONSTANT AT 20 DEGREES C
2ND ORDER ACID HYDROLYSIS RATE CONSTANT
AT 20 DEGREES C
TEMP. CORRECTION FACTOR FOR HYDROL.
NUMBER OF WATER COLUMN SEGMENTS
HENRYS CONSTANT (SET TO ZERO TO SKIP
VOLATILIZATION)
MOLECULAR WEIGHT OF TOXICANT
TEMPERATURE CORRECTION FACTOR FOR
VOLATILIZATION
CONCENTRATION OF TOXICANT
ATMOSPHERE
FLAG INDICATING A FLOWING
A NON-FLOWING ( 1.) WATER
MASS DEPENDENT ADSORPTION
FOR SOLIDS TYPE I
MASS DEPENDENT ADSORPTION COEFFICIENT
FOR SOLIDS TYPE I
DENSITY OF SOLIDS TYPE I
IN THE
( 0.) OR
BODY
EXPONENT
C
C
DIMENSION PART(4),RHO(3),PARTOX(3),BIOS(3)
EQUIVALENCE(TIME,!)
RRAL KD,KDT,KP,KPT,KD20,KP20,KHOH,KHQH20,KHN,KHN20,KHH,KHH20
REAL KH,KHT,KL,KLA,KLT,MOLWT
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
(CONST(1),KD20)
(CON5T(2),KDT),(CQNST(3),KP20),(CQNST(4),KPT)
(CONST(5),KHOH20),(CONST(6),KHN20),(CONST(7),KHH20)
(CONST(8),KHT),(CONST(9),XWSEG),(CONST(10),HENRY)
(CONST(li),MOLWT),(CONST(12),KLT)
(CONST(13),ATMOS),(CONST(14),IFLOW),(CONST(15),AIPTMP)
(CQNST(16),A1),(CONST(17),R1), (CONST (18 ), RHOI )
(CONST(19),A2),(CONST(20),B2),(CONST(21),RH02)
(CONST(22),A3),(CONST(23),B3),(CONST(24),RH03)
112
-------
FORTRAN IV
V02.5-2
Tu« 13-S6P-83 16:45)38
PAGE 004
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0036
0037
0038
0039
0040
0041
0042
DIMENSION PH(75),TEMP(75),BACW(75),8ACS1(75)
DIMENSION BACS2(75),BACS3(75)
DIMENSION DEPTH(75),TDEPTH(75),DWNSEG(75),AREA(75),WVEL(1)
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
EQUIVALENCE
(PARAM(1,1),DEPTH(1)),(PARAM(J,2),TDEPTH C1)J
(PARAM(1,3),DWNSEG(1)>,(PARAM(1,4),AREA(1))
(PARAM (1,5),PHU)),(PARAM(1,6),TEMP(1))
-------
FORTRAN IV
V02.5-2
Tue 13-SCP-83 16:45:38
PAGE 005
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0073
0074
0076
0078
0079
0081
0082
0083
0084
0085
0086
0087
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
C
C
C
c
XLAM2 = 4.
CDRAG * 0,0011
SEGMENT LOOP
no 100 i a I,NOSEG
STP20 a TEMP(I) - 20.
COMPUTE SEGMENT POROSITY
SLV = 0.0
DO 26 J=2,NOSYS
?6 SLV s SLV t C(J,I)*l.E-06/RHU(J-l)
PORE = 1. - SLV
SYSTEM 1
TOXICANT
C
C
COMPUTE PARTICIPATE AND DISSOLVED FRACTION FOR EACH SEGMENT
SIIMMP a 0,0
CALCULATE PARTITION COEFFICIENTS
PART(l) = 0.0
IF(C(2,I),EO.O.O) GO TO 501
PARTC2) a 81/(C(2,I)**A1)
501 IF(^OSYS.GE,3,ANn.CC3,I).GT.0,0) PARTC3) = B2/(C(3,I)**A2)
IF(NOSYS.Ge.4.AND.C(4,I).GT.O.O) PART(4) = B3/(C(4,I)**A3)
nn 30 J = 2,NOSYS
IF(SYSEXC(J),EQ.l) GO TO 30
SUMMP = SUMMP + C(J,I)*PART(J)
30 CONTINUE
FPART(1,I) = l./(PORE+SUMMP)
DISTOX = FPART(1,T)*C(1,I)
DO 35 0 = 2,NOSYS
FPAPT(J,I) = CCJ,I)*PART(J)/(PORE-»-SUMMP)
35 PARTOX(J-l) a FPART(J,I)*C(1,I)
114
-------
FORTRAN IV
V02.5-2
Tue 13-SCP-83 16:45138
PAGE 006
0088
0089
0090
0091
0092
0093
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0106
0107
0108
0109
0110
0112
0114
0116'
0117
0118
0119
0120
0122
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
BIODEGRADATION : SECOND ORDER KINETICS ASSUMED
KD s KD20*(KDT**STP20)
KP a KP20*(KPT**STP20)
A) DEGRADATION IN WATER
BIOW = KD*BACW(I)*DISTOX*PORE
B) DEGRADATION ON SEDIMENT
DO 40 J=2,NOSYS
JJ * J-l
IF(C(J,I).GT,0.)GO TO 39
BIQS(JJ) « 0.0
GO TO 40
39 BIOS(JJ) = KP*BACS1(I)*PAPTOX(JJ)*C(J,I)
40 CONTINUE
HYDROLYSIS : ALKALINE, NEUTRAL, AND ACID HYDROLYSIS CONSIDERED
COMPUTE HYDROGEN AND HYDROXIDE ION CONCENTRATION
PPHI = PH(I)
H a 1,/C10.**PPHI)
OH = 1.E-14/H
XXX ' KHT**(TEMP(I)-20.)
KHOH s KHQH20*XXX
KHN = KHN20*XXX
KHH = KriH20*XXX
KH s KHOH*OH t KHN t KHH*H
HYDROL = KH*DISTOX*PORE
PHOTOLYSIS: DIRECT PHOTOLYSIS ONLY
PHOTOL s DAYK(I)*DISTOX*POKE
VOLATILIZATION
VOLAT =0.0
IF(HENRY.EO.0.0) GO TO 47
IF (TDEPTH(I).NE.O.O) GO TU 47
IF (IFLOW,F.Q.I) GOTO 45
KK s 1
QSCAL = SCALEO(1,1)
KMAX = FIELDOU)
DO 931 K=1,KMAX
931 IF(JO(K).EQ.I.AND.IO(K).EQ.DWNSEG(I)) GO TO 933
WRITE(OUT,932) I
115
-------
FORTRAN IV
0123 932
V02.b-2
won 3l-Oct-83 19:57:17
PAGE 007
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
FOkMAT(//'s$s$ EKKOK IN SPECIFICATION OF DOWNSTREAM SEGMENT
2 FDR SEGMKNT',14,' $$$$')
CALL ABOR'JUB)
933 fin TO (44,43,43),IQOPT
43 KK = K
USF: CONVERSION FACTOR o.o8t>4 TO CONVERT FLOW FROM
10**b CUBIC FEET/DAY TO CFS
44 0 = (MQ(KK)*(T-NUT) + 8Q(K))*QSCAL/0.0864
CONVERSION FACTOR 0,3048 TO CONVERT VELOCITY FROM FPS
TO METF.HS/S
VEL = (U/AREAd)) * 0.3048
XKL = SQRT(D1FF*VEL/DEPTH(1)) * 86400,
XKGH = 100,*HENPY
KL = (1./C1./XKL + l./XKGH)) * KLT**STP20
GOTO 46
COMPUTE DENSITY (G/ML) AND VISCOSITY (M**2/S) OF AIR AND WATER
45 DENA = 0.001293/(1.+0,00367*AIRTMP)
OENW - \. ' B.8E-05*TEMP(I)
XNUA = (1.32 + 0.009*AIRTHP) * l.E-05
XNUW = (1U,**(1301./(998.333+R.1855*STP20+0.00585*STP20**2)-
$ 3.30233)/DF,iJW) * l.K-04
XKL = (DIFF/XNUW)**,666*SORT(CDRAG*DENA/DENW)*,74*WVEL(I)/XLAM2
XKG = (D1KF/XNUA)**.666*SORT(CORAG)*.74*1«VEL(I)/(;XLAM2)
XKGH a XKG*KENRY
KL = (!./(!./XKL + l./XKGH)) * KLT**STP20*86400.
46 KLA a KL/DEPTHCI)
VOLAT = KLA * (DISTOX*PORE - ATMOS/HENRY)
47 CONTINUE
KINETIC DERIVATIVE
2 + PHOTUL + VOLAT)
MULTIPLY HY VOLUME TO BE CONSISTANT WITH TRANSPORT CALCULATION
CPU,I) = CD(1,I)*BVOL(I)
C
C
C
SYSTEM 2
SOLIDS
C
C
116
-------
FORTRAN IV
V02.5-2
Fri 07-Qct-83 15;48S16
PAGE 008
0146
0147
0148
0149
0150
0152
0153
0154
0155
0156
015T
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
C
c
C
c
c
c
c
100
c
c
c
c
901
c
902
951
C
NO KINETICS FOR SOLIDS SYSTEMS
CD (2, I) * 0.0
CDC2,I) ' 0.0
cnc3,i)a o.o
CD(4,I)s 0.0
CONTINUE
DUMP RESULTS TO DISK AT SPECIFIED INTERVALS
IFUGQTQ .1*6, 1) GO TO 951
NUM a 11
CALL FILEOC(NUM)
DO 901 I*1,NOSEG
IRC1*IRECLC1)
DISTQXsC(l,I)*FPARTCl,I)*1000.0
PARTOX(1)=C(1,I)*FPART(2,I)*1000000,0/C(2,I)
TOTTOX a C(1,I)*1000,
WRITE (ll'IRCl) TOTTOX, DISTOX, PARTOX ( 1 )
MUM a 12
CALL FILEOC(NUM)
DO 902 I*1,NQSEG
IRC2sIRECL(2)
WRJTEC12'IPC2}CC2,n,CC3,I),C(4,I)
CONTINUE
RETURN
END
117
-------
APPENDIX 3
Test Program Input and Output
To test that the computer code for the exposure concentration component
of WASTOX is working correctly on the user's computer, a simple lake problem
is provided. Two segments, representing the epilimnion and hypolimnion, are
considered. An inflow and outflow of 1000 cfs occurs in the epilimniom.
7 3
Settling occurs at a flow rate of 10 cfs. An exchange of 10 ft /d is
inputted between the segments but is bypassed in the test run. The chemical
enters at a concentration of 10 yg/& (0.01 tng/J,). Adsorption, volatilization,
and photolysis are considered.
118
-------
INPUT
1
11122
KASTOX TEST RUN EXAMPLE 1-B
TIME VARIABLE
00
1 1
1 2
0,0
10.
*
1 2
1000.
*
1 2
2 2
1.0
1000.
*
2
#
1
*
1
*
1
1
DPTH
TEMP
DPTH
TEMP
KD20
KHN
ML«T
Al
0
1
27
1589.
1193.
577.
80.
0.01
40.
0.0028
0.0020
1
1.0
10.
1
0.01
1
100.
1.0
0
0
10.
20.
20.
20.
1000
0.0
0.00056
0.00016
0.0
TOX
SLD1
10000.
15.
0.1
0.0
100.
1
1 2
1.0
2000.
1.0
0
1
1 0
1.0
0
2
2 1
1.0
1
1.0
1
TDPT
irfVEL
TDPT
WVEL
KDT
KHN
KLT
HI
1673.
1195.
770.
IB.
80.
0.0028
0.0019
0.00043
0.00016
1.0
TOX
SLD1
10000.
360.
*****
*****
0.0 2 *****
*****
*****
*****
1.0 2
1000. 0
0.0
4.47
10.
1.0
1.0
1.0
1710.
1175.
945.
4.
1.
0.0026
0.0018
0.00035
0.00016
0.0
100.
*****
10. o
*****
*****
*****
DSEG 0.
BACK
DSEG 0.
HAC*
KP20
KHN 1 .
ATMS 0.
RH01 2.
*****
*****
1677.
938.
1007.
*****
*****
0.0025
0.00152
0.0003
*****
*****
*****
*****
*****
EXCHANGES *****
MATER COLUMN DISPERSION *****
NO DISPERSION *****
SCALE FACT3R
VOLUMES *****
SCALE FACTOR
FLOWS *****
HYDRONANMIC FLOW *****
1
SCALE
FACTOR
SETTLING OF SOLIDS *****
2
SCALE
FACTOR
BOUNDARY CONDITIONS *****
NO FORCING FUNCTIONS *
SCALE
SCALE
****
FACTOR
FACTOR
NO FORCING FUNCTIONS *****
0 AREA 10000.
BCS1
0 AREA 20000.
UCS1
KPT 1.0
0 WSEG 2.0
0 IFI.K 1.0
65
NO TIME FUNCTION *****
PHOTOLYSIS *****
1575. 1410.
571. 334.
1053. •)32.
QUATUM YIELD *****
PH
PH
KHOH
HNRY
ATMP
1279.
336.
597.
7.0
7.0
0.001
20.
1215.
434.
261.
LATITUDE LONGITUDE SEASON *****
0.0024 0,0023
0.00122 0.0010
0.00026 0.00023
INITIAL CONDITIONS OF
INITIAL CONDITIONS ***
0.0022
0.0021
O.OOOB2 0.00069
0.0002
0.00017
TOXICANT *****
**
STABILITY CRITERIA *****
PRINT INTERVAL *****
TIME STEP, TOTAL TIME
*****
TOT TOX DISSTOX PARTOX
1 1 2
2 1 2
3 1 2
SOLIDS
1 1 2
119
-------
OUTPUT
n
u.
u
u
I"
w
>•
w
H
(J
.J
o
o
o
o
u
I
Ul
a
o
33
I
u. o
4-
UJ
H
10
o
X
O
^
uu
< r>
U. ~^
o
Ul «c
120
-------
u
<
u.
u
o
10
U
II
10
II
U
-3
u
H
H
U
Ul
J
X
U
H
CO
h.
Ill
U
ta
o
K
I*.
X
o
z
o
10
U.
U
o
I
in
o
u.
Id
U.
3
O
I*
u.
o
•x
a
H «
O O
O
~3
U.
O «
a:
u.
o
u
UJ
fl
u
Q
U
V}
o
tx
Q
O •
U3 O
Hi
a.
o
O I
O I
U
' o
J
O I
o
o to
• X
o o
X
u
ac o
O M
[-. l-
JJ W
I-
-------
r- r- o t >4 • • •• O O O
fcl Ul <•< 00 00 + + +
^ tO^ VI a Q U uJ Ui
XX. 000
to «o uy nu ooo
v4 O Ol O CC* <^i *^CM^
X«t X«t <31«•« -• CM ^00 iMOO HUM
M« (O* -*-*+*t-(J*
(j o uj u: j3 tj a. to J
OCtO JH!O •SF 3T t^OOf-»OO i£i/>t.
3 O O 15 Cd^^OOlOO *M
U.U. W U.U. 'jj S-« «-» id O O U3 O O
o to ^> '-o towarooyoo
yj. tO« X >» L9OOUOO O-^O^-»
•zo zo w y> w • • u * • c>ooo
O^ I O "3 *S I U COOOtOOO + + * +
M ice ^» ito aoc UJUJCJw
k, u. ooaix oooo
< 'J < ^] U. U. 3 H II O II U 0000
a: o: tr a: u-u*i*.u* ooo»n
u. < H , tn >. 10 ** i»- xwu;xjju:H'N2:zo
3;?. 2CT -rT-T'«roo[/;a.*:Hx
^5 *tC3 o 'J a,ooxo3^^ O> ^,O^--«
^ *I UJ ^. < 'J H* i r- * _•_•**
Is °3 II 11 00^.0
,. S til • • •• OOOO
<5 S5 =0 oo ***»
^ 7 ^2 oooo
^ S3 - S- S - S£ S2. ?S°§
3S S <« »• £S ££ §§2§
«S „ o «*" u S o.o. o.o.o.o.
CD I tO *
g g H •• « •
S S O Z I- -«
2 O
-------
H
Ul
f.
< o o o o o o«o o^ooooooooooooooooooo
U I I I I I I I !*• I I I I I I I I I I I I I I I I I I
Ul M U) Ul Ul U) Ul UI Ul [21 la) U! U JJ Ul '.*1 '*) U) UJ 'oJ UJ U3 'J oJ 'J *1 'J Ul
O OOOOOOOOOOOOOOOOOOOO3>OOOOOO
OOOOOOOOOOO*NrMOOOOOOOOOOOOOO
Ul 0000'Oi/l^r*>!N^O>»r)r>*O'>*'7'O'»>^C'i)"OO1^-sOOO
z •••••••*•••••••••••••••»•••
•« oooooooooooooooooooo^oooor>o
o
O
;r z oooooooooooooooooooocoooooo o
O O O + + 4> + +. + + + + -*- + + + + + + + .t. + + + + + 4> + + <*> O
_3 t-H UUlUlUlial^lUl^lUUlUlUlUlaJUUlUlUlUlUlLtJUlUUJ'dJUJU •
x M of^'*>opfcino(^in'*>inirtoooooooor-«*>oooooo o
L> cor*«-*r«-r-»-ir^-^oiitT»r»oo-^*^o^J'r"-ot/^OinrNr^—«ooo a*
t-t <-«*4i-i*4r^v4«<4*«t*-*T^o^iDr4co*4^ o
H •• • O
X OOOOOOOOOOOOOOOOOOOOOOOOOOO J
M O
o
X
o
frl
CJ
u r*ocNir>r-of«*ior»o'**oooooooooooooooo u.
j ^ocoo*-*~-*i-«-*4f**f^'***'U>>£r-aD(ytOv4c>i'*'i^'tns£>r*>oD
123
-------
a
*4
H
u
J
a
Ht
H
<
oc
U
o
ct o
jj •
H O
J Z <£
< -. r»>
>
(X Z
Cx] O U-
Z t. l-l
f- U J
Z U 4
II f- !•
ac z a
a. « K
124
-------
ooooooooooooooooooooooooo
OOOOOOOOOOOOOOOOOOOOOOOOO
OOOOOOOOOOOOOOOOOOOOOOOOO
ooooooooooooooooooooooooo
X+-+- + + 4- + ++-4-4- + + + + + 4- + +. + 4-+- + + + +
OCxJUUIilUUtUblUUblUllilEjJUlUUUUJtalUULiJUlU
£^i*«9'«9'*^<*«j'«i<
ooooooooooooooooooooooooo
xooooooooooooooooooooooooo
WJ
X
10
Q
a:
u,
a,
I
ooooooooooooooooooooooooo
XOOOOOOOOOOOOOOOOOOOOOOOOO
Q + I I t I I I I I I I I I I I I I I I I I I I I I
HLUUUUCilUJUldUUUUlLiUJUIUUlUUuJ'alCUUUu}
tOO^^OOO'-<'-l<-t'-«'-''-i'-<»-'«-i'-i'^'-i»-<«-i«-<«-'«-»-^'-«'-<»-'
^^ ^^ ^* ^D OD 00 GO 00 GO GO GO 00 00 00 CD 00 GO GO GO GO GO GO 00 00> 00 GO
MOioaoaoacaoaoooaoooaoacaocoaoaoooaoaoacaoaoaoaoao
oo'*^'**»'^'*«i'^i^>'^'«*'«J'*J'^>'»'<«3t«*'»ii*«*'**'»'<*
................*..••»...
ooooooooooooooooooooooooo
O »-• •-« «•*•*»* «~f v-r «-«<-l«-l«>l«^<-tv-l<-l«-f<^>-l<-««-(v4>-l<^ ffft QQ gift eTrt pft Qp
ooioaoooaoaoaoaDOoaoaoaoacaoooaoaoaoaoaoaoaoGocoao
Ho<<<'4'«i'«9<<9i<'<'^-<9'^>4<«i<^r^<<«<«r<»<*'4'
ooooooooooooooooooooooooo
Ul
ooooooooooooooooooooooooo
125
-------
"O no r< cs CM rs rM :M M JN .-N :N ^s --N r>» rx CM CM ,-N CN tN CN rs rN .-><
OOOOOOOOOOOOOOOOOOOOOOOOO
UI
f-l
00
>*
U5
o
X
u.
ooooooooooooooooooooooooo
ooooooooooooooooooooooooo
(O ^h ^ ^H ^~ "^ ^ ^K 4* ^ ^" ^* ^ ^h ^ ^ 4* ^" *^ ^* ^ ^ 't* *^ ^~ ^*
^l '»^ LO t^ «J ^3 m tiJ U *T-^ CiJ -i3 f«j '«^ r^i f^ [j ^j f|j r^j c^^ 'jj £j ^O ro (^3
JOPJOOOOOOOOOOOOOOOOOOOOOOO
ooooooooooooooooooooooooo
ooooooooooooooooooooooooo
126
------- |