EPA
United States
Environmental Protection
Agency
Environmental Research
Laboratory
Athens GA 30613
EPA-600 3-83-005
March 1983
Research and Development
User's Manual for the
Chemical Transport and
Fate Model (TOXIWASP),
Version 1
-------
EPA-600/3-83-005
USER'S MANUAL FOR THE CHEMICAL TRANSPORT AND FATE MODEL TOXIWASP
VERSION 1
by
Robert B. Ambrose, Jr., Sam I. Hill, and Lee A. Mulkey
Environmental Research Laboratory
U.S. Environmental Protection Agency
Athens, Georgia 30613
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GEORGIA 30613
-------
DISCLAIMER
This document has been reviewed in accordance with the U.S. Environmental
Protection Agency's peer and administrative review policies and approved for
publication. Mention of trade names or commercial products constitute
endorsement or recommendation for use.
The TOXIWASP computer code has been tested against other computer
programs to verify its computational accuracy. Nevertheless, errors in the
code are possible. The U.S. Environmental Protection Agency assumes no
liability for either misuse of the model or for errors in the code. The user
should perform verification checks of the code before using it.
ii
-------
FOREWORD
As environmental controls become more costly to implement and the
penalties of judgment errors become more severe, environmental quality
management requires more efficient analytical tools based on greater
knowledge of the environmental phenomena to be managed. As part of this
Laboratory's research on the occurrence, movement, transformation, impact,
and control of environmental contaminants, the Technology Development and
Applications Branch develops management or engineering tools to help pollu-
tion control officials achieve water quality goals.
Concern about environmental exposure to synthetic organic compounds
has increased the need for techniques to predict the behavior of chemicals
entering natural waters as a result of the manufacture, use, and disposal
of commercial products. A number of mathematical models have been developed
to provide data on aquatic pollutant behavior from which exposure assess-
ments can be made. The modeling technique described in this manual permits
the user to examine the probability of exposure of receptor populations to
a range of chemical concentrations for a specified location in stratified
lakes and reservoirs, large rivers, estuaries, and coastal waters.
David W. Duttweiler
Director
Environmental Research Laboratory
Athens, Georgia
111
-------
ABSTRACT
This manual describes the dynamic model TOXIWASP for simulating the
transport and fate of toxic chemicals in water bodies. TOXIWASP combines the
kinetic structure adapted from the Exposure Analysis Modeling System (EXAMS)
with the transport framework provided by the Water Analysis Simulation Program
(WASP), along with simple sediment balance algorithms. TOXIWASP formulates
variable chemical degradation rates from chemical characteristics of a
compound and the environmental parameters of the aquatic system. These rates
combine calculated first order rates due to several processes, including
hydrolysis, biolysis, photolysis, oxidation, and volatilization. Sorption
onto sediments and onto biomass is calculated assuming local equilibrium,
using a chemical specific partition coefficient, and spatially varying
environmental organic carbon fractions.
TOXIWASP calculates total sediment and chemical concentrations explicitly
every time step for every segment, including surface water, subsurface water,
surface bed, and subsurface bed. Sediment concentrations are affected by
advection, dispersion, mass loading, settling, scour, and in the bed, by
burial and erosion. Chemical concentrations are affected by these same
processes, plus degradation, sediment-water dispersion, and percolation. No
lateral migration of chemical within the bed is allowed. Transport data must
be specified based on monitoring data or predictions from hydrodynamic
models.
TOXIWASP is designed to provide WASP users the capability for simple
dynamic simulations of chemicals, particularly pesticides, which enter the
aquatic environment in pulses that cannot be simulated with steady state
models. It can be used for cases requiring more dynamic transport and loading
capabilities than EXAMS, but less detailed and mechanistic sediment
predictions than SERATRA. Its toxic chemical simulation capabilities are
roughly equivalent to those in the Hydrologic Simulation Program in FORTRAN
(HSPF). Although it is less applicable to upland streams and one-dimensional
reservoirs than HSPF, TOXIWASP is more suited to stratified lakes and
reservoirs, large rivers, estuaries, and coastal waters.
This manual gives an overview of the theory pertaining to TOXIWASP, then
describes the program structure. Detailed data input instructions are given,
followed by operating instructions and a glossary.
This report covers the period March 1, 1981 to November 1, 1982, and work
was completed as of November 1, 1982.
iv
-------
CONTENTS
Disclaimer ii
Foreword iii
Abstract iv
Figures vii
Tables viii
Acknowledgments ix
1. Introduction 1
1.1 Background 1
1.2 Exposure Assessment 2
1.3 Overview of this Manual 3
2. Model Capability and Theory 5
2.1 Overview of TOXIWASP 5
2.2 Sediment Transport Processes 8
2.3 Chemical Partitioning (Sorption) Processes 10
2.4 Bed-Water Chemical Redistribution Processes 13
2.5 Bed Sedimentation and Erosion 17
2.6 Transformation Processes .... 22
2.7 Special Data Input Options 30
2.8 Special Simulation Control Options 33
2.9 Special Output Options 33
2.10 Summary of TOXIWASP Capabilities 34
3. Model Mainline and Subroutine Overview 37
3.1 Standard WASP Subroutines 37
3.2 TOXIWASP Kinetic Subroutines 43
4. Preparation of Data Input 49
4.1 Card Group A - Model Identification 52
4.2 Card Group B - Exchange Coefficients 53
4.3 Card Group C - Volumes 55
4.4 Card Group D - Flows 56
4.5 Card Group E - Boundary Conditions 59
4.6 Card Group F - Forcing Functions 62
4.7 Card Group G - Parameters 65
4.8 Card Group H - Constants 67
-------
CONTENTS (Continued)
4.9 Card Group I - Miscellaneous Time Functions 71
4.10 Card Group J - Initial Conditions , 73
4.11 Card Group K - Stability and Accuracy Criteria 73
4.12 Card Group L - Intermediate Print Control 74
4.13 Card Group M - Integration Control Information 74
4.14 Card Group N - Display Parameters 75
5, Operational Considerations 83
5.1 Acqisition Procedures 83
5.2 Installation Procedures 83
5.3 Testing Procedures 83
5.4 Machine Limitations 84
References 85
Appendix 1, Glossary 87
Appendix 2, TOXIWASP Subroutine Listings 119
Appendix 3, TOXIWASP Support Software 158
3.1 Listing of TOXICMP.CMD 158
3.2 Listing of TOXITKB.CMD 159
3.3 Listing of JCL to Compile and Link Edit 161
3.4 Listing of TOXIC.BIS 161
3.5 Listing of JCL to Run Simulations 162
3.6 Listing of DAILYCVT 164
3.7 Listing of HSPFCVT 167
3.8 Listing of SWRVCVT 169
3.9 Listing of USGSCVT 171
Appendix 4, Listing of Input Data File 175
vi
-------
FIGURES
Number Page
2-1 Examples of TOXIWASP Network Configurations 6
2-2 TOXIWASP Sediment Transport Processes 11
2-3 TOXIWASP Chemical Transport and Redistribution Processes . 15
2-4 TOXIWASP Sediment Burial 19
2-5 TOXIWASP Sediment Erosion 21
3-1 WASP Program Structure 38
3-2 Simplified WAS15-WASPB Flow Chart 42
3-3 TOXIWASPB Kinetic Subroutine Structure 44
vii
-------
TABLES
Number Page
2-1 Concentration Related Symbols Used in Mathematical
' Equations 9
2-2 Structure of Auxiliary Flow and Nonpoint Source Load
File 'DA1LY.DAT' 32
4-1 Summary of Card Groups 50
4-2 Pen Plot Example 82
viii
-------
ACKNOWLEDGMENTS
Staff at the Athens Environmental Research Laboratory provided critical
support for this project. Lawrence Burns provided useful discussions while
we were adapting EXAMS kinetics to the WASP framework, and later in testing
TOXIWASP against EXAMS (the final product, including any mistakes are, of
course, our responsibility). Bruce Bartell, with Computer Science Corpora-
tion, drafted the figures by computer. Annie Smith typed quick, accurate
drafts. Robert Ryans provided our technical editing. Sandra Ashe assisted
in the tedious typing, copying, and other preparation chores.
ix
-------
SECTION 1
INTRODUCTION
1.1 BACKGROUND
TOXIWASP is a dynamic model for simulating the transport and fate of
toxic chemicals in water bodies. Two state variables are simulated:
organic chemical and total sediment. TOXIWASP was created by adapting the
kinetic structure of the Exposure Analysis Modeling System (EXAMS) (Burns et
al., 1982) to the transport framework provided by the Water Analysis Simula-
tion Program (WASP) (DiToro et al., 1981), and adding simple sediment
balance algorithms along with special input and output software. The result
is a generalized chemical model that can be used for cases requiring more
dynamic transport and loading capabilities than EXAMS, but less detailed and
mechanistic sediment predictions than SERATRA (Onishi and Wise, 1982). The
TOXIWASP chemical simulation capabilities are roughly compatible with those
in the Hydrologic Simulation Program—FORTRAN (HSPF) (Johanson et al.,
1980). HSPF is most applicable to upland streams and one-dimensional
reservoirs, whereas TOXIWASP is more suited to stratified lakes and reser-
voirs, large rivers, estuaries, and coastal waters.
The main component of TOXIWASP is WASP, a dynamic compartment modeling
program for aquatic systems. Time-varying transport, including forcing and
boundary processes, advection, dispersion, and mass loading, are represented
in this program. Water quality processes are represented in a special kine-
tic subroutine that is either chosen from a library or written by the user.
WASP is structured to permit easy substitution of kinetic subroutines into
the overall package to form a particular model. The model provided a sound
basis for the rapid development of the TOXIWASP methodology. Indeed, ver-
sions of WASP have been used to examine PCB problems in the Great Lakes
(Thomann, 1982; Richardson, et al., 1982).
The kinetic component of TOXIWASP was modified from EXAMS. From chemi-
cal characteristics of a compound and the environmental parameters of the
system, EXAMS formulates a total transformation rate. This rate is based on
a simple addition of the pseudo-first order rate constants due to each pro-
cess. In TOXIWASP a kinetic derivative is calculated from this rate and
passed to standard WASP subroutines, where it is integrated along with
transport and loading derivatives to yield time-varying chemical concentra-
tions for a specified spatial network. EXAMS uses a sophisticated kinetic
structure that allows the study of five different ionic forms of a chemical,
several ways to calculate photolysis, and other capabilities. It also has
the user advantage of being interactive. In these respects, EXAMS was
simplified and reduced for use in the TOXIWASP application. TOXIWASP pro-
1
-------
vides the users of WASP with an expanded library of kinetic subroutines
employing state-of-the-art chemical transformation techniques.
TOXIWASP allows the user to perform simple dynamic simulations of
potentially toxic organic chemicals, particularly pesticides, that enter
the aquatic environment in pulses that cannot be simulated with steady-
state models. Dynamic simulations allow the consideration of pulse loads,
the prediction of peak events, and an estimate of time-varying chemical
exposure. TOXIWASP was created for use in exposure assessment.
1.2 EXPOSURE ASSESSMENT
Evidence of potentially harmful effects of pesticides and other toxic
organic chemicals on aquatic organisms has led to intensive efforts toward
environmental risk assessment for existing and new chemicals. The concept
of risk reflects the likelihood or probability of causing such effects.
Obviously, if an effect is induced, the organism affected must first have
been exposed to the toxicant for a time period and intensity sufficient to
inflict the damage. The study of this exposure, or exposure assessment, is
defined as a quantitative evaluation of the concentration of chemical toxi-
cants in contact with receptor populations within the various environmental
media as the toxicant is released, transported, and transformed into,
through and between environmental compartments.
The source of toxic chemicals in a river basin is an important factor
in selection of techniques to analyze the problem. Toxic chemical releases
may arise from both nonpoint sources such as agricultural fields from which
pesticides are transported in runoff and point sources such as chemical
manufacturing plants from which chemicals are discharged in effluent. Non-
point sources are characterized by highly variable loadings with rainfall-
runoff events often dominating the timing and magnitude of chemical inputs
to the receiving waters. Point sources are usually much less varied and are
often thought to provide steady inputs to the aquatic system. Clearly,
exposure assessments for toxic chemicals in the aquatic environment must
accommodate both the nonpoint or variable loading case and the point or
steady loading case.
Chemical concentrations in aquatic environments are also determined by
the transport and transformation processes that occur in river basins. As
chemicals are transported by water or sediment, they are subject to a number
of physical, chemical, and biological processes that can lead to either
transport or transformation. The influence of transport on cheimical beha-
vior is determined by the local environmental regime (flow rates, sediment
concentrations, etc.) and such chemical properties as the partition
coefficient. Transformation processes are determined by chemical properties
and by such environmental conditions as pH and temperature.
Exposure assessment as a part of risk evaluation implies the need to
estimate probability of exposure. Exposure can be viewed as the realization
of the joint probability of an ambient concentration and the existence of
-------
receptor populations. The probability of ambient concentrations can be
approximated by calculating their frequency of occurrence over a sufficient
period of time. A convenient way to represent frequency is the cumulative
distribution function (cdf). A cdf can be expressed as a graph of the
concentrations (on the abcissa) exceeded for various fractions of a time
period (on the ordinate). A cdf can be related to average expected return
periods for different concentrations. Given the cdf for a specific chemical
in a specific environment, one can, for example, say that once every 10
years the concentration of a specific chemical will exceed a specified
level. The concentration level can then be related to biological effects
information (e.g., LC5Q values) and a return interval for damage (or risk)
results. A cdf can readily be constructed given sufficient data in the form
of time-series chemical concentrations. Rarely are any such data available,
however, and when short time periods are sampled, the statistical reliabili-
ty of the resulting data is always suspect. An alternative to field sampling
is mathematical simulation with a dynamic model such as TOXIWASP.
If the user wishes, TOXIWASP produces a time-series of chemical
concentrations at a segment that reflects daily changes in flow and loading,
and seasonal changes in water temperature, light intensity, wind speed, pH,
and POH. TOXIWASP cannot simulate stochastic variability or variability be-
low the resolution of its time step, which is usually between one and twelve
hours. Nevertheless, the cdf constructed by TOXIWASP for a specified loca-
tion reflects the probability of exposure of local receptor populations to a
range of chemical concentrations. For populations that can avoid locally
high chemical concentrations, the problem of estimating exposure is more
difficult. The spatial variability of peak concentrations must be examined
to estimate the probability of avoidance. While TOXIWASP can provide the
spatial variability during peak concentration events, users must convert
this information to their own estimates of avoidance.
1.3 OVERVIEW OF THIS MANUAL
This manual documents the computer program TOXIWASP. It is meant to
provide a convenient reference for the users of TOXIWASP. It is not meant
to provide a tutorial on chemical simulation modeling. Nevertheless, an
example input data file Is provided in Appendix 4, and on the TOXIWASP dis-
tribution tape.
Following this introduction is a chapter on model capability and
theory. This gives a summary of the equations solved in TOXIWASP, and short
discussions of the sediment transport processes, the partitioning processes,
the redistribution processes, the transformation processes, special data in-
put options, control options, and output options.
Chapter 3 is an overview of the WASP and TOXIWASP subroutines. The
descriptions of WASP subroutines are only slightly modified from the
original WASP users manual.
Chapter 4 details the preparation of data input files. Data are
organized into "card groups" A through N in the same format as the original
-------
WASP users manual. Of course, "card group" can refer to a series of data
records on a computer file.
Chapter 5 discusses such operational considerations as how to acquire
TOXIWASP, how to install it on a computer and test it, and what its machine
limitations are.
The first appendix is a glossary of all variables in TOXIWASP. This
lists the type of variable (local, parameter, constant, time function), the
subroutines where it is used, its definition, and its units. Also in the
glossary is a list of abbreviations for commonly used units.
The second appendix gives source listings of the TOXIWASP subroutines.
The third appendix gives listings of TOXIWASP support software. The fourth
appendix gives a sample input dataset. All of these listings are available
as files on the distribution tape.
A note on units may be helpful. WASP uses mixed English and SI units.
For the general process equations in this document, we label dimensions
following the M.L.T. system where M stands for mass, L stands for length,
and T stands for time. Molarity is labelled by M, and moles by n. For
specific equations used in the TOXIWASP computer program, we label terms
with specific units. CHEMW, for example, has units of mgc/lw,, or milligrams
of chemical per liter of water. These specific units are defined in the
glossary, starting on page 87.
-------
SECTION 2
MODEL CAPABILITIES AND THEORY
All users of TOXIWASP should be familiar with and have access to the
WASP user's manual. WASP theory will only be considered where necessary to
explain differences between EXAMS and WASP. This section is intended to
provide an overview of the capabilities provided by the new toxics sub-
routines in TOXIWASP. This section is not intended to supplant the detailed
documentation of kinetic process theory available in the EXAMS manual. That
documentation is too large and exhaustive to be included here. The user who
has questions is referred to this excellent source.
2.1 OVERVIEW OF TOXIWASP
TOXIWASP is a dynamic chemical and sediment model that can be applied
to water quality problems in streams, lakes, reservoirs, estuaries, and
coastal waters. It uses the compartment modeling approach whereby segments
can be arranged in a one-, two-, or three-dimensional configuration, as in
Figure 2-1. Pollutant transport is based on user-specified flow and disper-
sive mixing between segments. All parameters and forcing functions may vary
in space and time. TOXIWASP calculates time-varying water quality concen-
trations using an explicit, backward difference numerical solution to the
mass flux form of the one dimensional advective-dispersive equation:
AM-i n AC
1 = E [Q. .C. + E. .A. .(—)..] + W. - KV. /n
At i=1 L^ID i ID 13 L 'ID D : u'
where M = constituent mass
C = constituent concentration, ML~~3
Q = water flow, L^T"!
E = longitudinal dispersion. L2T-1
A = cross-sectional area, L^
L = characteristic mixing length, L
W = mass loading, MT~1
K = kinetic degradation or transformation rate, ML~3x~l
V = segment volume L^
j = segment number
i = adjacent segments
ij = interface between segment j and adjacent segments i
TOXIWASP was created from WASP by adding a set of subroutines that
calculate chemical degradation (the last term in equation 1) and that
-------
Figure 2-1
Examples of TOXIWASP Network Configurations
Chemical and Sediment Loads
Flow
Chemical and
Sediment Loads
Flow
Chemical and
Sediment Loads
i
Flow
Pond
Advective Flow
Dispersive liixing
Segment Type
Surface Water
Subsurface Water
Surface Bed
Subsurface Bed
Pattern Scheme
Lake
-------
calculate special bed sediment transport processes such as sediment-water
dispersive exchange, pore water percolation, and net burial-scour. The two
differential equations solved by TOXIWASP for water-borne chemical and
sediment, written in the general concentration form, are:
3C1 3C1 9 ,P9C^ Wl K + q
TT = u 3T" + a£ (^T') ~V ~ 1 (2a)
r*. C* r\ /^ f"11 J^ f* 1 Uo *
fl ^V ^ *^X ^ ^^ ^
7 ~ H + ( E ' ) + r? ~^~ b y- rt i \
o
where C^ = concentration of chemical, ML
G£ = concentration of suspended sediment, ML
u = flow velocity of water, LT~1
W, = mass loading of chemical, MT""^
W£ = mass loading of sediment, MT
S, = net exchange of chemical with bed, ML"3 T~*
•
S2 - net exchange of sediment with bed, ML"3 T"1
x = longitudinal distance, L
t = time, T
Ms = mass of sediment
From chemical characteristics of a compound and the environmental
parameters of the system, TOXIWASP formulates a total transformation rate.
This rate is based on a simple addition of the pseudo-first order rates for
hydrolysis, oxidation, biolysis, volatilization, photolysis, and other
processes. Sorption onto sediments and onto biomass is calculated assuming
local equilibrium, using a chemical specific partition coefficient, and
spatially varying environmental organic carbon fractions.
TOXIWASP applies equations 2a and 2b, in the form of equation 1 to
calculate sediment and chemical concentrations for every segment, including
surface water, subsurface water, surface bed, and subsurface bed. Sediment
is assumed to be a conservative constituent that is advected and dispersed
among water segments and that can settle to and scour from bed segments.
In a simulation, the chemical undergoes first-order decay, based on
summation of several process rates, some of which are second-order. The
effective first order decay rate can vary with time, then, and is recalcu-
lated as often as necessary throughout a simulation. The chemical is ad-
vected and dispersed among water segments and can be transported to the bed
by sorption onto sediment and subsequent settling, or by dispersive mixing
between the top bed layer and the bottom water layer. Within the bed, the
chemical can migrate downward or upward through percolation and pore water
diffusion. No lateral migration of the chemical within the bed is allowed.
-------
Because the kinetic transformation processes are all expressed in terms
of first-order rates, large loadings may cause problems by changing the
characteristics of the systems. An example noted in the EXAMS manual is a
chemical subject to alkaline hydrolysis that kills the biological component,
making the system less alkaline, and thus reducing the alkaline hydrolysis
degradation rate. It is the responsibility of the user to anticipate these
problems and look for cases where they might happen. The task of keeping
track of all possible ecosystem variables is too complex, and would occupy
too much computer space for TOXIWASP. Therefore, feedback phenomena are not
considered.
TOXIWASP, like EXAMS, is based on the following assumptions.
1. All segments are well mixed.
2. Sorption is an essentially instantaneous process within each
segment (thus local equilibrium).
3. Chemical properties of the compound can be coupled with the
characteristics of the environment to formulate a pseudo-first-
order rate for each of the degradation processes using primarily
mechanistic equations.
4. Pseudo-first-order rates combine linearly to form a total
degradation rate for the chemical toxicant in each segment.
It should be noted that TOXIWASP, like WASP, requires the user to
specify the flow field. Flows can be based on measurements, simple con-
tinuity calculations, or on other hydrodynataic models. WASP does not check
the flow field for inconsistencies that may lead to mass balance errors.
The user should take care to check the specified flows for errors.
In the following development it is convenient to define concentration
related symbols as in Table 2-1.
2.2 SEDIMENT TRANSPORT PROCESSES
Sediment transport is an important process in hydraulics and river
mechanics, but its role in pollutant transport is not fully known. Several
investigators (Karickhoff et al., 1979; Rao and Davidson, 1980) have shown
that for most organic toxicants, particularly non-polar compounds, the sedi-
ment organic matter determines the extent of sorption and hence the impor-
tance of sorbed-pollutant transport. The particle size of the sediments is
influential, but Rao and Davidson (1980) have shown that the organic matter
content of each size fraction is the controlling factor. Organic matter
content and surface area relationships suggest that silt and clay size sedi-
ment fractions are most relevant in sediment transport. One outcome of this
rationale is the idea that, in river systems, wash load or suspended sedi-
ment is the important pollutant transport medium rather than bed load. If
bed load transport can be ignored, the problem remains to estimate suspended
sediment. In general, the stream transport capacity for wash load is in
excess, and the problem is one of estimating sediment source loading—
namely, watershed erosion. Thus, the key sediment-related term In equation
2b is the source term, ^2•
-------
Table 2.1. Concentration Related Symbols
Used In Mathematical Equations
Symbol
c,q
Cw
Cw
Cs
t
Cs
CP
t
cp
Cb
1
Cb
C2
s
1
s
sw
sb
0
IT
Definition
Concentration of chemical (bed or water).
Concentration of dissolved chemical in water.
Concentration of dissolved chemical in water.
Cw = V
Concentration of sorbed chemical in water.
Concentration of sorbed chemical in water.
Cs ' Cs/Sw
Concentration of dissolved chemical in bed (pore water).
Concentration of dissolved chemical in bed (pore water).
CP - V0
Concentration of sorbed chemical in bed.
i
Concentration of sorbed chemical in bed. Cb = Cft/S^
Concentration of sediment (bed or water).
Concentration of sediment (bed or water). S = Co^O
Concentration of sediment (bed or water). S = S/0
Concentration of sediment in water.
Concentration of sediment in bed.
Porosity or volume water per volume segment.
Partition coefficient of chemical on sediment.
Units
mgc/lT
mgc/lT
mgc/lw
mgc/lT
mgc/kgs
tngc/lf
mgc/lw
mgc/lT
mgc/kgs
rngs/lf
kgs/lT
kgs/lw
kgs/lT
kgs/lj
IW/IT
lw/kgs
-------
TOXIWASP treats sediment as a single size fraction that is advected and
dispersed in the water and that settles (and deposits) at a constant,
spatially-variable velocity. Sediment is returned from the bed to the water
column through a spatially-variable scour (or erosion) velocity. Thus the
net exchange of sediment with the bed, 82 in Equation 2b, is calculated as
WsSb "dsw
S2 ~ Lb Lw (3)
where W = scour (erosion) velocity of bed sediment, LT~*
—» 1
Wj = deposition velocity of suspended sediment, LT
Lb = depth of active bed sediment layer, L
Lw = depth of water layer, L
•
Nonzero values for 82 will lead to net sedimentation or erosion,, The
consequences of chemical burial or release are discussed in Section 2.5.
Settling and scour velocities are read in using the parameter WS(J).
Sediment transport processes are illustrated in Figure 2-2.
2.3 CHEMICAL PARTITIONING (SORPTION) PROCESSES
The interaction of dissolved organic toxicants with entrained sediment is
an important process that must be accounted for mathematically. If one
assumes sorption is kinetically a first order process, then a rate equation
similar to the other expressions can be written as:
Rs _
~s = k C ~ k. c r /. \
S s w d s (4)
where Rs = net rate of chemical transfer between dissolved and sorbed state,
ML
"3"1
!" (typically mgc/lT~hour)_
3~
_
k = rate constant for sorption, LM~ T (typically lw/kgs-hour)
kj = rate constant for desorption, T (typically per hour)
If one assumes Rs = 0 then
or _ ks r, (5)
(^ _ — _ _ ^
s kd w
Note that equation 5 is the linear form of the Freundlich isotherm with the
ratio ks/k
-------
Figure 2-2
TOXIWASP Sediment Transport Processes
r
V*
±*
Process
a) Mass Loading
b) Advection
c) Dispersion
d) Settling
e) Erosion
Segment Types
1,2 to Water
1,2 Water-Water
1,2 Water-Water
1,2 Water-Water or Water-Bed
3 Surface Bed-Water
11
-------
empirical expressions relating equilibrium coefficients to laboratory
measurements leading to fairly reliable means of estimating appropriate
values.
EXAMS assumes sorption and ionization to be essentially instantaneous
and considers them to be equilibrium processes. The model allows for three
sorption possibilities (dissolved, sediment sorbed, and biota sorbed) and
for four ionic species in addition to the unionized form. To conserve com-
puter space, TOXIWASP only considers the three sorption possibilities for
the unionized form of the chemical. As EXAMS notes, not all possible forms
of a chemical will exist. The specific fractions (1,2,3) of the chemical in
each form, ALPHA, are calculated for each compartment based on 12 segment
parameters, constants, and state variables:
1
ALPHA(l) = 1 + KP * S/FRW + KB * B/FRW = «1 (6)
_ KP * S/FRW _
ALPHA(2) = 1 + KP * S/FRW + KB * B/FRW = «2 (?)
_ KB * B/FRW _
ALPHA(3) = 1 + KP * S/FRW + KB * B/FRW = ot3 (8)
where KP = partition coefficient of the chemical on the sediment
(lw/kgs) = KOC * OCS(J) - IT
KB = partition coefficient of the chemical with biota (lw/kgb)
= KOC * OCB
S = concentration of sediment, kgB/l^
B = concentration of biomass, kgb/l-p
FRW = porosity of segment (water volume/total volume) -
-------
where CHEMW = chemical dissolved in water phase of segment, mgc/lw
CHEMS = chemical sorbed onto sediment phase of segment, mgc/kgs
CHEMB = chemical sorbed onto biological phase of segment, ugc/g|,iomass
(same as mgc/kgb)
2.4 BED-WATER CHEMICAL REDISTRIBUTION PROCESSES
An important, but poorly quantified, set of processes governs chemical
exchange between the bed and the water column. A chemical may sorb onto
suspended sediments and settle to the bed, or sorb directly to the bed sur-
when water concentrations are high. When water concentrations are low, the
chemical may desorb from the bed surface. Turnover of the bed surface due
to physical or biological processes should accelerate this sorption-desorp-
tion process. Diffusion of pore water into overlying water will also ex-
change dissolved chemicals. Percolation of pore water will move dissolved
chemical up through the bed into the water column; infiltration can leach
chemical down through the bed.
So even though bed load transport per se may be relatively unimportant
in pollutant transport, interaction between the stream bed and the overlying
water should not be ignored. Concentration gradients between the bed mate-
rial and the overlying water will result in sorption and desorption and thus
the stream bed can act either as a sink or a source of contaminants. For
the stream bed one can write the dissolved and particulate chemical conser-
vation equations:
s
wscb wdcs
- R ' K (16a)
___ .
3t ~ 3y 3y ' L8 LW s (16b)
where D^ = vertical dispersion coefficient for dissolved chemical,
I)bs = vertical dispersion coefficient for sorbed chemical, L^T~1
Up = velocity of net pore water movement into or out of the bed,
LT-1
Wj , Ws = deposition and scour velocity of sediment between bed and
water column, LT~1
Lw = depth of water compartment, L
LS = depth of active bed layer, L
y = vertical distance, L
Equation 16a accounts for diffusion-dispersion and pore water transport
of pollutants between the bed and the overlying water column. Equation 16b
accounts for sediment-bound transport by including the scour and deposition
terms, Wg and Wj. Estimation usually requires a rigorous treatment of
hydrodynamic forces and characterization of sediments (e.g., cohesiveness) .
Several approaches are available for this purpose, but many chemical expo-
sure models either omit these processes entirely or lump them into an over-
all calibrated dispersion or exchange coefficient.
13
-------
Comparing equations 16a and 16b with equation 2a reveals that S^, the
net exchange of chemical from the bed to the water column, is a summation of
several processes. These include pore water diffusion and percolation, bed
sediment dispersion (due to physical mixing and bioturbation), scour and
deposition, and direct sorption-desorption.
In TOX1WASP, the subroutines TOX1SETL and TOX1SEDW handle the redistri-
bution terms in equations 2a, 2b, 16a, and 16b. The processes are illus-
trated in Figure 2-3. The settling velocity W
-------
Figure 2-3
TOXIWASP Chemical Transport and Redistribution
Processes
Process
a) Mass Loading
b) Advection
c) Dispersion
d) Settling
e) Erosion
f) Pore Water Diffusion
g) Sediment Turnover
h) Percolation
i) Sedimentation
j) Volatilization
Segment Types
1,2 to Water
1,2 Water-Water
1,2 Water-Water
1,2 Water-Water or Water-Bed
3 Surface Bed-Water
3,4 Surface Bed-Water or Bed-Bed
3 Surface Bed-Water
4 Bed-Bed or Bed-Water
4 Bed-Bed
1 Surface Water
15
-------
The volumetric dispersion varies spatially, and is calculated by:
SCALE * Db * AS * FRW
DISPV= Li U8a)
where D^ = vertical pore water dispersion coefficient,
A = surface area, L
FRW = porosity, or volumetric fraction of pore water in the bed
Lc = characteristic mixing length, L
SCALE = units conversion factor
Because TOXIWASP and EXAMS assume local equilibrium between dissolved and
sorbed phases of a chemical in each segment, equation 17b can be rewritten:
EXCH = DISPV * 1 + gg * sl * [c; - KPB * c;i (18b)
where DSPSED is set equal to 1.0.
Equation 17a is in the same form as equation 18b, and so the HSPF-SERATRA
sorption-desorption rate can be equated with the EXAMS-TOX1WASP dispersion
rate by equating RSED * KT with DISPV * 1 + KPB * S* and solving:
1 KPB
SCALB * D * As 1 + KPB * S'
KT = BVOL * Lc * KPB * Sf (19)
where BVOL = EXAMS-TOX1WASP bed volume, L3
SCALB = unit conversion factor (=2.584 x 10~4 when E = m2/hr, A =
ft2, Lc = ft, BVOL = ft3 x 106).
RSED = BVOL * S' * FRW
For chemicals with high partition coefficients, the last term in equation 22
is approximately 1.0, simplifying the expression.
Although the basic sediment-water exchange mechanisms are mathemati-
cally equivalent, each model has its own elaborations. The exchange rates
are spatially variable in EXAMS and TOXIWASP, and these models include per-
colation of pore water. TOXIWASP adds spatially variable sediment turnover
ratios, spatially variable net bed scour rates, and settling of suspended
sediment. HSPF and SERATRA simulate three sediment size fractions. Scour-
deposition and settling are calculated as dynamic processes. SERATRA adds
multiple sediment layers, which can be added through deposition or sub-
tracted through erosion. HSPF allows spatially variable sorption-desorpLion
rates. All these mathematical formulations are limited by the lack of good,
thorough data necessary to validate the models on a site-specific basis.
TOXISEDW was inserted to calculate sediment-water exchange in a manner
similar to EXAMS. Bedload-water column interaction is a complex and only
partially understood problem. TOXIWASP handles the interaction in the
16
-------
following manner. When a fraction of the bed sediment compartment (sediment
plus pore water) is mixed into the water column, the sediment immediately
comes into equilibrium with the chemical in the overlying water compartment,
then is returned to the bed. The pore water is exchanged with an equal
amount of overlying water. The rate of exchange is based on D1SPV, the
volume exchange rate. DISPV is calculated for each sediment-water interface
"I" from the dispersion rate, BR(I), which is calculated from parameters in
Card Group B:
BR(I) = E(I)*A(I)/[(IL(I)+JL(I))/2.)*(5280**2)*1.E-06] (20)
DISPV(I) = CONVERSION FACTOR * BR(I) * FRW (21)
where 'BR(I) = volumetric dispersion rate for segment interface "I",
(MCF/day for water-water interface, mixed units for
sediment-water interface)
E(I) = dispersion coefficient for segment interface "I", L^T~1
(mi^/day for water-water interface, cm^/sec for
sediment-water interface)
A(l) = surface area of segment interface "I", L^ (ft2)
IL(I),JL(I) = characteristic mixing lengths for segments with respect to
interface "I", L (ft)
DISPV(I) = volumetric pore water dispersion for sediment-water interface
"I", L3T-1 (MCFw/day)
BR(I) is calculated from equation 20 in WASP2. For sediment-water inter-
faces, BR(I) is converted to DISPV(I) by equation 21 in TOXISEDW, then set
to 0. Actual sediment-water mass exchange is then calculated each time step
by equation 17b. The parameter DSPSED allows the user to control the rela-
tionship between pore water diffusion and sediment turnover. If DSPSED is
set to 0.0 for a segment, only pore water diffusion occurs. If DSPSED is
set to 1.0, the sediment turnover volume is equal to the pore water diffu-
sion volume.
2.5 BED SEDIMENTATION AND EROSION
The chemical redistribution processes in Equations 16 and 17 and Figure
2-3 operate in reference to a fixed bed. When sediment deposition (W^SW/LW)
exceeds sediment scour (WsSb/Ls), however, net deposition can leads to a
rising bed surface and burial of sorbed chemical. When sediment scour
exceeds sediment deposition, net erosion can lead to a falling bed surface
and release of previously buried sorbed chemical. In this section, the term
"sedimentation" will refer to net sedimentation or erosion.
A simple method for treating sedimentation is to define bed segment
locations in reference to the rising or falling bed surface. If the bed
surface rises at W cm per year, then buried chemical descends through the
bed at W cm per year in reference to the rising bed segments (all other bed
dispersion processes being ignored for the moment). There are two
complications in this approach. First, if the density of the bed increases
17
-------
with depth, then compression must be occuring, decreasing W and squeezing the
pore water and dissolved chemical upward. Second, the simulation time step
is extremely short in relation to the net sedimentation velocity. This
leads to rather severe "numerical dispersion," a form of computational
inaccuracy with in the bed.
To handle the first complication, TOXIWASP assumes two types of bed
sediment: an upper uncompacted layer with bulk density 85, and lower,
compacted layers with bulk density S£. Assuming that sedimentation is in
equilibrium with net deposition, all bed properties remain constant. The
sediment mass addition is:
' Wdsw WsSb
M= V ~ TT* vs (22)
where M = sediment mass added to bed, MT~1
Vs = volume of the upper layer, L , equal to AgLg
A = sediment surface area, L
The volume of sediment added to the upper layer is M/S^, while the
volume added to the lower layer is M/S. The difference in volumes is equal
to the pore water volume squeezed out through compaction. The rise in the
bed surface equals the volume added to the lower layer divided by the surface
area:
_ M li ^* ^
W - ASS£ - Wd Lw S£ - Ws S£ (23)
The amount of chemical lost from the upper bed layer due to sedimentation is:
v |£ = _WA S c' - WA -1 - 1) C1 (24)
s a t s&b sSfc P
where the first term represents sorbed chemical buried downward and the
second term represents dissolved chemical squeezed upward.
TOXIWASP implements a version of this simple approach that amounts to a
Lagrangian scheme for treating movement of chemical through the bed due to
sedimentation. The time step for sedimentation is calculated internally,
based on the scour and deposition rates. Whereas the simulation time step is
on the order of hours, the sedimentation time step is on the order of months
to years. The sedimentation time step varies with location and time and is
tailored to minimize numerical dispersion in sedimentation calculations.
For locations where sediment deposition exceeds scour, TOXIWASP responds
as in Figure 2-4. As sediment and sorbed chemical settle from the water
column, the top bed segment increases in volume, depth, chemical mass, and
sediment mass. Its density remains constant. When the sediment mass in the
top bed segment equals the initial sediment mass in the top two bed segments,
then sediment compression is triggered. At this time the top bed segment
depth (and volume) exceeds the initial top two bed segment depths (and
18
-------
Figure 2-4 TOXIWASP Sediment Burial
to
I
to + At
Segment
1
2
3
4
Time = t0
Depth Density Cone
dt 1.0 Cj(0)
d2 pz 0.0
da f>3 0.0
da Ps C.(0)
Time = t2
Depth Density Cone
dt-d3(2) 1.0 Ct(2)
d2(0)+d3£r Pz C2(2)
da p3 0.0
da P3 C4(0)
Time = t2 + At
Depth Density Cone
dt-d3 1.0 C!(2)
d2 p2 C2(2)
d3 P3 C2(2)g
d3 p3 0.0
+ Compaction: SVOL Compacted = Vol(3) £~ - 1
I I
Pore water volume SVOL squeezed
into water column.
19
-------
volumes) as indicated in Figure 2-4. The top bed segment is now compressed
into two segments.
The new top bed segment has the same depth, volume, and sediment mass
as the initial top bed segment. The new second bed segment has the same
depth, volume, and sediment mass as the initial second bed segment. (In
fact, these properties of lower bed segments always remain constant.) The
combined volume of the top two bed segments is now slightly less than the
volume of the top bed segment just before compression. This volume repre-
sents pore water squeezed into the water column during compression. Chemical
mass in the top two bed segments equals the chemical mass in the top bed
segment just before compression minus the dissolved chemical mass in the
pore water squeezed into the water column. Whereas chemical concentration
remains constant in the top bed segment during the compression step, concen-"
tration increases in the new second bed segment where compression actually
occurs.
Although this sedimentation algorithm minimizes numerical dispersion,
it does lead to discontinuous concentration histories in lower bed segments
as their locations are redefined. For "smoother" output with slightly higher
numerical dispersion, TOXIWASP can recalculate sedimentation ten times per
sedimentation time step. The user can adjust this ratio by altering the
value of "FAG" (Constant 53) from 0.1.
The compression step does not affect any of the properties of the lower
bed segments. In fact, even chemical concentration is unaffected at the
third bed segment and below. Immediately following compression,, however,
TOXIWASP renumbers the lower bed segments, dropping the old bottom segment
from the simulation. Chemical layers, then, are transferred downward intact
and finally lost through the bottom. The accumulated mass lost through the
bottom is saved and printed as variable BMASS(J), where "J" is the number of
the bottom segment. Compression and renumbering completes the TOXIWASP
sedimentation cycle (or time step).
For locations where sediment scour exceeds deposition, TOXIWASP responds
as in Figure 2-5. As sediment and sorbed chemical erode from the bed, the
top bed segment decreases in volume, depth, chemical mass, and sediment
mass. Its density remains constant. When the sediment mass in the top bed
layer equals zero, then segment renumbering is triggered. All the properties
of the remaining bed segments, including chemical concentration, remain
unaffected by renumbering. The new top bed segment, for example, has the
same depth, volume, sediment and chemical concentration as the old second
bed segment. A new bottom bed segment is created with the same physical
properties as the other bed segments. Its chemical concentration, however,
is zero. Renumbering and creation of a new bottom segment completes the
TOXIWASP erosion cycle (or time step).
As a consequence of how TOXIWASP treats sedimentation, certain con-
straints are imposed on the bed segment properties defined in the input data
set. The density (or sediment concentration) of a top bed segment (TYPEE
= 3) must be greater than or equal to the density of the lower bed segments
(TYPEE = 4) within a vertical stack. The volumes, depths, and densities of
20
-------
Figure 2-5 TOXIWASP Sediment Erosion
to
JZ_
1
t2 + At
Segment
1
2
3
4
Time = t0
Depth Density Cone
dj 1.0 0.0
^2 PZ C2(0)
d3 Ps C3(0)
d3 Ps Cj(0)
Time = t2
Depth Density Cone
dl+d2(0) 1.0 Ci(2)
0 Pz C2(0)
d3 P3 C3(0)
d3 p3 C^(0)
Time = t2 -f At
Depth Density Cone
dl+d2(0) 1.0 Ct(2)
d3 P3 C3(2)
d- f\ /O\
Q fs^ *sA.\.&)
d3 /?3 0.0
21
-------
lower bed segments must be constant within a vertical stack. If a smooth
representation of burial is desired, the volumes and depths of lower bed
segments should be significantly smaller than the top bed layer.
If lower bed segments are not included in a TOXIWASP network, then
chemical burial will not take place. The upper bed layer will increase in
volume and depth with settling, or decrease in volume and depth with scour.
TOXIWASP prevents both complete filling of the water column and complete
erosion of the bed.
2.6 TRANSFORMATION PROCESSES
The rate equation for the transformation term K in equations 1, 2a, 16,
and 17 cannot be specified in a general form. Depending on the process under
consideration, many variables may be influencing the rate, leading to a
multi-term and often non-linear relationship. The problem becomes tractible,
however, when it is assumed that the rates are first order with respect to the
pollutant concentration for given environmental conditions, leading to
n
K = z kjC (25)
j-l
where K = total transformation rate, ML~3x~l
k-: = pseudo-first order rate constant for the jth process, which can
vary in space and time, T~l
n = number of processes operating on chemical
This first-order assumption is generally accepted to be accurate for most
chemicals at environmental concentrations. The assumption is invalid at
concentrations near the solubility limit, however. EXAMS does not allow
simulations where predicted chemical concentrations exceed one half of the
solubility limit. TOXIWASP does not automatically abort such runs. If the
user sets CMAX(l) in Card Group K to zero, however, TOXIWASP will abort the
simulation if concentrations become high enough to violate the first-order
assumption. This is accomplished by setting CMAX(l) to one half of the
solubility or 10~5 M} whichever is less. In any case, the user should
analyze results and consider whether this and other assumptions are valid.
Only rarely does the simple mathematical form of equation 25 provide a
rigorous description of the kinetics of pollutant degradation in aquatic
systems. But if the influencing variables that lead to more complex
expressions are reasonably constant for a specific situation then the
influence can be lumped into the single value, k j. The problem of
calculating kinetic transformation then becomes one of calculating the rate
constants k^ from chemical properties and environmental characteristics.
Several investigators have reported processes that influence pesticides
(Mackay and Leinonen, 1975; Paris et al., 1975; Wolf et al., 1975; Wolf et
al. , 1977; Zepp et al., 1975). Smith et al. (1977) have identified the major
pathways that lead to transformation of organic pollutants in the aquatic
22
-------
environment. The role of these processes is receiving increasing attention
and laboratory protocols for related measurements have been developed (Mill
et al., 1982). The general processes identified include hydrolysis, photo-
lysis, oxidation, microbial degradation, and volatilization. TOXIWASP calcu-
lates the individual rate constants kj for each of these processes. Theory
and assumptions are summarized below.
2.6.1 Hydrolysis
Hydrolysis, or reaction of the chemical with water, is known to be a
major pathway for degradation of many toxic organics. Because hydrolysis can
be influenced by pH, the appropriate rate constant, k^, is written as
kh = ka[H+] + kn + kb[OH~J (26)
where k^ = net hydrolysis rate constant, T
ka, kfo = specific acid and base catalyzed rate constants, respectively,
L3n-lT-l
kn = neutral rate constant, T~*
Estimates for ka, kjj, and kn can be determined from laboratory studies
completed over a range of pH values (Mill et al., 1982). The remaining
unknowns in equation 26 can be determined from pH measurements taken in the
aquatic system under investigation or from data from simulations.
In TOXIWASP, hydrolysis by specific-acid-catalyzed, neutral, or
specific-base-catalyzed pathways is considered for dissolved, sorbed, and
biosorbed chemical. Nine second-order rate constants are read-in and
modified based on pH and pOH of the compartment. Hydrolysis, like oxidation,
has paired inputs (ENHG, etc.) that allow calculating or modifying the rates
based on the temperature-based Arrhenius function. For example if KNHG =
2.38E-4 and ENHG = 0.0, then the value of KNHG would be used. Alternately,
if values of ENHG = 20 and KNHG = 11.02 were used, at a temperature of 25°C,
KNHG = 2.36E-4/hr. The value, however, would then change if the water
temperature changed.
2.6.2 Photolysis
Photolysis rate constants are dependent upon both intensity of solar
radiation and structure of the compound (Smith et al., 1977; Zepp and Cline,
1977). This relationship can be written as
kp = 0pl£ALA (27)
where k = net photolysis rate constant, T
0p = reaction quantum yield, moles per einstein
e, = absorption coefficient of the toxicant at wavelength A, per cm per
mole
L, = parameter proportional to incident solar energy flux at wavelength
A, einsteins per square cm per day
23
-------
For most natural waters, adjustments must be made for light extinction (Mill
et al., 1982) as follows.
A
(28)
where a = attenuation coefficient of radiation of wavelength X in water, L~*
£^ = depth of water, L
TOXIWASP calculates photolysis for water compartments by means of the
subroutine PHOT01 from EXAMS. PHOT01 uses an average first order rate
constant measured for near surface waters during cloudless conditions:
3
kp = kpl [Lj Z 0pi a± (29)
1=1
where kpl = average first-order photolytic rate constant for water surface
during cloudless conditions in summer, day~~l
0pi = average reaction quantum yield for compound in form "i"
(dissolved, sorbed, biosorbed) , moles per einstein
[L] = fraction of cloudless summer surface light intensity in
segment, unitless:
-K £ D^
e f
= ( 1~e ) * (1 - .056 CLOUDG) * FACTOR * LIGHTN
Ke& Df
Ke = segment light extinction coefficient, per meter
Df = ratio of optical path length to vertical depth
CLOUDG = average cloud cover, tenths
.056 = cloud type reduction factor (altostratus to altocumulus)
FACTOR = latitude correction factor
LIGHTN = normalized time function for light, representing diurnal
or seasonal changes (0-1.0)
At the heart of this photolysis methodology are two principles. First, the
intensity of light in the water column decreases exponentially as a function
of depth according to the Beer-Lambert Law. The second principle involves
"reaction quantum yields," the efficiency of the reactions in the pollutant
due to absorption of electromagnetic radiation. The "reaction quantum yield1
is actually one of several such yields. Although a molecule activated by
light can undergo several transformations, we are solely interested in the
disappearance yield, that is, the change from parent compound to daughter
compounds.
24
-------
2.6.3 Oxidation
In aquatic systems, chemical oxidation of organic toxicants can be a
consequence of interactions between free radicals and the pollutants. Free
radicals can be formed as a result of photochemical reactions. Free radicals
that have received some attention in the literature include alkylperoxy
radicals, R02»; OH radicals; and singlet oxygen. The oxidation process is
modeled as a second order process:
ROX = kox 1R°2-] • c (30)
where R = net oxidation rate, ML T~*
[ROn.J = molar concentration of oxidant, nL
. k = second order oxidation rate constant, L n T
It should be noted that chemical oxidation is a result of a number of
sequential and competing reactions.
Because of the large number of alkylperoxy radicals that potentially
exist in the environment, it would be impossible to obtain estimates of kox
for each species. Mill et al. (1982) propose estimation of a rate coeffi-
cient using t-butyl hydroperoxide as a model oxidizing agent. They argue
that other alkylperoxides exhibit similar reactivities to within an order of
magnitude.
In addition to estimating a rate coefficient, an estimate of free radi-
cal concentrations must be made to completely define the expression for free
radical oxidation. Mill et al. (1982) report R02 concentrations on the
order of 10~9 M and OH concentrations on the order of 10~17 ^ for a limited
number of water bodies. Zepp et al. (1977) report an average value on the
order of 10~12 M for singlet oxygen in water bodies sampled. The source of
free radicals in natural waters is photolysis of naturally occurring organic
molecules. If a water body is turbid or very deep, free radicals are likely
to be generated only near the air-water interface, and consequently, chemical
oxidation will be relatively less important. In such cases, the concentra-
tions cited above are appropriate in only the near-surface zones of water
bodies.
Oxidation can be very important in an ecosystem. Like photolysis,
oxidation is driven by sunlight and is an important mechanism for pollutant
transformation. Through the interaction of sunlight, humic materials, and
dissolved oxygen, chemical oxidants can accelerate the transformation of
chemicals that may not be directly sensitive to sunlight. TOXIWASP handles
this by coupling second order rate constants for dissolved, sorbed, and
biosorbed chemical with a user entry of the oxidants in the ecosystem. It
corrects this on the basis of temperature. Several studies are cited in the
EXAMS manual in which the steady-state oxidant concentrations (mole/liter) in
natural waters are on the order of E-9 for peroxy radicals, E-13 for singlet
oxygen, and E-17 for hydroxyl radicals in near surface waters. Because these
25
-------
concentrations are presumably functions of solar flux, light extinction
effects may also need to be calculated. These can be calculated similarly to
the method used in TOXIPHOT.
2.6.4 Microbial Degradation or Biolysis
Rate expressions for microbial degradation in aquatic systems must
account for the concentration of microbes present, i.e., a second-order rate
expression usually written as
bac ~ bac bac (31)
where Rhac = net microbial degradation rate, ML~^T~^
= second-order bacterial rate constant, L M T
= concentration of the bacteria, ML
The use of kbac as defined in equation (31) follows from the Michaelis-Menten
equation as follows:
C (32)
Km
where v = specific bacterial growth rate,
,ro
ymax = maximum specific bacterial growth rate, T~
Kffl = half saturation constant, ML"
Bacterial growth is related to chemical substrate degradation with the yield
coefficient Y:
Rbac = I ycbac (33)
Y
where Y = yield coefficient, mass biomass produced per mass of chemical
degraded
If Km »C (usually the case at environmentally observed concentrations) then
we can define kbac as
kbac = Umax (34)
and equation 31 results. Practical utilization of equations 31-34 is
fraught with uncertainties and often one is left with the problem of
estimating kbacCbac as a single value for the system under investigation.
In TOXIWASP biolysis is a process rate based on equation 31, where six
second order rate constants are entered for dissolved, sorbed, and biosorbed
chemicals in the water and in the bed. TOXIWASP bases this calculation on
the concept that population size is an adequate measure of the microbial
community's capacity to degrade the toxicant. Temperature effects are
considered by use of the parameter BIOTMG, which is separated from
26
-------
environmental temperature to allow weighting to match observed data.
ACBACG, the proportion of cells actually degrading the toxicant, is another
means of weighting results.
The biolysis kinetics used in TOXIWASP are not as mathematically
complex as the Monod function or other well known techniques. The simpler
approach is desirable for several reasons. One is that less computer time is
needed to solve the first-order rate equations. Also, for all but very rare
cases, the toxicant is not going to be the predominant substrate for feeding
the microorganisms. In these rare cases, the treatment used here is less
acceptable. Yet another reason for simplification is the experimental error
introduced in measuring populations of microorganisms. In short, these
processes are not yet well defined. Thus, TOXIWASP allows several methods to
correct rates to reflect field data.
2.6.5 Volatilization
Another process that removes organic pollutants from the aquatic system
is volatilization. Volatilization is not a transformation process per se
because the mass of the parent compound is preserved but it can be lumped
into the transformation terms. Volatilization is usually modeled as
\ - kv (Cw - £g_) (35)
H
where Ry = net volatilization transfer rate, ML~^T~^
H = Henry's law constant
C = concentration of the chemical in the bulk gas phase, ML
kv = mass transfer coefficient divided by the average depth of the
water body, T~l
In most cases, the concentrations of organic toxicants in the atmosphere are
much lower than the partial pressures equilibrated with the organic concen-
tration in water. Consequently, equation 35 can be approximated by
Rv - kv Cw (36)
The value of kv, the mass transfer coefficient, depends on the
intensity of turbulence in a water body and in the overlying atmosphere.
Mackay and Leinonen (1975) have discussed conditions under which the value of
kv is primarily determined by the intensity of turbulence in the water. As
the Henry's law coefficient increases, the mass transfer coefficient tends to
be increasingly influenced by the intensity of turbulence in water. As the
Henry's law coefficient decreases, the value of the mass transfer coefficient
tends to be increasingly influenced by the intensity of atmospheric
turbulence.
27
-------
Because Henry's law coefficient generally increases with increasing
vapor pressure of a compound and generally decreases with increasing solubi-
lity of a compound, highly volatile low solubility compounds are most likely
to exhibit mass transfer limitations in water and relatively nonvolatile
high solubility compounds are more likely to exhibit mass transfer limita-
tions in the air. Volatilization is usually of relatively less magnitude in
lakes and reservoirs than in rivers and streams.
In cases where it is likely that the volatilization rate is regulated by
turbulence level in the water phase, estimates of volatilization can be
obtained from results of laboratory experiments. As discussed by Mill et al.
(1982), small flasks containing a solution of a pesticide dissolved in water
that have been stripped of oxygen can be shaken for specified periods of
time. The amount of pollutant lost and oxygen gained through volatilization
can be measured and the ratio of mass transfer coefficients for pollutants
and oxygen can be calculated. As shown by Tsivoglou and Wallace (1972),
this ratio should be constant irrespective of the turbulence in a water
body. Thus, if the reaeration coefficient for a receiving water body is
known or can be estimated and the ratio of the mass transfer coefficient for
the pollutant to reaeration coefficient has been measured, the pollutant
mass transfer coefficient can be estimated.
In TOXIWASP, volatilization is the first process considered when the
segment loop is entered. Only the unionized, dissolved molecular form of a
chemical in a water compartment with an air interface can be volatilized.
The volatilization calculation procedure used is based on the two-resistance
principle.
The two-resistance method assumes that two "stagnant films" are bounded
on either side by well mixed compartments. Concentration differences serve
as the driving force for the water layer diffusion. Pressure differences
drive the diffusion for the air layer. From mass balance considerations it
is obvious that the same mass must pass through both films, thus the two
resistances combine in series. There is actually yet another resistance
involved, the transport resistance between the two interfaces, but it is
assumed to be negligible. This may not be true in two cases: very turbulent
conditions and in the presence of surface active contaminants. Although this
two-resistance method, the Whitman model, is rather simplied in its assump-
tion of uniform layers, it has been shown to be as accurate as more complex
models. Laboratory studies of volatilization of organic chemicals confirm
the validity of the method as an accurate predictive tool (EXAMS, sec.
2.3.2).
The volatilization rate constant ky is calculated as follows:
1/Lw (37)
k =
v RESLIQ + RESGAS
__ 1 1 (38a)
RESLIQ - * or
28
-------
RESGAS = (38b)
H * WATA/18/MWTG
TKEL * 8.206E-5
K02L = reaeration velocity in segment, meters per hour
KVOG = measured liquid phase transport resistance, ratio
MWTG = molecular weight of compound
WAT = piston velocity of water vapor, meters per second
= .1857 + 11.36 * WIND
WIND = local, time varying wind speed, meters per second
TKEL = water temperature, degrees Kelvin
8.206E-5 = gas constant
While EXAMS requires that reaeration velocity K02G(J) be specified,
TOXIWASP calculates flow-induced reaeration in subroutine TOXIREOX based on
the Covar method (Covar, 1976). This method calculates reaeration as a
function of velocity and depth by one of three formulas, Owens, Churchill, or
O'Connor-Dobbins, respectively:
K20 = 27.6*AVVELEO-67/AVDEPE0.85 (39a)
K20 = 14.8*AVVELE0.969/AVDEPE°'673 (39b)
or K20 = 16.4*AVVELE°'5/AVDEPE°'5 (39c)
where K20 = reaeration velocity at 20°C, cm/hr
AVVELE = average segment velocity, ft/sec
AVDEPE = average segment depth, ft
The Owens formula is automatically selected for segments with depth less than
2 feet. For segments deeper than 2 feet, the O'Connor-Dobbins or Churchill
formula is selected based on a function of depth and velocity. Deeper,
slowly moving rivers select O'Connor-Dobbins, while moderately shallow,
faster moving streams select Churchill.
Whenever the volatilization rate is calculated during a simulation,
wind-induced reaeration is determined by
K20 = -0.46 * WINDSG + .136 * WINDSG2 (40)
where WINDSG = time-varying wind velocity, meters/sec
A minimum value of 2 cm/hr is imposed on K20. Windspeed affects reaeration,
then, above 6 meters/sec. The reaeration velocity used to compute
volatilization is either the flow-induced reaeration or the wind-induced
reaeration, whichever is larger. Segment temperatures are used to adjust K20
by the standard formula
K02L = 0.01 * K20 * 1.024(TCELG - 20.) (41)
where K02L = segment reaeration rate, meters/hour
TCELG = segment temperature, °C
29
-------
2.7 SPECIAL DATA INPUT OPTIONS
For some simulation problems, flows and mass loadings are quite
variable. Pesticides, in particular, typically enter streams in pulse
loadings associated with runoff events and increased stream flows. To
generate accurate frequency predictions in such situations, both flows and
loads should be specified with more precision than is convenient with the
standard WASP input options.
TOXIWASP provides a special option that causes the input of daily-
average flows and loadings from an auxiliary file. When the constant LOPT
is set to 1.0, TOXIWASP will read Unit 91 at the beginning of each new simu-
lation day. Unit 91 is an unformatted sequential file named 'DAILY.DAT1 in
the POP version. For the IBM version, Unit 91 is named in the JCL. Flows
from Unit 91 replace any flows that may have been read in Card Group D.
Mass loadings from Unit 91, however, are added to those read in Card Group
F. Loads from Card Group F are considered point source loads., while those
from Unit 91 are considered nonpoint source (NFS) loads.
Unit 91 can contain either flows, nonpoint source loads, or both. The
first four variables read are NWKS, MQOPT, MOQ, and SCALQ. NWKS is the
number of nonpoint source loads. If NWKS is 0, the NFS load input section
is bypassed. MOQ is the number of flows. If MOQ is 0, the flow input sec-
tion is bypassed. SCALQ is a scale factor to convert flows to cubic feet
per second (cfs). Daily NPS loads must already be in pounds per day. MQOPT,
a flow option, is read but not used.
If daily flows are to be read, the second set of variables on Unit 91
defines the segment flow pairs. MOQ segment flow pairs are read. Each flow
pair K consists of the upstream segment, JQ(K), and the downstream segment
IQ(K). Positive flows travel from JQ(K) to IQ(K). The Fortran Read
Statement in TOXIWASP is:
READ(91)((JQ(K),IQ(K),K = 1, MOQ)
If daily flows are not to be read, this set of variables should not appear on
Unit 91.
The rest of Unit 91 contains daily flows and/or NPS loads. At the
beginning of each new day, NPS loads of chemical and sediment are read for
each of NWKS input segments. These segments must be in the order defined by
Card Group F. The Fortran Read Statement in TOXIWASP is:
READ(91)((NPSWK(I,J),I = 1,2), J=l, NWKS)
where NPSWK (1,J) is NPS chemical load in Ib/day and NPSWK (2,,J) is NPS
sediment load in Ib/day. If NWKS is 0, this set of variables should not
appear on Unit 91.
30
-------
Next, flows are read for each of MOQS flow pairs. The flows must be in
the order defined by the second set of variables on Unit 91. The Fortran
Read Statement in TOXIWASP is:
READ(91)(BQ(K),K=1,MOQ)
where BQ(K) is the flow from JQ(K) to IQ(K) in cfs.
of variables should not appear on Unit 91.
If MOQS is 0, this set
The contents of Unit 91 are summarized on Table 2-2. Simple programs
can be written to produce this file. The following software for the PDF
computer is available to construct a proper Unit 91 file.
DAILYCVT - This program converts flows and loadings from user
specified files into an unformatted, sequential file for
input to TOXIWASP. The user is queried for the load file
name and format, and the flow file name and format. The
indicated files are read and data are written to the file
DAILY.DAT in the proper sequence.
HSPFCVT - This program converts an HSPF plot file to a formatted
WASP file with NFS loads. The user is queried for seg-
ment numbers and associated acreage for which the HSPF
unit loads apply. The user then enters the beginning
and ending date. The HSPF file is read, and data are
written to the file WASPLD.DAT. This file can be pre-
pared for input to TOXIWASP using DAILYCVT.
SWRVCVT - This program converts a SWRV plot file to a formatted
WASP file with NFS loads. It is equivalent to HSPFCVT.
The SWRV file is read, and data are written to the file
WASPLD.DAT. This file can be prepared for input to
TOXIWASP using DAILYCVT.
Some simulation problems require detailed, time-varying flows and
volumes. For these, WASP has been modified to accept hydrodynamic data
generated by the Stormwater Management Model (SWMM), receiving water block
(RECEIV), program SWFLOW (see Huber, et al., 1975). The network for SWFLOW
must be equivalent to the network for WASP (i.e., the SWFLOW junctions, or
nodes, are the same as the WASP segments). Hydrodynamic data from SWFLOW are
stored on Unit 92, and read by WASP through subroutine SWFLOW. To invoke
this option, both IVOPT and IQOPT in Card Groups C and D should be set to 4.
The water quality time step specified in the SWFLOW hydrodynamic simulation
is automatically used in the WASP simulation.
31
-------
Table 2-2. Structure of Auxiliary Flow and
Nonpoint Source Load File 'DAILY.DAT'
Condition Da)
0
MOQ > 0
Contents*+
NWKS,MQOPT,MOQ,SCALQ
JQ(1),IQ(1),JQ(2)>IQ(2),...,JQ(MOQ),IQ(MOQ)
NWKS >0
MOQ >0
NPSWK(1,1),NPSWK(2,1),NPSWK(1,2),...,NPSWK(2,NWKS)
BQ(1),BQ(2),...,BQ(MOQ)
NWKS >0
MOQ >0
NPSWK(1,1,NPSWK(2,1),NPSWK(1,2),...,NPSWK(2 ,NWKS)
BQ(1),BQ(2),...,BQ(MOQ)
NWKS > 0 TEND
MOQ > 0
NPSWK(1,1),NPSWK(2,1),NPSWK(1,2),...,NPSWK(2,NWKS)
BQ(1),BQ(2),...,BQ(MOQ)
* Data written sequentially.
+ The following variables are integers: NWKS, MQOPT,, MOQ, JQ, IQ.
All other variables are real (floating point) numbers.
32
-------
2.8 SPECIAL SIMULATION CONTROL
For seasonal or multiyear simulations, the environmental conditions
affecting chemical transformation rates change slowly relative to the
simulation time step. Consequently, it is inefficient to recalculate the
chemical transformation rates every time step, or even every day. TOXIWASP
recalculates the total transformation rate at constant intervals specified by
the user. The time interval between rate computations is given by the
constant TDINT. To investigate seasonal variability, the user may specify
TDINT of 15 or 30 days. To investigate diurnal variability, the user may
specify TDINT equal to a short time step of, say, 0.05 day. A TDINT value of
0 is automatically set to a very large number. Rates are calculated the
first time step only, and are held constant throughout the simulation.
The user can specify the total chemical transformation rates with the
spatially variable parameter TOTKG(J). If TOTKG(J) is positive for a seg-
ment, rate calculations are bypassed for that segment. Specified rates are
constant throughout the simulation. It is possible to specify rates for
some segments and have TOXIWASP calculate rates for the remaining segments.
If TOTKG(J) is 0 for a segment, TOXIWASP calculates a rate and stores it (as
a negative number) until the next recalculation TDINT days later.
WASP requires the user to specify simulation time steps, which can vary
throughout the simulation. Three problems can occur. If the time step is
too large, computational instability develops and halts the simulation. If
the time step is too small, numerical dispersion builds up and leads to
inaccuracies. A small time step is also inefficient and leads to high costs.
The optimal time step varies inversely with flow and dispersion, and directly
with segment size:
TRT = VOL (42)
TQ + TR
where TRT = optimal time step for segment, days
VOL = volume of segment, MCF
TQ = total flow through segment, MCF/day
TR = total dispersive exchange, MCF/day
For irregular networks with unsteady flow, TRT varies both spatially and with
time. If the constant DTOPT is set to 1.0, TOXIWASP calculates the optimal
time steps TRT for the network and selects the minimum to avoid instability.
TRTs are recalculated every time step for flow options 2 and 3, and every day
for the special flow input option through Unit 91.
2.9 SPECIAL OUTPUT OPTIONS
Standard output from WASP is limited to regular intervals, which must be
approximately ten days apart for a year's simulation. For problems with
variable mass loadings, this output can miss significant concentration
33
-------
events. TOXIWASP allows the user to have peak events printed out. Two
constants must be specified. The first is the reference concecitration CTRIG.
The second is the reference segment XJTR. Whenever concentrations rise above
CTRIG in segment XJTR, special event printout is triggered. Chemical
concentrations for all segments are printed out every three hours until the
event ends. If CTRIG is set to 0.0, this special event printout is skipped.
For many problems, statistical summaries are the most convenient form of
output. The use of the cumulative distribution function (cdf) to estimate
frequency for exposure assessment was discussed in Section 1.2. When a
reference segment XJTR is designated, TOXIWASP builds a file containing the
time history of chemical concentrations at that segment and the top benthic
segment below it. Total, dissolved, sediment-sorbed, and biosorbed chemical
concentrations are written to the file every three hours throughout the
simulation. Following the simulation, this file is processed by program
TOXIFREQ and subroutine TOXIORDR to obtain statistical summaries. Two tables
are printed. The first describes the frequency distribution with the mean,
standard deviation, skewness, kurtosis, and select percentiles. The second
gives the cumulative frequencies at intervals of 0.02.
2.10 SUMMARY OF TOXIWASP CAPABILITIES
The following finite difference equations summarize the capabilities of
TOXIWASP. They refer to segments 3 and 4 of the river shown in Figure 2-1.
The symbols used are defined in the appropriate sections of this chapter,
with subscripts added to indicate segment number. G£ 3> for example, revers
to sediment concentration (C2) in segment 3. Subscripts separated by a dash
refer to the interface between two segments. Qi-3 > f°r example, refers to
the flow between segments 1 and 3. Subscript "k" refers to the chemical
form: dissolved, sorbed, or biosorbed. ct^ .= , then, refers to the fraction
of chemical of form "k" in segment "j " , and «2 4 ^s tne sorbed fraction in
segment 4.
Solids in Segment 3 of River
= Qj_3 C£ i ~ 0.3-5 C2 3 Advection (2.1)
C2.1 - C2,3 C2>3- C2;>5
+ E^-3 A^_3 Ll-3 ~ E3-5 A3-5 L3-5 Dispersion
(2.1)
+ W£ 3 Mass Loading (2.1)
4 - Lw 3 Scour, Deposition
(2.2)
34
-------
Chemical in Segment 3 of River
\7 * = O C — O C
V3 dt ^1-3 1,1 W3,5 1,3 Advection (2.1)
A C1>3-Clt5
1-3 1-3 Lj_3 3-5 3-5 L3_5 Dispersion (2.1)
+ Wj 3 Mass Loading (2.1)
+ Ws,4G2,4 * Ci>4 "2>4
Lb>4 C2>4 Scour (2.2, 2.4)
,3 "2,3
3 C2 3 Deposition (2.2,
2.4)
+ Db A .Cl,4 al,4 - C1.3 al,3
3-4 3-4 I L_ J Pore Water
D A
- L_ Po
Diffusion (2.4)
Db A * f * A
3-4 3-4 s,4 04 Direct Sorptive
CT A a-> A ci -3 ao i ^ A Exchange with Bed
*[_iii__lti _ J-^3 2,3 t _4j Surface (2.4)
C2,4 C2,3 ^3
,C1^ aL* or C1.3 oci.3
P,4l 0^ 03 J Percolation or
Infiltration (2.4)
3
IE (ka k[H'f]3 + kn k + kb k[OH~]3) ak)3JC1 3 Hydrolysis (2.6.1)
k=l ' ' ' '
3
kpl[LJ3 [ £ 0p,k ak,3Jcl,3 Photolysis (2.6.2)
k=l
3
I E kox,k^R02-b ak,3^cl,3 Oxidation (2.6.3)
k=l
3
[ z kbac,ktcbach ak,3lcl,3 Bacterial Degrada-
k=l tion (2.6.4)
35
-------
-kv (wind, current, depth) al 3 GI 3 Volatilization
(2.6.5)
Chemical in Segment 4 of River
v dCl»4
4 dt = Deposition-Scour-Pore Water Diffusion
-Direct Sorptive Exchange with Bed Surface
-Percolation (or + Infiltration)
-Hydrolysis-Oxidation-Bacterial Degradation
Networks with multiple benthic segments may also transport chemical
through the process of sedimentation. For example, Segment 2 of the lake in
Figures 2-1 and 2-3 can be represented by the following finite difference
equation:
Lake Segment 2
dC 1 ?
V
2 dt = Deposition-Scour
-Pore Water Diffusion
-Direct Sorptive Exchange with Surface
-Percolation (or + Infiltration
-Hydrolysis-Oxidation
-Bacterial Degradation
wd ws G! 2 a2 2
(*T C2 1 ~ h. C2 2^ ^2 ^2~~2'~ Burial (2-5)
Wd C2)i ^s C2,2 C2>3
"^L, C9 o ~ L0 C9 o^ 2 ^C9 9 -1^ Pore Water Squeezed
1 2,3 2 2,3 2,2 (2>5)
36
-------
SECTION 3
WASP MAINLINE AND SUBROUTINE OVERVIEW
WASP is a modular program. Its many subroutines can be grouped into the
standard categories of "input," "process," "output," and "utility," as in
Figure 3-1.
The WASP mainline is really just a program control module. As such it
assigns the logical units (disk) for temporary scratch files and the dump-
save files generated in the kinetic subroutine, and controls the calling
sequence of the WASP subroutines.
The following subroutine descriptions include references (within paren-
theses), to Card (or Record) Groups. These groups of WASP input data are
described in Section 4, WASP Input Structure. Also, in the following
subroutine description, occasional reference to internal WASP program
variables will be made. A complete definition of these variables is pre-
sented in the Glossary Section entitled "WASP COMMON."
WASP1
WASPl reads the model identification and system by-pass options for the
user's model (Card Group A). WASPl also performs some variable initializa-
tion. Two variables of interest set by WASPl are NBCPSY and NWKPSY. NBCPSY
and NWKPSY are used, respectively, to indicate the maximum number of boundary
conditions and forcing functions for which WASP is dimensioned.
WASP2
WASP2, depending upon the input option selected by the user, reads Card
Group B, the time variable or constant exchange coefficients (or dispersion
coefficients, cross-sectional areas, and characteristic lengths with appro-
priate conversion to exchange coefficients). WASP2 also reads the exchange
by-pass options for each system.
WASP3
WASP3 reads the segment volumes (Card Group C).
WASP4
WASP4, depending upon the input option selected by the user, reads the
time variable or constant advective flows (Card Group D). WASP4 converts the
flows from cubic feet per second to million cubic feet per day. WASP4 also
reads the flow by-pass options for each system.
37
-------
Du 0 n
O O ?1
3 a <
s ;=! o
Et. fc CO
1 1
O
^ i
CO
0) M
d ^
^ K^
.2? o
DH
OH
^
-WMESS
c
c
• r
-t-
^
-+•
"
f
c
•«c
^
2
£
P
c/
<:
&
a
P
t
0
L)
H
J
H
H
J
J
H
H
li
H
^
-t
si
J)
H
0
U
»•
••
a
c
e
^ P
H E-
a 6
CO 0
1
-FMTER
a
o
i
i
OT 1 lO
CD i *-i
O CO
0 I <
1
1
-4-)
n
h~H
i-H
01
1
1-H
1
T— 1
CO
1
C\J
a,
CO
£
\
a
T"
C/
&
(
K
r
<
0
0
p
0
3 05
•3 CO
" 1
CD
.-. CO
Q) "3
6 OT
CD
U ^3 ^
H ? CO
3 .S
-------
WASPS
WASP5, depending upon the Input options selected by the user, reads the
time variable or constant boundary conditions for each system in the user's
model (Card Group E).
WASP6
WASP6, depending upon the input options selected by the user, reads the
time variable or constant forcing functions for each system in the user's
model (Card Group F).
WASP7
WASP? reads the kinetic constants, segment parameters, and kinetic
piecewise linear functions of time (Card Group G, H, and I, respectively).
WASPS
WASPS is used to update the piecewise linear functions of time, if any,
for exchange coefficients, advective flows, and the miscellaneous kinetic
functions. This means computing new slopes and intercepts, and setting a
variable to indicate the next simulation time that the functions are to be
updated. The following convention is used for the itn update.
slope = f(t)j+i - f(t)j
intercept
next update time
WAS8A
WAS8A is used to update the piecewise linear functions of time, if any,
for boundary conditions and forcing functions. This means computing new
slopes and intercepts for any system or state variable that requires an up-
date, and setting a variable to indicate the next simulation time that the
piecewise linear functions are to be updated. The same conventions used in
WASPS are used in WAS8A for computing slopes and intercepts.
WASP9
WASP9 reads the initial conditions or initial segment concentrations for
each system or state variable (Card Group J). WASP9 also reads the stability
and accuracy criteria (Card Group K) for each system.
39
-------
WAS 10
WAS10 reads the print control options (Card Group L), consisting of the
print interval and up to eight system-segment pairs for intermediate print-
out during the simulation.
WASH
WASH reads the integration control information (Card Group M). The
integration control information includes the integration step-size or sizes
to be used, the total simulation time, the starting time for the simulation
(if not zero), and whether negative solutions will be permitted.
WAS12 - WAI 2A
WAS12 and WA12A act together to complete the evaluation of the mass
balance Equation (2.16). Upon entry to WAS12, only the kinetic portion of
the mass balance equation has been evaluated, which for discussion purposes
will be noted as K V C . WAS12 then goes through the following steps.
a. Using the IQ and JQ vectors as drivers, WAS12 computes
or (V C-m) = K-V C-m + ZQ..C.m — EO. C-m
where Q may be computed as a function of time using the MQ (slopes)
and BQ (intercepts) vectors or Q may be a constant (using BQ vector)
b. Using the IR and JR vectors as drivers, WAS12 computes
or (ViCim) = K^C^ + SQjiCj111 - Z3±}!.C±m + ER^ (C^ - C/1)
where R may be computed as a function of time using the MR (slopes)
and BR (intercepts) vectors, or R may be a constant (using the BR
vector) .
c. Using the IWK vector as a driver, WAS12 computes
(ViC^) = (vicim) + v±m
or (V±C±m) = K^C^ + SQjiCi1" - ^QifcCi111 + £Ri-j(C.jm - C^) + W±m
where W may be computed as a function of time using the MWK (slopes)
and BWK (intercepts) vectors or W may be a constant (using the BWK
vector) .
40
-------
d. WAS12 completes the evaluation of the derivative, C^ by dividing
through by the volume and multiplying by the time scale factor
(SCALT)
Cj™ = SCALT x (V^111)/^
Note: C^ now has units M/L3/T - nominally mg/l/day.
WAS 13
WAS13 is used to print the intermediate system-segment pairs during the
simulation.
WAS14
WAS14 is used to adjust the integration stepsize as necessary, during
the simulation.
UNIT
TINIT is used to set the initial slopes and intercepts for any piecewise
linear functions of time if the user starts a simulation at any time other
than time equal to zero.
WAS15 or EULER
WAS15 or EULER is the heart of the integration procedure. It determines
the calling sequence for TINIT, WASPB (the user's kinetic subroutine), WAS12,
WA12A, WAS13, and WAS14. WAS15 performs a second order Runge-Kutta integra-
tion, while EULER performs a first-order Euler integration. A brief flow-
chart is presented in Figure 3-2.
WAS16
WAS16, dependent upon the user's input data, retrieves and prints the
state variables (or water quality concentrations) and any other variables of
interest (that the user computed in WASPB) from the auxiliary storage (disk
files) that were generated in WASPB during the simulation.
SWFLOW
SWFLOW reads hydraulic information created by SWMM (RECEIV block, program
SWFLOW). This subroutine is invoked when IQOPT in Card Group D is set to 4.
IVOPT in Card Group C must also be 4. The file is named "SWFLOW.DAT," and is
locted on unit 92.
41
-------
Figure 3-2
Simplified EULER/WAS15 Flow Chart
Call EULER^N
or WAS15 J
^
Calculate Kinetic Derivative
Calculate Transport, Load Derivative
Intermediate
Dumps
Calculate Kinetic Derivative
42
-------
WAS 17
WAS17, dependent upon the user's input data, retrieves and provides
printer plots for the state variables (or water quality concentrations) and
any other variables of interest (that the user computed in WASPB) from the
auxiliary storage (disk files) that were generated in WASPB during the simu-
lation.
WAS 18
WAS18, written especially for the DEC PDF 11/45, provides off-line
digital pen plotting capabilities. The pen plots generated are similar to
the printer plots available through WAS17, but in addition permit the user
to overplot his observed field data for model-data comparison.
WAS 19
WAS19, dependent upon the user's input data, retrieves and provides
printer plots of the spatial variation at selected times for the state vari-
ables (or water quality concentrations) and any other variables of interest
(that the user computed in WASPB) from the auxiliary storage (disk files)
that was generated in WASPB during the simulation.
Auxiliary plot subroutines called by WAS19 include SIR, PLOT, BLKPLN,
and ENCOD (for IBM version only).
Miscellaneous Subroutines
WASP also includes a number of subroutines that perform rather trivial,
but necessary, operations. These subroutines include SCALP, WMESS, WERR,
SETIA, SETRA and FMTER.
DEC PDP Subroutines
Two special subroutines were written for the DEC computer system. These
subroutines FILEOP and FILEOC were necessitated due to the way the DEC
operating 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.
TOXIWASP Kinetic Structure
The standard WASP kinetic subroutine shown in Figure 3-1 is actually a
set of subroutines in TOXIWASP. The structure of these subroutines is
illustrated in Figure 3-3. These WASPB subroutines combine chemical and
environmental parameters to first-order rates and, from chemical
concentrations passed by WASP, calculate kinetic derivatives. These
derivatives are passed back to WASP where they are integrated along with the
transport and loading derivatives every time step.
43
-------
£
^.
o
-p
rn
^
« rn ^
rH >•'•« -
(
(
I
X
3
— )
Q
><
0
Q
W
X
O
I-J
CO
o
Q
O
Ct,
X
o
EH
EH
O
w
o
EH
— .
E~*
0
EH
X
O
w
o
EH
44
-------
TOXIWASPB
TOXIWASPB serves as the main program for the kinetic portion of
TOXIWASP, calling other subroutines when appropriate. Initialization is
performed during the first time step by calling TOXINIT. As the simulation
progresses, the proper time- and space-variable environmental and chemical
characteristics are calculated, then passed to TOXIFORD. Kinetic derivatives
are calculated based on first-order rate constants returned from TOXIFORD.
These derivatives are then adjusted for settling, erosion, and percolation
by calling TOXISETL. Further adjustments in the derivatives due to pore
water mixing and sediment-water exchange are calculated by calling TOXISEDW.
Variables are periodically dumped to a save file and intermediate parameters
and results are printed by calling TOXIDUMP. Finally, daily flows and non-
point' source loads are read from an auxilliary file.
TOXINIT
TOXINIT is called during the first time step only to initialize
parameters and functions for the chemical simulation. A sequential file is
set up for reading daily flows and nonpoint source loads, and two tables are
printed. The top sediment layer for the special print segment is identified.
Effective partition coefficients for each segment are calculated from either
the organic carbon or octonol-water partition coefficient, the spatially
variable sediment organic carbon fractions, and the target organism organic
carbon content. For benthic segments, the units of biomass are adjusted and
porosity is calculated. The active bacterial population is calculated for
all segments. Initial and reference bed volumes are saved, values are ini-
tialized to 0, and CMAX(l) is set to assure the first-order decay assumption.
Finally, TOXINIT checks for proper segment alignment, with assigned numbers
increasing sequentially from water surface to bottom benthic segments.
TOXIFORD
TOXIFORD, a modification of the EXAMS subroutine FIRORD, calculates
total first order chemical transformation rates for each segment as the
simulation progresses. The total first-order rate for each segment is the
summation of the transformation (including degradation or transfer) rates
due to five processes: hydrolysis, oxidation, bacterial degradation, volati-
lization, and photolysis. These individual rates are calculated from ambient
environmental and chemical characteristics passed from TOXIWASPB.
Whenever total rates are to be recalculated during a simulation,
TOXIWASPB calls TOXIFORD once for each segment. TOXIFORD first calculates
the volatilization transfer rate of dissolved chemical from surface water
segments by calling TOXIVOLT. Next, the photolysis rate constant is calcu-
lated for water segments by calling TOXIPHOT. Then, for each species of the
chemical (dissolved, sediment-sorbed, and biosorbed), first-order transforma-
tion rates are calculated for photolysis, hydrolysis, oxidation, and bac-
terial degradation. Finally these individual rates are summed to give a
total first-order chemical disappearance rate, which is passed back to
TOXIWASPB.
45
-------
TOXIVOLT
When called by TOXIFORD, TOXIVOLT computes the volatilization transfer
rate constant for a surface water segment using a two-film model of movement
of toxicant across the air-water interface. Liquid phase resistance is
computed from oxygen reaeration rates modified by the chemical molecular
weight. Gas phase resistance is computed from wind speed and Henry's Law
constant, which is supplied by the user or calculated from vapor pressure and
solubility data. The volatilization rate is computed from the air and water
phase resistances and the depth of the surface water segment.
TOXIPHOT
When called by TOXIFORD, TOXIPHOT computes the photolysis transformation
rate constant for a water segment. TOXIPHOT accepts a measured (clear day)
photolysis rate constant at a specified reference latitude as input data.
Next, this rate constant is corrected for the latitude, cloud cover, and
light extinction in the water column. The rate is further modified by a
time-varying input function to approximate seasonal changes in incident light
intensity.
TOXIREOX
TOXIREOX is call by TOXIFORD for surface water segments to calculate
reaeration velocities. These are used in TOXIVOLT for calculating
volatilization rates. For the first time step only, TOXIREOX calculates
average flow-induced reaeration rates by the Covar method (Covar, 1976) using
average segment velocities and depths. In subsequent time steps, TOXIREOX
calculates time-varying wind-induced reaeration rates. The reaeration rate
returned to TOXIFORD is either the flow-induced rate or the wind-induced
rate, whichever is larger. TOXIREOX is slightly modified from the HSPF
subroutine OXREA.
TOXISETL
TOXISETL is called by TOXIWASPB for each segment and each time step to
calculate settling of chemical and suspended sediment from the water column,
erosion of chemical in the bed, and percolation of dissolved chemical
vertically through the bed. Concentration derivatives are adjusted within
TOXISETL. First, settling rates are calculated from spatially variable
settling velocities, segment depths, concentrations of sorbed chemical and
suspended sediment. Erosion is calculated from spatially variable yearly
depletion rates from the bed sediment surface, along with sediment density
and sorbed chemical concentration. Bed sediment segments are assumed to
maintain their physical characteristics, such as density, porosity, and
organic content. Vertical pore water percolation is calculated from spa-
tially variable vertical flow rates, along with porosity and dissolved con-
centrations. Positive flow is upward, and negative flow is downward. It
should be noted that downward percolation will eventually transport the
chemical out of the bottom benthic segment.
46
-------
TOXISEDW
TOXISEDW is called by TOXIWASPB for surface benthic segments each time
step to calculate dispersive exchanges of chemical between the bed and the
water column. Concentration derivatives are adjusted within TOXISEDW. The
two mechanisms are pore water diffusion and local surface sediment
equilibration with the water column. Pore water diffusion is calculated from
spatially variable diffusion coefficients, surface areas, characteristic
mixing lengths, and sediment porosity. Sorption-desorption of chemical
between overlying water and the benthic surface is calculated using sediment
turnover rates and chemical partition coefficients. Because local equili-
brium is assumed, sorptlon-desorption is controlled by the fraction of under-
lying sediment brought into contact with the overlying water per unit time.
This fraction is related to the pore water diffusion rate by a spatially
variable multiplier supplied by the user. For immobile, armored stream
reaches, the multiplier may be 0, whereas for reaches vigorously mixed by
physical or biological processes, the multiplier may be 1, with high pore
water dispersion coefficients as well.
TOXIDUMP
At specified intervals, TOXIWASPB calls TOXIDUMP to prepare or print
output from the simulation. Three kinds of output are handled. The first is
the standard WASP dump of select variables to a save file at fixed intervals.
Eight variables are dumped for "System 1": total chemical concentration
(mgc/l-jO , dissolved chemical concentration (mgc/lw), sediment-sorbed chemical
concentration (mgc/kgs), biosorbed chemical concentration (ugc/gb), chemical
fraction dissolved, chemical fraction sediment-sorbed, chemical mass in seg-
ment (kg), and chemical mass lost through volatilization or burial (kg).
Eight more variables are dumped for "System 2": total sediment concentration
(rags/If), segment depth (ft), the total transformation rate (per day), the
photolysis rate, the hydrolysis rate, the biodegradation rate, the oxidation
rate, and the volatilization rate (all specific rates in units of per hour).
The second kind of output is event-triggered. When the total chemical
concentration exceeds a reference value at a selected segment, total chemical
concentrations are printed for every segment every 3 hours until the
concentration again falls below the reference.
The third kind of output is the dump of chemical and sediment
concentrations at a selected water segment (and its top benthic segment) to a
save file every 3 hours. This file is processed by TOXIFREQ after the
simulation to yield statistical information, including a cumulative frequency
table. Concentrations analyzed include total chemical, dissolved chemical,
sediment-sorbed chemical, and biosorbed chemical.
TOXIFREQ
After completion of the simulation, TOXIFREQ processes a file of chemi-
cal concentrations at a water and a benthic segment, producing statistical
tables. A save file of concentrations written every 3 hours by TOXIDUMP is
47
-------
first ordered by TOXIORDR. Next, the statistics are computed, including
minimum, maximum, and mean; various percentiles; standard deviation;
skewness; and kurtosis. This table is printed, then a cumulative frequency
table is prepared, with concentrations corresponding to ascending even
probability values (0.0, 0.02, 0.04,..., 0.98, 1.00). This table is printed,
and a plot file is prepared.
48
-------
SECTION 4
PREPARATION OF DATA INPUT
This section describes the input required to run TOXIWASP. This section
is reproduced, mostly intact, from the WASP user's manual. Although the
general input structure allows many systems (or state variables), TOXIWASP
allows only two: organic chemical and sediment. Furthermore, whereas WASP
allows arbitrary segment numbering schemes, TOXIWASP requires that vertically
layered sections be numbered sequentially from surface water to bottom
benthic.
To arrange the input data into a logical format, the data cards required
are divided into 14 card groups, A through N. The card groups are briefly
summarized in Table 4-1. For each card group, a brief description of each
card is given to define the variables that appear within the group, and any
options that may be available. Depending upon the structure of the user's
model, a certain card group, or cards within a group, may not need to be
inputted. Where it is appropriate, the manual informs the user how to avoid
inputting unnecessary information.
To help the user visualize how the data input file is organized, a
sample dataset is provided as Appendix 4. This dataset describes the
constant addition of the chemical dibenzothiophene into the pond illustrated
in Figure 2-1 (page 6).
Provisions for handling a wide range of time and space scales are con-
tained in the data input structure via scale factors. These scale factors
facilitate the conversion of the user's time and space units to those con-
sistent with the WASP program. The standard units for the WASP input data
are detailed in this manual. Departure from these units necessitates the
use of appropriate scale factors. Scale factors may also be used to scale
input data in sensitivity analysis. For example, a user may wish to test
the effect of increased dispersion upon the model. Rather than altering all
the interfacial dispersion coefficients he may simply change the dispersion
coefficient scale factor to reflect the appropriate increase in dispersion
levels.
4.1 CARD GROUP A - MODEL IDENTIFICATION AND SYSTEM BYPASS OPTIONS
The variables that appear on each card are as follows.
4.1.1 Model Identification Numbers
49
-------
TABLE 4-1. SUMMARY OF CARD GROUPS
Card Group
A. Model Identification and System Bypass Options
1. Model Identification Numbers
2. Title Card
3. Simulation Option
4. System Bypass Option
B. Exchange Coefficients
1. Number of Exchange Coefficients, Input Option Number
2. Scale Factor
3. Exchange Coefficients
4. Exchange Bypass Options
C. Segment Volumes
1. Number of Volumes, Input Option Number
2. Scale Factor
3. Volumes
D. Flow
1. Number of Flows, Input Option Number
2. Scale Factor
3. Flows
4. Flow Bypass Option
E. Boundary Conditions
1. Number of Boundary Conditions, Input Option Number
2. Scale Factor
3. Boundary Conditions
Cards 1-3 are inputted for each system of the model
F. Forcing Functions
1. Number of Forcing Functions, Input Option Number
2. Scale Factor
3. Forcing Functions
Cards 1-3 are inputted for each system of the model
G. Parameters (Environmental Characteristics)
1. Number of Parameters
2. Scale Factors
3. Parameters
H. Constants (Chemical Characteristics and Special Options)
1. Number of Constants
2. Constants
50
-------
TABLE 4-1. (CONT.)
Card Group
I. Miscellaneous Time Functions (Environmental Variability)
1. Number of Time Functions
2. Functions Name, Number of Breaks in Function
3. Time Function
Cards 2 and 3 are inputted for each time function required by
the model
J. Initial Conditions
1. Initial Conditions for each system of the model, i.e.,
chemical and sediment
K. Stability and Accuracy Criteria
1. Stability Criteria
2. Accuracy Criteria
L. Intermediate Print Control
1. Print Interval
2. Display Compartments
M. Integration Control Information
1. Integration Option
2. Time Warp Scale Factor
3. Integration Interval and Total Time
N. Display Parameters
1. Variable Names
2. Dump Parameters
3. Printer Plot Parameters (Time History) Cards 1 and 2 are
read for each system; etc.
4. Printer Plot Parameters (Spatial Profile)
5. Pen Plot Parameters
51
-------
Card Group A.
Model Identification and System Bypass Options
5 10 15 20 25 30 35
MODEL ISER IRUN NOSEG NOSYS LISTG LISTC
FORMAT (715)
MODEL = model designation.
ISER = series designation.
IRUN = run number.
NOSEG = number of model segments.
NOSYS = number of systems, here always = 2.
LISTG = 0, print input data for exchange coefficients, volumes,
flows, and boundary conditions on the principal output
device.
= 1, suppress printing of input data for exchange coeffi-
cients, volumes, flows, and boundary conditions.
LISTC = 0, print input data for forcing functions, segment para-
meters, constants, miscellaneous time functions, and initial
conditions on the principal output device.
= 1, suppress printing of input for forcing functions, segment
parameters, constants, miscellaneous time functions, and
initial conditions.
MODEL, ISER, IRUN, although not actually used by the WASP program, can
assist the user in maintaining a log of computer simulations.
4.1.2 Title
80
VERIFICATION OF AUGUST 1973 RIVER SURVEY, REACH 1
FORMAT (20A4)
Card column 1-80 contain any information the user feels would be helpful
in describing the run and identifying the output for later reference.
4.1.3 Simulation Option
Presently WASP permits only time variable simulations, therefore include
the following card.
1 24
TIME VARIABLE SIMULATION
FORMAT (6A4)
52
-------
4.1.4 Systems Bypass Options
2 4
SYSBY(l) SYSBY(2) ... SYSBY(NOSYS)
FORMAT (1912)
SYSBY(K) = 0, perform the kinetic and transport phenomena associated
with system K (numerically integrate the differential
equations).
= 1, bypass all kinetic and transport phenomena associated
with system K (concentrations read as initial conditions
for system K apply throughout simulation period).
4.2 CARD GROUP B - EXCHANGE COEFFICIENTS
Exchange coefficients in TOXIWASP may be inputted in only one of WASP's
six formats—coefficients calculated from inputted dispersion coefficients
and accompanying cross-sectional areas and characteristic lengths.
4.2.1 Data Input Option Number; Number of Exchange Coefficients
5 10
IROPT NOR
FORMAT(2I5)
Data input options:
IROPT = 4, constant exchange coefficients calculated from the
dispersion, coefficient, cross-sectional area, and
characteristic lengths specified for each interface.
NOR = number of exchange coefficients.
If no exchange coefficients are to be read, set NOR equal to zero, and
continue with Card Group C.
4.2.2 Scale Factor for Exchange Coefficients
10
SCALR
FORMAT (E10.3)
SCALR = scale factor for exchange coefficients. This can be used
for units conversion or for sensitivity analysis.
53
-------
4.2.3 Exchange Coefficients
The data input format is determined by the option selected.
Option 4
Each card in this package contains the information to calculate the
exchange coefficients for two interfaces. The number of dispersion co-
efficients is equal to NOR. The information on each card is described
below:
10
E(K)
50
T7/V'_|_T ^
I-* \ Jx> J- J
20
A(K)
60
A(K+1)
25
IL(K)
65
IL(K+1)
30
JL(K)
70
JL(K+1)
35
IR(K)
75
TD f VJ-1 ^
J- R \ R. i J- /
40
JR(K)
80
JR(K+1)
E(K)
FORMAT (2(2F10.0, 2F5.0, 215))
= dispersion coefficient for the interface between segment
IR(K) and JR(K) in square miles/day. If exchange is
between bed and water column, the units are square
centimeters/sec.
= the interfacial cross-sectional area between segments 1R(K)
and JR(K), in square feet.
= the length of segment IR(K), with respect to the IL(K)-
JL(K) interface, in feet.
the length of segment JR(K) in relation to the IR(K)-
JR(K) 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 segment with which it is
exchanging.
IR(K),JR(K) = segments between which exchange takes place. NOTE: for
exchange only, order is not important — if a segment ex-
changes with a boundary, the boundary is specified as zero.
4.2.4 Exchange Bypass Options
2 4
A(K)
IL(K)
JL(K)
RBY(l) RBY(2) RBY(NOSYS)
FORMAT (1912)
RBY(K) = 0, exchange phenomena occurs in system K.
= 1, bypass exchange phenomena for system K (effectively set
for all exchange coefficients equal to zero for system K).
54
-------
4.3 CARD GROUP C - VOLUMES
4.3.1 Data Input Option Number; Number of Volumes
5 10
IVOPT NOV
FORMAT (215)
Data input options:
IVOPT = 1, constant volumes.
= 4, volumes from RECEIV hydrodynamic model read in each time
step through unit 92 (SWFLOW.DAT). IQOPT in Card Group D must
also be 4.
NOV = number of volumes; normally NOV is equal to NOSEG, the
number of segments, but for some special input structures,
NOV need not equal NOSEG.
4.3.2 Scale Factor for Volumes
10
SCALV
FORMAT (E10.3)
SCALV = scale factor for volumes; volumes are normally expressed in
units of million cubic feet. If other units are necessi-
tated by alteration in the space scale, SCALV should contain
the appropriate conversion factor; if normal units are em-
ployed, SCALV =1.0.
4.3.3 Volumes
The data input format is determined by the option selected. If IVOPT =
4, these volumes are not entered.
Option 1
Each card in this package contains the volume information for eight seg-
ments. The number of volumes is equal to NOV. The information on each card
is described below.
10 20 30 70 80^
VOL(K) VOL(K+1) VOL(K+2) .... VOL(K+6) VOL(K+7)
FORMAT (8F10.0)
VOL(K) = volumes of segment K, in million cubic feet. The volumes
are to be input consecutively, starting with segment 1, and
ending with segment NOV.
55
-------
4.4 CARD GROUP D - FLOWS
4.4.1 Data Input Option Number; Number of Flows
5 10
IQOPT NOQ
FORMAT (215)
Data Input Options:
IQOPT = 1, constant flows.
= 2, all flows proportional to one piecewise linear approxi-
mation.
= 3, each flow is represented by its own piecewise linear
approximation.
= 4, flows from RECEIV hydrodynamic model read in each time step
through unit 92 (SWFLOW.DAT).
NOQ = number of flows.
If no flows are to be inputted, set NOQ to zero, and go to Card Group E.
4.4.2 Scale Factor for Flows
10
SCALQ
FORMAT (E10.3)
SCALQ = scale factor for flows, flows are normally read in cubic
feet per second.
4.4.3 Flows
The data input format is determined by the option selected. If IQOPT =
4, no flows are read (skip to 4.4.4).
Option 1—
Each card in this package contains the flow information for four inter-
faces, the number of flow specifications is equal to NOQ. The information on
each card is described below.
10 15 20 30 35 40
BQ(K) JQ(K) IQ(K) BQ(K+1) JQ(K+1)
50 55 60 70 75 80
BQ(K+2) JQ(K+2) IQ(K+2) BQ(K+3) JQ(K+3) IQ(K+3)'
FORMAT (4(F10.0, 215)
56
-------
BQ(K) = flow between segment JQ(K) and IQ(K) in cfs. AESOP conven-
tion is: if the flow value is positive, then flow is from
segment JQ(K) to IQ(K).
JQ(K) = upstream segment.
1Q(K) = downstream segment.
If flow is from a segment to a boundary, then IQ(K) is set equal to
zero; if a flow is from a boundary to a segment, then JQ(K) is set equal to
zero.
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 Flows
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 information on each card is described below:
10 15 20 30 35 40
BQ(K) JQ(K) IQ(K) BQ(K+1) JQ(K+1)
50 55 60 70 75 80
BQ(K+2) JQ(K+2) IQ(K+2) BQ(K+3) JQ(K+3) IQ(K+3)
FORMAT (4(F10.0, 215))
BQ(K) = ratio of the flow between segments JQ(K) and IQ(K) to the
piecewise linear flow approximation.
JQ(K) = upstream segment.
IQ(K) = downstream segment.
If flow is from a segment to a boundary, then JQ(K) is set equal to
zero; if a flow is from a boundary to a segment, then IQ(K) is set equal to
zero.
Sub-Package II - Piecewise Linear Flow—The number of breaks required to
describe the piecewise linear approximation is followed by a time series
describing the piecewise linear flow approximation. Each time series element
consists of two parts; a flow and a time. The input is as follows.
10
NOBRK
FORMAT (15)
57
-------
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 ? 80
_ _. _ _
QT(K) T(K) QT(K+1) T(K+1) QT(K+2) T(K+2) QT(KH-3) T(K+3)
FORMAT (8F10.0)
QT(K) = value of the piecewise linear approximation at time T(K) , in
cubic feet per second.
T(K) = time in days, if the length of the simulation exceeds
T(NOBRK), the broken line function will repeat itself,
starting at time T(l), i.e., the approximation is assumed to
be periodic, with period equal to T(NOBRK).
Option 3 —
Each flow is defined by a package of cards consisting of two sub-pack-
ages. The first sub-package identifies the two segments between which the
flow occurs, and the number of values comprising the piecewise linear flow
approximation. The second sub-package defines the piecewise linear approxi-
mation which describes the flow. The input is as follows.
Sub-Package I —
5 10 15
JQ(K) IQ(K) NOBRK
FORMAT (315)
JQ(K) = upstream segment, flow from segment JQ(K) to IQ(K), assuming
positive flow.
IQ(K) = downstream segment flow from segment IQ(K) , 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,
and all breaks must occur at the same time relative to one
another.
Sub-Package II —
Sub-package II is a time series describing the piecewise linear
approximation. Each time series element consists of two parts: a flow and a
time. The input is as follows.
10 20 _ 30 _ 40 _ 50 _ 60 _ 70 _ 80
QT(K) T(K) QT(K+1) T(K+1) QT(K+2) T(K+2) QT(K+3) T(K+3)
FORMAT (8F10.0)
QT(K) = value of the piecewise linear flow approximation at time T(K)
in cubic feet per second.
58
-------
T(K) = time in days, if the length of the simulation exceeds T(NOBRK)
the broken line function will repeat itself, starting at time,
T(l). All break times must agree for all flows, T(2) must be
the same for all flows, etc.
4.4.4 Flow Bypass Options
The flow bypass options permit the flow transport to be set equal to
zero in one or more systems, while maintaining the flow regime (as defined
by one of the above options) in the remaining systems.
QBY(l) QBY(2) .... QBY(ISYS)
FORMAT (1912)
QBY(K) = 0, flow transport occurs in system K.
= 1, bypass the flow transport for system K (effectively set all
flows equal to zero in system K).
4.5 CARD GROUP E - BOUNDARY CONDITION
All input is read NOSYS times; once for each system of the model.
4.5.1 Data Input Option Number; Number of Boundary Conditions
5 10
IBCOP(K) NOBC(K)
FORMAT (215)
Data Input Options:
IBCOP(K) = 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(K) = number of boundary conditions used for system (K).
If no boundary conditions are to be inputted, set NOBC(K) equal to zero,
and continue with the next system, or go to the next card group.
4.5.2 Scale Factor for Boundary Conditions
10
SCALE
FORMAT (E10.3)
59
-------
SCALE = scale factor for boundary conditions. Boundary conditions are
normally expressed as milligrams per liter (mg/1), or parts
per million parts (ppm).
4.5.3 Boundary Conditions
The data input format is determined by the option selected.
Option 1—
10 15 25 30 40
BBC(K) IBC(K) BBC(K+1) 1BC(K+1) BBC(K+2)
45 55 60 70 7_5
IBC(K+2) BBC(K+3) IBC(K+3) BBC(K+4) IBC(K+4)
FORMAT (5(F10.0, 15))
BBC(K) = boundary condition of segment IBC(K) in mg/1.
IBC(K) = segment number to which boundary condition BBC(K) is to be
applied.
Option 2—
The card package consists of two sub-packages. Sub-package I contains
the boundary condition data, while sub-package II contains a detailed speci-
fication of the piecewise linear approximation to which all the 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 information on each card is described
below.
10 15 25 30 40
BBC(K) IBC(K) BBC(K+1) IBC(K+1) BBC(K+2)
45 55 60 70 7_5
IBC(K+2) BBC(K+3) IBC(K+3) BBC(K+4) IBC(K+4)
FORMAT (5(F10.0, 15))
BBC(K) = ratio of the boundary condition for segment IBC(K) to the
piecewise linear approximation.
IBC(K) = segment number.
Sub-Package II—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 concentration, and a time. The input is as follows.
60
-------
NOBRK
FORMAT (15)
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 70 80
BCT(K) T(K) BCT(K+1) T(K+1) BCT(K+2) T(K+2) BCT(K+3) T(K+3)
FORMAT (8F10.0)
BCT(K) = value of the broken line approximation at time T(K) in mg/1.
T(K) = 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 as-
sumed to be periodic, with period equal to T(NOBRK).
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 piece-
wise linear approximation. The second sub-package defines the piecewise
linear approximation which describes the boundary conditions. All boundary
conditions within a system must have the same number of breaks. The input is
as follows.
Sub-Package I—
5 10
IBC(K) NOBRK(K)
FORMAT (215)
IBC(K) = boundary segment number.
NOBRK(K) = number of values and times used to describe the broken line
approximation. The number of breaks must be equal for all
boundary conditions within a system.
Sub-Package II—The segment number and the number of breaks required to
describe the broken line approximation is followed by a time series describ-
ing the broken line approximation. Each time series element consists of two
parts: a boundary 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 input is as follows.
61
-------
10 20 30 40 50 60 70 8J3
BCT(K) T(K) BCT(K+1) T(K+1) BCT(K+2) T(K+2) BCT(K+3) T(K+3)
FORMAT (8F10.0)
BCT(K) = value of the boundary approximation at time T(K) in milligrams
per liter,
T(K) = 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 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 boundary
conditions, T(2) must be the same for all boundary conditions,
etc.
4.6 CARD GROUP F - FORCING FUNCTIONS
All input is read NOSYS times, once for each system of the model.
4.6.1 Data Input Option Number; Number of Forcing Functions
5 10
IWKOP(ISYS) NOWK(ISYS)
FORMAT (215)
Data Input Options:
IWKOP(ISYS) = 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(ISYS) = number of forcing functions used for system IS?S. NOTE:
forcing functions may also be considered as sources (loads)
or sinks of a water quality constituent. If no forcing
functions are to be inputted, set NOWK(ISYS) to zero, and
continue with next system or go to next card group.
4.6.2 Scale Factor for Forcing Functions
10
SCALW
FORMAT (E10.3)
SCALW = scale factor for forcing functions. Forcing functions are
normally read as pounds per day.
62
-------
4.6.3 Forcing Functions
The data input format is determined by the option selected.
Option I
10 15 25 30 40
BWK(K) IWK(K) BWK(K+1) IWK(K+1) BWK(K+2)
45
55
60
70
75
BWK(K)
IWK(K)
Option 2
IWK(K+2) BWK(K+3) IWK(K+3) BWK(K+4) IWK(K+4)
FORMAT (5(F10.0, 15))
= forcing function of segment IWK(K), in pounds per day.
= segment number to which forcing function BWK(K) is to be
applied.
The card package consists of two sub-packages. Sub-package I contains
the forcing function data, while sub-package II contains a detailed specifi-
cation of the piecewise linear approximation to which all the forcing func-
tions 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. The information on each card is described below:
10
15
25
30
40
BWK(K) IWK(K) BWK(K+1) IWK(K-fl) BWR(K+2)
45
55
60
70
75
BWK(K)
IWK(K)
IWK(K+2) BWK(K+3) IWK(K+3) BWK(K+4) IWK(K+4)
FORMAT (5(F10.0, 15))
= ratio of the forcing function for segment IWK(K) to the
piecewise linear approximation.
= segment number to which forcing function BWK(K) is to be
applied.
Sub-Package II—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 input is as follows.
NOBRK
FORMAT (15)
63
-------
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 70_ 8£
WKT(K) T(K) WKT(K+1) T(K+1) WKT(K+2) T(K+2) WKT(K+3) T(K+3)
FORMAT (8F10.0)
WKT(K) = value of the forcing function at time T(K), in pounds per day.
T(K) = time in days, if the length of the simulation exceeds
T(NOBRK), the forcing function approximation is repeated,
starting at 1(1), i.e., the approximation is assumed to be
periodic, with period equal to T(NOBRK).
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 approximation which
describes the forcing function. The input is as follows.
Sub-Package I—
5 10
IWK(K) NOBRK(K)
FORMAT (215)
IWK(K) = segment number which has forcing function BWK(K).
NOBRK(K) = number of breaks used to describe the forcing function
approximation. The number of breaks must be equal for all
forcing functions within a system.
Sub-Package II—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 input is as
follows:
10 20 30 40 50 60 70 80
WKT(K) T(K) WKT(K+1) T(K) WKT(K+2) T(K+2) WKT(K+3) T(K+3)
FORMAT (8F10.0)
WKT(K) = value of the forcing function at time T(K), in pounds per day.
-------
T(K) = time in days, if the length of the simulation exceeds
T(NOBRK), the approximation is repeated, starting at T(l),
i.e., 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.
4.7 CARD GROUP G - PARAMETERS
The definition of the parameters will vary, depending upon the structure
and kinetics of the systems comprising each model. The input format however
is constant and is detailed below.
4.7.1 Number of Parameters
NOPAM
FORMAT (15)
NOPAM = number of parameters required by TOXIWASP = 18.
4.7.2 Scale Factors for Parameters
10 20 30
SCALP(l) SCALP(2) SCALP(3) ... SCALP(NOPAM)
FORMAT (8E10.3)
SCALP(K) = scale factor for parameter group K.
4.7.3 Segment Parameters
Segment parameters (Numbers on top line are column numbers for input.
*.* denotes values for parameters on left.)
5
TEMPM
BACTO
OXRAD
WS
TEMPM
15 20
* . *DEPTH
*.*ACBAC
*.*OCS
*.*CMPET
*.*DEPTH
30 35 45 50 60 65 75
*.*VELOC *.*WINDS *.*TYPEE *.*
*.*BIOMS *.*BIOTM *.* POHG *.*
*.*PCTWA *.*DSPSD *.* PHG *.*
*.*TOTKG *.*
FORMAT (5(A5,F10.0))
65
-------
Note that only 18 parameters are input here. DISPV is calculated
internally and occupies two spaces for each segment. Although contained in
the glossary section, a brief description of parameters and units follows.
Fortran
Parameter Variable
Description
Average temperature for segment.
1 TEMPM(J)
2 DEPTHG(J) Depth of segment.
3 VELOC(J) Average velocity of water in
segment.
4 WINDG(J) Average wind velocity 10 cm
above the water surface (sur-
face water segments only).
5 TYPEE(J) Flag designating segment type:
1 = surface water, 2 = subsurface
water, 3 = surface bed, 4 =
subsurface bed.
6 BACTOG(J) Bacterial population density in
segment.
ACBACG(J) Proportion of bacterial popula-
tion that actively degrades
chemical.
BIOMAS(J) Total actively sorbing biomass
in segment.
9 BIOTMG(J) Biotemperature in segment.
10 POHG(J) Hydroxide ion activity in
segment.
11 OXRADG(J) Molar concentration of
environmental oxidants in
segment.
12 OCS(J) Organic carbon content of sedi-
ments as fraction of dry weight.
Units
degrees C
feet
feet per second
meters per second
cells per milli-
liter (water)
cells per 100 grams
dry weight (bed)
dimensionless ratio
mg (dry weight) per
liter (water)
grams (dry weight)
per square meter
(bed)
degrees C
pOH units
moles per liter
dimensionless
66
-------
13
PCTWA(J)
14
15
16
DSPSED(J)
PHG(J)
WS(J)
17
18
CMPETG(J)
TOTKG(J)
Percent water in benthic sedi-
ments, expressed as fresh/dry
weight; all values must be
greater than or equal to 1.
Fraction of sediment volume
that mixes.
dimensionless
dimensionless
Hydrogen ion activity in segment. pH units
Spatially variable parameter
denoting settling rate of
suspended sediment in water
column segments (types 1 and 2),
erosion rate of surface bed
sediment in surface bed segment
(type 3), or percolation of pore
water through subsurface bed
segments (type 4).
Single-valued zenith light
extinction coefficients (for
water segments only).
Total first order decay rates
calculated externally. If
equal to zero, the program
evaluates all separate processes
and calculates a combined total
first order decay rate internally.
meters per day
(1, 2)
cm per year (3)
cubic feet per
second (4)
per meter
per day
4.8 CARD GROUP H - CONSTANTS
4.8.1 Number of Constants
NCONS
FORMAT (15)
NCONS = number of constants required by the model
66.
67
-------
4.8.2 Constants
5
15 20
30 35
45 50
60 65
75
EBHG1
ENHG3
KAHG2
KNHG1
EOXG3
KBWG2
KB SGI
QTBS3
KVOG
CLOUD
QUAN3
LOPT
*.*EBHG2
*.*EAHG1
*.*KAHG3
*.*KNHG2
*.*KOXG1
*.*KBWG3
*.*KBSG2
*.* KOC
*.* SOLG
*.* LATG
*.* XJTR
*.*
*.*EBHG3
*.*EAHG2
*.*KBHG1
*.*KNHG3
*.*KOXG2
*.*QTBWl
*.*KBSG3
*.* KOW
*.* MWTG
*.*ESOLG
FAC
*.*DFACG
*.*CTRIG
*.*ENHG1
*.*EAHG3
*.*KBHG2
*.*EOXG1
*.*KOXG3
* . *QTBW2
*.*QTBS1
*.* OCB
*.*HENRY
*.*EVPRG
*.* KDPG
*.*QUAN1
*.*DTOPT
*.*ENHG2
*.*KAHG1
*.*KBHG3
* . *EOXG2
*.*KBWG1
*.*QTBW3
*.*QTBS2
*.*
*.*VAPRG
* . *EHENG
*.*RFLAT
* . *QUAN2
*.*TDINT
*.*
*.*
*.*
*.*
*.*
*.*
*.*
*.*
*.*
*.*
*.*
*.*
Note that there are additional elements in the constant array but that
they are internally calculated and not read-in. Numbers 1, 2, and 3 refer to
dissolved, sediment sorbed, and biosorbed chemical, respectively. KNHG2, for
example, is the neutral hydrolysis rate for the sediment-sorbed chemical
fraction. Although contained in the glossary section, a brief description of
constants and units follows.
Fortran
Constant Variable
1,2,3
4,5,6
EBHG(I,1)
ENHG(I.l)
Description
Arrhenius activation energy of
specific-base-catalyzed hydro-
lysis of the toxicant.
Arrhenius activation energy of
neutral hydrolysis of the
toxicant.
Units
kcal/gram mole
kcal/gram mole
68
-------
7,8,9 EAHG(I,1)
10,11,12 KAHG(I.l)
13,14,15 KBHG(I.l)
16,17,18 KNHG(I.l)
19,20,21 EOXG(I,1)
22,23,24 KOXG(I,1)
25,26,27 KBACWG(I,1)
28,29,30 QTBAWG(I,1)
31,32,33 KBACSG(I,1)
34,35,36 QTBASG(I,1)
37
38
KOC
KOW
Arrhenius activation energy of
specific-acid-catalyzed hydro-
lysis of the toxicant.
Second-order rate constants for
specific-acid-catalyzed hydro-
lysis of chemical.
Second-order rate constants for
specific-base-catalyzed hydro-
lysis of chemical.
Rate constants for neutral hydro-
lysis of organic chemical.
Arrhenius activation energy of
oxidative transformation of the
chemical.
Second-order rate constants for
oxidative transformation of
toxicant.
Second-order rate constants for
water column bacterial biolysis
of the organic chemical.
Q-10 values for bacterial trans-
formation rate in the water
column. Q—10 is the increase
in the second-order rate constant
resulting from a 10 degree C
temperature increase.
Second-order rate constants for
benthic sediment bacterial
biolysis of the organic.
Q-10 values for bacterial trans-
formation of organic chemical in
benthic sediments. The Q-10 is
the increase in the second-order
rate constant resulting from a 10
degree C temperature increase.
Organic carbon partition
coefficient.
Octanol water partition
coefficient
kcal/gram mole
per mole[H+]
per hour
per mole [OH-]
per hour
per hour
kcal/gram mole
liter per mole
environmental
oxidant per
hour
ml/cell-hour
dimensionless
ml/cell-hour
dimensionless
lw/kg organic
carbon
lw/loct
69
-------
39
43
44
45
46
47
48
49
50
53
54
55
56
OCB
MWTG
HENRYG
VAPRG
KVOG
SOLG
ESOLG
EVPRG
EHENG
FAC
KDPG
RFLATG
CLOUDG
dimensiionless
grams per mole
Atmosphere-
cubic meters
per mole
torr
dimensionless
Organic carbon content of the
compartment biomass as a fraction
of dry weight.
The molecular weight of the
chemical.
Henry's Law constant of the
toxicant.
Vapor pressure of compound.
Measured experimental value for
(volatilization) liquid-phase
transport resistance, expressed
as a ratio to the reaeration rate.
Aqueous solubility of toxicant
chemical species.
Exponential term for describing
solubility of the toxicant as a
function of temperature (see SOLG)
Molar heat of vaporization for
vapor pressure described as a
function of temperature (See
VAPRG).
Constant used to compute Henry's
Law constants for volatilization
as a function of environmental
temperatures (TCELG). When EHENG
is non-zero, the Henry's Law
constant is computed as follows:
log HENRY = HENRYG-((1000.*EHENG)/(4.58*(TCELG+273.15)))
Multiplication factor for sedimen- dimensionless
tation time step. Recommended
0.1.
mg/lT
kcal/gram mole
kcal/gram mole
kcal/gram mole
A near-surface photolytic rate
constant for the chemical.
Reference latitude for correspond-
ing direct photolysis rate constant
KDPG.
Average cloudiness in tenths of
full sky cover.
per hour
degrees and deci-
mal fraction (e.g.,
40.72)
dimensionless, range
of 0.0 to 10.0
70
-------
57
58
LATG
DFACG
Geographic latitude of ecosystem.
Distribution function (ratio of
optical path length to vertical
depth).
59,60,61 QUANTG(I,1) Reaction quantum yield in photo-
lytic transformation of chemical.
62
63
64
XJTR
CTRIG
DTOPT
Reference segment for triggering
event and frequency output. A
value of zero will disable it.
Trigger concentration that de-
fines a peak event.
Option to optimize time step,
if set to 1.
degrees and
tenths (e.g.
37.2)
dimensionless
dimensionless
dimensionless
mgc/lT
dimensionless
65
66
TD1NT
LOPT
Time interval between recalcula-
tion of decay rates.
Option to read in flows and/or
nonpoint source loads on a daily
basis:
days
dimensionless
0 = skip daily read option
1 = read sequential tape
containing daily flows
and/or loads
4.9 CARD GROUP I - MISCELLANEOUS TIME FUNCTIONS
The definition of the miscellaneous piecewise linear time functions will
vary depending upon the structure and the kinetics of the systems comprising
each model. The input format, however, is constant and is detailed below.
4.9.1 Number of Time Functions
5
NFUNC
FORMAT (15)
NFUNC
number of time functions required by the model = 5.
If no time functions are to be inputted, set NFUNC equal to zero, and go
to Card Group K.
71
-------
4.9.2 Time Functions
The following package of cards is required for each time function. The
first sub-package defines the function name and the number of breaks in the
time function. The second sub-package contains a detailed specification of
the piecewise linear time function.
Sub-Package _!_
5 10
ANAME(K) NOBRK(K)
FORMAT (A5, 15)
ANAME(K) = an optional one to five alpha-numeric character descriptive
name for the time function K.
NOBRK(K) = number of breaks used to describe the time function K.
Sub-Package II
10 20 30 40
VALT(K) T(K) VALT(K+1) T(K+1)
50 60 70 80
VALT(K+2) T(K+2) VALT(K+3) T(K+3)
FORMAT (8F10.0)
VALT(K) = value of the function at time T(K).
T(K) = 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.
The five time-variable functions used in TOXIWASP are listed below.
These functions will usually describe seasonal variability for long
simulations. They can be set up, however, to describe diurnal variability
over short simulations if the time step is sufficiently short (1-3 hours).
Time Fortran
Function Variable Description Units
1 TEMPN Normalized temperature dimensionless
2 WINDN Normalized wind speed dimemsionless
3 PHN Normalized pH dimemsionless
4 POHN Normalized pOH dimensionless
5 LIGHTN Normalized light intensity dimensionless
72
-------
4.10 CARD GROUP J - INITIAL CONDITIONS
The initial conditions are the segment concentrations for the state
variables at time zero (or the start of the simulation).
4.10.1 Initial Conditions
5 15 20 3£
ANAME(K) C(ISYS,K) ANAME(K+1) C(ISYS,K+1)
35 45 50 6£
ANAME(K+2) C(lSYS,K+2) ANAME(K+3) C(lSYS,K+3)
FORMAT (4(A5, F10.0))
ANAME(K) = an optional one to five alpha-numeric character descriptive
name for the initial condition in segment K of system ISYS.
C(ISYS,K) = initial condition in segment K of system ISYS in the
appropriate units (normally mg/1 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 the concentrations being read from segment 1 through
NOSEG within a system packet.
4.11 CARD GROUP K - STABILITY AND ACCURACY CRITERIA
4.11.1 Stability Criteria
10 20 30 80
CMAX(K) CMAX(K+1) CMAX(K+2) .... CMAX(NOSYS)
FORMAT (8F10.0)
CMAX(K) = stability criteria for system K, i.e., the maximum concentra-
tion (normal units mg/1 or ppm) for system K, which if exceeded
by any segments in system K indicates that the numerical inte-
gration procedure has become unstable. If unstability occurs
an appropriate message is printed and the integration procedure
is terminated and a call is made to the display subroutines.
If the user sets CMAX (1) to zero, TOXIWASP will reset it to one half of
the solubility or 10~5 M, whichever is less. This would abort runs for which
the first-order rate assumption becomes invalid.
4.11.2 Accuracy Criteria
10 20 30 8JD
CMIN(K) CMIN(K+1) CMIN(K+2) .... CMIN(NOSYS)
FORMAT (8F10.0)
73
-------
CMIN(K) = originally WASP read the accuracy criteria for system K; i.e.,
for time variable simulations the minimum concentration (normal
units mg/1 or ppm) that governs integration step-size control,
if the user chooses option 1. As it is not recommended that
the user use this option due to numerical difficulties, just
set CMIN = 0.0 for each system.
4.12 CARD GROUP L - INTERMEDIATE PRINT CONTROL
4.12.1 Print Interval
FORMAT (F10.0)
PRINT = print interval in days for time-variable applications. NOTE:
The maximum number of print outs = total time/print interval +
1 (for time zero) must be equal to or less than 41.
4.12.2 Compartments (System - Segment) to be Displayed
ISYS(l) ISEG(l) ISYS(2) ISEG(2) ... ISYS(8) ISEG(8)
FORMAT (1613)
ISYS(K), = system, segment combinations user wishes to have displayed
ISEG(K) during simulation - user may select a maximum of 8. NOTE: All
system-segment concentrations as well as other miscellaneous
calculations may be displayed at the end of the simulation, see
Card Group N.
4.13 CARD GROUP M - INTEGRATION CONTROL INFORMATION
4.13.1 Integration Option - Negative Solution Option
2 4
INTYP NEGSLN
FORMAT (212)
INTYP = 1, user wishes the WASP program to determine the integration
step size (based upon its own accuracy criteria). It is re-
commended that the user not use this option.
= 2, the user will supply the integration step sizes to be used
by WASP. It is recommended that the user utilize this option.
NEGSLN = 0, a user wishes to restrict integration to the positive plane
only - this is the normal option selected.
74
-------
= 1, user will permit the integration procedure to go negative -
used for special applications (ex., DO deficit, pH -
alkalinity).
4.13.2 Time Warp Scale Factor - Starting Simulation Time (TZERO)
10 20
SCALT TZERO
FORMAT (2E10.4)
SCALT = time warp scale factor - allows the user to completely change
the time scale of his simulation via just one card. This scale
factor changes all times employed in piecewise linear functions
or piecewise linear approximations for volumes, exchanges,
flows, etc., as well as print out time.
TZERO = prototype time for start of simulation, usually equal to zero,
but user may start at time other than zero (used to initialize
any of the piecewise linear time functions).
4.13.3 Number of Integration Step Sizes
5
NOSTEP
FORMAT (15)
NOSTEP = number of integration step sizes to be used in the simulation.
4.13.4 Integration Step Size History
10 20 30 40
DT(K) TIME(K) DT(K+1) TIME(K+1)
FORMAT (8F10.0)
DT(K)
TIME(K)
integration step size (normal units-days).
time until which step size DT(K) will be used, then switching
to DT(K+1) until TIME(K+1).
4.14 CARD GROUP N - DISPLAY PARAMETERS
The cards presented for this group are those required by the normal WASP
display package. A special subroutine that permits off-line digital pen
plotting capability is available for the DEC POP system. The input data re-
quired are described here.
The following two sub-groups are read for each system, starting with
system 1 and running through system NOSYS. For each variable-segment com-
bination chosen a time history of the segment will be displayed (dumped).
75
-------
4.14.1 Variable Names
8 16 80
ANAME(l) ANAME(2) .... ANAME(IO)
FORMAT (10A8)
ANAME(K) = a one to eight alpha-numeric character descriptive name for
display variable K. The order of these names is determined
via the appropriate disk file WRITE in the users kinetic sub-
routine. NOTE: In TOXIWASP, the following variables are
dumped: for system 1, CHEM, CHEM1, CHEMS, CHEMB, ALPHA(l),
ALPHA(2), XMASS, and BMASS; for system 2, SED, DEPTHG, TOTKL,
PHOTKL, HYDRKL, BIOLKL, OXIDKL, VOLKL.
4.14.2 Variable Number, Segment Numbers
369 27
VARNO SEG(K) SEG(K+1) .... SEG(K+7)
FORMAT (913)
VARNO = the position of the desired display variable in the WRITE file
statement in the kinetic subroutine (see previous note).
VARNO = 1 to 8.
SEG(K) = segment number to be displayed. NOTE: Order of display
unimportant, i.e. , need not be sequential.
A blank card terminates display for system, ISYS. Then another Variable
Name card, followed by Variable Number, Segment Number card(s) Is read until
system NOSYS has been read, then the plot cards will be read.
4.14.3 Printer Plot Display Cards
The following cards are read for each system, starting with system 1,
and running through NOSYS. The printer plot display sub-group requires 3
cards per plot, and because two plots are formulated per page of printout,
the user should input an even number of plots.
3a. Number of Segments and Variable Number for This Plot
2 4
NSPLT VARNO
FORMAT (212)
NSPLT = number of segments to be plotted (maximum of five.).
VARNO = the position of the desired variable to be plotted, in the
WRITE file statement in the kinetic subroutine.
76
-------
3b. Plotting Scales
10 20
PMIN PMAX
FORMAT (2F10.0)
PMIN, = minimum and maximum values, respectively, to be used for
PMAX this plot.
3c. Segment to be Plotted
36 9
SEG(l) SEG(2) .... SEG(INSPLT)
FORMAT (513)
SEG(K) = segment numbers to be plotted (a maximum of 5 segments per plot
is allowed).
A blank card terminates the plotting for system, ISYS.
At this point, unless the user has written additional display routines
(such as the PDP graphics), WASP input is finished, and the user should end
his deck with the appropriate end of data set indicator.
4.14.4 Plot Display Cards (Spatial Profile)
Cards 4a and 4b are read-in once, and apply to all spatial plots. There
are two types of spatial plots. The first type plots predicted variables,
and is controlled by Cards 4c and 4d. The second type plots observed data,
and is controlled by Cards 4e, 4f, and 4g. Any number of plots desired can
be produced by repeating card combinations 4c-4d and/or 4e-4g.
4a. Spatial Scale
5 10
RM1 RM2
FORMAT (2F5.0)
RM1, RM2 = minimum and maximum river mile values, respectively, to be used
for all spatial plots.
4b. Segment River Miles to be Plotted
5 10 15 20 25 30 35_
SEG(K) RM(K) SEG(K+1) RM(K+1) SEG(K+2) RM(K+2) SEG(K+3)
40 45 50 55 60 65^
RM(K+3) SEG(K+4) RM(K+4) SEG(K+5) RM(K+5) SEG(K+6)
70 75 80
RM(K+6) SEG(K+7) RM(K+7)
FORMAT (8(15, F5.0))
77
-------
SEG(K) = segment number to be plotted.
RM(K) = river mile value for SEG(K).
A maximum of "NOSEG" combinations of SEG-RM pairs are allowed.
For each spatial plot, Cards 4c-4d, or 4e-4g are read.
4c. Predicted Variable Plot Control Information
30
10
15
20
25
70
MXTIM IVAR YSTR YSTP SYSOPT OVRLAY
TITL1
FORMAT (215, 2F5.0, 215, 40A1)
MXTIM = number of time selections to be included on this plot (maximum
of 5).
IVAR = the position of the desired variable to be plotted in the WRITE
file statement in the kinetic subroutine.
YSTR,YSTP = minimum and maximum values, respectively, to be used for the
Y-axis of this plot.
SYSOPT = system number of the desired variable to be plotted.
OVRLAY = flag to cause this plot to be overlaid with the following
plots.
= 0, causes this plot to be printed alone (or with preceeding
plot, if OVRLAY on the preceding plot cards is set: to 1).
= 1, causes this plot to be overlaid on the following plot.
(Note, although any number of plots can be overlaid, we suggest
a maximum of three; YSTR and YSTP values should be compatible
for overlaid plots.)
TITL1 = title for plot; when overlaying plots, the first two titles and
the last title will be printed.
4d. Time Selections and Characters for Predicted Variable Plot
5 10
TIM(l) TIM(2)
27
15 20
TIM(3) TIM(4)
28 29
SYMTAB(2) SYMTAB(3) SYMTAB(4)
25 26
TIM(5) SYMTAB(l)
30
SYMTAB(5)
FORMAT (5F5.0, 5A1)
TIM(K) = time selections for this plot (1-MXTIM).
SYMTAB(K) = plot symbol associated with time TIM(K).
78
-------
4e. Observed Data Plot Control Information
10
15
20
25
30
70
71
FLAG IUNIT YSTR YSTP NOOBS OVRLAY TITL1 OBSSYM
FORMAT (215, 2F5.0, 215, 40A1, 1A1)
FLAG = flag to indicate observed data.
= 99999, indicates observed data are to be plotted.
IUNIT = unit device number where observed data are to be found (default
=5; optional unit numbers are 82-89).
YSTR,YSTP = minimum and maximum values, respectively, to be used for the
Y-axis of this plot.
NOOBS = number of observed data points for this plot.
OVRLAY = flag to cause this plot to be overlaid with the following plot:
= 0, causes this plot to be printed alone (or with preceding
plot, if OVRLAY on the preceding plot cards is 1).
= 1, causes this plot to be overlaid on the following plot.
TITL1 = title for this plot.
OBSSYM = plot symbol associated with observed data for this plot.
If "IUNIT" on Card 4e equals zero or five, Card 4f is read and Card 4g
is skipped. If "IUNIT" equals 82-89, then Card 4f is skipped and Card 4g is
read.
4f. River Mile - Observed Data Values
10
20
30
40
RIVMIL(K) VALUE(K) RIVMIL(K-fl) VALUE(K+1)
50 60 70 80
RIVMIL(K+2) VALUE(K+2) RIVMIL(K+3) VALUE(K+3)
FORMAT (8F10.0)
RIVMIL(K) = river mile location for observed data point "K".
VALUE(K) = observed value of variable at RIVMIL(K).
If "IUNIT" on Card 4e equals zero or five, Card 4g is skipped.
79
-------
4g. Format Specification for Data on Auxiliary Input File "IUNIT"
1 8£
FMT
FMT
FORMAT (20A4)
= format specification for observed river mile - observed data
values on auxiliary input file "IUNIT" (specified on Card 4e).
Must begin and end with parentheses, and contain valid formats,
such as (2F5.0), (16F5.0), or (F5.0/F5.0)
4.14.5 Pen Plot Display Cards
The following cards permit the user to obtain plots of WASP results on a
digital pen plotter. The pen plot software has been written such that six
"variable vs. time" are generated per 11 inch x 11 inch frame. Therefore,
the user must supply the following input data cards in multiples of six.
The user should note that unlike the Dump and Printer Plot Display Cards
which are organized by system, the Pen Plot Display Cards permit the user to
mix the various system outputs on the same frame.
5a. System, Segment, Variable, and Units for This Plot
5 10 15 16 29
SYS
SEG
IVAR
UNITS
SYS
SEG
IVAR
FORMAT (315, A14)
= system number of the desired variable to be plotted.
segment number of the desired variable to be plotted.
= the position of the desired variable, to be plotted, in the
WRITE file statement in the kinetic subroutine.
UNITS = an alphanumeric descriptor, which describes the units of the
variable to be plotted.
5b. Plotting Scales
10 20
YMIN
YMAX
FORMAT (2F10.0)
YMIN, YMAX = minimum and maximum values, respectively, to be used for this
plot.
5c. Field Data
The pen plot subroutine, WAS18, h*s been written so as to permit the
user to overplot theory with observed field data. Basically the input data
to be supplied for overplotting field data consists of four pieces of
80
-------
information: a survey number, the time (relative to time zero of the WASP
simulation) at which the data was collected, the value (mean value if several
samples were collected and are to be aggregated) of the particular water
quality parameter to be plotted, and the standard deviation of the data if
applicable. The input data is read as follows.
5 10 20 30 60 70 80
SURVEY SURVYT(l) Y(l) SD(1) SURVYTQ) Y(3) SD(3T
FORMAT (15, 3(F5.0, 2F10.0))
SURVEY = survey number for which the data were collected. Must be a
number from 1 to 5. In actuality SURVEY selects the
appropriate plotting symbol for plotting the field data.
SURVYT(I) = time, relative to time zero of the WASP simulation, at which
data were collected. Nominal units are days.
Y(I) = the value for the water quality parameter to be plotted. If a
number of samples were collected and are aggregated together
than Y(l) is the aggregated sample mean.
SD(I) = the standard deviation of the aggregated sample mean (if
applicable).
The following description, used together with Table 4-2, should demon-
strate to the user how to prepare data input for the plots. The user will be
plotting phytoplankton and secchi depth results for segment 1 through 3 of
his model. The kinetic subroutine was written such that the phytoplankton
were written to the first system file as the first variable (SYS = 1, IVAR =
1) and, the secchi depths were written to the second system file as the
fourth variable (SYS = 2, IVAR =4). Cards 1, 7, 16, 19, 22, and 27 specify
that phytoplankton and secchi depth are to be plotted for segments 1, 2, and
3, respectively. Cards 2, 8, 17, 20, 23 and 28 select the plotting scales
to be used for the appropriate variables. Cards 3 through 5 contain the
field data (times, means, and standard deviations) to be overplotted with
the theory for segment 1. Card 6 indicates that no more survey data was
collected for segment 1. Cards 9 through 14 contain the field data, from
three different surveys, to be overplotted with segment 2 theory. Note that
the data for the third survey are just grab samples (no associated standard
deviations). Card 15 indicates the end of survey data for segment 2. Card
17 contains grab sample data collected on the third survey; Card 18 indicates
no more survey data for segment 3.
No secchi depth field data was collected for segments 1 and 3; there-
fore the cards following cards 20 and 28 are blank. Cards 24 and 25 contain
single point measurements of secchi depth collected during survey 2. Card 26
indicates the end of survey data. Card 30 indicates that no more plots are
required.
81
-------
w
s
w
H
Z
w
PH
.
1
a
i
M
(U
|
w
o
<}• r~H
\O \^ *.O l-O "*^"
O O O O CD
in co co Kf -r-inr^-*
mrHCOO inrHCMOCMOrH
i-H rH
rH O CM O
omorH omorHooo
rHCOl-H rHcOrHrHCM««O
rH CO rH CO CM CO
i-t rHiHrH rH iHrHrHCSlCMCO
I-H CM co *d* in ^o r^ oo cjN o rH CM co -^
rH rH rH rH rH
CT> in
OO rH
o in
m m
O OO
CO CM
r- m
rH -
-------
SECTION 5
OPERATIONAL CONSIDERATIONS
This chapter describes how to obtain the computer program TOXIWASP,
how to install it on a DEC POP mini computer or an IBM mainframe, how to
test the program with sample datasets, and what machine limitations limit
the program.
5.1 ACQUISITION PROCEDURES
To obtain the program TOXIWASP along with sample datasets and support
software, write to:
Center for Water Quality Modeling
Environmental Research Laboratory
U.S. Environmental Protection Agency
College Station Road
Athens, GA 30613
Specify which version of the program you desire (POP or IBM) and the system
on which you will install the program. A nine-track magnetic tape will be
mailed to you. Please copy the contents and return the tape.
5.2. INSTALLATION PROCEDURES
Among the datasets on the magnetic tape are the standard WASP
subroutines and the special TOXIWASP subroutines. These must be compiled
and linked into a task image. This is accomplished on the PDF IAS
operating system by running the command file "TOXICMP.CMD," which is listed
in Appendix 3.2. If the compilation suceeds, then linkage is automatically
attempted with the command file "TOXITKB.CMD," which is listed in Appendix
3.2. The Job Control Language (JCL) necessary to accomplish these two
tasks on an IBM 370 OS/MVS operating system is listed in Appendix 3.3.
5.3. TESTING PROCEDURES
Once TOXIWASP is installed, the sample input datasets should be run
and compared with the sample output datasets to verify that the program is
calculating correctly. Data sets describing a sample pond, river, and
oligotrophic lake polluted with various chemicals are available on the
tape. Simulations are run in a batch mode. To perform a simulation on the
83
-------
PDF, submit the batch input sequence "TOXIC.BIS," which is listed in
Appendix 3.4. This BIS file is set up to run the test pond. To run the
test river, change the line
$COPY POND.INP WASP.INP
to
$COPY RIVER.INP WASP.INP
and resubmit "TOXIC.BIS." The test lake is titled LAKE.INP. Output
written to FOR006.DAT should be compared to POND.OUT, RIVER.OUT, and
LAKE.OUT, respectively.
JCL necessary to run simulations for the three test cases on an IBM
370 is listed in Appendix 3.5.
5.4. MACHINE LIMITATIONS
Currently, TOXIWASP is set up for the following configurations.
POP 11/70 Hardware IBM 370 Hardware
IAS operating system OS/MVS operating system
FORTRAN IV FORTRAN IV+
25 parameters 25 parameters
85 constants 85 constants
50 segments 100 segments
2 systems 2 systems
The PDP 11/70 computer utilizing an IAS operating system allocates a
32k word (64k byte) user area for execution of programs. TOXIWASP occupies
at least 31k words of memory. Any enlargement of this program may result
in an overflow of the user area.
A compromise in the size of the program arrays (regulated by the
PARAMETER statements at the beginning of the common area) will result in
more free area to add other applications. Care must be taken to assure
that any changes made to the code are completely researched and tested.
84
-------
REFERENCES
Burns, L. A., D. M. Cline, and R. R. Lassiter. 1982. Exposure Analysis
Modeling System (EXAMS): User Manual and System Documentation. U.S.
Environmental Protection Agency, Athens, GA. EPA-600/3-82-023.
Covar, A. P. 1976. "Selecting the Proper Reaeration Coefficient for Use
in Water Quality Models." Presented at the U.S. EPA Conference on
Environmental Simulation and Modeling, April 19-22, 1976.
DiToro, D. M. , J. J. Fitzpatrick, and R. V. Thomann. 1982. Water Quality
Analysis Simulation Program (WASP) and Model Verification Program
(MVP) - Documentation. Hydroscience, Inc., Westwood, NY, for U.S.
Environmental Protection Agency, Duluth, MN, Contract No. 68-01-3872.
Huber, W. C. , Heaney, J. P., Medina, M. A., Peltz, W. A., Sheikh, H. ,
and Smith, G. F., Storm Water Management Model User's Manual, Version
II, University of Florida, Gainesville, FL, for U.S. Environmental
Protection Agency, Cincinnati, OH, EPA-670/2-75-017 , March 1975.
Johanson, R. C. , J. C. Imhoff, and H. H. Davis, Jr. 1980. Users Manual
for Hydrological Simulation Program — FORTRAN (HSPF). U.S.
Environmental Protection Agency, Athens, GA. EPA-600/9-80-015.
Karickhoff, S. W. , D. S. Brown, and T. A. Scott. 1979. Sorption of
hydrophohic pollutants on natural sediments. Water Research, 13:
241-248.
Mackay, D. and P. J. Leinonen. 1975. Rate of evaporation of low
solubility contaminants from water bodies to atmosphere. Environ.
Sci. Technol.,
Mill, T. , W. R. Mabey, P. C. Bomberger, T. W. Chou, D. G. Hendry, and J.
H. Smith. 1982. Laboratory Protocols for Evaluating the Fate of
Organic Chemicals in Air and Water. U.S. Environmental Protection
Agency, Athens, GA. EPA-600/3-82-0220 (In press).
Onishi, Y. , S. M. Brown, A. R. Olsen, M. A. Parkhurst, S. E. Wise, and W.
H. Walters. 1982. Methodology for Overland and Instream Migration
and Risk Assessment of Pesticides. U.S. Environmental Protection
Agency, Athens, GA. EPA-600/3-82-0240 (In press).
Onishi, Y. and S. E. Wise. 1982. User's Manual for the Instream Sediment-
Contaminant Transport Model, SERATRA. U.S. Environmental Protection
Agency, Athens, GA. EPA-600/3-82-0550 (In press).
85
-------
Paris, D. F., D. L. Lewis, J. T. Barnett, and G. L. Baughman. 1975.
Microbial Degradation and Accumulation of Pesticides in Aquatic
System. U.S. Environmental Protection Agency, Washington, DC.
EPA-660/3-75-007.
Rao, P. S. C. and J. M. Davidson. 1980. Estimation of Pesticide
Retention and Transformation Parameters Required in Nonpoint Source
Pollution Models. In: Environmental Impact of Nonpoint Source
Pollution, M. R. Overcash and J. M. Davidson, Eds. Ann Arbor Science,
Ann Arbor, MI. pp. 23-67.
Richardson, W. L., V. E. Smith, and R. Wethington. 1982. Dynamic Mass
Balance of PCB and Suspended Solids in Saginaw Bay—A Case Study.
U.S. Environmental Protection Agency, Grosse lie, MI, In Press.
Smith, J. H. , W. R. Mabey, N. Bohonos, B. R. Holt, S. S. Lee,, T-W. Chou,
D. C. Bomberger, and T. Mill. 1977. Environmental Pathways of
Selected Chemicals in Freshwater Systems. Part I: Background and
Experimental Procedures. U.S. Environmental Protection Agency,
Athens, GA. EPA-600/7-77-113.
Thomann, R. V. 1982. Development of a Physio-Chemical Model of Toxic
Substances in the Great Lakes. Manhattan College, NY, for U.S.
Environmental Protection Agency, Grosse He, MI, Cooperative Agreement
No. CR105916010, In Press.
Tsivoglou, E. E. and J. R. Wallace. 1972. Characterization of Stream
Reaeration Capacity. U.S. Environmental Protection Agency,
Washington, DC. EPA-R3-72-012.
Wolfe, N. L., R. G. Zepp, G. L. Baughman, and J. A. Gordon. 1975.
Kinetic Investigation of Malathion Degradation in Water. Bull.
Environ. Contain. Toxicol. , 13:707-713.
Wolfe, N. L., R. G. Zepp, G. L. Baughman, and D. M. Cline. 1977.
Kinetics of Chemical Degradation of Malathion in Water. Environ. Sci.
Technol., 11:88-93.
Zepp, R. G., N. L. Wolfe, J. A. Gordon, and G. L. Baughman. 1975.
Dynamics of 2,4-D Esters in Surface Waters: Hydrolysis, Photolysis,
and Vaporization. Environ. Sci. Technol., £:144-150.
Zepp, R. G. and D. M. Cline. 1977. Rates of Direct Photolysis in
Aquatic Environment. Environ. Sci. Technol., 11:359-366.
86
-------
APPENDIX 1:
GLOSSARY
PART 1: GLOSSARY ABBREVIATIONS
L - Local Variable
P - Parameter (i.e. depends on segment)
C - Constant
T - Time function
P* - Parameter from main program
SE
SD
PO
WA
VO
DU
IN
FI
RE
TOXISETL
TOXISEDW
TOXIPHOT
TOXIWASPB
TOXIVOLT
TOXIDUMP
TOXINIT
TOXIFORD
TOXIREOX
PART 2: UNIT ABBREVIATIONS
°C - degrees Celcius
cm - centimeters
ft - feet
g - grams
hr - hours
kcal - kilocalories
kgb - kilograms biomass (dry weight)
kgs - kilograms sediment (dry weight)
L - length
loct - liters of octonol
lw - liters of water in segment
IT - liters total volume of segment
M - mass
M - molarity, moles per liter
m - meters
MCF - million cubic feet total volume
MCFW - million cubic feet of water in segment
mi - miles
mg - milligrams
mgb - milligrams biomass (dry weight)
mgc - milligrams chemical
mgs - milligrams sediment
ml - milliliters
n - moles chemical
sec - seconds
T - time
ug - micrograms
ugc - micrograms chemical
87
-------
PART 3: VARIABLES AND CONSTANTS IN TOXIWASP SUBROUTINES
A(J) p*
Cross sectional area between exchanging sediment and water
compartments, input in subroutine WASP2.
Units: square feet
ACBACG(J) P(WA,IN,FI)
Proportion of bacterial population that actively degrades
toxicant. If the biolysis rate constants are not based on natural
mixed bacterial populations, the total bacterial populations
(BACTOG) given for each compartment can be modified via ACBACG to
give the size of the population that is actively degrading the
toxicant.
Units: dimensionless ratio, nominal range 0.0 to 1.0
ACBACL L(FI)
Active bacterial population for segment being considered. Equal to
BACTOG(J).
Units: Water column compartments: cells per milliliter
Benthic compartments: cells per milliliter pore water
ALPHA (I) L(WA,FI,DU,SE,SD)
The values of ALPHA are distribution coefficients (fraction of
total concentration of toxicant (Y) present as a particular
species/form configuration of the molecule) for each ecosystem
compartment. ALPHA vector represents the partitioning of each
species among three physical forms (dissolved, sediment-sorbed,
biosorbed). ALPHAs are calculated internally from the parameters,
constants, and state variables describing the segment, including
the predicted concentration, the sediment, the biomass, the
octonol-water partition coefficient, etc. In the output these
ALPHAs are designated by DISSF, SEDF and BIOLF.
ALPHA(1)
Fraction of toxicant present as the neutral molecule (SH2)
dissolved in the water phase of the compartment.
ALPHA(2)
Fraction present as neutral molecule sorbed with sediment phase
of compartment.
ALPHA(3)
Fraction present as neutral molecule sorbed with compartment
biomass.
88
-------
ALPH1M L(WA,SD)
Fraction of chemical dissolved in the water phase of the segment
immediately above the bed. Equal to ALPHA(l) for that segment.
Transferred to subroutine TOXISEDW for calculation of bed-water
column mixing of chemical.
ALPH2M L(WA,SD)
Fraction of chemical sorbed onto sediment phase of the segment
immediately above the bed. Equal to ALPHA(2) for that segment.
Transferred to subroutine TOXISEDW for calculation of bed-water
column mixing of chemical.
BACTOG(J) P(WA,IN,FI)
Bacterial population density in each compartment.
Units: Water column compartments: cells per milliliter
Benthic compartments: cells per 100 grams dry weight of
sediment
Internally, BACTOG(J) is multiplied by ACBACG(J) to give active
population density. For the sediment, BACTOG is converted to cells
per milliliter pore water. The conversion used is:
BACTOG(J) = BACTOG(J) * SED/FRW(lE08)
where:
_ cells cells
T00~g 100 * 1000 mgs
FRW = —^ = —— * 1000
SED(mg/ml) = (l/1000)SED(mg/lT)/FRW
BETA L(VO)
Same as ALPHA(l) for water compartments that intersect surface.
BIOFAC L(WA)
Intermediate variable for determining fraction of chemical sorbed
onto biological phase of a segment.
BIOMAS(J) P(WA,IN,FI)
Total actively sorbing biomass in each ecosystem compartment. This
parameter is used in computation of sorption of the toxicant with
plant/animal material in the ecosystem compartments. The parameter
is interpreted differently for the water column versus the benthic
compartments. For a water column compartment, total biomass must
be expressed as milligram (dry weight) per liter of water in the
compartment, and it includes all biomass subject to biosorptive
exchange with that water. In the case of benthic compartments,
89
-------
BIOMAS is the total biomass of the benthic infauna and other
components in grams (dry weight) per square meter of bottom.
Note that in this simplification from EXAMS, movable biomass (e.g.
plankton) is not distinguished from stationary biomass (e.g.,
roots). This is a potential source of error in systems having high
biotic content.
Units: Water column compartments: mg (dry weight) per liter
Benthic compartments: grams (dry weight) per square meter
BIOLKL L(WA,FI,DU)
Total pseudo-first-order degradation rate constant (per hour) for
bacterial biolysis.
Units: per hour
BIOTMG(J) P(WA,FI)
Biotemperature in each ecosystem compartment, i.e. tempeirature to
be used in conjunction with Q-10 expressions for biolysis rate
constants. This parameter is separated from the physical
temperature input data (TCELG) in order that the input data can
reflect Q-10 averaging of an observed temperature time-series.
Units: degrees C.
BMASS(J) L(WA,FI,SE)
Mass of chemical lost from the network during the simulation.
Chemical can be lost by volatilization through a surface water
segment or burial through a bottom bed segment.
Units: kg
BOTLIT L(PO,IN)
Light level at bottom of compartment
Units: dimensionless, fraction of incident light level
BURY L(SE)
Net burial or erosion rate of top benthic segment. Equal to the
time variable depth of sediment settling in from the overlying
water column minus the time-constant depth of sediment eroding.
Both the settling and erosion rates are entered through the
spatially-variable parameter WS (J).
Units: meters per year
BVOLO(J) L(WA,FI,SE)
Volume of the surface bed segment at time 0. The volume of the
surface bed segment is reset to BVOLO during the compaction cycle.
Units: MCF
90
-------
CBB, CBW, CSB, CSW, CTB, CTW, C1B, C1W L(WA,DU)
CBB: Concentration of chemical sorbed onto biological phase of
bed segment JSTR. Values transferred to TOXIDU every 3
hours for saving on the statistical file.
Units: ug/g
CBW: Concentration of chemical sorbed onto biological phase of
water segment JTR. Values transferred to TOXIDU every 3
hours for saving on the statistical file.
Units: ug/g
CSB: Concentration of chemical sorbed onto sediment phase of
bed segment JSTR. Values transferred to TOXIDU every 3
hours for saving on the statistical file.
Units: mgc/kgs
CSW: Concentration of chemical sorbed onto sediment phase of
water segment JTR. Values transferred to TOXIDU every 3
hours for saving on the statistical file.
Units: mgc/kgs
CTB: Total concentration of chemical in bed segment JSTR.
Values transferred to TOXIDU every 3 hours for saving on
the statistical file.
Units:
CTW: Total concentration of chemical in water segment JTR.
Values transferred to TOXIDU every 3 hours for saving on
the statistical file.
Units: ragc/l-p
ClB: Concentration of chemical dissolved in water phase of bed
segment JSTR. Values transferred to TOXIDU every 3 hours
for saving on the statistical file.
Units: mgc/1
w
Ciw: Concentration of chemical dissolved in water phase at
water segment JTR. Values transferred to TOXIDU every 3
hours for saving on the statistical file.
Units: mgc/lw
91
-------
CHEM L(WA,SE,SD)
Chemical concentration in segment, equivalent to C(J,1)
Units: mg/lf
CHEM1, CHEM2, CHEM3, CHEMB, CHEMS, CHEMW L(DU)
CHEM1: Chemical dissolved in water phase of current segment.
Units: mgc/l-p
CHEM2: Chemical sorbed onto sediment phase in current segment.
Units: rage/I^
CHEMS: Chemical sorbed onto biological phase in current segment.
Units: mgc/l^
CHEMB: Chemical sorbed onto biological phase in current segment.
Units: ugc/gbiomass
CHEMS: Chemical sorbed onto sediment phase in current segment.
Units: mgc/kgs
CHEMW: Chemical dissolved in water phase of current segment.
Units: mgc/lw
CHEM1S, CHEM2S, CHEMSS L(SD)
CHEM1S: Chemical concentration dissolved in pore water in bed
segment
Units: mgc/lw
CHEM2S: Chemical concentration sorbed on sediment in bed segment
Units: mgc/li
CHEMSS: Chemical concentration sorbed on sediment in bed segment
Units: mgc/kgs
CHEM1W, CHEM2W, CHEMSW L(SD)
CHEM1W: Chemical concentration dissolved.in water segment above
bed segment
Units:
92
-------
CHEM2W: Chemical concentration sorbed on sediment in water segment
above bed segment
Units: mgc/lj
CHEMSW: Chemical concentration on sorbed sediment in water segment
above bed segment
Units: mgc/kgs
CLOUDG C(WA,FI,PO)
Average cloudiness in tenths of full sky cover.
Units: dimensionless, range of 0.0 to 10.0
CMAX(9) L(WA,IN)
The concentration of a 10" 5 molar solution of chemical, or half
the chemical solubility, whichever is less. CMAX(9) is used to
abort the simulation whenever dissolved chemical concentration
rises above this limit. This prevents violation of the model's
first-order kinetics assumption. This check is activated only when
CMAX(l) is set to 0.
Units: mgc/lw
CMAX(IO) L(IN)
Half the chemical solubility.
Units: mgc/lw
CMPETG(J): P(WA.PO)
Single-valued zenith light extinction coefficients for water
columns, dummy variable for benthic compartments.
Units: per meter
This coefficient is based on external calculations described in
EXAMS. The user is referred to the EXAMS manual (2.3.3.2.2.) and
Fig. 2.18.
COND L(VO)
Inverse of addition of series resistances of gas and liquid
interfaces.
Units: meters per hour
CTRIG C(WA,IN,DU)
Trigger concentration that defines a peak event. When the chemical
concentration in segment JTR rises above CTRIG, a peak event is
flagged. Concentrations for all segments are printed out every 3
hours until the concentration falls below CTRIG and the event
93
-------
ends. This option is designed to catch high transient concentrations
that would be missed by the regular WASP dumps. If less frequent
peak printout is desired, statement 68 in TOXIDUMP can be changed
from PNEXT = 3.0 to, say PNEXT = 8.0. This would cause printouts
every 8 hours during peak events. A trigger concentration of 0
will disable event printouts.
Units: mgc/1^
BELT L(WA)
Intermediate value of simulation time step used for maximizing time
step to minimize numerical dispersion and simulation cost. Ranges
between 0.01 and 0.50 days.
Units: days
DENOM L(WA)
Intermediate variable for calculating fractions of chemical
dissolved, sorbed onto sediment, and sorbed onto biomass.
DEPTH L(WA,SE)
Depth of segment being considered.
Units: feet
DEPTHG(J) P(WA,IN,FI,VO,PO,SE)
Depths of segments.
Units: feet
DEPTHM L(WA,PO,SE)
Depth of segment being considered.
Units: meters
DFACG C(WA,PO)
Distribution function (ratio of optical path length to vertical
depth).
Units: dimensionless, values constrained to range of 1.0 to 2.0
Default value 1.19, explained in TOXIPHOT.
DISP L(SW)
Volumetric dispersion between a bed and an overlying water segment,
or between two vertically adjacent bed segments. DISP is brought
from WASP2 using BR(I), and is corrected internally for its
resulting mixed units.
Units: MCF per second times square cm per square mile
94
-------
DISPV(J) L(WA,SD)
Dispersive exchange volumes between sediment and water
compartments. Water-water exchanges are calculated in WASP2.
Sediment-water exchanges are calculated in TOXISEDW using
characteristic lengths, areas, and dispersion fromWASPZ, along
with porosity and other factors from TOXISEDW.
Units: million cubic feet per day
DSPSED(J) P(WA,SD)
Fraction of sediment that mixes. 1.0 is equivalent to full bed
sediment dispersion, 0.0 is equal to pore water diffusion only.
Units: dimensionless
DTOPT C(WA)
Option to optimize time step. If set to 1, program computes
maximum time step that preserves numerical stability throughout
network. Both flow volumes and dispersive exchange volumes are
kept less than or equal to segment volumes. Time step can vary
between 0.01 days and 0.5 days.
E(J) P*
Sediment-water dispersion coefficent input in WASP2. E is a
composite of direct sorption to the sediment surface, mixing of the
sediments by benthic animals, stirring by demersal fishes, etc.
Units: square centimeters per second
EAHG(I.l) C(FI)
Arrhenius activation energy of specific-acid-catalyzed hydrolysis
of the toxicant.
Units: kcal/gram mole
EBHG(I.l) C(FI)
Arrhenius activation energy of specific-base-catalyzed hydrolysis
of the toxicant.
Units: kcal/gram mole
EHENG C(VO)
Constant used to compute Henry's Law constants for volatilization
as a function of environmental temperatures (TCELG). When EHENG is
non-zero, the Henry's Law constant is computed as follows.
log HENRY = HENRYG-((1000.*EHENG)/(4.58*(TCELG+273.15)))
Units: kcal/gram mole
95
-------
ENHG(I,1) C(FI)
Arrhenius activation energy of neutral hydrolysis of the toxicant.
Units: kcal/gram mole
EOXG(I.l) C(FI)
Arrhenius activation energy of oxidative transformation of the
toxicant.
Units: kcal/gram mole
ESOLG C(VO)
Exponential term for describing solubility of the toxicant as a
function of temperature (see SOLG).
Units: kcal/gram mole
EVPRG C(VO)
Molar heat of vaporization for vapor pressure described as a
function of temperature (See VAPRG).
Units: kcal/gram mole
EXDO L(SD)
Chemical sorption rate from lower water column segment to top bed
segment.
Units: MCF per day times ragc/l-p
EXUP L(SD)
Chemical desorption rate from upper bed segment to lower water
column segment.
Units: MCF per day times mgc/lj
FAC C(IN,SE)
Multiplication factor for sedimentation time step. Set FAC to 1.0
to eliminate numerical dispersion in the bed. For smoother output
with some numerical dispersion, set FAC to 0.1.
FACTOR L(PO)
Latitude correction factor by which the photolysis rate is adjusted
from the rate at the reference latitude.
FRW(J) L(WA,IN,SE,SD)
Porosity, or water fraction in sediment on a volumetric basis.
Calculated from the parameter PCTWA(J).
Units: lw/lx
96
-------
FVOL L(SE)
The fractional volume of the surface bed segment that is to be
compacted into the volume of the second bed segment during the
compaction cycle. The difference between FVOL and the volume of
the second bed segment is SVOL, the pore water squeezed out.
Units: MCF
HENRYG C(FI,VO)
Henry's Law constant of the toxicant. If parameter EHENG is
non-zero, HENRYG is used as the pre-exponential factor in computing
the Henry's Law constant as a function of environmental temperature
(TCELG).
Units: Atmosphere-cubic meters per mole
HENRYL L(VO)
Local value of HENRYG. If HENRYG is zero, model will calculate
value.
Units: Atmosphere-cubic meters per mole
HPLUS L(FI)
Intermediate calculation in obtaining temporally averaged
concentration of hydronium ions.
HPLUS = 10**(-PHG)
HYDRKL L(FI,DU)
Total pseudo-first-order rate constant (per hr) for hydrolytic
transformations of the toxicant in each compatment.
Units: per hour
HYDRX L(FI)
Intermediate calculation in obtaining temporally averaged
concentration of hydroxide ions.
HYDRX = 10**(-POHG)
II, 12 L(SD)
Local variables representing segments involved in exchange.
(See "IR(I), JR(I)" card group B)
ICHK L(WA)
Integer flag denoting day on which last set of daily flows and
nonpoint source loads were read from auxiliary file. Saved in
COMMON as TIMCHK.
ICOUNT L(DU)
Number of entries into auxiliary file for statistical analysis.
Integer value of TCOUNT.
97
-------
ICOL L(IN)
Subscript indicating column in table of daily flows from auxiliary
file (1-10).
IL(J) P*
Characteristic length for dispersive exchange, input in WASP2.
Units: feet
INDEXS L(WA,FI)
INDEXS is an internally calculated flag designating a benthic
compartment. INDEXS is 1 for benthic compartments (TYPEE = 3 or
4), 0 for water compartments (TYPEE = 1 or 2).
1NDEXW L(WA,FI)
INDEXW is an internally calculated flag designating a water column
compartment. It takes a value of 1 for water compartments (TYPEE =
1 or 2) and a value of 0 for benthic compartments (TYPEE = 3 or 4).
IOPT L(DU)
Flag designating type of output desired in TOXIDU.
ITCHK L(WA,IN)
Integer flag denoting whether a full day has passed since the last
set of flows and nonpoint source loads were read from auxiliary
file.
ITIMEC L(WA)
Integer flag denoting the current simulation day. The integer
value of TIME.
ITYPE L(IN,SE)
Flag designating segment type. The integer value of TSfPEE(J).
IZERO L(PO)
Light intensity at top of lower level water compartment.
0.0 to 1.0.
Units: dimensionless
JROW L(IN)
Subscript indicating row in table of daily flows from auxiliary
file.
JSTR L(WA,IN,DU)
Top bed segment below water segment FTR. Concentrations in JSTR
are printed out during peak "events," and are saved every 3 hours
on an auxiliary statistical file.
98
-------
JTR L(WA,DU;
Segment used to trigger event printouts when total chemical
concentration exceeds CTRIG. Saved in COMMON as XJTR.
Concentrations in JTR are printed out during peak events, and are
saved every 3 hours on an auxiliary statistical file.
JL(J) P*
Characteristic length for dispersive exchange input in WASP2.
Units: feet
KAHG(I.l). C(FI)
Second-order rate constants for specific-acid-catalyzed hydrolysis
of toxicant. If the corresponding entry in the Arrhenius
activation energy matrix (EAHG) for this reaction is zero, the
value entered in KAHG is taken as the second-order rate constant.
If the corresponding entry in the activation energy matrix (EAHG)
is non-zero, the value entered in matrix KAHG is interpreted as the
base-10 logarithm of the frequency factor in an Arrhenius function
for the reaction, and local values (KAHL) of the second-order rate
constant are computed as a function of temperature (TCELG) in each
system compartment.
Units: per mole[H+]per hour
KAHL L(FI)
Local value of KAHG(I,1) corrected for temperature.
Units: per mole [H+J per hour
KB L(WA,FI,IN)
Effective biomass partition coefficient. Calculated internally by
multiplying OCB and KOW.
Units: lw/kgb
KBACSG(I,1) C(FI)
Second-order rate constants for benthic sediment bacterial biolysis
of the organic toxicant. If the corresponding entry in the Q—10
matrix (QTBASG) for this process is zero, the number entered in
matrix KBACSG is taken as the second-order rate constant. If the
corresponding entry in the Q-10 matrix is non-zero, the value
entered in matrix KBACSG is interpreted as the numerical value of
the second-order rate constant at 20 degrees C. , and local values
(KBACSL) of the rate constant are computed as a function of
temperature (BIOTMG) in each ecosystem compartment.
Units: ml/cell-hour
99
-------
KBACSL L(WA,FI)
Local value of KBACSG(I,1), corrected for temperature.
Units: ml/cell-hour
KBACWG(I,1) C(FI)
Second-order rate constants for water column bacterial biolysis of
the organic toxicant. If the corresponding entry in the Q-10
matrix (QTBAWG) for this process is zero, the number entered in
matrix KBACWG is taken as the second-order rate constant. If the
corresponding entry in the Q-10 matrix is non-zero, the value
entered in matrix KBACWG is interpreted as the numerical value
KBACWL of the rate constant are computed as a function of
temperature (BIOTMG) in each ecosystem compartment.
Units: ml/cell-hour
KBACWL L(WA,FI)
Local value of KBACWG(J), corrected for temperature.
Units: mg/cell-hour
KBHG(I.l) C(FI)
Second-order rate constants for specific-base-catalyzed hyrolysis
of toxicant. If the corresponding entry in the Arrhenius
activation energy matrix (EBHG) for this reaction is zero, the
value entered in KBHG is taken as the second-order rate constant.
If the corresponding entry in the activation energy matrix (EBHG)
is non-zero, the value entered in matrix KBHG is interpreted as the
base-10 logarithm of the frequency factor in an Arrhenius function
for the reaction, and local values (KBHL) of the second-order rate
constant are computed as a function of temperature (TCELG) in each
system compartment.
Units: per mole[OH-] per hour
KBHL L(WA,FI)
Local value of KBHG(I,1), corrected for temperature.
Units: per mole [OH-] per hour
KDPG C(FI,PO)
A near-surface photolytic rate constant for the toxicant. The
value of KDPG represents the outcome of an experiment conducted in
natural sunlight. The constant is a temporally averaged (e.g. over
whole days, seasons, etc.) first-order photolytic transformation
rate constant pertaining to cloudless conditions at some reference
latitude RFLATG.
Units: per hour
100
-------
KDPL L(FI,PO)
Locally adjusted value of KDPG returned from TOXIPHOT.
KNHG(I.l) C(FI)
Rate constants for neutral hydrolysis of organic toxicant. If the
corresponding entry in the Arrhenius activation energy matrix
(ENHG) for this reaction is zero, the value entered in KNHG is
taken as the rate constant. If the corresponding entry in the
activation energy matrix (ENHG) is non-zero, the value entered in
matrix KNHG is interpreted as the base-10 logarithm of the
frequency factor in an Arrhenius function for the reaction, and
local values (KNHL) of the rate constant are computed as a function
of temperature (TCELG) in each system compartment.
Units: per hour
KNHL L(FI)
Local value of KNHG(I,1), corrected for temperature.
Units: per hour
K02G(J) L(WA,IN,VO)
Reaeration parameter at 20 degrees C in each ecosystem
compartment. Calculated from segment depths and velocities, and
time-varying wind.
Units: centimeter per hour
K02L L(WA,VO)
Reaeration parameters in each compartment (m/hr) after temperature
adjustment and units conversion.
KOC C(WA,IN)
Organic carbon partition coefficient. Value of KOC read in or if
equal to zero, calculated from KOW. Multiplication of KOC by the
fractional organic carbon content (OCS) of each system sediment
yields the partition coefficient for sorption of unionized (SH2)
compound to the sediment.
Units: lw/kg organic carbon
KOW C(WA,IN)
Octanol water partition coefficient Value of KOW read in, or if
equal to zero, calculated from KOC.
Units: (mgc/loct)-(mgc/lw), or lw/loct
101
-------
KOXG(I,1) C(FI)
Second-order rate constants for oxidative transformation of
toxicant. If the corresponding entry in the Arrhenius activation
energy matrix (EOXG) for this reaction is zero, the value entered in
KOXG is taken as the second-order rate constant. If the
corresponding entry in the activation energy matrix (EOXG) is
non-zero, the value entered in matrix KOXG is interpreted as the
base-10 logarithm of the frequency factor in an Arrhenius: function
for the reaction, and local values (KOXL) of the second order
constants are computed as function of temperature (TCELG) in each
system compartment.
Units: liter per mole environmental oxidant per hour
KOXL L(FI)
Local value of KOXG(I,1), corrected for temperature.
Units: liter per mole environmental oxidant per hour
KP(J) L(WA,IN,SE)
Effective sediment partition coefficient. Internally calculated as
the product of KOC and OCS(J), and saved as Parameter 12.
Units: lw/kgs
KT L(IN)
Daily counter for reading nonpoint source loads from auxiliary tape
and printing table (1-1000).
Units: days
KVOG C(VO)
Measured experimental value for (volatilization) liquid-phase
transport resistance, expressed as a ratio to the reaeration rate.
Units: dimensionless ratio
K20 L(FO,RE)
The computed reaeration rate at 20?oOC. Used to calculate the
chemical volatilization rate.
Units: per day
LATG C(PO)
Geographic latitude of ecosystem.
Units: degrees and tenths (e.g. 37.2)
LIGHTL L(PO)
The average light intensity in the current compartment, as a
fraction of the near-surface light intensity (taken as 1.0 or 100%).
102
-------
LIGHTN T(WA,PO)
Normalized light time function, transferred through COMMON as
constant 78 to TOXIPHOT. There it adjusts the average photolysis
rate for seasonal light variability.
LOPT C(WA,IN)
Option to read in flows and/or nonpoint source loads on a daily
basis:
0 = skip daily read option
1 = read sequential tape containing daily flows and/or loads
MOQ L(WA,IN)
Number of flow pairs read from auxilliary file daily. MOQ is read
from auxilliary file by TOXINIT and passed to TOXIWASPB through
COMMON as the floating point variable MOQS.
MOQS L(WA,IN)
Number of flow pairs read from auxilliary file daily. MOQS is set
equal to MOQ in TOXINIT and saved in COMMON as constant 83.
MQOPT L(IN)
Flow option read from auxilliary file. Not used.
MTYPE L(IN)
Flag designating type of segment immediately above current segment.
Integer value of TYPEE(J-l).
MWTG C(VO)
The molecular weight of the toxicant.
Units: grams per mole
NCOL L(IN)
Variable indicating the number of column entries in the last row of
daily flows from auxilliary file.
NOWKS L(WA,IN)
Number of nonpoint source loads read from auxilliary file daily.
NOWKS is read from auxilliary file by TOXINIT and passed to
TOXIWASPB through COMMON as constant 84.
NPSWK(I,J) L(WA,IN)
Nonpoint source load J for constituent I (chemical or sediment),
read from auxilliary file daily. The segment into which load J
discharges is defined in Card Group F: Forcing Functions.
Units: million pounds per day
NWKS L(WA,IN)
Number of nonpoint source loads read from auxilliary file daily.
NWKS is the integer value of NOWKS.
103
-------
OCB C(WA,IN)
Organic carbon content of the compartment biomass as a fraction of
dry weight. Coupled to KOW to generate biomass partition
coefficient.
Units: dimensionless
OCS(J) P(WA,IN)
Organic carbon content of sediments as fraction of dry weight.
Parameter is coupled to KOC to generate the sediment partition
coefficient as a function of a property of the sediment.
Units: dimensionless
OXIDKL L(FI,DU)
Pseudo-first-order rate constants for oxidative transformation of
toxicant.
Units: per hour
OXRADG(J) P(FI)
Molar concentration of environmental oxidants (e.g. peroxy radicals)
in each ecosystem compartment.
Units: moles per liter
PCTWA(J) P(WA,IN)
Percent water in bottom sediments of benthic compartments. PCTWA(J)
should be expressed as the conventional soil-science variable
(fresh/dry weight); all values must be greater than or equal to 1.
Units: dimensionless
PERC L(SE)
Pore water percolation, calculated from the spatially-variable
parameter WS(J) for type 4 (subsurface benthic) segments.
Units: MCF per day
PERCMS L(SE)
Mass transport rate of dissolved chemical in pore water.
Units: MCF per day times mgc/lw
PH L(WA,FI)
Hydrogen ion activity of segment being considered. Local value of
PHG(J).
Units: pH units
104
-------
PHG(J) P(WA)
Hydrogen ion activity. The negative value of the power to which 10
is raised in order to obtain the temporally averaged concentration
of hydronium ions [H30+] in gram-molecules per liter.
Units: pH units
PHOTKL L(FI,DU)
Pseudo-first-order rate constant for photolytic transformation of
the toxicant.
Units: per hour
PHN T(WA)
Normalized time function for PH. Adjusts PHG(J) for seasonal
variability.
POHN T(WA)
Normalized time function for POH. Adjusts POH'XJ) for seasonal
variability.
PNEXT L(WA,DU)
Triggers event printout every 3 hours during peak event period.
Local value of POHG(J).
POH L(WA,FI)
Hydroxide ion activity of segment being considered. Local value of
POHG(J).
Units: POH units
POHG(J) P(WA)
Hydroxide ion activity. The negative value of the power to which 10
is raised in order to obtain the temporally averaged concentration
of hydroxide [OH-] ions in gram-molecules per liter.
Units: pOH units
POROS L(SD)
Porosity of top bed segment, used in computing dispersive exchange
between bed and water column. For computing dispersive exchange
between two bed segments, POROS is the mean porosity.
Units: IW/IT
PTIME L(DU)
Rounded-off value of TIME for output.
Units: days
QT(I) L(IN)
Variable for temporary storage of daily flows from auxilliary file.
105
-------
QTBASG(I,1) C(FI)
Q-10 values for bacterial transformation (c.f. KBACSG) of organic
toxicant in benthic sediments. The Q-10 is the increase in the
second-order rate constant resulting from a 10 degree C temperature
increase.
Units: dimensionless
QTBAWG(I.l) C(FI)
Q-10 values for bacterial transformation (c.f. KBACWG) of chemical
in the water column of the system. Q-10 is the increase in the
second-order rate constant resulting from a 10 degree C temperature
increase.
Units: dimensionless
QUANTG(l.l) C(FI)
Reaction quantum yield in photolytic transformation of chemical.
The quantum yield is the fraction of total quanta absorbed by the
toxicant resulting in transformations. Separate values are provided
for each molecular configuration of the toxicant in order to make
assumptions concerning their relative reactivities readily available
to the user.
Units: dimensionless
RATEK L(VO,PO)
Internally calculated rate returned from TOXIVOLT or TOXIPHOT to
TOXIFORD for use as VOLKL or KDPL.
RESGAS L(VO)
Gas film resistance to volatilization.
Units: hours per meter
RESLIQ L(VO)
Liquid film resistance to volatilization.
Units: hours per meter
RFLATG C(WA)
Reference latitude for corresponding direct photolysis rate constant
(c.f. KDPG).
Units: degrees and decimal fraction (e.g. 40.72)
RVOL(J) L(IN,SE)
Reference volume at which the sediment compaction cycle is initiated
for surface bed segment "J". RVOL(J) is computed in TOXINTT for
each surface bed segment and is stored in location RVOL(J-l)
[equivalenced to VVOL(J-l)]. RVOL(J) is the surface bed volume that
will hold the mass of sediment contained in the upper two bed
106
-------
segments at time 0. RVOL(J) will be greater than the combined
volume of the upper two bed segments at time 0 if the density of the
second bed segment is greater than the top bed segment.
Units: MCF
SCALQ L(IN)
Scale factor to convert daily flows from auxiliary file to cubic
feet per second.
SED L(WA,DU,SE,SD)
Sediment concentration of segment being considered.
Units: mgs/l-p
SEDCOL L(SD)
Sediment concentration after conversion to internal units (kg/liter
of water) in sediment compartment.
Units: kg/lw
SEDFAC L(WA)
Intermediate variable for determining fraction of chemical sorbed
onto sediment phase of a segment.
SEDFL L(SD)
Rate of sediment mixture to the surface of the bed, allowing
sorption or desorption with the overlying water. Related to the
dispersive exchange of water by the spatially-varying parameter
DSPSED(J).
Units: MCF (water) per day times mgs/lw
SETINS, SETOUS L(SE)
Sediment settling into or out of compartment. These variables are
transferred to TOXIWASPD to hold the values until the next call to
TOXISETL.
Units: rags/lf-day
SETINC, SETOUC L(SE)
Chemical settling into or out of a compartment. Chemical can only
settle when adsorbed to sediment.These variables are transferred to
TOXIWASPD to hold the values until the next call to TOXISETL.
Units: mgc/lf-day
SEDW L(SD)
Sediment concentration in water compartment above bed.
Units: kg/lT
107
-------
SMASS L(SE)
Mass of dissolved chemical in the pore water squeezed from the
surface bed segment to the water column during the compaction
cycle.
Units: MCFW times mgc/lw
SOLG C(VO)
Aqueous solubility of toxicant chemical species. If the
corresponding value in ESOLG (c.f.) is zero, SOLG is interpreted as
an aqueous solubility in mg/liter. If ESOLG is non-zero, SOLG is
used in an equation describing the molar solubility of the toxicant
species as a function of environmental temperature (TCELG), i.e.
SOLL(mg/l) = 1000.*MWTG*10.**J(SOLG-(1000.*ESOLG/(2.303*R*(TCELG+
273.15)))) Solubilities are used (inter alia) to limit the
permissible external loadings of the toxicant on the system to
values that generate final residual concentrations .LE. 50% of
aqueous solubility (or 1.E-5M). This constraint is imposed in order
to help ensure that the assumption of linear sorption isotherms is
not seriously violated.
Units: mg/lj
SOLL L(VO)
Temperature-corrected aqueous solubility of chemical
SVOL L(SE)
Volume of pore water squeezed from the surface bed segment to the
water column during the compaction cycle. SVOL is the difference in
volume between the surface bed segment before compaction and the top
two bed segments after compaction.
Units: MCF
TCELG,TKEL L(Wa,FI,VO)
Product of segment temperature TEMPM(J) and temperature function
TEMPN. Immediately converted to Kelvin temperature (TKEL) from
Celsius.
TCOUNT L(IN,DU)
Number of entries into auxilliary file for statistical analysis.
Stored in COMMON as constant 73.
TDINT C(WA.DU)
Time interval between recalculation of decay rates.
Units: days
TEMPM(J) P(WA)
Average temperature for compartment.
Units: degrees C
108
-------
TEMPN T(WA)
Value of normalized time function describing temperature. Multiples
TEMPM to give TCELG,
Unitless
TIMCHK L(WA,IN)
Day on which last set of daily flows and nonpoint source loads were
read from auxiliary file. Saved in COMMON as constant 82.
Units: day
TMARK L(WA,IN)
Day on which rate constants will be recalculated. Incremented
throughout simulation by TDINT.
Units: day
TMASS L(SE)
Mass of chemical transferred to the next lower bed segment during
the compaction cycle. TMASS is first set to the mass of chemical in
the top bed layer that is compacted into the new second bed layer
(FVOL*CHEM). Once the new concentration of the second bed layer is
determined, TMASS is set to the mass of chemical in the second bed
layer that is buried into the new third bed layer. This is
continued to the bottom bed segment, where TMASS is added to
BMASS(J), the amount of chemical mass lost from the network by
burial below segment J.
Units: MCF times mgc/lT
TMPM L(SE)
Temporarily held value of TMASS, the amount of chemical transferred
to the next lower bed segment during the compaction cycle.
Units: MCF times mgc/lx
TOTKG(J) P(WA,IN)
Total first order decay rates calculated externally. If equal to
zero, the program evaluates all separate processes and calculates a
combined total first order decay rate, TOTKL.
Units: per day
TOTKL L(WA,FI,DU)
Value of total first-order decay rate for compartment, derived from
TOTKG(J).
Units: per day
TQ L(WA)
Temporary value for flow or exchange between two segments, used for maxi-
mizing time step to minimize numerical dispersion and simulation cost.
109
-------
Units: million cubic feet per day
TRT L(WA)
Test value of time step that would eliminate numerical dispersion in
a segment.
Units: days
TYPEE(J) P(WA,IN,FI,PO,SE,SC)
Numerical code designating segment types used to define ecosystem.
Available types: 1 = Epilimnion, 2 = Hypolimnion, 3 = Benthic
(active), and 4 = Benthic (buried).
Units: dimensionless
TYP1,TYP2 L(SD)
Internal variables describing type (e.g. TYP1 = TYPEE(Il)) of
interchanging compartments.
V1,V2,...,V8 L(WA,DU)
Dummy variables used in argument of calls to TOXIDU.
VAPRG C(VO)
Vapor pressure of compound. Used to compute Henry's Law constant if
the latter input datum (HENRYG) is zero, but VAPRG is non-zero.
Units: torr
VAPRL is the converted value to atmospheres. L(VO)
VELOC(J) P(WA,IN)
Average velocity of water in the compartment
Units: feet per second
VOL L(WA)
Equal to BVOL(J), the volume of the compartment.
Units: million cubic feet
VOLKL L(DU,FI)
Pseudo-first-order rate constants for volatilization losses from
surface water compartments.
Units: per hour
VVOL(J) L(WA,FI,SE)
A dummy array to hold parameters describing the varying surface bed
segment volume.
WAT L(VO)
Piston velocity term for water vapor.
110
-------
Units: meters per second
WATFL L(SD)
Dispersive exchange of water between top bed segment and water
column, or between two vertically adjacent bed segments.
Units: million cubic feet per day
WIND L(WA,VO)
Local, time-varying wind speed. Equal to the average wind velocity
over a segment WINDG(J) times the time function WINDN.
Units: meters per second
WINDG(J) P(WA)
Average wind velocity at a reference height of 10 cm above the water
surface. Parameter is used to compute a piston velocity for water
vapor (Liss 1973, Deep-Sea Research 20:221) in subroutine VOLAT.
Units: meters per second
WINDN T(WA)
Normalized wind speed as function of time.
Units: unitless
WS(J) P(SE)
Spatially variable parameter denoting settling rate of suspended
sediment in water column segments (types 1 and 2), erosion rate of
surface bed sediment in surface bed segment (type 3), or percolation
of pore water through subsurface bed segments (type 4).
Units: Segment types 1,2: meters per day
Segment type 3: cm per year
Segment type A: cubic feet per second
XJSTR L(WA,IN,DU)
Top bed segment below water segment JTR. Equal to JSTR, this
variable is saved in COMMON as constant 75.
XJTR C(WA,IN,DU)
Reference segment for triggering event and frequency output. A
value of zero will disable it. A non-zero value will cause printout
of all segment concentrations when the concentration at XJTR is
greater than CTRIG. Also, a time history of concentrations (every 3
hours) at XJTR and its associated benthic segment will be stored for
later frequency analysis.
XMASS L(DU)
Mass of chemical in segment.
Units: kg
111
-------
XXX L(PO)
Exponential light disappearance term.
XTES L(PO)
Maximum light disappearance term. After this value (87.4) no light
is present.
XTST L(IN,DU)
Test variable to determine whether to set up event-triggered
printout. If either CTR1G or JTR is zero, then event-triggered
printout is disabled.
PART 4: VARIABLES IN WASP COMMON
The following list defines the variables contained in blank COMMON.
Blank COMMON is used by WASP as the vehicle to pass information from
subroutine to subroutine within the program. The R, I, and * con-
tained within parentheses after the variable name indicate, respec-
tively, whether the variable is a REAL (floating point), or INTEGER
(fixed point), and whether the variable is read as input data.
Variable Name Definition
IN(I) Device number for reading input data.
OUT(I) Device number for printer output.
NOSYS(I*) Number of systems or water quality constituents in
the user's model.
ISYS(I) System currently having its derivatives evaluated.
ISEG(I) Segment currently having its derivatives evaluated.
ISIM(I*) Simulation type - currently only time variable is
permitted.
LISTG(I*) User selected option to print exchange coefficient,
segment volume, advective flow and boundary condi-
tion input data.
LISTC(I*) User selected option to print forcing function
(waste load), kinetic constants, segment parameters,
and miscellaneous kinetic time functions, and ini-
tial condition input data.
INITB(I) Internal program indicator which permits the user to
perform initialization or to execute special code
upon initial entry to the WASPB kinetic subroutine.
Initially equal to zero, INITB must be reset by the
user in WASPB.
112
-------
IPRNT(I)
IDUMP(I*)
IDISK(I)
IREC(I)
MXDMP(I)
IDFRC(I)
NBCPSY(I)
NWKPSY(I)
SYSBY(I*)
RBY(I*)
QBY(I*)
Not currently used.
System - segment combinations to be printed out
during the integration procedure.
When checked by the user in the kinetic subroutine,
WASPB, IDISK acts as internal program indicator which
informs the user when a print interval has been reached,
permitting the user to write the current state variables
or segment concentrations to auxiliary storage (disk).
Normally IDISK equals zero, but at a print interval it
is externally set to one; must be reset by the user
before exiting from WASPB.
Internal counter used to keep track of the number of
print intervals generated during the course of the
simulation.
Blocking factor or the maximum number of variables
saved per segment at each print interval.
Used only in the DEC-POP version as the record address
pointers for the direct access dump files. Not needed
for the IBM 370 version because sequential files are
used.
Maximum number of boundary conditions permitted per
system; set for a particular WASP configuration in
subroutine WASP1.
Maximum number of forcing functions (waste loads)
permitted per system; set for a particular WASP con-
figuration in subroutine WASP1.
User selected system by-pass indicators. If a user
wishes he may choose to by-pass computations for a
particular system (or systems), ISYS, for a simula-
tion run by setting SYSBY (ISYS) appropriately.
User selected exchange transport by-pass indicators.
If a user wishes he may by-pass exchange transport for
a particular system, ISYS, by setting RBY(ISYS)
appropriately. (Example: if a user had incorporated
rooted aquatic plants in his model he would not wish to
have them "disperse").
User selected advective transport indicators. If a
user wishes he may by-pass advective transport for a
particular system, ISYS, by setting QBY(ISYS) appro-
priately. (Example: if a user had incorporated
rooted aquatic plants in his model he would not wish
to have them transported via flow).
113
-------
NEGSLN(I*)
TIME(R)
DT(R)
TZERO(R*)
SCALT(R*)
TEND(R)
PRNT(R*)
OMEGA(R)
ITCHCK(R)
MXITER(R)
C(R)
CD(R)
CMAX(R*)
Indicates whether the user has chosen to permit WASP
to compute negative water quality concentrations
(Example: permit negative D.O. deficit, i.e.,
supersaturation). NEGSLN normally equals zero, but
will equal one, if user choses to permit negative
solutions.
Current simulation time.
Current integration time step.
User selectable time for start of simulation. If
for example a user's input data for a model was set
up such that time zero was January 1, a user may
skip computations for January and February and start
March 1 by setting TZERO (on input) to 59.
Time scale factor. The nominal unit for time in
WASP in days. A user may choose to run his simula-
tion in hours by inputting SCALT to be 0.041667 (or
1/24). WASP will then interpret any time specifica-
tions (such as the print interval, integration step
sizes and total simulation time and the time breaks
for any piecewise linear functions) to be hours
rather than days.
Ending time for use of the current integration step
size. For single integration stepsize input this
will be the total simulation time. For multiple
integration stepsize histories, when TIME equals
TEND, a new integration step size will be chosen and
TEND reset.
Print interval.
Not used in current version of WASP.
Not used in current version of WASP.
Not used in current version of WASP.
State variable or water quality concentration array.
Derivative array.
Stability criteria vector. The vector contains the
maximum allowable segment concentration for each
system. If any segment exceeds the stability
criteria for any system (usually because the inte-
gration stepsize is too large) the simulation is
terminated.
114
-------
CMIN(R*) Not used in current version of WASP (although user
must include in his input data check).
PARAM(R*) Segment parameters for use in the WASPB kinetic sub-
routine.
CONST(R*) Constants for use in the WASPB kinetic subroutine.
EVOL(R*) Segment volumes (nominal input units MCF, nominal
internal units MFC).
MVOL(R) Not used in current version of WASP.
The value for any time variable exchange coefficient, advective flow,
boundary condition, forcing function or time-variable function utilized
by the WASPB kinetic subroutine is computed using a form of the
following equation
VAL = M*TIME + B
where
VAL is the desired value at time = TIME
M is the slope of the piecewise linear function used to approximate
the exchange coefficient, flow, etc.
B is the intercept of the piecewise linear function.
BR(R) Exchange coefficient intercepts.
BQ(R) Advective flow intercepts.
BBC(R) Boundary condition intercepts.
BWK(R) Forcing function intercepts.
BFUNC(R) Intercepts for the time variable functions required
for the WASPS kinetic subroutine.
MR(R) Exchange coefficient slopes.
MQ(R) Advective flow slopes.
MBC(R) Boundary condition slopes.
MWK(R) Forcing function slopes.
MFUNC(R) Slopes for the time-variable functions required for
the WASPB kinetic subroutine.
Note; If any of the above are time invariant (i.e., constant in time), then
the "B" vector (array) will contain the time invariant value of exchange,
flow, etc., and the "M" vector (array) will contain zero slope.
115
-------
IRCI*), JR(I*)
, JQCI*)
IBC(I*)
IWK(I*)
IVOPT(I*j
NOV(I*)
IROPT(I*)
NOR(I*)
IQOPT(I*)
NOQ(I*)
IBCOP(I*)
NOBC(I*)
IWKOP(I*)
NOWK(I*)
Contain the segment numbers between which exchange
is to take place.
Contain the segment numbers between which advective
flow is to take place. If the advective flow is
positive then JQ will contain the upstream segment
number (from which flow is leaving) and IQ will con-
tain the downstream segment number (to which flow
will go). If, however, the advective flow is nega-
tive then JQ will be considered the downstream seg-
ment (flow to) and IQ will be considered the up-
stream segment (flow from).
Contains the segment numbers for which boundary con-
ditions have been specified.
Contains the segment numbers for which forcing func-
tions have been specified (i.e., receiving water
segments for waste loads).
User selected volume input option. Currently WASP
only permits time-invariant or constant volumes
(IVOPT=1).
Number of volumes (normally NOV equals NOSEG).
User selected exchange coefficient input option.
IROPT flags the exchange coefficients as constant in
time (IROPT=1) or time-variable (IROPT=2,3).
Number of exchange coefficients read.
User selected advective flow input option. IQOPT
flags the flows as constant in time (IQOPT=1) or
time-variable (IQOPT=2,3).
Number of advective flows read.
User selected boundary condition input options for
each system. IBCOP(ISYS) flags the boundary condi-
tions for system ISYS as being constant in time
(IBCOP(ISYS)=1) or time-variable (IBCOP(ISYS)=2,3).
Number of boundary conditions read for each system.
User selected forcing functions input option for
each system. IWKOP(ISYS) flags the forcing func-
tions for system ISYS as being constant in time
(IWKOP(ISYS)=1) or time-variable (IWKOP(ISYS)=2,3).
Number of forcing functions read for each system.
116
-------
NOPAM(I*)
NCONS(I*)
NFUNC(I*)
NVOLT(R)
NRT(R)
NQT(R)
NBCT(R)
NWKT(R)
NFUNT(R)
ITIMV(I)
Number of segment parameters (for use in the WASPB
kinetic subroutine) read.
Number of constants (for use in the WASPB kinetic
subroutine) read.
Number of time variable functions (for use in the
WASPB kinetic subroutine) read.
Not used in the current version of WASP.
Used if the exchange coefficients are time-variable
(approximated by a piecewise linear functions of
time). NRT will contain the time at which the next
break in the piecewise linear functions will occur,
at which point it will be necessary to obtain new
slopes (MR) and intercepts (BR).
Used if the advective flows are time-variable (ap-
proximated by a piecewise linear functions of time).
NQT will contain the time at which the next break in
the piecewise linear functions will occur, at which
point it will be necessary to obtain new slopes (MQ),
and intercepts (BQ).
Used if the boundary conditions for a system are
time-variable (approximated by a piecewise linear
functions of time). NBCT(ISYS) will contain the
time at which the next break in the piecewise linear
functions for system 1SYS, will occur, at which
point it will be necessary to obtain new slopes
(MBC) and intercepts (BBC) for system ISYS.
Used if the forcing functions for a system are time-
variable (approximated by a piecewise linear func-
tion of time). NWKT(ISYS) will contain the time at
which the next break in the piecewise linear func-
tions, for system ISYS, will occur, at which point
it will be necessary to obtain new slopes (MWK) and
intercepts (BWK) for system ISYS.
Used if time variable functions (approximated as
piecewise linear functions of time) have been read
for use in the WASPB kinetic subroutine. NFUNT will
contain the time at which the next break in the
piecewise linear functions will occur, at which
point it will be necessary to obtain new slopes
(MFUNC) and intercepts (BFUNC).
Not used in current version of WASP.
117
-------
ITIMR(I)
ITIMQ(I)
ITIMP(I)
ITIMB(I)
ITIMR(I)
ITIMQ(I)
ITIMP(I)
ITIMB(I)
ITIMW(I)
Used as a breakpoint counter for obtaining correct
slope and intercept values for the time-variable ex-
change coefficients.
Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable advec-
tive flows.
Used as a breakpoint counter for obtaining correct
slope and intercept values for the time-variable
functions required by the WASPB kinetic subroutine.
Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable bound-
ary conditions.
Used as a breakpoint counter for obtaining correct
slope and intercept values for the time-variable ex-
change coefficients.
Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable advec-
tive flows.
Used as a breakpoint counter for obtaining correct
slope and intercept values for the time-variable
functions required by the WASPB kinetic subroutine.
Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable bound-
ary conditions.
Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable forcing
functions.
118
-------
APPENDIX 2:
TOXIWASP Subroutine Listings
£*AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA*AAAAAAAAAAAAAAAAAA
c
c
C TOXIWASPD
C
C LAST REVISED 10/13/82
C*AAAA**AAAAAAAAAAA*AA*AAAAAAAAA*AAAAAAAAAAAAAAAAAAAAAAA*AAAAAAAAAAAAAAAA*A
C
SUBROUTINE WASPB
C
C THIS SUBROUTINE COMBINES CHEMICAL AND ENVIRONMENTAL PARAMETERS
C TO FIRST ORDER RATES, AND CALCULATES KINETIC DERIVATIVES.
C IT IS BASED ON THE SUBROUTINE FIRORD IN THE EXAMS MODEL BY DR. L.
C A.BURNS. INCLUDED ARE HYDROLYSIS, BIOLYSIS, OXIDATION,PHOTOLYSIS 1,
C AND VOLATILIZATION.
p__ —___________________ _ ________________________________________________
C TOXIWASPD REDUCES THE KINETICS OF THE TOXICANT TO PSEUDO-FIRST-
C ORDER FORM BY COUPLING CHARACTERISTICS OF THE TOXICANT AND THE
C ENVIRONMENT. IF THE INPUT DATA PERMITS, THE PROPERTIES OF THE
C TOXICANT ARE ADJUSTED TO REFLECT THE EFFECT OF TEMPERATURE IN EACH
C PHYSICAL ELEMENT OF THE SYSTEM. THE PSEUDO-FIRST-ORDER REACTION RATE
C CONSTANTS ARE COMPUTED FROM THE INTERACTIVE EFFECTS OF PARTITIONING,
C TEMPERATURE, AND OTHER CHARACTERISTICS OF THE ENVIRONMENT FOR EACH
C SPECIES. FOR SIMPLICITY, ONLY THE NEUTRAL ION FORM IS CONSIDERED HERE.
C TOXIWASPD CALLS THE FOLLOWING SUBROUTINES:
C INIT (TOXINIT)
C TOXIDU (TOXIDUMP)
C FIRORD (TOXIFORD)
C VOLAT (TOXIVOLT)
C PHOT01 (TOXIPHOT)
C REOX (TOXIREOX)
C SETTLE (TOXISETL)
C SEDWAT (TOXISEDW)
C
C
C SYSTEMS
C
C 1 CHEMICAL, MG/L
C 2 SEDIMENT, MG/L
C
c
119
-------
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
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
CONSTANTS
HYDROLYSIS
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
EBHG(l.l)-
EBHG(2,1)
EBHG(3,1)
ENHG(1,1)-
ENHG(2,1)
ENHG(3,D
EAHG(1,1)-
EAHG(2,1)
EAHG(3,1)
KAHG(l.l)-
KAHG(2,1)
KAHG(3,1)
KBHG(1,D-
KBHG(2,1)
KBHG(3,1)
KNHG(l.l)-
KNHG(2,1)
KNHG(3,1)
ARRHENIUS ACTIVATION ENERGY FOR BASE HYDRO.
N DENOTES NEUTRAL HYDROLYSIS
A DENOTES ACID HYDROLYSIS
2ND ORD. RATE CONST. FOR ACID HYD.
B DENOTES BASE HYDROLYSIS
N DENOTES NEUTRAL HYDROLYSIS
OXIDATION
19 EOXG(l.l)- ARHENIUS ACT. ENERGY OF OXI. TRANSFORMATION
20 EOXG(2,1)
21 EOXG(3,D
22 KOXG(1,1)- 2ND ORD. RATE CONST. FOR OXI. TRANS.
23 KOXG(2,1)
24 KOXG(3,1)
BIOLYSIS
25 KBACWG(1,1)- 2ND ORD. RATE FOR BACT. BIOL. IN WATER
26 KBACWG(2,1)
27 KBACWG(3,1)
28 QTBAWG(1,D- CHANGE IN KBACWG DUE TO TEMP CHANGE OF 10 C
29 QTBAWG(2,1)
30 QTBAWG(3,1)
31 KBACSGd.D- 2ND ORD. RATE FOR BACT. BIOL. IN SEDIMENT
32 KBACSG(2,1)
33 KBACSG(3,1)
34 QTBASGd.D- CHANGE IN KBACSG DUE TO TEMP CHANGE OF 10 C
35 QTBASG(2,1)
36 QTBASG(3,1)
PARTITIONING
120
-------
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
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
37 KOC
38 KOW
39 OCB
40 (ALPHA(l))
41 (ALPHA(2))
42 (ALPHA(3))
VOLATILIZATION
43 MWTG
44 HENRYG
45 VAPRG
46 KVOG
47 SOLG
48 ESOLG
49 EVPRG
50 EHENG
51 (WIND)
52 (K02L)
53
PHOTOLYSIS
54 KDPG
55 RFLATG
56 CLOUDG
57 LATG
58 DFACG
59 QUANTG1
60 QUANTG2
61 QUANTG3
PART. COEFF.-ORG. CARBON
PART . COEFF . -OCTANOL-WATER
ORG. CARBON FRACTION OF BIOMASS
FRACTION OF CHEMICAL DISSOLVED (CALCULATED)
FRACTION OF CHEMICAL SORBED TO SED (CALCULATED)
FRACTION OF CHEMICAL BIO-SORBED (CALCULATED)
MOLECULAR WEIGHT OF TOXICANT
HENRY'S LAW OF THE TOXICANT
VAPOR PRESSURE OF THE TOXICANT TORR
CONST. EVAP. FOR LIQ. PHASE RESIS.
AQUEOUS SOLUBILITY
SOLUBILITY EXPONENT
MOLAR HEAT OF VAPORIZATION
CONSTANT TO COMPUTE HENRY'S LAW CONSTANT
LOCAL TIME-VARYING WIND (CALCULATED)
LOCAL TIME-VARYING REAERATION VELOC (CALCULATED)
PHOTOLYTIC RATE CONSTANT(1/HOUR)
REF. LAT. FOR KDPG
AVG CLOUD. IN TENTHS SKY COVER
GEOGRAPHIC LAT OF ECOSYSTEM
OPTICAL PATH TO VERTICAL DEPTH
REACTION QUANTUM YIELD IN PHOTO. TRANS.
SPECIAL PRINT AND CONTROL OPTIONS
62 XJTR
63 CTRIG
64 DTOPT
65 TDINT
66 LOPT
INTERNALLY SAVED
67 BURY
68 KDPL
69 KB
70 TEMPN
71 TKEL
72 PNEXT
SEG FOR TRIGGERING EVENT PRINTS AND FREQ ANALYSIS
CONCETRATIONS DEFINING EVENTS
OPTION TO OPTIMIZE TIME STEP
TIME BETWEEN NEW DECAY CALCULATIONS
OPTION TO READ IN DAILY FLOWS AND/OR LOADS
0 = SKIP DAILY READ OPTION
1 = READ DAILY FLOWS AND/OR NPS LOADS
CONSTANTS
LOCAL BURIAL / SCOUR RATE
LOCALLY CALCULATED PHOT. RATE CONSTANT
PARTITION COEFF FOR TARGET BIOMASS
TIME-VARYING NORMALIZED1 TEMPERATURE
LOCAL TEMPERATURE, DEC K
COUNTER FOR NEXT EVENT PRINT CHECK
73 TCOUNT
121
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C 1=1,2,3
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
c
c
c
c
c
c
74
75
76
77
78
79
80
81
82
83
84
85
TMARK
XJSTR
PH
POH
LIGHTN
INDEXW
INDEXS
BOTLIT
TIMCHK
MOQS
NOWKS
TMASS
BENTHIC SEGMENT FOR EVENT PRINTS AND FREQ
HYDROGEN ION ACTIVITY
HYDROXYL ION ACTIVITY
NORMALIZED LIGHT INTENSITY
INDEX=1. FOR WATER, =0. FOR BED SEGMENT
INDEX=0. FOR WATER, =1. FOR BED SEGMENT
LIGHT AT BOTTOM OF SEGMENT
CHECK FOR NEW DAY
NUMBER OF DAILY FLOWS
NUMBER OF NPS LOADS
CHEMICAL MASS BURIED
ANALYSI
INDICATES DISSOLVED FORM, SEDIMENT SORBED, BIO-SORBED
PARAMETERS
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
22
23
TEMPM(J)
DEPTHG(J)
VELOC(J)
WINDG(J)
TYPEE(J)
BACTOG(J)
ACBACG(J)
BIOMAS(J)
BIOTMG(J)
POHG(J)
OXRADG(J)
OCS(J)
PCTWA(J)
DSPSED(J)
PHG(J)
WS(J)
CMPETG(J)
TOTKG(J)
DISPV(J)
NPSWK(J)
VVOL(J)
BMASS(J)
MAX SEG TEMP DEC C
AVERAGE SEG DEPTH FT
AVERAGE SEG CURRENT FT/ SEC
AVG WIND VEL. AT 10 CM M/S
SEG TYPE INDICATOR
1 SURFACE WATER
2 SUBSURFACE WATER
3 TOP BENTHIC
4 SUBSURFACE BENTHIC
BACTERIAL POPULATION
FRACTION OF BACT POP THAT IS ACTIVE
TARGET BIOMASS, MG/L (WATER) , G/SQ.M (BED)
BIOMASS INTERNAL TEMPERATURE, DEC C
HYDROXYL ION ACTIVITY
MOLAR CONCENTR. OF ENV. OXI. N/L
ORGANIC CARBON FRACTION OF SEDIMENT
SEDIMENT FRESH WEIGHT TO DRY WEIGHT
SEDIMENT TURNOVER DISPERSION RATIO
HYDROGEN ION ACTIVITY
SETTLING, BURIAL, OR PERCOLATION RATES
TYPEE=1 WS=SETTLING, M/DAY
TYPEE=2 WS=SETTLING, M/DAY
TYPEE=3 WS=SCOUR, CM/YEAR
TYPEE=4 WS=PERCOLATION, CFS
SINGLE VALUE ZENITH LIGHT EXTINC. COEFF.
PRE-CALCULATED TOTAL DEC. RATE, 1/HR
INTERNALLY CALCULATED SED-WATER DISPERSION
EXTERNALLY READ NPS LOADS
INTERNALLY SAVED VARIABLE VOLUME PARAMETERS
(INCLUDES BVOLO(J),RVOL(J))
INTERNALLY SAVED PHYSICAL MASS LOSS,, KG
122
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
™
c
c
14
TIME
1
2
3
4
5
VOLKG(J)
FUNCTIONS
TEMPN
WINDN
PHN
POHN
LIGHTN
INTERNALLY SAVED VOLATILIZATION RATE CONST
NORMALIZED TEMPERATURE
NORMALIZED WIND SPEED
NORMALIZED PH
NORMALIZED POH
NORMALIZED LIGHT INTENSITY
C
c
c
c
c
c
INCLUDE 'TOXICMN.COM'
REAL KAHL,KNHL,KBHL,KOXL,KBACWL,KBACSL,MWTG
REAL KAHG, KNHG,KBHG,KOXG,KBACWG,KBACSG,K02G,K02L,IL
REAL 1NDEXW,INDEXS,KOC,KOW,KP,KB, KVOG,TKEL,JL
REAL LATG.KDPG.KDPL,LIGHTN,LOPT,NPSWK(SY,WK),MOQS,NOWKS
DIMENSION EBHG(3,1),ENHG(3,1),EAHG(3,1),QUANTG(3,1),
KAHG(3,1),KBHG(3,1),KNHG(3,1),EOXG(3,1),KOXG(3,1),KBACWG(3,1),
KBACSG(3,1),QTBAWG(3,1),QTBASG(3,1),TEMPM(SG),DEPTHG(SG),
VELOC(SG),ALPHA(3),BIOTMG(SG),K02G(SG),PHG(SG),WS(SG),
OXRADG(SG),WINDG(SG),TYPEE(SG),TOTKG(SG),ACBACG(SG)
DIMENSION POHG(SG),OCS(SG), BACTOG(SG),BIOMAS(SG),CMPETG(SG)
DIMENSION DISPV(SG), KP(SG),PCTWA(SG),FRW(SG),DSPSED(SG)
DIMENSION VVOL(SG),BVOLO(SG),RVOL(SG),BMASS(SG),VOLKG(SG)
CONSTANTS
HYDROLYSIS
EQUIVALENCE (CONST(1),EBHG(1,1)),(CONST(2),
EBHG(2,1)),(CONST(3),EBHG(3,1)),(CONST(4),ENHG(1,1)),
(CONST(5),ENHG(2,1)),(CONST(6),ENHG(3,1)),(CONST(7),
EAHG(1,1)),(CONST(8),EAHG(2,1)),(CONST(9),EAHG(3,1))
EQUIVALENCE
(CONST(10),KAHG(!,!)),(CONST(11),KAHG(2,1)),(CONST(12),
KAHG(3,1)),(CONST(13),KBHG(1,1)),(CONST(14),KBHG(2,1)),
(CONST(15),KBHG(3,1)),(CONST(16),KNHG(1,1)),
(CONST(17),KNHG(2,1)),CCONST(18),KNHG(3,1))
OXIDATION
EQUIVALENCE (CONST(19),EOXG(1,1)),
(CONST(20),EOXG(2,1)),(CONST(21),EOXG(3,1)),(CONST(22),
123
-------
KOXG(1,1)),(CONST(23),KOXG(2,1)),(CONST(24),KOXG(3,1))
C
C BIOLYSIS
EQUIVALENCE
(CONST(25),KBACWG(1,1)),(CONST(26),KBACWG(2,1)),(CONST(27),
KBACWGC 3,1)), (CONST(28),QTBAWG(1,1)),(CONST(29),QTBAWG(2,1)),
(CONST(30),QTBAWG(3,1)),(CONST(31),KBACSG(1,1)),(CONST(32),
KBACSG(2,1)),(CONST(33),KBACSG(3,1)),(CONST(34),QTBASG(1,1)),
(CONST(35),QTBASG(2,1)),(CONST(36),QTBASG(3,1))
C
C PARTITIONING
EQUIVALENCE (CONST(37),KOC), (CONST(38),KOW),(CONST(39),OCB),
(CONST(40),ALPHA(1)),(CONST(4l),ALPHA(2)),(CONST(42),ALPHA(3))
C
C VOLATILIZATION
EQUIVALENCE
(CONST(43),MWTG),(CONST(44).HENRYG),(CONST(45),VAPRG),
(CONST(46),KVOG),(CONST(47),SOLG),(CONST(48).ESOLG),
(CONST(49),EVPRG),(CONST(50),EHENG),(CONST(51),WIND),
(CONST(52),K02L),(CONST(53),DUMMY6)
C
C PHOTOLYSIS
EQUIVALENCE (CONST(54),KDPG),
(CONST(55),RFLATG),(CONST(56),CLOUDG),(CONST(57),LATG),
(CONST(58),DFACG),(CONST(59),QUANTG(1,1)),
(CONST(60),QUANTG(2,1)),(CONST(61),QUANTG(3,1))
C
C SPECIAL PRINT OPTIONS
EQUIVALENCE (CONST(62),XJTR),(CONST(63),CTRIG),(CONST(64).DTOPT),
(CONST(65),TDINT),(CONST(66),LOPT)
C
C INTERNALLY SAVED CONSTANTS
EQUIVALENCE (CONST(67),BURY),(CONST(68),KDPL),
(CONST(69),KB),(CONST(70),TEMPN),(CONST(71),TKEL),
(CONST(72),PNEXT),(CONST(73),TCOUNT),
(CONST(74),TMARK),(CONST(75),XJSTR)
EQUIVALENCE (CONST(76),PH),(CONST(77),POH),(CONST(78),LIGHTN) ,
(CONST(79),INDEXW),(CONST(80),INDEXS),(CONST(81).BOTLIT),
(CONST(82),TIMCHK),(CONST(83),MOQS),(CONST(84),NOWKS),,
(CONST(85),TMASS)
C
C PARAMETERS
C
EQUIVALENCE (PARAM( 1,1),TEMPM(1)),(PARAM(1,2),DEPTHG(1)),
(PARAM(1,3),VELOC(1)),(PARAM(1,4),WINDG(1)),(PARAM(1,5),
TYPEE(1)),(PARAMC1,6),BACTOG(1)),(PARAM(1,7),ACBACG(1)),
(PARAM(1,8),BIOMAS(1)),(PARAM(1,9),BIOTMG(1)),
(PARAM(1,10),POHG(1)),(PARAM(1,11),OXRADG(1))
EQUIVALENCE
(PARAM(1,12),OCS(1)),(PARAMC1,13),PCTWA(1)),
CPARAM(i,u),DSPSED(i)),(PARAM(i,15),PHG(i)),
124
-------
(PARAM(1,17),CMPETG(1)),(PARAM(1,16),WS(1)),
(PARAM( 1 , 18 ) ,TOTKG( 1 ) ) , (PARAM( 1 , 19 ) ,DISPV( 1 ) )
EQUIVALENCE (PARAM(1 ,22) ,VVOL(1)) ,(VVOL(1) ,BVOLO(1 )) ,
(VVOL(1),RVOL(1))
EQUIVALENCE (PARAM( 1,14) , VOLKG( 1 ) ) , (PARAM( 1 ,23) ,BMASS( 1 ) )
EQUIVALENCE (PARAM(1 ,12) ,KP(1)) ,(PARAM(1 ,13) ,FRW(1) )
EQUIVALENCE (PARAM( 1 , 20) ,NPSWK( 1,1)) , (PARAM( 1,3) ,K02G( 1 ) )
C
C
C
C PERFORM CALCULATIONS FIRST TIME STEP
C
IF(INITB.EQ.O) CALL INIT
JTR=XJTR
JSTR=XJSTR
C
C EVALUATE PIECEWISE LINEAR FUNCTIONS OF TIME
C
IF(NFUNT.GT.TIME) GO TO 25
CALL WASP8(MFUNC,BFUNC,NFUNC,4,ITIMF,NFUNT,73)
ITIMF = ITIMF + 1
25 CONTINUE
TEMPN = MFUNC(l) * (TIME-NFUNT) + BFUNC(l)
WINDN = MFUNC(2) * (TIME-NFUNT) + BFUNC(2)
PHN = MFUNC(3) * (TIME-NFUNT) + BFUNC(3)
POHN = MFUNC(4) * (TIME-NFUNT) + BFUNC(4)
LIGHTN = MFUNC(5) * (TIME-NFUNT) + BFUNC(5)
C
C CHECK FOR PRINT INTERVAL
C
IF(IDISK.NE.O) CALL TOXIDU( J,3,V1 ,V2 ,V3,V4 ,V5,V6,V7 ,V8)
C
IF(TIME.LT.TMARK) GO TO 30
TMARK=TMARK+TDINT
DO 27 J=1,NOSEG
IF(TOTKG(J).LT.O.) TOTKG(J)=0.0
27 CONTINUE
C
30 CONTINUE
C
C COMPUTE CURRENT SEGMENT CONDITIONS
C
DO 500 J=1,NOSEG
IF(C(l,J).LT.l.E-37) C(1,J) = 0.
CHEM=C(1,J)
SED=C(2,J)
125
-------
VOL=BVOL(J)
DEPTH=DEPTHG(J)
DEPTHM=DEPTH/3.281
TCELG=TEMPN*TEMPM(J)
TKEL=TCELG+273.15
PH=PHN*PHG(J)
POH=PHN*POHG(J)
WIND=WINDN*WINDG(J)
INDEXS=0.
INDEXW=1.
IF(TYPEE(J).LE.2.) GO TO 51
INDEXS=1.
INDEXW=0.
51 CONTINUE
SEDFAC = KP(J) * SED * 1.E-06/FRW(J)
BIOFAC= KB*BIOMAS(J)*l.E-06/FRW(J)
DENOM= 1.0 +SEDFAC +BIOFAC
ALPHA(1)= 1.0/DENOM
ALPHA(2)= SEDFAC/DENOM
ALPHA(3)= BIOFAC/DENOM
CHEM1 = CHEM*ALPHA(1)
C
C CHECK SOLUBILITY LIMIT FOR DISS CHEM
C
IF(CHEM1.LT.CMAX(9)) GO TO 60
TEND = TIME
WRITE(6,59) CHEM1,CMAX(9),J
59 FORMAT(1HO/////15X,'DISSOLVED CHEMICAL CONCENTRATION OF1,
E10.3,' EXCEEDS LIMIT OF',E10.3,' IN SEGMENT1,13//
5X,'THIS LIMIT WAS SET TO ASSURE FIRST ORDER KINETICS. TO',
' BYPASS THIS CHECK, CHANGE CMAX(l) FROM 0.0 TO A LARGE',
' NUMBER.'/////25X,'THIS SIMULATION IS ABORTED!')
60 CONTINUE
C
C CALCULATE CURRENT SEGMENT DEGRADATION RATES
C
TOTKL = TOTKG(J) * 24.
IF(TOTKL.EQ.O.) CALL FIRORD(J,TOTKL,PHOTKL,HYDRKL,BIOLKL,,
OXIDKL,VOLKL)
IF(TOTKL.LT.O.) TOTKL=-TOTKL
C
CD(1,J) = -TOTKL * CHEM * VOL
C
C ADJUST DERIVATIVES FOR SETTLING, EROSION/BURIAL, AND PERCOLATION
C
CALL SETTLE(J,SETINS,SETOUS,SETINC,SETOUC)
C
C CALCULATE PORE WATER MIXING, INCLUDING SEDIMENT WATER EXCHANGE
C
IF(TYPEE(J).GE.3.)CALL SEDWAT(J,ALPH1M,ALPH2M)
126
-------
ALPH1M= ALPHA(l)
ALPH2M= ALPHA(2)
C
C SAVE CHEMICAL MASS VOLATILIZED FROM WATER
C
IF(TYPEE(J).EQ.l.) BMASS(J)=BMASS(J)+VOLKG(J)*CHEM*VOL*DT*28.317
C
C PRINT VARIABLES AND DECAY RATES
C
IF(IDISK.NE.O) CALL TOXIDU(J,1,CHEM,SED,TOTKL,PHOTKL,HYDRKL,
BIOLKL,OXIDKL,VOLKL)
C
C SAVE CONCENTRATIONS FOR STAT FILES
C
IF (J.NE.JTR) GO TO 80
CTW= CHEM
C1W= ALPHA(1)*CHEM
CSW= l.E06*ALPHA(2)*CHEM/SED
CBW= l.E06*ALPHA(3)*CHEM/BIOMAS(J)
GO TO 85
80 CONTINUE
IF (J.NE.JSTR) GO TO 85
CTB= CHEM
C1B= ALPHA(1)*CHEM
CSB= l.E06*ALPHA(2)*CHEM/SED
CBB= l.E06*ALPHA(3)*CHEM/BIOMAS(J)
85 CONTINUE
500 CONTINUE
IDISK = 0
C
C EVENT TRIGGERED PRINTOUT
C
IF (JTR.NE.O) CALL TOXIDU(J,2,CTW,C1W,CSW,CBW,CTB,C1B,CSB,CBB)
C
C READ AND ADD DAILY FLOWS AND NPS LOADS FROM AUXILIARY FILE
C
IF(LOPT.EQ.O) GO TO 600
C CHECK FOR NEW DAY
ICHK=TIMCHK
ITIMEC=TIME
ITCHK=ITIMEC-ICHK
IF(ITCHK.LE.O) GO TO 700
TIMCHK = TIME
NWKS = NOWKS
DO 540 J=1,NWKS
DO 530 I=1,NOSYS
BWK(I,J)=BWK(I,J)-NPSWK(I,J)
530 CONTINUE
540 CONTINUE
IF(NOWKS.GT.O.) READ(91) ((NPSWK(I.J),1=1,NOSYS),J=1,NWKS)
IF(MOQS.GT.O.) READ(91) (BQ(K),K=1,NOQ)
127
-------
DO 560 J=l,NWKS
DO 550 I=1,NOSYS
BWK(I,J)=BWK(I,J)+NPSWK(I,J)
550 CONTINUE
560 CONTINUE
600 CONTINUE
C
C OPTIMIZE TIME STEP
C
IF(DTOPT.EQ.O.) GO TO 700
DELT=0.5
C CHECK FLOWS
DO 650 K=l,NOQ
I-IQ(K)
J=JQ(K)
IF(I.EQ.O) I=J
IF(J.EQ.O) J=I
TQ=BQ(K) + MQ(K)*(TIME-NQT)
IF(TQ.EQ.O.) GO TO 650
TRT=0.5*(BVOL(I)+BVOL(J))/TQ
IF(TRT.LT.DELT) DELT=TRT
650 CONTINUE
C CHECK DISPERSIVE EXCHANGES
DO 660 K=l,NOR
I=IR(K)
J=JR(K)
IF(I.EQ.O) I=J
IF(J.EQ.O) J=I
TQ=BR(K) + MR(K)*(TIME-NQT)
IF(TQ.EQ.O.) GO TO 660
TRT=0.5*(BVOL(I)+BVOL(J))/TQ
IF(TRT.LT.DELT) DELT=TRT
660 CONTINUE
IF(DELT.LT.O.Ol) DELT=0.01
DT=DELT
700 CONTINUE
C
C
INITB=3
RETURN
END
C
C**************************************************************************
128
-------
c
c
C TOXICMN
C
C
C
c***********************************************************
c
c
C PDF 11/70 VERSION OF WASP
C
C THIS VERSION OF WASP IS CONFIGURED FOR A MAXIMUM OF
C 2 SYSTEMS AND 50 SEGMENTS
C
C
C***********************************************************
C
C THE FOLLOWING PARAMETER STATEMENTS ESTABLISH PARAMETERS
C FOR SYSTEMS AND SEGMENT CHANGES. CHANGES HERE WILL
C BE EVIDENT IN ALL PARTS OF THE PROGRAM WHERE THESE PARA-
C METERS ARE USED.
C
C SY = NUMBER OF SYSTEMS.
PARAMETER SY=2
C
C SG = NUMBER OF SEGMENTS.
PARAMETER SG=50
C
C S2 IS ALWAYS SET TO 2 * SG IE. SG=50, S2=100
PARAMETER 32=100
C
C BC = MAX. NO. OF BOUNDARY CONCENTRATIONS PER SYSTEM
PARAMETER BC=30
C
C WK = MAX. NO. OF FORCING FUNCTIONS PER SYSTEM (WTH'S)
PARAMETER WK=30
C
C++++++++ PARAMETER SUBSTITUTIONS FOLLOW ++++++++
C
REAL NVOLT,NRT,NQT,NBCT(SY),NWKT(SY),NFUNT
REAL MVOL(SG),MR(S2),MQ(S2),MBC(SY,BC),MWK(SY,WK),MFUNC(10)
INTEGER OUT,SYSBY(SY),RBY(SY),QBY(SY)
COMMON IN,ICRD,OUT,NOSYS,NOSEG,ISYS,ISEG,ISIM,LISTG,LISTC
COMMON INITB,IPRNT,IDUMP(8,2),IDISK,IREC,MXDMP,IDFRC(SY)
COMMON NBCPSY,NWKPSY,SYSBY,RBY,QBY,NEGSLN
COMMON TIME,DT,TZERO,SCALT,TEND,PRNT
COMMON OMEGA,ITCHCK.MXITER
COMMON C(SY,SG),CD(SY,SG),CMAX(10),CMIN(SY)
COMMON PARAM(SG,23),CONST(85)
COMMON BVOL(SG),BR(S2),BQ(S2),BBC(SY,BC),BWK(SY,WK),BFUNC(10)
129
-------
COMMON MVOL,MR,MQ,MBC,MWK,MFUNC
COMMON IR(S2) , JR(S2) ,IQ(S2) ,JQ(S2) ,IBC(SY,BC) ,IWK(SY,,WK)
COMMON IVOPT,NOV,IROPT,NOR,IQOPT,NOQ,IBCOP(SY),NOBC(SY)
COMMON IWKOP(SY) ,NOWK(SY) ,NOPAM,NCONS,NFUNC
COMMON NVOLT,NRT,NQT,NBCT,NWKT,NFUNT
COMMON ITIMV,ITIMR,ITIMQ,ITIMF,ITIMB(SY),ITIMW(SY)
COMMON /PDP/ MXSYS.MXSEG
c***************************************************************************
C
C TOXINIT
C
C LAST REVISED 03/10/83
C***************************************************************************
C
SUBROUTINE INIT
C
INCLUDE 'TOXICMN.COM'
C
DIMENSION OCS(SG) ,KP(SG) ,TYPEE(SG) ,PCTWA(SG) ,FRW(SG) ,
BIOMAS(SG) ,DEPTHG(SG) ,TOTKG(SG) .BACTOG(SG) ,
ACBACG(SG),QT(10),VELOC(SG)
DIMENSION VVOL(SG),BVOLO(SG),RVOL(SG)
C
REAL KP,KOC,KOW,KB,LOPT,NPSWK(SY,WK),MOQS,NOWKS,MWTG
C
EQUIVALENCE (CONST(37) ,KOC) , (CONST(38) ,KOW) , (CONST(39) ,OCB) ,
(CONST(62) ,XJTR) ,(CONST(63) ,CTRIG) , (CONST(69) ,KB) ,
(CONST(73) .TCOUNT) ,(CONST(74) .TMARK) ,(CONST(75) ,XJSTR) ,
(CONST(66) ,LOPT) , (CONST(81 ) ,BOTLIT) , (CONST(82) ,TIMCHK) ,
(CONST(83) ,MOQS) , (CONST(84) ,NOWKS) , (CONST(65) ,TDINT) ,
(CONST(43),MWTG),(CONST(47),SOLG)
EQUIVALENCE (CONST(53) .DUMMY6)
C
EQUIVALENCE (PARAM( 1,2) ,DEPTHG( 1 ) ) , (PARAM( 1,5) ,TYPEE( 1 ) .) ,
( PARAM ( 1 , 8 ) , B I OMAS ( 1 ) ) , ( PARAM ( 1 , 1 2 ) , 0 C S ( 1 ) ) ,
(PARAM( 1,13) ,PCTWA( 1 ) ) , (PARAM( 1 , 18) ,TOTKG( 1 ) ) ,
(PARAM( 1 ,6) ,BACTOG( 1 ) ) , (PARAM( 1 ,7) ,ACBACG( 1 ) ) ,
(PARAM(1 ,3) ,VELOC( 1 ) )
C
EQUIVALENCE (PARAM( 1 ,22) ,VVOL( 1 ) ) , (VVOL( 1 ) ,BVOLO( 1 ) ) ,
(VVOL(1),RVOL(1))
C
EQUIVALENCE (PARAM(1 ,12) ,KP(1) ) ,(PARAM(1 ,13) ,FRW(1))
EQUIVALENCE (PARAM( 1 ,20) ,NPSWK( 1,1))
C
INITB=1
130
-------
MXDMP=8
REWIND 2
TCOUNT=0.
TMARK=0.
IF(TDINT.EQ.O.O) TDINT=1.E06
C
C
C SET UP FOR READING OF NEW FLOWS AND LOADINGS
C
TIMCHK = -30000.
IP(LOPT.EQ.O.) GO TO 100
OPEN(UNIT=91,NAME='DAILY.DAT',TYPE='OLD',
FORM='UNFORMATTED',ACCESS='SEQUENTIAL')
TIMCHK = TIME
C
C READ UNIT 91 FOR CONTENT
READ(91) NWKS,MQOPT,MOQ,SCALQ
NOWKS=FLOAT(NWKS)
MOQS=FLOAT(MOQ)
IF(MOQ.EQ.O) GO TO 20
C
C READ UNIT 91 FOR FLOWS
IQOPT=1
WRITE(6,2001) SCALQ
2001 FORMAT(1H1/38X,'DAILY FLOWS READ FROM AUXILIARY FILE',
' "DAILY.DAT" '/35X,'FLOWS WILL BE MULTIPLIED BY ',E10.4,
' TO CONVERT TO CFS'/5IX,'UPSTREAM FLOWS FOLLOW :'//
IX,IOC DAY FLOW ')/lX,10(' ')/)
READ(91) «JQ(K),IQ(K)),K=1,MOQ)
JROW=0
5 CONTINUE
JROW=JROW+1
DO 10 ICOL=1,10
IF(NOWKS.GT.O.) READ(91,END=15)
((NPSWK(I,J),1=1,NOSYS),J=1,NWKS)
READ(91) (BQ(K),K=1,MOQ)
QT(ICOL)=BQ(1)/(SCALQ*0.0846)
DO 8 K=1,MOQ
MQ(K)=0.0
8 CONTINUE
10 CONTINUE
WRITE(6,2003) (((JROW-1)*10+ICOL),QT(ICOL),ICOL=1,10)
2003 FORMAKIH ,10(14,F8.1,IX))
GO TO 5
15 CONTINUE
NCOL=ICOL-1
IF(NCOL.LE.O) GO TO 18
WRITE(6,2003) (((JROW-1)*10+ICOL),QT(ICOL),ICOL=1,NCOL)
18 CONTINUE
REWIND 91
READ(91) NWKS.MQOPT.MOQ,SCALQ
131
-------
READ(91) ((JQ(K),IQ(K)),K=1,MOQ)
WRITE(6,2004) ((K,JQ(K),IQ(K)),K=1,MOQ)
2004 FORMAT(1HO//55X,'FLOW PAIRS FOLLOW :'//10X,10(10(313,3X)/))
20 CONTINUE
IF(NWKS.EQ.O) GO TO 30
C
C READ UNIT 91 FOR NFS LOADS
WRITE(6,1001) (IWK(1,I),I=1,5)
1001 FORMAT(1H1/35X,'DAILY LOADS (LB/DAY) READ FROM AUXILIARY',
' FILE "DAILY.DAT"'/40X,'NONZERO LOADS FOR THE',
1 FIRST FIVE SEGMENTS FOLLOW:1//
6X,5(3X,5(1H*),' SEGMENT',13,IX,5(1H*))/
IX,'DAY',2X,5(4X,'CHEMICAL SEDIMENT ')/
IX,' ',2X,5(3X,10(1H_),2X,10(1H_))/)
DO 25 KT=1,1000
READ(91,END=30) ((NPSWK(I,J),1=1,NOSYS),J=1,NWKS)
IF(MOQS.GT.O.) READ(91) (BQ(K),K=1,MOQ)
X=NPSWK(1,1)+NPSWK(1,2)
IF(X.EQ.O.) GO TO 25
WRITE(6,1002) KT,((NPSWK(I,J),I=1,2),J=1,5)
1002 FORMAT(1X,I3,2X,5(3X,E10.4,2X,E10.4))
25 CONTINUE
30 CONTINUE
C
C REPOSITION UNIT 91
REWIND 91
READ(91) NWKS,MQOPT,MOQ,SCALQ
IF(MOQ.EQ.O) GO TO 35
READ(91) ((JQ(K),IQ(K)),K=1,MOQ)
35 CONTINUE
DO 50 J=1,NWKS
DO 40 I=1,NOSYS
NPSWK(I,J)=0.0
40
50
100
C
C
C
CONTINUE
CONTINUE
CONTINUE
DETERMINE SED
PRINT SEGMENT
DO 205 1=1,10
JSTR=XJTR+I
IF(TYPEE(JSTR).EQ.3.)GO TO 206
205 CONTINUE
206 CONTINUE
XJSTR=JSTR
IF(XJTR.EQ.O)XJSTR=0.
C
C CALCULATE EFFECTIVE PARTITION COEFF.(LOG KOC+.09+.91*LOG KOW)
C
IF(KOC.EQ.O.O)KOC=KOW* 0.41
IF(KOW.EQ.O.O)KOW=KOC/0.41
132
-------
DO 207 J=1,NOSEG
KP(J) = KOC*OCS(J)
207 CONTINUE
KB=OCB*KOW
XTST=XJTR*CTRIG
IF(XTST.NE.O.)WRITE(6,3073)CTRIG,XJTR
3073 FORMAT(1HI,20X,'EVENT TRIGGERED PRINTOUT...CHEM EXCEEDS',
F7.3,2X,'IN SEGMENT',F4.0//13X,'TIME',7(4X,'SEG CONC'))
C
C IF THE SEGMENT IS BENTHIC, THE UNITS OF BIOMASS MUST BE
C CONVERTED FROM G/SQ.M TO MG/L, AND FRW MUST BE CALCULATED.
C FRW IS THE FRACTION WATER VOLUME TO TOTAL SEGMENT VOLUME.
C
DO 309 J=l,NOSEG
IF(BIOMAS(J).EQ.O.) BIOMAS(J)=1.E-06
IF(TYPEE(J).LE.2.) GO TO 308
FRW(J)=1.E-06*(PCTWA(J)-1.)*C(2,J)
BIOMAS(J)=BIOMAS(J)/(DEPTHG(J)*3.281)
GO TO 309
308 CONTINUE
FRW(J) =1.0
309 CONTINUE
C
C COMPUTE ACTIVE BACTERIAL POPULATION FOR BIOLYSIS. IF THE
C SEGMENT IS BENTHIC, CONVERT INPUT UNITS FROM 'CELLS/100 GRAMS
C DRY WEIGHT SEDIMENT' TO 'ACTIVE CELLS/ML PORE WATER' :
C SED(MG/ML) = (1/1000)*SED(MG/L)/FRW(J).
C
DO 320 J=1,NOSEG
BACTOG(J) = BACTOG(J)*ACBACG(J)
IF(TYPEE(J).GE.3.) BACTOG(J)=BACTOG(J)*C(2,J)/
(FRW(J)*100.*1000.*1000.)
320 CONTINUE
C
C SET INITIAL AND REFERENCE VOLUMES FOR TOP BED SEGMENTS
C
DO 330 J=1,NOSEG
IF(TYPEE(J).NE.3.) GO TO 330
M = J-l
N = J+l
BVOLO(J) = BVOL(J)
FAC = DUMMY6
RVOL(M) = BVOL(J) + FAC*BVOL(N)*C(2,N)/C(2,J)
C (RVOL FOR SEG "J" IS STORED IN LOG "M")
IF(TYPEE(N).NE.4.) RVOL(M)=1.E+10
330 CONTINUE
C
C ZERO SLOPES FOR FLOWS AND EXCHANGES
C
DO 350 K=1,NOQ
MQ(K)=0.0
133
-------
350 CONTINUE
DO 360 K=1,NOR
MR(K)=0.0
360 CONTINUE
C
C SET CMAX TO ASSURE FIRST ORDER
C
CMAX(9) = l.E+10
IF(CMAX(1).GT.O.) GO TO 400
CMAX(9) = l.E-05 * MWTG * 1000.
CMAX(10)= 0.5 * SOLG
IF(CMAX(10).LT.CMAX(9)) CMAX(9)=CMAX(10)
CMAX(l) = CMAX(9) * 1000.
WRITE(6,399) CMAX(9)
399 FORMAT(1H1/////10X,'WARNING : CMAX FOR DISSOLVED CHEMICAL',
1 HAS BEEN SET TO',E10.3,' TO ASSURE FIRST ORDER DECAY',
' ASSUMPTION'//)
400 CONTINUE
C
C CHECK FOR PROPER SEGMENT ALLIGNMENT.
C
M=0
J=l
NBED=0
ITYPE=TYPEE(J)
MTYPE=0
IF(ITYPE.NE.l) GO TO 500
DO 450 J=2,NOSEG
M=J-1
ITYPE=TYPEE(J)
MTYPE=TYPEE(M)
GO TO (410,420,430,440),ITYPE
410 CONTINUE
IF(MTYPE.LT.3) NBED=NBED+1
GO TO 450
420 CONTINUE
IF(MTYPE.GT.2) GO TO 500
GO TO 450
430 CONTINUE
IF(MTYPE.GE.3) GO TO 500
GO TO 450
440 CONTINUE
IF(MTYPE.LE.2) GO TO 500
450 CONTINUE
IF(NBED.GT.O)WRITE(6,459)NBED
459 FORMAT(1H,5X,'WARNING:',13,'WATER',
'SEGMENTS LACK BENTHIC SEGMENTS')
RETURN
500 CONTINUE
WRITE(6,2005) M.MTYPE,J.ITYPE
2005 FORMAT(1HO////40X,'SEGMENT ALLIGNMENT PROBLEM:'/40X,'SEGMENT' ,
134
-------
13,' IS TYPE',13,' WHILE SEGMENT',13,' IS TYPE',13/
40X,'EXECUTION DISCONTINUED')
STOP
END
£**************************************************************************
c**************************************************************************
C
C TOXIDUMP
C
C LAST REVISED 03/02/83
C**************************************************************************
C
SUBROUTINE TOXIDU(J,IOPT,VI,V2,V3,V4,V5,V6,V7,V8)
C
INCLUDE 'WSPCMN.F4P'
C
DIMENSION ALPHA(3), BIOMAS(SG),FRW(SG),DEPTHG(SG)
DIMENSION BMASS(SG)
C
EQUIVALENCE(CONST(72),PNEXT),(CONST(73),TCOUNT),
(CONST(62),XJTR),(CONST(75),XJSTR),(CONST(63),CTRIG),
(CONST(64),TPO),(CONST(65),TDINT),(CONST(40),ALPHA(1))
C
EQUIVALENCE (PARAM(1,8),BIOMAS(1)),(PARAM(1,13),FRW(1)),
(PARAM(1,2),DEPTHG(1))
C
EQUIVALENCE (PARAM(1,23),BMASS(1))
C
C
GO TO (100,200,50),IOPT
C
50 CONTINUE
C
C SET UP DUMPS
C
PTIME = TIME + .0001*TIME
C
C COMPUTE RECORD ADDRESS COUNTERS
IDFRC(l) = MXSEG * (IREC - 1) + 1
IDFRC(2) = IDFRC(l)
C
C CALL DISK BUFFER OPEN / CLOSE SUBROUTINE
CALL FILEOC(IO)
WRITE(IO'IREC) PTIME
RETURN
C
100 CONTINUE
C
C DUMP VARIABLES TO FILE
135
-------
VOL =BVOL(J)
CHEM=V1
SED=V2
TOTKL=V3
PHOTKL=V4
HYDRKL=V5
BIOLKL=V6
OXIDKL=V7
VOLKL=V8
CHEM1= ALPHA(1)*CHEM
CHEM2= ALPHA(2)*CHEM
CHEM3= ALPHA(3)*CHEM
CHEMW= CHEM1/FRW(J)
CHEMS= l.E+06 * CHEM2/SED
CHEMB= l.E+06 * CHEM3/BIOMAS(J)
XMASS= CHEM*VOL*28.317
CALL FILEOC(ll)
WRITE(11'IDFRC(1))CHEM,CHEMW,CHEMS,CHEMB,
ALPHA(1),ALPHA(2),XMASS,BMASS(J)
CALL FILEOC(12)
WRITE(12'IDFRC(2))SED,DEPTHG(J),TOTKL,PHOTKL,
HYDRKL,BIOLKL,OXIDKL,VOLKL
DO 77 J = 1, NOSEG
IF(TOTKG(J).LT.O.)TOTKG(J) = 0.0
77 CONTINUE
RETURN
C
200 CONTINUE
C
C EVENT-TRIGGERED PRINTOUT
C
CTW=V1
C1W=V2
CSW=V3
CBW=V4
CTB=V5
C1B=V6
CSB=V7
CBB=V8
PNEXT = PNEXT - DT*24.
IF (PNEXT.GT.O.) RETURN
PNEXT=3.
C
C DUMP VARIABLES TO STAT FILES
C
TCOUNT = TCOUNT +1.0
ICOUNT = TCOUNT
JTR=XJTR
JSTR=XJSTR
WRITE(2)TZERO,TIME,ICOUNT,JTR,JSTR
136
-------
WRITE(2)ICOUNT,TIME,CTW,C1W,CSW,CBW
WRITE(2)ICOUNT,TIME,CTB,C1B,CSB,CBB
XTST=XJTR*CTRIG
IF(XTST.EQ.O.) RETURN
IF(INITB.EQ.l) WRITE(6,473) CTRIG.XJTR
473 FORMAT(1H1,20X,'EVENT TRIGGERED PRINTOUT...CHEM EXCEEDS',
F7.3,2X,'IN SEGMENT',F4.0//13X,'TIME',7(4X,'SEG CONC'))
IF(C(1,JTR).GE.CTRIG) WRITE(6,470)TIME,(I,C(1,1),1=1,NOSEG)
470 FORMAT(1HO,9X,F10.3,7(I5,E10.3),15(/20X,7(I5,E10.3)))
RETURN
END
C
C********A*************************************A***************************
c**************************************************************************
C
c
C TOXIFORD
C
C LAST REVISED 09/24/82
C
C
SUBROUTINE FIRORD( J , TOTKL , PHOTKL , HYDRKL , BIOLKL , OXIDKL , VOLKL)
C
C
C
REAL KAHL,KNHL,KBHL,KOXL,KBACWL,KBACSL,MWTG
REAL KAHG, KNHG,KBHG,KOXG,KBACWG,KBACSG,K02G,K02L,K20
REAL INDEXW,INDEXS,KOC,KOW,KP,KB, KVOG.TKEL, JL,IL
REAL LATG,KDPG,KDPL,LIGHTN,LOPT
C
INCLUDE 'WSPCMN.F4P'
C
DIMENSION EBHG(3 , 1 ) ,ENHG(3 , 1 ) ,EAHG(3 , 1 ) ,QUANTG(3 , 1 ) ,
KAHG(3 , 1 ) ,KBHG(3 , 1 ) ,KNHG( 3,1) ,EOXG(3 , 1 ) ,KOXG( 3,1) ,KBACWG( 3 , 1 ) ,
KBACSG(3 ,1 ) ,QTBAWG(3 ,1) ,QTBASG(3,1) ,TEMPM(SG) ,DEPTHG(SG) ,
VELOC(SG),ALPHA(3),BIOTMG(SG),K02G(SG),PHG(SG),WS(SG),
OXRADG(SG) ,WINDG(SG) ,TYPEE(SG) ,TOTKG(SG) ,ACBACG(SG)
DIMENSION OCS(SG), BACTOG(SG) ,BIOMAS(SG) ,POHG(SG) ,CMPETG(SG)
DIMENSION DISPV(SG), KP(SG) ,PCTWA(SG) ,FRW(SG) »DSPSED(SG)
DIMENSION VOLKG(SG)
C
C CONSTANTS
C
C HYDROLYSIS
EQUIVALENCE (CONST(1 ) ,EBHG(1 ,1)) ,(CONST(2) ,
137
-------
EBHG(2,1)),(CONST(3),EBHG(3,1)),(CONST(4),ENHG(1,1)),
(CONST(5),ENHG(2,1)),(CONST(6),ENHG(3,1)),(CONST(7),
EAHG(1,1)),(CONST(8),EAHG(2,1)),(CONST(9),EAHG(3,1))
EQUIVALENCE
(CONST(10),KAHG(1,1)),(CONST(11),KAHG(2,1)),(CONST(12),
KAHG(3,1)),(CONST(L3),KBHG(1,1)),(CONST(14),KBHG(2,1)),
(CONST(15),KBHG(3,1)),(CONST(16),KNHG(1,1)),
(CONST(17),KNHG(2,1)),(CONST(18),KNHG(3,1))
C
C OXIDATION
EQUIVALENCE (CONST(19),EOXG(1,1)),
(CONST(20),EOXG(2,1)),(CONST(21),EOXG(3,1)),(CONST(22),
KOXG(1,1)),(CONST(23),KOXG(2,1)),(CONST(24),KOXG(3,1))
C
C BIOLYSIS
EQUIVALENCE
(CONST(25),KBACWG(1,1)),(CONST(26),KBACWG(2,1)),(CONST(27),
KBACWG(3,1)),(CONST(28),QTBAWG(1,1)),(CONST(29),QTBAWG(2,1)),
(CONST(30),QTBAWG(3,1)),(CONST(31),KBACSG(1,1)),(CONST(32),
KBACSG(2,1)),(CONST(33),KBACSG(3,1)),(CONST(34),QTBASG(1,1)),
(CONST(35),QTBASG(2,1)),(CONST(36),QTBASG(3,1))
C
C PARTITIONING
EQUIVALENCE (CONST(37),KOC), (CONST(38),KOW),(CONST(39),OCB),
(CONST(40),ALPHA(1)),(CONST(41),ALPHA(2)),(CONST(42),ALPHA(3))
C
C VOLATILIZATION
EQUIVALENCE
(CONST(43) ,MWTG) , (CONST(44) .HENRYG) ,(CONST(45) .VAPRG.),
(CONST(46),KVOG),(CONST(47),SOLG),(CONST(48),ESOLG),
(CONST(49),EVPRG),(CONST(SO).EHENG),(CONST(51),WIND),
(CONST(52),K02L),(CONST(53),DUMMY6)
C
C PHOTOLYSIS
EQUIVALENCE (CONST(54),KDPG),
(CONST(55),RFLATG),(CONST(56),CLOUDG),(CONST(57),LATG),
(CONST(58),DFACG),(CONST(59),QUANTG(1,1)),
(CONST(60),QUANTG(2,1)),(CONST(61),QUANTG(3,1))
C
C SPECIAL PRINT OPTIONS
EQUIVALENCE (CONST(62),XJTR),(CONST(63),CTRIG),(CONST(64),DTOPT),
(CONST(65),TDINT)
C
C INTERNALLY SAVED CONSTANTS
EQUIVALENCE (CONST(66),LOPT),(CONST(67),BURY),(CONST(68),KDPL),
(CONST(69),KB),(CONST(70),TEMPN),(CONST(71),TKEL),
(CONST(72),PNEXT),(CONST(73),TCOUNT),
(CONST(74),TMARK),(CONST(75),XJSTR)
EQUIVALENCE (CONST(76),PH),(CONST(77),POH),(CONST(78).LIGHTN),
(CONST(79),INDEXW),(CONST(80).INDEXS),(CONST(81),BOTLIT),
(CONST(82),TIMCHK)
138
-------
c
C PARAMETERS
C
EQUIVALENCE (PARAM(1,1),TEMPM(1)),(PARAM(1,2),DEPTHG(1)),
(PARAM(1,3),VELOC(1)),(PARAM(1,4),WINDG(1)),(PARAM(1,5),
TYPEE(l)),(PARAMO,6),BACTOG(l)),(PARAM(1,7),ACBACG(l)),
(PARAM(1,8),BIOMAS(1)),(PARAM(1,9),BIOTMG(1)),
(PARAM(1,10),POHG(1)),(PARAM(1,11),OXRADG(1))
EQUIVALENCE
(PARAM(1,12),OCS(1)),(PARAM(1,13),PCTWA(1)),
(PARAM(1,14),DSPSED(1)),(PARAM(1,15),PHG(1)),
(PARAM(1,17),CMPETG(1)),(PARAM(1,16),WS(1)),
(PARAM(1,18),TOTKG(1)),(PARAM(1,19),DISPV(1))
C
EQUIVALENCE (PARAM(1,12),KP(1)),(PARAM(1,13),FRW(1))
,(PARAM(1,3),K02G(1))
C
EQUIVALENCE (PARAM(1,14),VOLKG(1))
C
C
C*—————-.———————————.——————.———————————__ _ _ _ ___ _ _ __ __—»_ » _-__ « —.—___ __ _____
c
C ZERO OUT PROCESS RATE CONSTANTS:
C
HYDRKL=0.
OXIDKL=0.
BIOLKL=0.
VOLKL=0.
PHOTKL=0.
C*****************
C VOLATILIZATION:
C*****************
C ONLY THE DISSOLVED SH2 MOLECULE IN SURFACE WATERS CAN VOLATILIZE:
IF(TYPEE(J).GE.2.) GO TO 170
C CALCULATE CONST REAERATION FROM DEPTH AND VELOC FIRST TIME STEP
IF(INITB.GT.2) GO TO 150
CALL REOX(J,0.,K20)
K02G(J) = K20
150 CONTINUE
C CHECK FOR INFINITE RESISTANCE TO TRANSPORT (E.G. ICE COVER):
IF (WIND.EQ.O.) GO TO 170
C CHECK ON DATA AVAILABILITY:
IF(MWTG.EQ.O. .OR. (HENRYG.EQ.O. .AND. VAPRG.EQ.O.)) GO TO 170
C CALCULATE LOCAL REAERATION
K02L = K02G(J)
CALL REOX(J,WIND,K20)
IF(K20.GT.K02L) K02L=K20
C CALL 2-FILM VOLATILIZATION SUBROUTINE, RETURN 1ST-ORDER RATE CONST.:
CALL VOLAT(J,VOLKL)
VOLKG(J) = VOLKL*24.
170 CONTINUE
139
-------
c
C COMPUTE ACTIVE BACTERIAL POPULATIONS FOR ****-BIOLYSIS-***:
ACBACL = BACTOG(J)
C COMPUTE ACID/BASE CONCENTRATIONS FOR ***-HYDROLYSIS-***:
HPLUS=10.**(-PH)
HYDRX=10.**(-POH)
C INITIALIZE COUNTER FOR COMPUTING INDEX FOR ALPHA VECTOR:
II=-1
C NEUTRAL FORM OF COMPOUND ONLY
K=l
C
C PHOTOLYSIS :
C*************
C COMPUTE RATE CONSTANT, ADJUST BY QUANTUM YIELD AND DISTRIBUTION
C COEFF(ALPHA)
C CASE 1. SEDIMENT COMPARTMENT WHERE NO LIGHT PRESENT
IF(TYPEE(J).GE.3.) GO TO 180
C CASE 2. DIRECT LOAD BY USER OF RATE CONSTANT @40N LATITUDE
C THIS SUBROUTINE WILL ADJUST IT FOR LAT, TIME OF DAY, AND LAT
C BASED ON THE ZENITH LIGHT EXTINCTION COEFF.
IF(KDPG.GT.O) CALL PHOT01(KDPL,J)
180 CONTINUE
c -----------------------------------------------------------------------
c -----------------------------------------------------------------------
c
C 3 FORMS OF EACH SPECIES (DISSOLVED, SORBED WITH SEDIMENT, BIOSORBED)
C
DO 190 1=1,3
C
C************
C PHOTOLYSIS:
C************
IF (TYPEE(J).GE.3.) KDPL = 0.
C ADD CURRENT CONTRIBUTION TO RATE CONSTANT
PHOTKL = PHOTKL + ALPHA(I) * KDPL * (QUANTG(I ,K) )
C
(]** A* ********
C HYDROLYSIS:
C* ***********
C LOAD THE GLOBAL (INPUT) VALUE OF THE RATE CONSTANTS AND ADJUST
C FOR ENVIRONMENTAL TEMPERATURE, ETC.:
C ACID HYDROLYSIS:
KAHL=KAHG(I,K)
IF(EAHG(I,K).NE.O.) KAHL=10.**(KAHG(I,K)-( (1000.*
EAHG(I,K))/(4.58*TKEL)))
C NEUTRAL HYDROLYSIS:
KNHL=KNHG(I,K)
IF(ENHG(I,K).NE.O.) KNHL=10.**(KNHG(I,K)-( (1000.*
ENHG(I,K))/(4.58*TKEL)))
140
-------
C ALKALINE HYDROLYSIS:
KBHL=KBHG(I,K)
IF(EBHG(I,K).NE.O.) KBHL=10.**(KBHG(I,K)-( (1000.*
EBHG(I,K))/(4.58*TKEL)))
C
C ADD CURRENT CONTRIBUTION TO TOTAL RATE CONSTANT:
HYDRKL=HYDRKL+ALPHA( I ) * ( KAHL*HPLUS+KNHL+KBHL*HYDRX )
C
C***********
C OXIDATION:
C***********
KOXL=KOXG(I,K)
IF(EOXG(I,K).NE.O.) KOXL=10.**(KOXG(I,K)-( (1000.*
EOXG(I,K))/(4.58*TKEL)))
C ADD CURRENT CONTRIBUTION TO TOTAL RATE CONSTANT:
OXIDKL=OXIDKL+ALPHA(I)*KOXL*OXRADG(J)
C
C BIOLYSIS:
C**********
C LOAD GLOBAL (INPUT) VALUE OF RATE CONSTANTS:
KBACWL=KBACWG( I ,K)
KBACSL=KBACSG( I ,K)
IF(QTBAWG(I,K).NE.O.) KBACWL=KBACWG(I,K)*
QTBAWG(I,K)**((TEMPN*BIOTMG(J)-20.)/10.)
IF(QTBASG(I ,K) .NE .0 . ) KBACSL=KBACSG( I ,K)*
QTBASG(I,K)**((TEMPN*BIOTMG(J)-20.)/10.)
C ADD CURRENT CONTRIBUTION TO TOTAL RATE CONSTANTS:
BIOLKL=BIOLKL+ALPHA( I ) *ACBACL
*(KBACWL*INDEXW + KBACSL*INDEXS)
C
190 CONTINUE
C END OF FORMS (DISSOLVED, SORBED, BIOSORBED) LOOP
f^t _ •_—•_-_•_—— ^_ „_,,.• _.•-— -g--^ -_-^ _—-,_-_—.— __ -^-rr.— — T— TT —— — -_— r __ .»_— . — ._ -—-_-..„— -^f-g-.—-.—— .-— ,— — — • -~ T--._ —
C
C INCREMENT ALPHA INDEX COUNTER:
200 CONTINUE
11=11+1
210 CONTINUE
C
C COMPUTE TOTAL FIRST-ORDER RATE COEFFICIENT:
C
TOTKL=( VOLKL+BIOLKL+OXIDKL+HYDRKL+PHOTKL) *24 .
TOTKG(J)=-TOTKL/24.
C
RETURN
END
C
Q**************************************************************************
141
-------
C***************************************************A****A*A***************
C
C
C TOXIVOLT
C
C LAST REVISED 08/05/82
C**********************A*****************A*********************************
C
SUBROUTINE VOLAT(K,RATEK)
("'__________________ __ __________ _ ___ __
C VOLAT. VOLAT COMPUTES VOLATILIZATION RATE CONSTANTS USING A TWO-FILM
C MODEL OF MOVEMENT OF TOXICANT ACROSS THE AIR-WATER INTERFACE.
C
C BETA IS THE FRACTION OF THE TOXICANT PRESENT AS THE NEUTRAL DISSOLVED
C MOLECULE, RATEK IS THE FIRST-ORDER RATE COEFFICIENT TO BE COMPUTED.
C K IS THE NUMBER OF THE CURRENT COMPARTMENT.
C
C LOCAL VARIABLES FOR THIS SUBROUTINE:
C K02L, BETA,COND,HENRYL,RATEK,RESGAS,RESLIQ,SOLL,TKEL,VAPRL,WAT
C
REAL MWTG, K02L,KVOG
C
INCLUDE 'WSPCMN.F4P'
C
DIMENSION WINDG(SG), DEPTHG(SG)
EQUIVALENCE (CONST(40),BETA),(CONST(43),MWTG),(CONST(44) .HENRYG),
(CONST(45),VAPRG), (CONST(46), KVOG),
(CONST(47),SOLG), (CONST(48),ESOLG),
(CONST(49), EVPRG), (CONST(50),EHENG),
(CONST(51),WIND),(CONST(52),K02L),(CONST(71),,TKEL)
EQUIVALENCE (PARAM(1,2), DEPTHG(l)), (PARAM(1,4),WINDG(1))
O_ _____ _ __________ _ _____ — ________________ _ _____..________________..__--___-
C EVALUATE OXYGEN TRANSFER RATE FROM INPUT DATA:
C REAERATION IN CM/HR @ 20 DEG. C. IN INPUT DATA, HERE CONVERTED TO
C UNITS OF METERS/HOUR AND ADJUSTED FOR LOCAL TEMPERATURE
K02L=(K02L/100.)*(1.024**((TKEL-273.15)-20.))
C — —
C LIQUID PHASE RESISTANCE: USE EXPERIMENTAL VALUE IF AVAILABLE,
C OTHERWISE ESTIMATE FROM MOLECULAR WEIGHT:
IF(KVOG.GT.O.) GO TO 100
C COMPUTE LIQUID PHASE RESISTANCE:
RESLIQ=1./(K02L*SQRT(32./MWTG))
GO TO 150
100 CONTINUE
RESLIQ=1./(K02L*KVOG)
C
150 CONTINUE
C COMPUTE GAS PHASE RESISTANCE:
142
-------
*-l _______—— — -._____ — — — — _______._.__.____
IF(HENRYG.GT.O.) GO TO 500
C
C COMPUTE HENRY'S LAW CONSTANT FROM VAPOR PRESSURE AND SOLUBILITY:
C CONVERT VAPOR PRESSURE TO ATMOSPHERES:
VAPRL=VAPRG/760.
C
C TEMPERATURE CORRECTION:
IF(EVPRG.NE.O.) VAPRL=(10.**(VAPRG-((1000.*EVPRG)/(4.58*TKEL)
2 )))/760.
C
C COMPUTE SOLUBILITY AS MOLES PER CUBIC METER:
C
SOLL=SOLG
C TEMPERATURE CORRECTION:
IF(ESOLG.NE.O.) SOLL=(10.**(SOLG-((1000.*ESOLG)/
2 (4.58*TKEL))))*MWTG*1000.
C
C CONVERT MG/L TO MOLES/CU. M.
SOLL=SOLL/MWTG
C
HENRYL=VAPRL/SOLL
C
GO TO 600
500 CONTINUE
HENRYL=HENRYG
IF(EHENG.NE.O.) HENRYL=10.**(HENRYG-((1000.*EHENG)/(4.58*TKEL)
D)
600 CONTINUE
C
C COMPUTE WATER VAPOR TRANSPORT TERM A FUNCTION OF WINDSPEED
C (LISS 1973. DEEP-SEA RESEARCH 20:221)
WAT=.1857 + 11.36 * WIND
C
C GAS PHASE RESISTANCE:
RESGAS = I./ ((WAT*HENRYL*SQRT(18./MWTG))/(TKEL*8.206E-5))
C (GAS CONSTANT IS 8.206E-5, MOL WT OF WATER IS 18.)
P ——— — ——— — — — — — — — — — — — —— ——— — -.—— — —— — — — __—_ ——_ — _ — — — — — __ —______«_
C COMPUTE OVERALL FIRST ORDER VOLATILIZATION COEFFICIENT:
COND=1./(RESLIO>RESGAS)
RATEK=(COND*BETA/(DEPTHG(K)/3.2 8))
C THE ADDITION TO THE ABOVE EXPRESSION ACCOUNTS FOR A FT TO METER
C CONVERSION.
Q——— — — — — — ———— — — — — — —— — — — — — — — — — — — — —— — _—_—_____ — — .«. — — — _—______
RETURN
END
143
-------
C*****************************************************************A**********
c
c
C TOXIREOX
C
C LAST REVISED 08/05/82
C****************************************************************************
C
SUBROUTINE REOX(J,WINDSG,K20)
C
INCLUDE 'WSPCMN.F4P'
C
REAL K20
C
DIMENSION DEPTHG(SG),VELOC(SG)
C
EQUIVALENCE (PARAM(1,3),VELOC(1)),(PARAM(1,2),DEPTHG(1))
C calculate oxygen reaeration
C
C
C
C DATA CFOREA/1.0/
CFOREA=1.0
C
C
C
AVDEPE = DEPTHG(J)
AVVELE = VELOC(J)
TEMPSG = 20.
C
IF(WINDSG.EQ.O.) GO TO 30
C this reach/res is a lake or reservoir
C
C calculate reaeration coefficient based on windspeed, surface
C area, and volume; WINDSG is windspeed in m/sec;
C WINDF IS IN CM/HR
C
IF (WINDSG.LE.6.0) GO TO 10
C determine windspeed factor empirically
WINDF= WINDSG *(-.46 + .136*WINDSG)
GO TO 20
10 CONTINUE
C assign value of 2.0 as windspeed factor
WINDF= 2.0
20 CONTINUE
C
C calculate reaeration coefficient
K20= ( .032808*WINDF*CFOREA/AVDEPE)
C CONVERT 1/HR TO CM/HR
K20 = K20*DEPTHG(J)/0.03281
C
144
-------
RETURN
C
30 CONTINUE
IF(AVVELE.GT.O.) GO TO 40
K20 = 0.0
RETURN
40 CONTINUE
C calculate reaeration coefficient for free-flowing reach
C
C calculate reaeration coefficient as a power function of
C average hydraulic depth and velocity; determine exponents
C to depth and velocity terms and assign value to REAR
C
IF (AVDEPE.GT.2.0) GO TO 50
C use Owen's formulation for reaeration
REAK= .906
EXPREV= 0.67
EXPRED= -1.85
GO TO 100
C
50 CONTINUE
C calculate transition depth; transition depth determines
C which method of calculation is used given the current
C velocity
C
IF (AVVELE.GE.1.7) GO TO 60
TRANDP= 0.0
GO TO 70
60 CONTINUE
TRANDP= .4263*(AVVELE**2.9135)
70 CONTINUE
C
DIF= AVDEPE- TRANDP
IF (DIP.GT.0.0) GO TO 80
C use Churchill's formulation for reaeration
REAK= .484
EXPREV= .969
EXPRED= -1.673
GO TO 90
C
80 CONTINUE
C use O'Connor-Dobbins formulation for reaeration
REAK= .538
EXPREV= 0.5
EXPRED= -1.5
C
90 CONTINUE
C
100 CONTINUE
C
C calculate reaeration coefficient
145
-------
K20 = REAK*(AVVELE**EXPREV)*(AVDEPE**EXPRED)
C CONVERT 1/HR TO CM/HR
K20 = K20*DEPTHG(J)/.03281
C
130 CONTINUE
C
C
RETURN
END
C****A**A*******AA******A*A*********AA**AA**************A*******A*A******A*
C
C
C TOXIPHOT
C
C
C*AA*****AA*********AAA*********************AAAA*A**AAAAAAAAA**************
C
SUBROUTINE PHOT01(RATEK.J)
c
C PHOT01 ACCEPTS A MEASURED (CLEAR DAY) PHOTOLYSIS RATE CONSTANT
C AT A SPECIFIED REFERENCE LATITUDE AS INPUT DATA. PHOT01 CORRECTS THIS
C RATE CONSTANT FOR DEVIATION OF THE LATITUDE OF THE ECOSYSTEM FROM THE
C REFERENCE LATITUDE, EFFECTS OF CLOUD COVER, AND LIGHT EXTINCTION IN
C THE WATER COLUMN BASED ON A SINGLE-VALUED ZENITH LIGHT EXTINCTION
C COEFFICIENT ENTERED AS PART OF THE CANONICAL ENVIRONMENTAL DATA.
C THE RATE IS FURTHER MODIFIED BY THE TIME-VARYING FACTOR 'LIGHTN1
C TO APPROXIMATE SEASONAL CHANGES IN INCIDENT LIGHT INTENSITY.
REAL LATG, KDPG, KDPL, IZERO, LIGHTL, LIGHTN, LOPT
C
C INTERNAL VARIABLES:LIGHTL,IZERO,BOTLIT,XXX,XTES
C
C
INCLUDE 'WSPCMN.F4P'
C
C
DIMENSION TYPEE(SG), DEPTHG(SG), CMPETG(SG)
EQUIVALENCE (CONST(54),KDPG), (CONST(55),RFLATG),
(CONST(56),CLOUDG), (CONST(57),LATG),
(CONST(58),DFACG),(CONST(66),LOPT),
(CONST(78),LIGHTN),(CONST(81),BOTLIT),
(CONST(82),TIMCHK)
EQUIVALENCE (PARAM(1,2), DEPTHG(l)), (PARAM(1,5),TYPEE(1)),
(PARAM(1,17), CMPETG(l))
DEPTHM = DEPTHG(J) / 3.281
146
-------
c
C EXPONENTIAL LIGHT DISAPPEARANCE PRODUCES UNDERFLOWS IF THE EXPONENT
C IS GREATER THAN 87.4 (ON MACHINES WITH RANGE E-38).
XTES=87.4
C
C THE DISTRIBUTION FACTOR (DFACG) IS THE OPTICAL PATH IN THE COMPARTMENT
C AS A RATIO TO THE DEPTH. IT IS DIFFICULT TO COMPUTE, BUT A PROBABLE
C BEST VALUE IS 1.19 (HUTCHINSON, TREATISE LIMNOLOGY). HOWEVER, IN
C THE PRESENCE OF A LARGE CONCENTRATION OF SCATTERING PARTICLES, IT
C MAY APPROACH OR REACH 2.0. IN ORDER TO ENSURE THAT AN IMPROPER VALUE
C IS NOT LOADED AND USED IN COMPUTATIONS, THE INPUT DFACG IS CHECKED
C AND SET TO 1.19 IF THE INPUT IS INVALID.
IF (DFACG.LT.l. .OR. DFACG.GT.2.) DFACG=1.19
C COMPUTE AVERAGE LIGHT LEVEL (AS FRACTION OF INCIDENT LIGHT) FROM
C THE COMPOSITE ZENITH LIGHT EXTINCTION COEFFICIENT ENTERED BY USER:
C
C INITIAL SEGMENT IS NECESSARILY A WATER COLUMN WITHOUT AN OVERLYING
C WATER MASS:
IF (TYPEE(J).GT.l.) GO TO 160
C
C COMPUTE LIGHT LEVEL AT BOTTOM OF COMPARTMENT:
C
XXX = CMPETG(J)*DEPTHM*DFACG
IF(XXX.GT.XTES) BOTLIT=0.
IF(XXX.LE.XTES) BOTLIT=EXP(-XXX)
LIGHTL=(1.-BOTLIT)/(CMPETG(J)*DEPTHM*DFACG)
GO TO 200
160 CONTINUE
C
C WATER COLUMN WITH AN OVERLYING WATER MASS: BOTLIT CONTAINS THE
C (RELATIVE) LIGHT INTENSITY AT THE TOP OF THE CURRENT COMPARTMENT:
IZERO=BOTLIT
XXX=CMPETG(J)*DEPTHM*DFACG
IF(XXX.GT.XTES) BOTLIT=0.
IF(XXX.LE.XTES) BOTLIT=IZERO*EXP(-XXX)
LIGHTL=(IZERO-BOTLIT)/(CMPETG(J)*DEPTHM*DFACG)
C
200 CONTINUE
C
C COMPUTE LATITUDE CORRECTION FACTOR:
C
FACTOR=(191696.65+87054.63*COS(0.0349*LATG))/(191696.65+
. 87054.63*COS(0.0349*RFLATG))
C
C ADJUST RATE FOR CLOUDINESS, LATITUDE, LIGHT EXTINCTION, AND VARIABILITY
C
RATEK=KDPG*LIGHTL*(1.-0.056*CLOUDG)*FACTOR*LIGHTN
C
C
RETURN
END
147
-------
c**************************************************************************
c
C TOXISEDW
C
C LAST REVISED 03/09/83
C**************************************************************************
C
SUBROUTINE SEDWAT (J,ALPH1M,ALPH2M)
C
INCLUDE 'WSPCMN.F4P'
REAL KP
DIMENSION DISPV(SG),ALPHA(3),TYPEE(SG),FRW(SG),DSPSED(SG),KP(SG)
EQUIVALENCE (CONST(40),ALPHA(1))
EQUIVALENCE (PARAM(1,13),FRW(1)),(PARAM(1,14),DSPSED(1)),
(PARAM(1,5),TYPEE(1)),(PARAM(1,19),DISPV(1)),
(PARAM(1,12),KP(1))
C
p_____ _______ ___«________,___________ ______________.»._ — __________ — _________________
C
IF (INITB.GT.l) GO TO 300
INITB=2
C
C FIND PORE WATER EXCHANGES AND CALCULATE DISPV
C
DO 200 1=1,NOR
II=NOR-I+1
I2=JR(II)
IF(I1.EQ.O.OR.I2.EQ.O)GO TO 200
TYP1=TYPEE(I1)
TYP2=TYPEE(I2)
DISP=BR(II)
IF(TYP1.LE.2.) GO TO 105
POROS=FRW(I1)
IF(TYP1.EQ.4.) POROS=(FRW(ll)+FRW(I2))/2.
DISPV(I1)=3.336E-06*DISP*POROS
BR(II)=0.0
IF(TYP1.NE.3.) DSPSED(I1)=0.0
105 CONTINUE
IF(TYP2.LE.2.) GO TO 200
POROS=FRW(I2)
IF(TYP2.EQ.4.) POROS=(FRW(ll)+FRW(I2))/2.
DISPV(I2)=3.336E-06*DISP*POROS
BR(II) = 0.0
IF(TYP2.NE.3.) DSPSED(I2)=0.0
200 CONTINUE
C BR COMES FROM WASP2 IN UNITS OF (MCF/SEC)*(CM/MI)**2. THE 3.336E-06
C PERFORMS THE CONVERSION TO THE DESIRED UNITS OF MCF/DAY.
300 CONTINUE
VOLPW=BVOL(J)*FRW(J)
148
-------
M=J-1
CHEM=C(1,J)
SED=C(2,J)
CHEMW=C(1,M)
SEDW=C(2,M)/1.E06
SEDCOL=SED*1.E-06/FRW(J)
CHEM1W=C(1,M)*ALPH1M/FRW(M)
CHEM2W=C(1,M)*ALPH2M
CHEMSW=CHEM2W/SEDW
CHEM1S=C(1,J)*ALPHA(1)/FRW(J)
CHEM2 S=C(1,J)*ALPHA( 2)
CHEMSS=1.E06*CHEM2S/SED
C
WATFL=DISPV(J)*DT
IF(WATFL.GT.VOLPW) WATFL=VOLPW
WATFL=WATFL/DT
SEDFL=WATFL*SEDCOL*DSPSED(J)
C
EXUP = WATFL * CHEM1S + SEDFL * CHEMSS
EXDO = WATFL * CHEM1W + SEDFL * CHEMSW*KP(J)/KP(M)
C
CD(1,M) = CD(1,M) + EXUP - EXDO
CD(1,J) = CD(1,J) - EXUP + EXDO
C
RETURN
END
C
£**************************************************************************
149
-------
C*****************A********************************************************
c
C TOXISETL
C
C LAST REVISED 01/12/82
C**************************************************************************
C
SUBROUTINE SETTLE(J,SETINS,SETOUS,SETINC,SETOUC)
C
Q
INCLUDE 'WSPCMN.F4P1
REAL KP
DIMENSION DEPTHG(SG),TYPEE(SG),WS(SG),FRW(SG),ALPHA(3),KP(SG)
DIMENSION VVOL(SG),BVOLO(SG),RVOL(SG),BMASS(SG)
EQUIVALENCE (CONST(67),BURY),(CONST(40),ALPHA(1)),
(CONST(41),ALPHA(2)),(CONST(85),TMASS)
(CONST(53),DUMMY6)
(PARAM(1,5),TYPEE(1)),(PARAM(1,16),WS(1)),
(PARAM(1,13),FRW(1)), (PARAM(1,2), DEPTHG(l)),
(PARAM(1,12),KP(1))
EQUIVALENCE (PARAM( 1,22) ,VVOL( 1) ) ,(VVOL( 1) ,BVOLO( 1).),
EQUIVALENCE
EQUIVALENCE
(VVOL(1),RVOL(1)),(PARAM(1,23),BMASS(1))
C
C
CHEM=C(1,J)
SED =C(2,J)
VOL =BVOL(J)
DEPTH=DEPTHG(J)
DEPTHM=DEPTHG(J)/3.281
SA=VOL/DEPTH
ITYPE=TYPEE(J)
M=J-1
N=J+1
NN=J+2
FAC=DUMMY6
GO TO (100,200,300,400), ITYPE
C**************************************************************************
C SURFACE WATER SEGMENT
C**************************************************************************
100 CONTINUE
SETINSO.
SETINC=0.
GO TO 250
C**************************************************************************
C SUBSURFACE WATER SEGMENT
C**************************************************************************
200 CONTINUE
SETINS=SETOUS
SETINC=SETOUC
150
-------
c
C ALL WATER SEGMENTS
C
250 CONTINUE
S ETOUS=SED*WS(J)*VOL/DEPTHM
SETOUC=CHEM*ALPHA(2)*WS(J)*VOL/DEPTHM
C (MASS RATE UNITS ARE MCF*MG/L PER DAY)
CD(1,J)=CD(1,J)+SETINC-SETOUC
CD(2,J)=SETINS-SETOUS
C CHECK FOR INFILTRATION
IF(TYPEE(NN).NE.4) RETURN
PERC=WS(NN)*1.£-06*86400.
IF(PERC.GE.O.) RETURN
PERCMS=PERC*ALPHA(1)*CHEM/FRW(J)
CD(1,J)=CD(1,J) + PERCMS
CD(1,N)=CD(1,N) - PERCMS
RETURN
£**************************************************************************
C TOP BENTHIC SEGMENT
£**************************************************************************
300 CONTINUE
VOLW=BVOL(M)
SETINS=SETOUS
SETINC=SETOUC
SETOUS=SED*WS(J)/(365.*100.)*VOL/DEPTHM
C (MASS RATE OF SEDIMENT SCOUR)
SETOUC=CHEM*ALPHA(2)*SETOUS/SED
C (MASS RATE OF SORBED CHEMICAL SCOUR)
CD(1,J)=CD(1,J)+SETINC-SETOUC
CD(2,M)=CD(2,M)+SETOUS
CD(1,M)=CD(1,M)+SETOUC
BURY*(SETINS-SETOUS)*DEPTHM/(VOL* SED)
C (NET BURIAL IN M/DAY)
EDEPTH=3.208*BURY*DT
EVOL=VOL*EDEPTH/DEPTH
IF(BURY.LT.O.) GO TO 350
C
C SEGMENT IS ACCRETING
C
IF(EVOL.GT.VOLW) GO TO 375
DEPTHG(J)=DEPTHG(J)+EDEPTH
DEPTHG(M)=DEPTHG(M)-EDEPTH
BVOL(J)=BVOL(J)+EVOL
BVOL(M)=BVOL(M)-EVOL
C
C CHECK FOR LAYER BURIAL
C
TMASS=0.
IF(BVOL(J).LT.RVOL(M)) GO TO 375
C (RVOL FOR SEG "J" IS STORED AT LOG "M")
FVOL = BVOL(J) - BVOLO(J)
C (FVOL = FAC*BVOL(N)*SED(N)/SED(J))
151
-------
SVOL = FVOL - FAC*BVOL(N)
C (VOL PORE WATER SQUEEZED OUT TO WATER COLUMN)
C (SVOL = FAC * BVOL(N) * (SED(N)/SED(J) - 1)
DEPTHG(M) = DEPTHG(M) + SVOL/SA
DEPTHG(J) = DEPTHG(J) - FVOL/SA
BVOL(M) = BVOL(M) + SVOL
BVOL(J) = BVOL(J) - FVOL
TMASS = FVOL * CHEM
SMASS = SVOL * ALPHA(l) * CHEM/FRW(J)
C(1,M) = C(1,M) + SMASS/BVOL(M)
TMASS = TMASS - SMASS
C (MASS CHEM BURIED TO LOWER LAYER, IN MCF*MG/L)
IF(TMASS.EQ.O.) TMASS=1.E-12
GO TO 375
C
C SEGMENT IS ERODING
C
350 CONTINUE
EDEPTH=-EDEPTH
EVOL=-EVOL
IF(EVOL.GT.VOL) GO TO 370
DEPTHG(J)=DEPTHG(J)-EDEPTH
DEPTHG(M)=DEPTHG(M)+EDEPTH
BVOL(J)=BVOL(J)-EVOL
BVOL(M)=BVOL( M)+EVOL
GO TO 375
C
C LAYER ERODED
C
370 CONTINUE
IF(TYPEE(N).NE.4.) GO TO 375
DEPTHG(J) = DEPTHG(N)
DEPTHG(M) = DEPTHG(M) + EDEPTH
BVOL(J) = BVOL(N)
BVOL(M) = BVOL(M) + EVOL
C(1,J) = C(1,N)
C(2,J) = C(2,N)
FRW(J) = FRW(N)
KP(J) = KP(N)
TMASS = -1.
375 CONTINUE
C
C CHECK FOR PERCOLATION/INFILTRATION
C
IF(TYPEE(N).NE.A.) RETURN
PERC=WS(N)*1.£-06*86400.
PERCMS=PERC*ALPHA(1)*CHEM/FRW(J)
IF(PERCMS.LT.O.) GO TO 380
CD(1,J)=CD(1,J) - PERCMS
CD(1,M)=CD(1,M) + PERCMS
RETURN
380 CONTINUE
152
-------
CD(1,J)=CD(1,J) + PERCMS
CD(1,N)=CD(1,N) - PERCMS
RETURN
C************************************************************************
C LOWER BENTHIC SEGMENT
c************************************************************************
400 CONTINUE
SETINCO.
SETINS=0.
SETOUS=0.
C
C CHECK FOR PERCOLATION/INFILTRATION
C
PERC=WS(J)*1. £-06*86400.
PERCMS=PERC*ALPHA( 1 ) *CHEM/FRW( J )
IF ( PERCMS. LT.O.) GO TO 480
CD(1 , J)=CD( 1 , J)-PERCMS
CD(1 ,M)=CD(1 ,M)+PERCMS
GO TO 500
480 CONTINUE
CD( 1 , J)=CD( 1 , JHPERCMS
IF(TYPEE(N) .EQ.4.) CD(1 ,N)=CD(1 ,N)-PERCMS
500 CONTINUE
C
C CHECK FOR BURIAL OR EROSION OF A LAYER
C
IF(TMASS.EQ.O.) RETURN
IF(TMASS.LT.O.) GO TO 525
C
C BURIAL OF LAYER
C
TMPM = FAC * C(1,J) * BVOL(J)
C(1,J) = C(1,J) + (TMASS-TMPM)/BVOL(J)
TMASS = TMPM
IF(TYPEE(N).NE.4.) BMASS(J)=BMASS(J)+TMASS*28.317
C (28.317 KG PER MCF*MG/L)
IF(TMASS.EQ.O.) TMASS=1.E-12
RETURN
525 CONTINUE
C
C EROSION OF LAYER
C
IF(TYPEE(N).EQ.4.) RETURN
C(1,J) - 0.0
TMASS = 0.0
RETURN
C
C
END
C
C**************************************************************************
153
-------
C***********************************************A************** * ***********
c
c
C TOXIFEIEQ
C
C LAST REVISED 07/30/82
C**************************************************************************
C
PROGRAM FREQ
C
C
P__
COMMON/STAORD/Z(10000)
C
DIMENSION CC(4),CUM(8,50),DIST(10),CUMZ(8)
REAL NAME(4)
C
DATA NAME/4HCTOT,4HCDIS,4HCSED,4HCBIO/
DATA DIST/.01,.05,.10,.30,.50,.70,.90,.95,.99,1.007
C
REWIND 2
C
C
10 CONTINUE
READ(2,END=25) TZERO,TIME,ICOUNT,JTR.JSTR
READ(2,END=25) ICOUNT,TIME,(CC(K),K=1,4)
READ(2,END=25) ICOUNT,TIME,(CC(K),K=1,4)
GO TO 10
25 CONTINUE
REWIND 2
C
C
IP = 1
LP = ICOUNT
NUMC = 4
NJSTAT = 2
NU = 2
IPL = 0
C
WRITE(6,600) ICOUNT,TZERO,TIME,(DIST(K),K=1,10)
C
DO 116 I1=1,NUMC
DO 114 JJ=1,NJSTAT
IPL = IPL+1
J = JTR
IF(JJ.EQ.2) J=JSTR
REWIND NU
III = JJ
KK = 0
DO 106 I=IP,LP
154
-------
READ(NU) TZERO,TIME,ICOUNT,JTR,JSTR
KK = KK + 1
DO 102 L=1,II1
READ (NU) ICYC.TIME, (CC(K), K=l,NUMC)
102 CONTINUE
Z(KK) = CC(II)
L2 = NJSTAT - JJ
IF (L2.EQ.O) GO TO 106
DO 104 M=1,L2
READ (NU) ICYC.TIME, (CC(K), K=l,NUMC)
104 CONTINUE
106 CONTINUE
CALL ORDER (IP, LP)
C
C COMPUTE STATISTICS USING ORDERED DATA
C
C MIN MAX AVG
ND = LP - IP + 1
ZAVG = 0.
DO 108 N=1,ND
ZAVG = ZAVG + Z(N)
108 CONTINUE
ZAVG = ZAVG / FLOAT(ND)
ZMIN = Z(l)
CUMZ(IPL) = Z(l)
C PERCENTILES
C PROB = FLOAT(ND) / 10.
DO 110 K=l,10
C IDIST = FLOAT(K) * PROB
IDIST = DIST(K) * FLOAT(ND)
IF(IDIST.EQ.O) IDIST=1
CUM(IPL,K) = Z(IDIST)
110 CONTINUE
C SKEWNESS + KURTOSIS
SDIFF2 = 0.
SDIFF3 = 0.
SDIFF4 = 0.
DO 112 N=1,ND
DIFF = Z(N) - ZAVG
SDIFF2 = SDIFF2 + DIFF**2
SDIFF3 = SDIFF3 + DIFF**3
SDIFF4 = SDIFF4 + DIFF**4
112 CONTINUE
XND = 1. / FLOAT(ND)
XNU = FLOAT(ND) - 1.0
SDEV = SQRT(SDIFF2/XNU)
SKEW = (XND * SDIFF3) / ((XND * SDIFF2)**!.5)
CURT = ((XND * SDIFF4) / (((XND * SDIFF2)**2)) -3)
155
-------
C OUTPUT STATISTIC TABLES
C
WRITE (6,602) NAME(II),J,ZAVG,SDEV,SKEW,ZMIN,
(CUM(IPL,K),K=1,10)
C
C PREPARE PROBABILITY TABLES
C
DO 150 K=l,50
PDIST=FLOAT(K)/50.
IDIST=PDIST*FLOAT(ND)
CUM(IPL,K)=Z(IDIST)
150 CONTINUE
114 CONTINUE
116 CONTINUE
REWIND NU
C
C OUTPUT CUM PROBABILITY TABLES
C
WRITE(6,611)((NAME(I),JTR,NAME(1),JSTR),I=1,4)
WRITE(6,612)(CUMZ(I),I=1,8)
DO 220 K=l,50
CUMPR=0.02*FLOAT(K)
WRITE(6,613)((CUM(I,K),CUMPR),I=1,8)
220 CONTINUE
C
C SET UP PLOTS
C
IFIL=3
WRITE(IFIL,612) (CUMZ(I),1=1,8)
DO 250 K=l,50
CUMPR=0.02*FLOAT(K)
WRITE(IFIL,613) ((CUM(l,K).CUMPR),1=1,8)
250 CONTINUE
C
C FORMAT STATEMENTS
C
600 FORMAT (1H //20X,83(1H*),/,1X,19(1H*),' STATISTICAL ANALYSIS OF',
*I5,IX,'CHEMICAL CONCENTRATIONS FROM DAY',F5.1,' TO DAY',F6.1,2X,
*22(1H*),/,20X,83(1H*),//,34X,34(1H*),' CUMULATIVE PROBABILITIES
*',34(1H*),/,1X,'VAR SG MEAN ST DEV SKEWNESS 0.00',
*10F9.2,/,1X,128(1H-),/ )
602 FORMAT (A4,12,3E9.2.11E9.2/)
604 FORMAT (1H1)
611 FORMAT(1H1/53X,'CUMULATIVE PROBABILITY'///
2X,8(1X,A4,I4,2X,'PROB',1X)/
3X,8(14(1H-),2X))
612 FORMATUH ,1X,8(E10.3,' 0.00',IX))
613 FORMAT(1H ,1X,8(E10.3,F5.2,1X))
STOP
END
C
c**************************************************************************
156
-------
c**************************************************************************
c
c
C TOXIORDR
C
C
C**************************************************************************
C
SUBROUTINE ORDER (IP, LP)
COMMON/STAORD/ Z(10000)
C
C ORDER THE X/ARRAY
C
NM = LP - IP
C FORWARD PASS
DO 104 11=1,NM
NCHNG = 0
DO 100 K=1,NM
IF (Z(K).LE.Z(K+1)) GO TO 100
TEMP = Z(K)
Z(K) = Z(K+1)
Z(K+1) = TEMP
NCHNG = NCHNG + 1
100 CONTINUE
IF (NCHNG.EQ.O) RETURN
C BACKWARD PASS
NCHNG = 0
DO 102 N = 1,NM
K = NM - N + 1
IF (Z(K).LE.Z(K+1)) GO TO 102
TEMP = Z(K)
Z(K) = Z(K+1)
Z(K+1) = TEMP
NCHNG = NCHNG + 1
102 CONTINUE
IF (NCHNG.EQ.O) RETURN
104 CONTINUE
RETURN
END
C
Q**************************************************************************
157
-------
APPENDIX 3:
TOXIWASP Support Software
Appendix 3.1 Listing of TOXICMP.CMD
SON WARNING CONT
SSET DEF DBO:
SDEL W/SPCMN.F4P;*
SDEL OUTCMN.F4P;*
SCOPY DBl:[55,116]
SCOP* DBl:[55,116]
SDEL *.OBJ;»
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
$FO/F4P DBl:[55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
$FO/F4P DBl:[55,11
$FO/F4P DBl: [55,11
$FO/F4P DBl: [55,11
SFO/F4P DBl: [55,11
$FO/F4P DBl:[55,11
$FO/F4P DBl:[55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,ll
SFO/F4P DBl:[55,11
$FO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
SFO/F4P DBl:[55,11
SFO/F4P DBl: [55,11
$FO/F4P DBl:[55,11
2X50.CDM WSPCMN.F4P
OUTCMM.F4P OUTCMN.F4P
6]WASPMAI
6]*'ASP8
6]WAS8A
6]WERR
6JWMESS
6JSCALP
6]SETRA
61SETIA
6JFILEDC
61FMTER
61WAS12
61EULER
6]WAS13
6DWAS14
6]WASP1
6]WASP2
6]WASP3
6]WASP4
6]WASPS
61WASP6
6]WASP7
6]WASP9
6]WAS10
6JWAS11
63WAS16
6JWAS17
63WAS19
6JOUTPUT
6JTINIT
6]WA12A
61SWFLOW
N
$FO/F4P/OBJ:WASPB DBl
SFO/F4P DB1:TOXIVOLT
DBl:TOXISETL
DBI:TOXISEDW
DBl:TOXIPHOT
DB1:TOXINIT
DBl:TOXIDUMP
DBI:TOXIFDRD
DBl:TOXIREOX
DBl:TOXIFREQ
DB1:TOXIORDR
TOXIWASPD
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
SFO/F4P
$9DBl:TOXITKB
!$DEL *.OBJ;»
158
-------
Appendix 3.2 Listing of TOXITKB.CMD, FREQTKB.CMD and TOXIODL.ODL
SON EARNING CONT
SSET DEF DBO:
SDEL TOXIFREQ.TSK;*
SLINK/OPT/READ TDXIFREQ TDXIORDR
ACTF1L=6
SON WARNING CONT
SSET DEF DBO:
SDEL TQXIWASP.TSK;*
SLINK/READ/OPT/OVER:DBI:TOXIODL/rASK:TOXIWASP
MAXBUF=1520
ACTFIL=5
UNITS=92
/
SSET DEF DBl:
SMESSAGE/TERMINALtTTl4
159
-------
CD
Jt
I
Q
CX
O
X
o
I
a.
X
o
f-l
I (O
h-t M S>
CL «t F~< ^
J I— X »• 3
J I
I X
X O
O hJ
tO (X
CO fl
M X
X O
O f-
H •>
I I
*-> _/ cQ O
x o i a. »
I O OL CO
I
I CM
ixu[ri
«•
3
ex
CX H
CX (X CX H O
H H H CJ U.
U O U bu •
CL, Cb U. • ••
• • • •• cO
tNj CO
X >* C>J SJ 3
^
I
H
_3
O
>
l-l CE
X J
o v.
f-> in
- CJ
I-H tr.
_3 CL
I ^f
H It
O ^i
X «
a "
l-l «-4
X ^
a ••
H o
v-> CO
* Q
CX CX
H H
O O
[j[, [^
.. o
6-1 •• z;
CO l-l il
160
-------
Appendix 3.3 Listing of JCL to Compile and Link Edit
//EPAXXXCL JOB (ACCTWASPM(JSR,MUSR), 'TOXIWASP CMPLNK ' , PRTY
// EXEC FTG1CL
//FORT.SYSIN UD DSN=CN.EPAXXX.ACCT.TOXIWASP.FTN,DlSP=SHR
//LKED.SYSLMOD DD UNIT=,SPACE=,DCB=,
// DSN=CN.EPAXXX.ACCT.LOAD(TOXIWASP),DISP=SHR
Appendix 3.4 Listing of TOXIC.BIS
SJOB ROBERT TOXICS 9999
SON WARNING CONT
$SET DEF DBO:
$DEL WASP.INP;*
SDEL *.XYZ;*
$DEL DBOtFQROOb.DAT;*
SDEL FOR002.DAT;*
SCOPY DBltPOND.lNP WASP.INP
SASS DBO; 6
SASS DBl: 91
$ASS DBO: 5
SRUN TOXIWASP
SCOPY DBO:FOR006.DAT POND.OUT
SPRINT/DELETE DBOrFOR006.DAT
SDELETE *.XYZ;*
SRUN TOXIFREQ
SPRINT DBO:FOR006.DAT
SEOJ
161
-------
ME=2
ON=400K
w
a
o
•H
-P
a
.9
CO
u
o
CP
c
•H
-p
U]
•H
in
n
X
•H
0)
a
13
Cx3
CC
Z
r>
cc
CL
co
<:
X
>-i
X
O
cc
X
CO
II
CL
to
•— O
CC .J
CO •
3 H
O
u
«f
o ^
ON UJ
f-i
ro U
II _3
uj u
t-l IS] O
Q M *
» tO 2
O ^ U
3 Z
O
«*
cc
UJ
co
3
II
cc
O
^Ji
a.
U
co
O
•*
cc
U
to
fi
u
J
UJ
Q
1.9 U 13
»-t E-II-I
f-iOUJ6-«
z^r>_3Z
OO-liJO
M
UJ
l!)
M ui
» CO
u
o
cc
to
ra
3T 0.
CL to
CO o
U
< X
2 U)
X
X 4
X CL
< U
a u
u to
X X
X. X.
X
X
X
«t »
CL «t
u u
• H
z n
CJ O
II CO
Z X
to co
o
a
o a
o
*-t
CO O
Ht O
•JU.
CL vo
hi O
H H
to u.
X X.
co *-* i
•• II
m a, i
m co
»-i t-t
II O
J «.,
CJ fi
[£) CO
CC I-H
j ac i
~u
•rf >-
CO t-i
(b ft
II l«
f II
u, z
CJ CO
u o
cc
s-*Q
II O
CO
u -«
o o
o
u,
o
J «-* UJ
a u J
> U] U
-tsl Q
<* CL
II to
II CC -i
Ld .J -<
CJ •• CL
UJ J II U3
u o u a
[»3 > N «.
a »M 3e
CttiiJCt
- co uj -
oOJX
WZTHZ'^Jf-f
•JOCCUOiXUZ
UCJLiJOOlAlCiO
Q » CO Ul »5/1 »O
^ [»j rD Q bu ;r> w •.
U CO II •> CO II E-i LL)
f-.JJCttJ.JJXfcJCO
UJ CC UJ h- (X I jj J J
•j «.couj «.cociia.
U ^ II J ^ II O «,
OOhJUO,J - ^->
••LOOQlOOSO
3: »> «-«.>CUfM
Ldl
F-i
U
O _3
>* tx]
CC O
U •-
CO Ui
=> f-t
II UJ
CC J
u w
CO Q
II -
^-> O
o •*
M CC
H U
z co
O 3
O II
». CC
U3 Li)
CO CO
*3 II
O U3 «-O
> Z ^ >
II W
II O O II H* O
J • CN Ul O CN
O O II LO O
O Z CM O II CM
CN ^-- »^ 01 a, >-*
^11 »•!*> IO •>
II CL o II t-t o
EJ CO J t-l "O ISJ «. OO
M a:
u cn
u co
u cn
H to
u co «. u no
H
U
CJ
>-> II
o CL
(N CO
ro h-(
II C
UJ -
M M
M CO
CO t-(
IS
CL *
- II uB >. II
L0 -
u « bl O >« U ' ' - - - -
cc in cc co cc co cc
to cc
II ^ Q
H II Q
H4 CO
to
II
II
CO
to
O II
o e-i
cc
«^ Q
II O
CO
CO CC CO
II ^ O II
M II O E-i
l-l CO hH
CJ
u
cc
H
CO
to
O II
Q e-<
CL
k. U- CO IJL. CO CO
II Ul •> H UI -
5" it *t 5. M rt
b. Z O b Z O
CJ CO CO O CO tO
t£> O >i tO O ><
CC CO CC CO
^ Q II —' Q II
II 0 E-i II O El
CO t-l OD 1-4
CO
U. t-i
II ufl
7 II
u, z
CJ CO
bJ O
cc
^ Q
H O
CO
O ^
Q O
O
U.
CO
—i O
«_> Tj<
•>. *
O II
•* U3
vD tsi
<—' t—I
II CO
U bC
CJ J
<• cc
CL »
co b,
- u
«* 3F
O U.
CO CJ
>i U
to cc
II v^
H II
l-l CD
Z U
X, X
X. X
162
-------
u:
a.
s-
II
u;
CX
IO
UJ O O
»O *™< ^3*
LU Z oJ
l-i O tO
J •> II
f.> r.i fy
O CO U
«._] tO
X Ct II
z ^ o
•^ ID i>
II «• *
CLt CN '""*
tO «-> O
l-l •-' OO
Q * II
bJ 00 N3
<£ II tO
Z til bC
a t; ij
S «* CD
O
in
«. r-i o U3
•X.
o
H Ct
• X U3 ••
>« O J ^
Lxl O
o in
o
ct.
w
ii
ct
to
II
< 3T
. Ci tO
H H
u u
o o
O II
^* Cxi
."* w-t
ua m
2- CN
V— ' \M*
II -
a o
to »
o n
Ci S*"*
II
X
u.
u
o o o
o o o
u, u. u.
—' CN 00
&* & CT*
f-i E-i H
II O
CD
CJ ^
D O
O
u.
1/1
o
ft
163
-------
Appendix 3.6 Listing of DAILYCVT
c ---------------------------------------------------------------------
c
C PROGRAM DAILYCVT
C
C LAST REVISED 08/27/82
f1— — . — — _-_.-.«_______-___ — __.__ — _. -.__ — -.__.-._-.____________ — — — ___ — _______»___-__,__________
C THIS PROGRAM CONVERTS FLOWS AND LOADINGS FROM USER SPECIFIED
C FILES INTO UNFORMATED FORM FOR READING BY TOXIWASP. THE
C USER IS QUERIED FOR THE LOAD FILE NAME FORMAT AND NUMBER OF
C SYSTEMS. THE INDICATED FILES ARE THEN READ AND DATA WRITTEN
C TO THE FILE DAILY.DAT ON THE DEFAULT DISK.
BYTE FILNM(30),FILNM2(30),FORMAT(80)
REAL LOADS(10,30)
VIRTUAL QT(10,1850)
DIMENSION IQ(10),JQ(10)
SCALQ=0.0
NQOPT=0.0
NOSYS=1
C
c
—
c
GET OPTIONS
WRITE(5,*) ' ENTER LOAD OPTION AND FLOW OPTION'
WR1TE(5,*) ' 1 = DO LOADS (OR FLOWS)'
WRITE(5,*) ' 0 = SKIP LOADS (OR FLOWS)'
READ(5,*) NOWK,NOQ
IF(NOWK.EQ.O) GO TO 190
C GET FILE NAME FOR LOADS
f _________________________ _ _ ___«.______.«___
C
1 CONTINUE
WRITE(5,*) 'ENTER NAME OF INPUT LOADINGS FILE'
READ(5,10) NC.FILNM
10 FORMAT (Q, 3 OA1)
NC = NC + 1
IF(NC.GT.30) GO TO 90
. FILNM(NC) = 0
OPEN(UNIT=1 ,NAME=FILNM,TYPE='OLD' ,ERR=90)
NC=0
GO TO 100
C
C FILE NAME ERROR
C
90 CONTINUE
WRITE(5,*) 'FILE NAME SPECIFICATION ERROR, OR FILE NOT FOUND'
164
-------
GO TO 1
C
c ------------------------------------
C GET FORMAT OF FILE
c ------------------------------------
C
100 CONTINUE
WRITE(5,*) 'ENTER FORMAT OF LOADS FILE IN PARENTHESIS. (23X,3E14.5) '
READ(5,110) FORMAT
110 FORMAT(SOAl)
C
c ------------------------------------
C GET NUMBER OF SYSTEMS
c ------------------------------------
C
115 CONTINUE
WRITE(5,*) 'ENTER NUMBER OF SYSTEMS AND LOAD SEGMENTS'
READ(5,*) NOSYS.NOWK
IF(NOSYS.LE.IO) GO TO 190
WRITE(5,*) 'SYSTEM LIMIT CURRENTLY SET AT 10 ... REENTER'
GO TO 115
C
C
C
c ------------------------------------
C GET FILE NAME FOR FLOWS
C
190 CONTINUE
IF(NOQ.EQ.O) GO TO 400
WRITE(5,*) 'ENTER NAME OF INPUT FLOW FILE'
READ(5,10) NC.FILNM2
NC = NC + 1
IF(NC.GT.30) GO TO 90
FILNM2(NC) = 0
OPEN(UNIT=3,NAME=FILNM2,TYPE='OLD' ,ERR=290)
GO TO 300
C
C FILE NAME ERROR
C
290 CONTINUE
WRITE(5,*) 'FILE NAME SPECIFICATION ERROR, OR FILE NOT FOUND'
GO TO 190
C
C READ FLOW FILES
c -------------------------------------
C
300 CONTINUE
READ(3,3001,ERR=850) NQOPT.NOQ.SCALQ
3001 FORMAT(2I5/E10.4)
IF(NOQ.LE.O) GO TO 400
165
-------
DO 380 K=1,NOQ
READ(3,3003) JQ(K) ,IQ(K) , NDAYS
3003 FORMAT(3I5)
READ(3 ,3005) ( (QT(K,N) ,XN) ,N=1 , NDAYS)
3005 FORMAT(100(4(E10.3,F10.2)/))
DO 375 N=l , NDAYS
QT(K,N)=QT(K,N)*SCALQ*0.0846
375 CONTINUE
380 CONTINUE
400 CONTINUE
C
P — — ————————————————————————— ________
C READ FILES AND CONVERT
OPEN(UNIT=2,NAME='DAILY.DAT' ,TYPE='NEW , FORM= ' UNFORMATTED ' ,
ACCESS= ' SEQUENTIAL ')
C
C
WRITE(2) NOWK,NQOPT,NOQ,SCALQ
IF(NOQ.GT.O) WRITE(2) ((JQ(K) ,IQ(K) ) ,K=1 ,NOQ)
N=0
700 CONTINUE
N=N+1
IF(NOWK.EQ.O) GO TO 710
READ (1, FORMAT, ERR=800,END=900) ((LOADS(I, J) ,1=1 ,NOSYS) , J=l ,NOWK)
WRITE(2) ((LOADS(I.J) ,1=1 ,NOSYS) ,J=1 ,NOWK)
710 CONTINUE
IF(NOQ.GT.O) WRITE(2) (QT(K,N) ,K=1 ,NOQ)
IF (N.GE. NDAYS) N=0
NEND=N+NOWK
IF(NEND.EQ.O) GO TO 900
GO TO 700
C
C READ ERROR
C
800 CONTINUE
WRITE(5,*) 'ERROR IN LOAD INPUT FILE FORMAT - RUN ABORTED1
GO TO 900
850 CONTINUE
WRITE(5,*) 'ERROR IN FLOW INPUT FILE FORMAT - RUN ABORTED'
GO TO 900
C
C END OF INPUT
C
900 CONTINUE
CLOSE (UNIT=1)
CLOSE (UNIT=2)
CLOSE (UNIT=3)
STOP
END
166
-------
Appendix 3.7 Listing of HSPFCVT
CA*AAAAAAA*AAAAAAAAAAA*AAAAAAAA*AA*AAAAAAAAAA**A*AAAAAAA**AAA*AAAAAAAA*AAAA*A
c
C PROGRAM HSPFCVT
C
C LAST REVISED 08/27/82
CAAAAAAAAAAAAAAAAAA**AA*AAA*A*AAAA*AAAAAAAAAAAAAAAAAAAAAA*AAAAAAA**AAAA*AA**A
C
C CONVERTS AN HSPF PLOT FILE TO A WASP LOAD FILE
C
C**AA*********A*******************A********************A*********************
DIMENSION IWK(30),ACRES(30),FACTOR(10),LOAD(10),INLOAD(10)
INTEGER YEAR,MONTH,DAY,DUR,BEGYR,BEGMO,BEGDA,ENDYR,ENDMO,ENDDA
REAL INLOAD.LOAD
C
WRITE(5,1005)
1005 FORMAT(1H1//////////////5X,'PLEASE ENTER ',
.'NUMBER OF SYSTEMS AND NUMBER OF LOADING SEGMENTS'//////////)
READ(5,*) NOSYS.NOWK
DO 25 I=1,NOSYS
WRITE(5,1021) I
1021 FORMAT(1HO/5X,'SYSTEM',13,' MULTIPLICATION FACTOR = ')
READ(5,*) FACTOR(I)
25 CONTINUE
FACTOR(2)=FACTOR(2)*2000.
WRITE(5,1026)
1026 FORMAT(1HO///5X,'FACTOR(2) HAS BEEN MULTIPLIED BY 2000',
' TO CONVERT TONS/ACRE TO LBS/ACRE'////)
DO 50 K=1,NOWK
WRITE(5,1031)
1031 FORMAT(1HO/5X,'SEGMENT NUMBER AND ACREAGE = ')
READ(5,*) IWK(K),ACRES(K)
50 CONTINUE
DO 75 1=1,28
READ(4,1074) DUMMY
1074 FORMAT(80A1)
7 5 CONTINUE
WRITE(5,1081)
1081 FORMAT(1HO///5X,'PLEASE ENTER BEGINNING YEAR, MONTH, AND DAY'/)
READ(5,*) BEGYR,BEGMO,BEGDA
WRITE(5,1083)
1083 FORMAT(1HO/5X,'NOW ENTER ENDING YEAR, MONTH AND DAY'/)
RE AD (5,*) ENDYR,ENDMO,ENDDA
C
C
NRECR=0
NRECW=0
ISTARTO
IENDO
100 CONTINUE
READ(4,1100,END=300) ALPHA,YEAR,MONTH,DAY,DUR,IDUM,
167
-------
(LOAD(I),I=1,NOSYS)
1100 FORMAT(A4,I6,4I3,3E14.5/5E14.5)
NRECR=NRECR+1
IF(ISTART.EQ.I) GO TO 105
IF(YEAR.GE.BEGYR.AND.MONTH.GE.BEGMO.AND.DAY.GE.BEGDA) ISTART=1
IF(ISTART.NE.l) GO TO 100
105 CONTINUE
IF(YEAR.LT.ENDYR) GO TO 110
IF(MONTH.LT.ENDMO) GO TO 110
IF(DAY.GE.ENDDA) IEND=1
110 CONTINUE
DO 120 K=1,NOWK
DO 115 I=1,NOSYS
INLOAD(I)=LOAD(I)*FACTOR(I)*ACRE S(K)
IF(INLOAD(I).LT.O.) INLOAD(I)=0.0
115 CONTINUE
WRITE(6,1200) ALPHA,YEAR,MONTH,DAY,DUR,IWK(K),
(INLOAD(I),I=1,NOSYS)
1200 FORMAT(1X,A4,I6,4I3,3E14.5/5E14.5/5E14.5)
120 CONTINUE
NRECW=NRECW+1
IF(IEND.EQ.l) GO TO 300
GO TO 100
300 CONTINUE
WRITE(5,1110) NRECR,NRECW,NOWK
WRITE(5,1120)
1110 FORMAT(1H1///////////////20X,'PROCESSING COMPLETE'/////
15X,I5,' UNIT LOADS READ'/
15X,I5,' SETS OF LOADS WRITTEN FOR EACH OF',13,' SEGMENTS'/
/////20X,'FORMAT IS 23X,3E14.5/5E14.5/5E14.5'/)
1120 FORMAT(1HO//20X,'TO CONVERT THESE LOADS TO WASP INPUT FORMAT,'/
20X,'RUN "DAILYCVT" USING INTERMEDIATE FILE "WASPLD.DAT"'/)
STOP
END
168
-------
Appendix 3.8 Listing of SWRVCVT
C*AAAAA**AA***********A***AA*AA***AA**AA*AAA***A***A**AA********AAAA**AAA**AA
C
C PROGRAM SWRVCVT
C
C LAST REVISED 03/15/82
CA**AAAAA***AAAAAAA****AAA**AA******AA***A****A*********A*AAAA*******AA*AA*A*
C
C CONVERTS AN SWRV PLOT FILE TO A WASP LOAD FILE
C
QA************A*************A**********A**A***************A***AA*************
DIMENSION IWK(30),ACRES(30),FACTOR(10),LOAD(10),INLOAD(10)
INTEGER YEAR,MONTH,DAY,DUR,BEGYR,BEGMO,BEGDA,ENDYR,ENDMO,ENDDA
REAL INLOAD.LOAD
C
WRITE(5,1005)
1005 FORMAT(1H1//////////////5X,'PLEASE ENTER ',
.'NUMBER OF SYSTEMS AND NUMBER OF LOADING SEGMENTS'//////////)
READ(5,*) NOSYS,NOWK
DO 25 I=1,NOSYS
WRITE(5,1021) I
1021 FORMAT(1HO/5X,'SYSTEM',13,' MULTIPLICATION FACTOR = ')
READ(5,*) FACTOR(I)
25 CONTINUE
FACTOR(2)=FACTOR(2)*2000.
WRITE(5,1026)
1026 FORMAT(1HO///5X,'FACTOR(2) HAS BEEN MULTIPLIED BY 2000',
1 TO CONVERT TONS/ACRE TO LBS/ACRE'////)
DO 50 K=l,NOWK
WRITE(5,1031)
1031 FORMAT(1HO/5X,'SEGMENT NUMBER AND ACREAGE = ')
READ(5,*) IWK(K),ACRES(K)
50 CONTINUE
DO 75 1=1,3
READ(4,1074) DUMMY
1074 FORMAT(80A1)
75 CONTINUE
WRITE(5,1081)
1081 FORMAT(1HO///5X,'PLEASE ENTER BEGINNING YEAR AND JULIAN DAY1/)
READ(5,*) BEGYR.BEGDA *
WRITE(5,1083)
1083 FORMAT(1HO/5X,'NOW ENTER ENDING YEAR AND JULIAN DAY'/)
READ(5,*) ENDYR.ENDDA
C
C
NRECR=0
NRECW=0
ISTART=0
IEND=0
169
-------
100 CONTINUE
READ(4,1100,END=300) ALPHA,YEAR,MONTH,DAY,DUR.IDUM,
(LOAD(I),I=1,NOSYS)
1100 FORMAT(A4,I6,4I3,2E14.5)
NRECR=NRECR+1
IF(ISTART.EQ.l) GO TO 105
IF(YEAR.GE.BEGYR.AND.DAY.GE.BEGDA) ISTART=1
IF(ISTART.NE.l) GO TO 100
105 CONTINUE
IF(YEAR.LT.ENDYR) GO TO 110
IF(DAY.GE.ENDDA) IEND=1
110 CONTINUE
DO 120 K=1,NOWK
DO 115 I=1,NOSYS
INLOAD(I)=LOAD(I)*FACTOR(I)*ACRES(K)
IF(LOAD(I).LE.0.00005) INLOAD(I)=0.0
115 CONTINUE
WRITE(6,1200) ALPHA,YEAR,MONTH,DAY,DUR,IWK(K),
(INLOAD(I),I=1,NOSYS)
1200 FORMAT(1X,A4,I6,4I3,3E14.5/5E14.5/5E14.5)
120 CONTINUE
NRECW=NRECW+1
IF(IEND.EQ.l) GO TO 300
GO TO 100
300 CONTINUE
WRITE(5,1110) NRECR,NRECW,NOWK
WRITE(5,1120)
1110 FORMAT(1H1///////////////20X,'PROCESSING COMPLETE'/////
15X.I5,' UNIT LOADS READ'/
15X.I5,' SETS OF LOADS WRITTEN FOR EACH OF',13,' SEGMENTS'/
/////20X,1FORMAT IS 23X.3E14.5/5E14.5/5E14.5'/)
1120 FORMAT(1HO//20X,'TO CONVERT THESE LOADS TO WASP INPUT FORMAT,'/
20X,'RUN "DAILYCVT" USING INTERMEDIATE FILE "WASPLD.DAT"'/)
STOP
END
170
-------
Appendix 3.9 Listing of USGSCVT
C*******************A********************************************************
C
C PROGRAM USGSCVT
C
C LAST REVISED 03/18/82
C****************************************************************************
C
C CONVERTS A USGS FLOW FILE TO WASP "DAILY" FORMAT
C
C****************************************************************************
C
VIRTUAL IQ(10),JQ(10),QT(10,1850),A(20)
C
OPEN(UNIT=1,NAME='CSU.DATI,TYPE='OLD')
OPEN(UNIT=2,NAME='WASPFL.DAT',TYPE='NEW')
C
WRITE(5,7001)
7001 FORMAT(1H1////////////////////////10X,'THIS PROGRAM WILL READ A',
. ' USGS FLOW FILE,'/10X,'THEN COPY AND INTERPOLATE DAILY FLOWS',
. ' FOR A WASP REACH'/)
WRITE(5,7003)
7003 FORMAT(1HO/10X,'PLEASE ENTER THE NUMBER OF SEGMENTS IN THE ',
. 'REACH'/)
READ(5,*) NSEG
WRITE(5,7005)
7005 FORMAT(1H /lOX.'NOW ENTEK A SCALE FACTOR TO MULTIPLY USGS FLOWS'/)
READ(5,*) SCALQ
C
C CALCULATE REACH PARAMETERS
C
ICHK=0
NOQ=NSEG+1
NQI=NSEG-1
XNQI=FLOAT(NQI)
XNSEG=FLOAT(NSEG)
IQOPT=1
WRITE(5,7007)
7007 FORMAT(1HO/10X,'PLEASE ENTER ADJOINING SEGMENT PAIRS,FROM UP-',
. ' TO DOWNSTREAM'/10X,'(BOUNDARIES ARE REPRESENTED BY A "O")'/)
DO 5 1=1,NOQ
WRITE(5,7009) I
7009 FORMAT(1H ,14X,'SEGMENTS IN PAIR',13,' ARE :')
READ(5,*) JQ(I),IQ(I)
5 CONTINUE
C
C READ USGS FILE FOR POSITION
C
WRITE(5,7011)
171
-------
7011 FORMAT(1HO/IOX,'WHAT YEAR AND MONTH DOES THE WASP FLOW RECORD',
. ' START?')
READ(5,*) IYF.IMF
WRITE(5,7013)
7013 FORMAT(1H ,9X,'WHAT YEAR AND MONTH DOES THE WASP FLOW RECORD',
. 'END?')
READ(5,*) IYL.IML
WRITE(5,7015)
7015 FORMATC1H /lOX.'WE WILL NOW SEARCH FOR THE FIRST STATION.')
C
10 CONTINUE
READ(1,8001,END=99) A
8001 FORMAT(20A4)
READ(1,8001,END=99) A
WRITE(5,7022) A
7022 FORMAT(1H ,20X,'STATION DESCRIPTION :'/1X,20A4/10X,
. 'IS THIS THE CORRECT STATION?'/15X,'1=YES, 2=NO, 3=QUITV)
READ(5,*) NSTA
GO TO (200,100,99) NSTA
C
99 CONTINUE
WRITE(5,7023)
7023 FORMAT(1HO//10X,'END OF RECORD ENCOUNTERED... PROGRAM ENDS.1)
STOP
C
C SEARCH FOR NEXT STATION IN FILE
100 CONTINUE
READ(1,8007,END=99) TEST,YEAR
IF(TEST.NE.O.OR.YEAR.NE.1826.) GO TO 100
READ(1,8001,END=99) A
GO TO 10
C
C FOUND CORRECT STATION
200 CONTINUE
WRITE(5,7025)
7025 FORMAT(1H /lOX.'IS THIS THE UP- OR DOWNSTREAM STATION?'/
10X,'1=UPSTREAM, 2=DOWNSTREAM'/)
READ(5,*) NQ
IF(NQ.GE.2) NQ=NOQ
C
C LOCATE FIRST YEAR
210 CONTINUE
READ(1,8007) TEST,YEAR
8007 FORMAT(F6.0,F4.0)
IF(TEST.NE.O.) GO TO 210
NYR=YEAR
IF(NYR.EQ.1826) GO TO 10
IF(NYR.LT.IYF) GO TO 210
WRITE(5,7027) NYR
7027 FORMAT(1H ,19X,'YEAR1,15,' LOCATED')
172
-------
c
C LOCATE FIRST MONTH
IF(IMF.EQ.l) GO TO 300
MM1=IMF-1
220 CONTINUE
READ(1,8009) MO.NDA
8009 FORMAT(I2,8X,I2)
READ(1,8011) (Q,N=1,NDA)
8011 FORMAT(4(8F10.0/))
IF(MO.LT.MMl) GO TO 220
WRITE(5,7029) IMF
7029 FORMAT(1H ,9X,'FILE POSITIONED AT MONTH',13)
C
C READ AND STORE FLOW FILE
C
300 CONTINUE
NDAYS=0
IM1=IMF
IM2=12
IYR=IYF
310 CONTINUE
IF(IYR.EQ.IYL) IM2=IML
DO 325 IMO=IM1,IM2
READ(1,8013) MO,NDA
8013 FORMAT(12,8X,12)
MDAY=NDAYS+1
NDAYS=NDAYS+NDA
READ(1,8015) ((QT(NQ,N)),N=MDAY,NDAYS)
8015 FORMAT(4(8F10.0/))
325 CONTINUE
IM1=1
IYR=IYR+1
IF(IYR.GT.IYL) GO TO 330
READ(1,8001,END=99) A
GO TO 310
330 CONTINUE
IYR=IYR-1
WRITE(5,7031) NQ,IMF,IYF,IM2,IYR
7031 FORMAT(1HO//10X,'FLOWS STORED FOR SEGMENT PAIR1,13,', FROM1,
. I3,'/',I4,' THROUGH1,13,V',14/)
ICHK=ICHK+1
IF(ICHK.GE.2) GO TO 400
WRITE(5,7033)
7033 FORMAT(1HO/10X,'WE WILL NOW SEARCH FOR THE SECOND STATION.')
GO TO 100
C
C INTERPOLATE FLOWS FOR NETWORK
C
400 CONTINUE
DO 450 N=1,NDAYS
173
-------
DELTQ=(QT(NOQ,N)-QT(1,N))/XNSEG
DO 425 K=2,NSEG
XK=FLOAT(K-1)
QT(K,N)=QT(1,N) + XK*DELTQ
425 CONTINUE
450 CONTINUE
C
C WRITE FLOWS IN WASP FORMAT
C
WRITE(2,3001) IQOPT,NOQ,SCALQ
3001 FORMAT(1H ,I4,I5/E10.4)
DO 500 K=1,NOQ
WRITE(2,3003) JQ(K),IQ(K),NDAYS
3003 FORMAT(1H ,14,215)
WRITE(2,3005) ((QT(K,N).FLOAT(N-l)),N=1,NDAYS)
3005 FORMAT(100(4(1X,E9.4,F9.2,1X)/))
500 CONTINUE
WRITE(5,7035)
7035 FORMAT(1HO//10X,'FLOWS SUCCESSFULLY WRITTEN TO WASP FLOW FILE.1)
STOP
END
PDS>
174
-------
APENDIX 4:
Listing of Input Data File
1
1
0
0
A: MODEL OPTIONS
TEST POND DIBENZOTHIOPHENE RUN
TIME VARIABLE SIMULATION
0 0
4 3
B: EXCHANGE COEFF
2.78
2.847E-05 107638. .164 .164 1
O.OOOE-05 107638. .082 .082 3
0 0
1 4
l.E-6
714285. 17857. 8928.5 8928.5
2 O.OOOE-05 107638. .164 .082
4
C: VOLUMES
2
2
0 0
2
2
2
2
1
0
1
18
2
1.0
0.10 0 1 0.10 1 0
1.0 0 1.0 1200
1
1.0
0.00 1
1.0 0 1.0 1200
1
1.0
2000. 1
1.0 0 1.0 1200
1
0.1
.0804 1
1
1.0
0. 1
1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0
D: FLOWS
E: BOUNDARIES
F: FORCING FUNCTIONS
G: SEGMENT PARAMETERS
1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0
175
-------
TEMPM
BACTO
OXRAD
WS
TEMPM
BACTO
OXRAD
SCOUR
TEMPM
BACTO
OXRAD
PERC
TEMPM
BACTO
OXRAD
PERC
66
EBHG1
ENHG3
KAHG2
KNHG1
EOXG3
KBWG2
KB SGI
QTBS3
DUMMY
KVOG
DUMMY
CLOUD
QUAN3
LOPT
5
TEMPN
WINDN
PHN
POHN
LIGHT
CHEM
SED
1 1
2 0
1.0 1.0
15.0DEPTH
1 .E05ACBAC
l.E-9 OCS
2.00 CMPET
15.0DEPTH
2.E07ACBAC
1.0 OCS
0.5CMPET
15.0DEPTH
2.E06ACBAC
0.0 OCS
O.OOCMPET
15.0DEPTH
2.E05ACBAC
0.0 OCS
O.OOCMPET
OO.OEBHG2
O.OEAHG1
O.OKAHG3
OO.OOKNHG2
O.OKOXG1
O.OKBWG3
5.3E-5KBSG2
0.000 KOC
O.ODUMMY
0.0 SOLG
DUMMY
4.00 LATG
0.0 XJTR
0.0
2
1.00 .000
2
1.00 .000
2
1.00 .000
2
1.00 .000
2
1.00 .000
0.002217 2:
200. 2:
00. 3000000.
0. 0.
20.
12131
6 . 5 6VE LOG 0 . OOWI NDS
l.OBIOMS 12.9BIOTM
0.01PCTWA 1.37DSPSD
2.0TOTKG .000000
0.164VELOC O.OWINDS
l.OBIOMS 50.0BIOTM
0.01PCTWA 1.37DSPSD
2.0TOTKG .000000
0.082VELOC O.OWINDS
l.OBIOMS l.OBIOTM
0.01PCTWA 1.10DSPSD
2.0TOTKG .000000
0.082VELOC O.OWINDS
l.OBIOMS 0.1BIOTM
0.01PCTWA 1.10DSPSD
2.0TOTKG .000000
O.OEBHG3 O.OOENHG1
O.OEAHG2 O.OEAHG3
O.OKBHG1 OO.OOKBHG2
O.OKNHG3 O.OEOXG1
O.OOOKOXG2 O.OOKOXG3
O.OQTBW1 2.0QTBW2
O.OKBSG3 O.OQTBS1
1.38E+05 KOW 0.0 OCB
MWTG 184.3HENRY
1.1 IE SOLG O.OEVPRG
FAC 0.1KDPG
40.0DFACG 1.19QUAN1
O.OCTRIG O.OODTOPT
1.00 5000.
1.00 5000.
1.00 5000.
1.00 5000.
1.00 5000.
0.000 3: 0.0000 4:
1350400.0 3: 2000000.0 4:
0. 0.
0.7TYPEE
15.0 POHG
0.0 PHG
O.OTYPEE
15.0 POHG
1 . 0 PHG
O.OTYPEE
15.0 POHG
PHG
O.OTYPEE
15.0 POHG
PHG
H:
O.OOENHG2
O.OKAHG1
O.OKBHG3
O.OEOXG2
O.OOKBWG1
2.0QTBW3
2.0QTBS2
0.003715DUMMY
O.OOVAPRG
OO.OOEHENG
0.732RFLAT
5.E-04QUAN2
O.OTDINT
I: TIME
1.0
6.0
8.0
3.0
8.0
6.0
4.0
8.0
6.0
4.0
8.0
6.0
CONSTANTS
0.0
0.0
0.0
0.0
5.3E-7
2.0
0.0
0.002
0.0
40.0
0.0
20.
FUNCTIONS
0.00500 J: INITIAL CONG.
2000000.0
K:
STABILITY
L: INT. PRINT CONTROL
421222324
M: INTEGRATION CONTROL
1.0 000.0
176
-------
CHEM
1
2
3
4
5
6
7
8
1
1
1
1
1
1
1
1
0.25
1. 0.25 800.0 0.5 6400.
CHEM1 CHEMS CHEMB DISSF SEDF MASS BMASS
2
2
2
2
2
2
2
2
3 4
3 4
3 4
3 4
3 4
3 4
3 4
3 4
JED DEPTH TOTKL PHOTKL HYDRKL BIOLKL OXIDKL VOLKL
4
4
4
4
4
4
4
4
4
4
1
1
2
3
4
6
8
1
1
1
1
1
1
1
1
1
1
1
1
1
1
3
4
7
7
8
8
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
0.
3
0.
3
0.
3
0.
3
0.
3
0.
3
0.
3
0.
3
0.
3
0.
3
3 4
3 4
3 4
3 4
3 4
3 4
N: DISPLA'
0.01
4
0.1
4
0.5
4
10.0
4
10.0
4
10.0
4
100.
4
10.
4
100.
4
10.
4
N:DISPLAY TABL
0. 1500000.0
177
-------
1
3 1
0. 4000000.0
234
1 2
0. 10.0
1
3 2
0. 1.0
234
4 3
0. 1.0
1234
4 3
0. 0.1
1234
4 3
0. .01
1234
4 3
0. .001
1234
4 6
0. .05
1234
1 8
0.0 0.01
1
END OF N CARD GROUP
*US GOVERNMENT PRINTING OFFICE 19B3 - 659-095/1908 178
------- |