xvEPA
United States
Environmental Protection
Agency
Research and Development
Environmental Research
Laboratory
Duluth MN 55804
EPA-600/3-81-044
May 1983
Documentation for
Water Quality
Analysis Simulation
Program (WASP) and
Model Verification
Program (MVP)
-------
EPA-600/3-8I-044
May 1983
DOCUMENTATION FOR
WATER QUALITY ANALYSIS SIMULATION PROGRAM (WASP)
AND
MODEL VERIFICATION PROGRAM (MVP)
Dominic M. DiToro*
James J. Fitzpatrick*
Robert V. Thomann*
Hydroscience, Inc., Westwood, N.J. 0767b
*(Presently at Hydroqual Inc., Mahwah, New Jersey 07430)
In Partial Fulfillment of
Contract No. 68-01-3872
Project Officer
William L. Richardson
Large Lakes Research Station
U.S. Environmental Protection Agency
Grosse He, Michigan 48138
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
DULUTH, MINNESOTA 55804
US. Environmental Protection Agency
Reg'on V. >J':--;.-y
230 South D"-::.\:'S\\ Street
Chicago, Illinois 6C504 . ,
-------
NOTICE
This document has been reviewed in accordance with
U.S. Environmental Protection Agency policy and
approved for publication. Mention of trade names
or commercial products does not constitute endorse-
ment or recommendation for use.
U,S. Environmental Protection Agency
11
-------
FOREWORD
One of the primary goals of the U.S. EPA's research program on Great
Lakes is to develop methodologies for simulating and predicting transport
and fate of contaminants in aquatic systems. Although methodologies have
specific research applications to the Great Lakes, they may be useful for
other water bodies. The computer programs documented in this report have
been used by researchers at Manhattan College and the EPA Large Lakes
Research Station for research ranging from transport of conservative tracers
like chloride'to complex interactions and transport of phytoplankton and
nutrients, and finally of toxic substances with suspended solids.
In documenting these programs, our intent is to allow researchers and
engineers to develop new theories and apply existing models to their prob-
lem. Operating complex, dynamic models should be approached with caution,
however. An experienced Fortran programmer should be employed to operate
and modify the computer programs to fit the computer and application con-
fronted. Experienced environmental scientists and engineers and preferably
a team consisting of biologists, limnologists, and hydrodynamicists would
ideally be involved in a program to develop and apply models and interpret
results. Modeling research should dovetail with surveillance and experi-
mental research which provides calibration and verification data and esti-
mates of model process rates.
No claims are made that the programs are applicable to every problem not
that they are error free. Procedures for obtaining copies of model codes
can be obtained by writing the project officer or contacting the EPA Water
Quality Modeling Center, ERL-Athens.
William L. Richardson
Environmental Scientist
Large Lakes Research Station
9311 Groh Road
Grosse He, Michigan 48138
-------
ABSTRACT
A generalized water quality modeling program and a model verification
analysis program have been developed that have application to a wide variety
of water resource management problems. The Water Quality Analysis Simula-
tion Program, WASP, is based on the flexible compartment approach. It may
be applied to water bodies in a one, two, or three-dimensional configura-
tion, and kinetic models may be structured to include linear and non-linear
reactions. The user may choose, via input options, to employ constant or
time-variable transport and kinetic processes, as well as point and non-
point waste discharges. The Model Verification Program, MVP, may be used
as an indicator of "goodness of fit" or adequacy of the model as a repre-
sentation of the real world.
To date, WASP has been applied to over twenty water resource management
problems. These applications have included one, two and three-dimensional
configurations and a number of different physical, chemical and biological
modeling frameworks, such as BOD-DO, eutrophication, and toxic substances.
A user's manual and program listings have been prepared. The user's
manual was oriented towards the system analyst, whose responsibility it
would be to design, develop and debug new kinetic models for end users, as
well as the end user who must prepare the data input to the program.
iv
-------
CONTENTS
Foreword iii
Abstract iv
Figures vii
Tables ix
Acknowledgments x
1. Summary and Conclusions 1
2. Recommendations 3
3. Introduction 4
Background 4
WASP 4
Applications 5
Program overview 28
References 32
4. WASP Theory 34
Mass balance equations 34
Finite differences 35
Stability and numerical error 41
Computational procedure 43
References 43
5. WASP Program Logic 46
Introduction 46
Time variable functions 47
Kinetic data 47
Units 49
WASP mainline and subroutine overview 51
WASP common 57
Writing (or revising) a WASPB kinetic subroutine 63
Example WASPB kinetic subroutine 64
References 69
6. WASP Input Structure 73
Introduction 73
Input data 73
7. MVP Theory 107
Introduction 107
Statistical theory 108
References 114
8. MVP Program Logic 115
Introduction 115
MVP mainline and subroutine overview 116
MVP common 118
MVP analysis overview, aggregation 121
-------
9. MVP Input Structure 123
Introduction 123
MVP input 124
Appendices
A. WASP Listings 137
B. MVP Listings 138
C. Major Program Modifications 139
VI
-------
FIGURES
Number Page
1 Segmentation map for the Western Delta-Suisun Bay
study 10
2 Systems diagram for the Western Delta-Suisun Bay study . . 11
3 Model calibration, observed vs. computed 12
4 Spatial scales used in the Lake Ontario analysis 13
5 Major physical features included in the LAKE! model .... 15
6 Systems diagram - LAKE1 model 16
7 Systems diagram - updated lake kinetics 17
8 LAKE 1A model calibration results (1972) 18
9 Long term chlorophyll 'a' calibration, epilimnion 19
10 Regression analyses and relative error analyses of
LAKE 1A model results vs. field data (chlorophyll) ... 21
11 Median relative error analyses of LAKE1A model for
several state-variables 22
12 Student's "t" test analyses of LAKE1 and LAKE1A model ... 23
13 Acid spill problem - lab data vs. model output 24
14 Spatial river pH response to a 72 minute spill of 96
percent sulfuric acid at 0.8 gpm 25
15 Temporal river pH response to a 72 minute spill of 96
percent sulfuric acid at 0.8 gpm
16 Schematic of major features in PCB distribution through-
out food chain
17 Total PCB and suspended solids calibration profiles
(upper Hudson River) 29
vn
-------
Number ffL'-
18 Land side simulator - WASP systems diagram ........
19 WASP/MVP program overview 31
20 Coordinate system and infinitesimal volume ........ ,Vi
21 Finite-difference approximations ..... 36
22 Finite-difference grid 37
23 Completely mixed finite segments 39
24 Cross-sectional areas and characteristic lengths 41
25 Second order runge - kutta method 4'1
26 Piecewise linear functions 48
27 Simplified WAS15 - WASPB flow charts 56
28 Simplified BOD-DO WASPB kinetic subroutine (Version 1) . . 66
29 Simplified BOD-DO WASPB kinetic subroutine (Version 2) . . 70
30 Determination of verification score V 109
31 Possible cases in regression between calculated and
observed values 113
32 Simplified MVP2A flow chart 117
vi n
-------
TABLES
Number Page
1 Partial List of WASP Applications 6
2 WASP Availability 32
3 Piecewise Linear Approximation for Figure 26(b) 47
4 Kinetic Term Units 50
5 Summary of WASP Input Card Groups 74
6 Pen Plot Example 106
7 Summary of MVP Input Card Groups 123
C-l WASP Scalar Variables 140
C-2 WASP Arrays 144
IX
-------
ACKNOWLEDGMENTS
The authors are especially grateful to William L. Richardson for his
suggestions and assistance in implementing WASP and MVP on the EPA computer
systems, and for his hospitality during their visits to Grosse lie. The
assistance provided by Dennis Parrott and Steve Golus in implementing WASP
and MVP is also acknowledged.
Appreciation is given to the following Hydroscience staff members who
contributed to this project: Richard Askins, for developing MVP; Richard
Winfield, for his assistance in implementing the LAKE 1 kinetics in WASP
and the verification data sets in MVP; Eric Ringewald, for assistance in
implementing WASP and MVP on the EPA COMNET system, Kathleen Whartenby, for
typing this document.
Recognition is also given Joan Nies for assisting in editing this docu-
mentation report and to Debra L. Caudill for retyping and formatting the
draft for final publication.
-------
SECTION 1
SUMMARY AND CONCLUSIONS
This report provides detailed documentation for the Water Quality
Analysis Simulation Program (WASP) and the Model Verification Program (MVP).
These two programs provide a generalized computational framework to be used
as part of the decision making process in water resource management prob-
lems.
This report is addressed to the needs of the systems analyst, whose re-
sponsibility it will be to design, develop and debug new kinetic models, as
well as the end user, whose main task it might be to implement an already
developed kinetic subroutine for a new river, lake, estuarine or ocean en-
vironment. As such, the report will detail the input formats needed to run
the WASP and MVP programs for the end user, and provide a complete descrip-
tion of the program logic of WASP and a detailed procedure to following for
the programming of new kinetic routines for the systems analyst.
CONCLUSIONS
WASP has proved itself to be a very versatile program, capable of study-
ing time variable or steady state, one, two or three dimensional, linear or
non-linear kinetic water quality problems. To date WASP has been employed
in over twenty modeling applications that have included river, lake, estu-
arine and ocean environments and that have investigated dissolved oxygen,
bacterial, eutrophication and toxic substance problem contexts. MVP pro-
vides a statistical framework to assist a user in determining whether a
model developed utilizing WASP is a "reasonable" representation of the real
world. MVP permits the user to make a detailed comparison of the model to
observed field data, so as to made a judgment concerning the adequacy of the
model. Naturally, the firmer the faith in the adequacy of the model the
more likely it will be able to play a useful role in the water resource
management decision making process.
Model Limitation
WASP does not compute hydrodynamics. As a result, the transport me-
chanisms, both advective and dispersive, must be specified by the user. The
user may choose to generate the transport input data from a separate program
capable of performing hydrodynamic computations; or by calibrating the model
against conservative substances such as salinity, electrical conductivity,
or dye from a dye tracer study.
1
-------
WASP cannot solve directly for steady state applications. Instead the
program must numerically integrate the differential equations involved in a
model structure until steady state is achieved, a very CPU time consumina
approach.
-------
SECTION 2
RECOMMENDATIONS
1. The Water Quality Analysis Simulation Program (WASP) should be put
to use in as many applications, both research and real problems, as pos-
sible. The evaluation of WASP and the development of new kinetic structures
are important to demonstrate the use of WASP as an acceptable tool in
addressing water quality modeling problems.
2. The Model Verification Program (MVP) can be used as one methodology
for judging model "accuracy" during the calibration and verification stages
of model development.
3. Due to the generality of the WASP program, the user should be knowl-
edgeable in water quality modeling techniques, and in particular knowledge-
able in the associated physical, chemical and biological principles that are
to be incorporated in the modeling framework. The ultimate responsibility
for a model's success lies entirely with the modeler, as WASP just provides
a computational framework for the development of a model.
4. The model developer, including both the system analyst and the pro-
grammer member of an administrator/engineering/programmer team, should be
knowledgeable in FORTRAN programming and operating systems interfacing.
5. The end user should review this report in some detail to fully
understand the principal underlying assumptions of the WASP/MVP package.
-------
SECTION 3
INTRODUCTION
BACKGROUND
The application of mathematical modeling techniques to water quality
problems has proved to be a powerful tool in water resource management. As
a diagnostic tool, it permits the abstraction of a highly complex real
world. Realizing that one can never fully expect to detail all the physical
phenomena that comprise our natural world, one attempts to identify and in-
clude only the phenomena, be they natural or man-made, that are relevant to
the water quality problem under consideration. As a predictive tool, mathe-
matical modeling permits the forecasting and evaluation of the effects of
changes in the surrounding environment on water quality. Although engineer-
ing insight and political and socio-economic concerns play important roles
in water resource management, some water quality problems are of such a
highly complex nature that the predictive capability of mathematical models
provides the only real means for screening the myriad number of management
alternatives.
It is important for a computer program that is to serve as the basis for
the mathematical modeler to be very general in nature. The program should
be flexible enough to provide the modeler with the mechanisms to describe
the kinetic process and the inputs to these processes, as well as the trans-
port processes and the geophysical morphology or setting, that go into the
framework of the model. Transport processes, basically hydrodynamic in
nature, include advection, turbulent diffusion, and, when spatial averaging
is included, dispersion. Kinetic (or reactive) processes are the sources
and sinks which act upon a particular water quality parameter and may be
physical, chemical or biological: for example, sedimentation and floccula-
tion of organics, the assimilative capacity of a water body to receive an
acid waste discharge and the predator-prey relationship of zooplankton-
phytoplankton.
WASP
While many of the water quality modeling programs available today,
DOSAG-I (Texas Water Quality Board, 1970), ES001 (EPA/Hydroscience, 1970),
QUAL-I (Masch, 1971), QUAL-II (WRE, 1974), the RECEIV block of SWMM (EPA/-
WRE), DEM (WRE), and the Hydrocomp Simulation Model (Hydrocomp), provide
flexibility in some of the aforementioned areas, no one program provides the
increased flexibility afforded by the Water Quality Analysis Simulation Pro-
-------
gram (WASP) developed by Hydroscience in 1970. WASP permits the modeler to
structure one, two, and three-dimensional models; allows the specification
of time-variable exchange coefficients, advective flows, waste loads and
water quality boundary conditions; permits the structuring of the kinetic
processes, within the larger modeling framework, without having to write or
rewrite large sections of computer code. Although WASP's multi-dimension-
ality and time-variable input capabilities are strong points, it is probably
the ease with which one may develop new kinetic or reactive structures that
is WASP's main strength. However, WASP's generality requires an additional
measure of judgment and insight on the part of the modeler. The kinetic and
transport structures are not "hard wired" in WASP (i.e., the equations are
not "fixed" and "buried" in the code). Therefore, the burden is on the
modeler (perhaps together with a programmer) to write the applicable kine-
tic equations (or use those already implemented) for a given problem con-
text.
APPLICATIONS
WASP has been used in many modeling frameworks since 1970. Table 1 is
a partial listing of the applications made to date of the WASP program. As
shown, a variety of problem contexts, varying from the more traditional D.O.
and bacterial problems, to eutrophication problems and the fate of hazardous
substances, have been studied. Spatial scales have ranged from meters to
100 km and time scales have ranged from minute to minute simulation of a
discharge of sulfuric acid to a 10 year long year to year analysis of phyto-
plankton and nutrients in Lake. Ontario. Kinetic structures have encompassed
simple linear kinetics, interactive linear kinetics, and a wide range of
non-linear interactive kinetic frameworks. One of the first applications of
WASP was in modeling eutrophication or algal population dynamics in the
Western Delta-Suisun Bay region of San Francisco Bay (1,2). The model seg-
mentation and the systems diagram of the model are shown in Figures 1 and 2
respectively. The model incorporated time variable (monthly averaged)
flows, nutrient waste loadings, time-variable boundary conditions, and spa-
tially, as well as temporally, resident extinction coefficients. In addi-
tion, the bulk dispersion coefficients were adjusted as a function of net
Delta outflow. Figure 3 presents some model calibration results for the
principal state variables.
The basic WASP structure has also been applied to eutrophication
analyses of the Great Lakes, specifically:
a. Lake Ontario, including the Rochester Embayment
b. Lake Erie
c. Lake Huron, including Saginaw Bay
d. Lake Michigan
A variety of spatial scales have been utilized in the framework as illu-
strated in Figure 4, where for Lake Ontario, the range in longitudinal scale
was from 10-100 km2 (Rochester embayment model) to 13,000 km2 (the whole
lake model). Details on these analyses are given in a series of reports on
Great Lakes eutrophication models (3,4,5,6,7).
-------
.* v.
r
- . .-
g
as
o
a*
i"t
Q.
OO
3
LL.
O
OO
^^
O
t i
1
e^
O
i (
1
CL.
O-
e^
U_
O
oo
»_H
1
<
f (
t
Q-
1
LU
_J
CO
cz£
1
oo
\x
S-
e
cu
a:
cu
cu
E re
i- O
1 OO
00
cu
JO
f& (&
Q.-1-
LJ re
c >
r
s- cu
D_ 4-1
re
4_)
OO
00
CU
S_ CU f
CU 4-1 .Q
r"> ( re re
E O -M -
=1 00 S-
2: re
=>
oo
s_ re 4-1
CU -r- C
_Q V- 4-1 CU
E 0 re E
3 Q. O
Z. CO CU
CO
>,
41
r
i re
re c.
r- O
4-* /
re oo
Q. C
oo cu
_E
Q
re
cu
S-
^^
r
re
o
r
Q.
re
^
1 Cn
! O
i
cu
CJ3
E
01
r- t4-
J3 0
O
i- S- +->
Q. O E
M- CU O
fO C7) dXO^
E C 0 r
r- !- i 1
cn-t-> cu Q-
r- 4-> > CO
j_ CU CU <
O to "O 3
^
CU
cu
3
1
NX
CU
cu
3
,
O
S-
4_>
ss
» «t
i re
i CJ
>^»i
-C i
Q.-1-
o co
S-
o «
i C
j: cu
c_5 en
r
r
CTi
oo
"re
oo c >^
4-1 O r
X -r- r
cu oo re
4-> C S-
c cu cu
o E «-»
cj -^ re
C3 I
O) CM
o
i- 0
0. re u
CU 00
C 1 S- -r-
o re ec o
^ 4_) ^
4-> i >i re
re cu re s-
O Q CO Lu
.c c c c
Q s_ zs re
o cu oo co
s- +-> <- >»
4-> oo 3 1*- re
3 CU CO O CO
UJ 3
_^
CU
cu
3
\x
CU
CU
3
,
O
S-
i OO
^
" O
i -G
i Q.
>> 01
-C 0
Q..C
0 Q-
S-
O "
^ cu
o en
CO
^O
oo
03
^
0
00
CU
^
r
Q
CM
>^
S-
re
^
i »
00
LLJ
O
re
E
O
o
0.
r
O)
T3
0
s:
i
cu
re
i
J*£ S-
cu re
cu cu
3 >-
1 T3 1
cu re re
cu cu
3 >-
\ **
o c
O I CU
hsi oo en
0 0
«-C S-
i Q. 4-J
r !
^) «« ^^
.C C
Q. O >
o 4-1 oo re
s_ _^: =3 <_>
o c s. >-
re o '
-G r- .C -r-
CJ Q. Q-CO
00
1
CM
're >,
C i
o
r- re
00 O
C !-
O) 4-1
E S-
i- cu
Q >
o
f
S-
co re
LU 4->
^ c
^ o
cu
I -^
«a: re
LLJ 1
as
r
O)
o
s:
00
cu
re
i
-^ S.
cu re
cu cu
3 >-
1 TD 1
-^ c s-
cu re re
cu cu
3 >-
1 A
O C
O 1 CU
rvi oo en
0 0
** -C ^_
i D- 4J
^ -r
>^ «« 2^
^r c:
Q. O >
o 4-J oo re
S- -* 3 O
o c s_ -r-
r re o '
-C ^ !-
C_) Q. Q.OO
oo
f*«^
^
re
c
o
^«
oo
c-
Ol
E
Q
oo
o
r
S-
re
4.)
e~
O
O)
\x
re
t
-S^ S-
cu re
cu cu
3 -o >-
\x re i_
cu re
cu cu
3 >-
1 n
O C
O 1 O)
I^sJ t/} O^
0 0
i Q. 4-1
r .|
>} *» "Z.
c~ c
Q. O «
o 4-> oo re
s- -i<: ^ o
o c s- -
i re o '
-C i -C '-
C_) Q. Q.CO
CO
CM
^
re
c
o
r
00
c~
CU
E
r-
Q
OO
S- 4-1
CU C
4-1 CU
00 E
CU >,
^ re
(J JO
o E
rr- UJ
-------
" ""
*
*I
o
CJ
*
r J
UJ
1
GO
a:
1
S-
O)
_Q
E
i.
O)
1
^^
Ol
E
1
r
re
Q.
r
(J
C
r
S-
D.
O)
4_>
tj re
o -i->
oo
i
re
i
o re
Q.
oo
i
re
f
4_>
re
Q.
oo
CO
V
i-
re
c
Ol
o:
OJ
i
re
{j
oo
CO
O)
r
O
re
r
e
re
>
Ol
re
4-1
oo
CO
Ol
_a
re
r
re
>
CO
4_>
c~
01
S
Ol
oo
+j'
r
r
re
C"
O
i
CO
01
E
Q
re
-
1 £ 1
.*: re s_
oi re
Ol O)
3 >-
|
Q.
O £Z
O 1 O)
r*-j co o)
0 0
-C 5-
i Q_ 4-1
r !
>, "Z.
JC C
Q. 0 -
o 4-1 co re
i. .* ^ u
0 C S_ -
r re o '
.C i .C -i-
cj CL CLOO
co
LT>
t-1
O
ft
CO
X
O)
J *
c-
o
CJ
E
Ol
J3
O
S-
CL
c
o
r-
^_>
re
o
i
c~
Q.
O
i-
UJ
r r
re re
c c >i
O O r-
co co re
c c s-
0) 01 0)
E E-M
i -r re
Q Q _1
CO CM
i*x
re
CO
c~
0 5
s- re
3 C
T" «r~
O)
oi re
^; oo
re
i
ir 1
-a o
CD re
CO S_ TO
O) O)
CO -!-> -r-
O) C c>_
-o »- v-
3 CO r
r -(-> C Q.
o c o E
C 01 -i- -r-
i i E 4-> OO
>s
TQ
re
OJ
4-)
OO
1
O
S
4-*
r
Z
A
r
r
>i
c"
CL
O
S-
o
^-*
5
LO
CT)
, ,
re re
c c
o o
CO CO
cr c:
OI O)
E E
o o
CO CM
>j
re
CO
O)
Ol NX
i- re
i- 0)
uj a.
re
O) co
.* Ol
re s=
_J CJ
c
o
i
i *
re
o
r- CO
-C O
CL-I-
0 4->
S- Ol
4-> C
3 *l
0)^;
Oi
re
4_>
OO
CO
^3
S-
O
-C
CL
CO
O
.c
a.
n
c:
en
>,
^^
re
s-
0)
re
i
i
0
S-
4-J
r- CO
Z^ ^5
i-
- 0
r- .C
Q.
>> CO
-E O
f^ C~
O Q-
O >
i C
-C O)
CJ en
r^
r
,
re >,
0 ^
t- re
CO O
£T !
Ol +J
II
CM
>^
C
o
1
o
S_ CJ
O)
> O)
i- O)
D: co
CO
>, O)
r- c:
C 01
r- h-
s_
1
1
0
4->
r- CO
z. ^
S-
« 0
i JC
i CL
>, CO
-c 0
CL.C
o a.
o «
sz
_E O)
CJ cn
00
^3-
,_
re
c~
O
r-
CO
c~
Ol
E
r
Q
CvJ
C"
O
4-*
to
cn
c
r
>
r
1
OJ
_Sf
fO
_J
>>J
-o
fC
cu
+J
oo
,
o
J_
4-*
1^
2:
r,
T
r
>,
a.
O
s-
o
n
o
CM
1
o
CM
r_
re
c~
O
CO
C"
ai
E
Q
oo
r-
O)
4->
re
4-1
OO
ft
co
~3
S-
o re
-C -i-
CL S-
CO Ol
O 4-*
c- C)
a_ re
CO
A
f~ A
Ol (Z)
cnS
CLOO
Q.
r
CO
CO
i
CO
CO
r
s:
s_
<*>
CL
=>
O
CM
CO
1
^
0
CL
re
Oi
C
r-
21
-------
i-
01
r"i
E
3
S-
J2
i
z.
oo
vx
c
re
E
O)
C£.
1)
<1J I
c r3
r- (J
1 00
00
oj
>^
i _Q
re re
O.T-
r- S-
o re
c~ ^>
1
S- QJ
Q- +->
re
4_>
oo
00
QJ
QJ ^~
+-> .a
tt_ re re
0 !-><-
oo s-
re
i 00
re +->
i- C
t+- +-> O>
O «3 E
ex a
00 !-
re to
Q. C
OO QJ
E
r
Q
re
QJ
£_
C£
^
re
o
i
JT
Q.
re
s_
Ol
s
C3
"O
O) 00
4-* 4-*
o c
o re
S- i
Q.
QJ O
T3 T-
3 4->
i re
CJ 3
c cr
i i re
S-
o
re
i
^_
^
o
^
8
V
cri
CO
oo
4_>
X
QJ
c
o
I
jQ
O
i-
Q_
re
'si
QJ
4^
Q.-1-
0-OL
1
re
r
O 3
4-^ E
r
"O 00
QJ
O QJ
re *o
^
re
o
i
>.
re
Q
re
r
s-
QJ
O
re
CQ
«\
8
s
Q
LO
Lf)
"re
c:
O
r
00
c
QJ
_E
Q
CO
S_
0
S-
re
j^
^_
o *
> co
o
S CM
z.
i-
o
re
s-
^
o
re
re
i
s_
OJ
4J
O
re
CQ
«
8
^
Q
CO
ps.
to
re
c~
O
r
00
f~
QJ
_E
Q
i
S-
QJ
>
i
OL
QJ
QJ
\x
^
re
^
^
SI
00
4_>
X
Q)
C
O
o
E
QJ
o
O
s-
Q.
QJ
O
re
00
ri
3
OO
O
r-"
X
o
1
>^
"O QJ
re 4->
QJ re
4-J | ^
OO OO
00 00
T3 "CQ
r- 00 O
r CQ O_
0 0
00 Q- QJ
.{_>
^3 *^3 re
QJ QJ '
73 > 3
C i O
QJ O -r-
a. oo +J
oo oo s_
3 -r- re
oo en Q.
co
r
CM
"re
c
o
r
00
QJ
£
Q
r-^
>
z.
\
^_
OJ CQ
> O
i- a.
QC
l-
c o
O J3
oo i.
~o re
=3 rc
re
i
QJ QJ
-t-> *->
3 3
C C
s: IE
00
X3 UJ
i- Q
i -0 Q
0 C
oo re QJ QJ
_4_) f~
*u ~o re re
QJ QJ i T3
"O > 3 C
C i O T-
QJ 0 - _J
Q. 00 +->
oo oo i.. TD
3 -I- re c
oo Q d. re
Lf)
Lf)
r~
^ r^
o ^
r TO
00 U
C !-
QJ +J
r- QJ
Q >
^«.
a
c:
re
UJ
Q
Q
1 QJ
c:
>, re
s- -a
s_ c
re !-
3 _l
o-
i
QJ QJ
^ ^
C C
1 'r
s: 2:
r-
r
3
^~
QJ
QJ
-t-1 re
re Q.
c
O 1
s_ E
re 3
O -r-
- s-
CQ JO
CM
co
r
re
£= >>
0 r
r r
oo re
C S-
QJ Q)
P "re
CTi 1
CM
n
n
r-
a.
oo
rc
a.
i
4-
QJ
.>
o;
i
r- O
re
s- QJ
+-> 3 OJ
oo -a c
-o c r^
QJ O O
-o -r- o
3 +J E
r re s_
o o QJ
i 1 14- -i->
i
QJ QJ
3 3
C C
s s:
QJ
oo
re
zz
T3
r
0
^
O
CM
're >,
C i
O i
r- (O
00 O
c: T-
ai +->
Sfc
Q 5>
00
1
r^
re
oo
O QJ
CL-P
oo oo
i- re
Q 3
C T3
re ,-
QJ O
O ef.
o
-------
1 *
*--
8
UJ
1
CO
h-
t/l
^vf
£_
fO
gz
d)
o:
O)
0) r
£EI ^o
r- O
1 00
CO
il CD
Q. -t->
ro
+J
OO
CO
OJ
S- o '
i~i 14. ro ro
E 0 4J -
3 00 S_
z. ro
=>
i CO
^_ f^ -l~*
cu - c:
r"i »[ +j ^j
e o ro E
rs Q. CT
2: oo oj
00
§
r ro
ro c
i- O
ro co
a. e
OO QJ
E
Q
ro
OJ
i_
^^
ro
o
i
.E
Q.
ro
i_
CT
0
QJ
cu
to
"Q
1
r^
O CO
oo ro
1 CD
CO
i
u_
o
r"
J3
0
S_
O)
ro
c
ct
CO
S- 1
O O)
+-> S-
o
ro «
O) co
s- 01
i
^J S-
O) <1J
S_ CO r
O O
c >,
i -^ O
-O O)
ro -4^
d) ro
-4-^ -ft-*
OO OO
s
n
0)
ro
S-
4-^
CO
o
00
OJ
o
CM
S-
O)
4 '
i
£
CT>
c~
r
r
s^
0
r
j_
h-
o
i
3
CT
p
I
I O>
CO
co ro
ro .c
CD Q.
^1
ro
O)
^_)
oo
c
O)
o
S-
1
"^
8
A
Q
81.
to
'
^~
c:
O)
CO
>^
X
o
Ol
s-
3
Q-
-------
to
I
to
QJ
T3
CD
2
o;
s-
o
cu
CD
oi
CO
O)
S-
3
O1
10
-------
E O
O -
to
-Q
OJ
"O
S-
0)
-(->
1/1
CD
O)
S-
£
HE
S-
cn
(0
O)
O)
CD
11
-------
0
120
u
? 1*0
t-
240
300
360
COdPUTID
A-OMMVIO
j I
«0
120
180
240
300
20 40 60 80 100
CHLOROPHYLL (/ig/l)
360
0.3 0.6 0.9 1.2 1.5
NO,-N (mg/l)
60
120
180
240
300
360
60
120
ISO
240
300
5 10 IS 20
SILICA (mg/l)
360
29 0
3 6 9 12 15
DISSOLVED OXYGEN (mg/l)
Figure 3. Model calibration observed vs. computed.
12
-------
MODEL
DESIGNATION
HORIZONTAL
NUMBER OF SCALE (km2)
SEGMENTS EPI LIMN ION
SEGMENTS
LAKE
LAKE 3
LAKE 3
(AGGREGATED)
67
13,000
200-1000
6000-13,000
ROCHESTER
EMBAYMENT
72
10-100
Figure 4. Spacial scales used in the Lake Ontario analysis.
13
-------
As an illustration of the application of WASP to lake models, the
earlier work in Lake Ontario can be considered. This work utilized a whole
lake analysis, vertically segmented at the average thermocline depth. This
model, designated the LAKE1 model to indicate its vertical one-dimensional
nature, extended the earlier work on eutrophication of estuaries to a large
lake such as Lake Ontario.
Figure 5 shows the major physical features of the LAKE1 model. The
principal physical features are a) horizontal transport, b) time variable
vertical exchange between the epilimnion and the hypolimnion to permit
conditions of a vertically mixed lake or a vertically stratified lake and c)
vertical setting of phytoplankton and other particulate nutrient forms. The
systems diagram for the LAKE1 model is shown in Figure 6 where eight vari-
ables were implemented: the nitrogen and phosphorus variables, chlorophyll
and the two zooplankton groups. This kinetic structure was later expanded
in the Lake Erie work (6), and also subsequently incorporated in further
work on Lake Ontario. The updated Lake Erie and Lake Ontario kinetic struc-
ture is shown in Figure 7. This expansion included the addition of silica
as a state variable (in two forms) and the division of the phytoplankton
variable into "Diatom" and "Non-Diatom". In each case, the coupling between
the variables is both linear and nonlinear and the parameters may be speci-
fied as time dependent straight line functions. A typical calibration using
the updated kinetics, LAKE1A, is shown in Figure 8.
The WASP structure proved particularly useful in all of the Great Lakes
work by permitting the ready expansion of models in the spatial dimension
as well as expansion of the kinetic structures. In each application for
each lake system, the principal effort was in deciding on the spatial con-
figuration to be used, together with the appropriate kinetics. Once the
decision was made, WASP permitted the formulation of a "new" lake model by
requiring only a minimum amount of effort to write the specific kinetics for
insertion into the program and in the preparation of any new input on trans-
port and dispersion.
WASP has also been used ab part of studies on Lake Ontario in a long
term, year to year analysis of lake behavior (8). Such an analysis provides
a basis for estimating lake responses under different external nutrient
loadings. In this mode, the seasonal dynamics were computed for each year
and the output after a two-year run became the initial conditions for a
subsequent two-year run. Output from each two-year run is saved on a per-
manent file for computations using the Model Verification Program, MVP.
In the multi-year analysis, various kinetic schemes were employed to
arrive at the best overall model that represented the long-term data. The
WASP framework proved particularly useful in this type of study. Figure 9
shows the comparison between the observed chlorophyll data in the epilim-
nion for 1967-76 and the calculated values using the LAKE1A kinetics. It
should be noted that the results shown in Figure 9 were calculated by
setting initial conditions in 1966 and continuing the calculation to 1976
without any resetting of conditions. Actual year to year temperature data
and loads were used.
14
-------
CMVI/WMHfMTAL INPUTS
NIAGARA MlVf M
TftWUTAMIU MUNMWAl
INOUCTItlAl WAfTE*
not/Ml
p»l«MTIIIl
(UOMTb
L*YtTM
UONTIXTINCnON.
PAMMfTIW
fP/LHINION
VERTICAL EXCHANGE » THAWfOWT 17 M
HYPOUMHON
i i i i 1
3) BENTHOS
II" II Illl Illl Mill III I Illl Illl Mil Illl Ml"
Figure 5. Major physical features included in the LAKE1 model
15
-------
5-.
So
cu
T3
O
s-
cn
03
r
-a
co
QJ
co
cu
CD
16
-------
fHOSfMOtus srecics
SILICA SftClfS
Stdinwntflion
UNAVAILABLE
PHOSPHORUS
AVAILABLE
PHOSPHORUS
HITHOOEH SffCICS
Figure 7. Systems diagram - updated lake kinetics.
17
-------
M
W
e
o
K
HJ
Z
1
K
|_0
^
t i i
5
6
1 1
o
o
o
o
oo
d
o
*.
d
CO
<1J
NOlXNVtdOOZ
NOlXNVTdOOZ
(I/Z0!S 6uu)
aivonis
I I I I
re
o
ai
T3
oo
01
o
o
o
o
ID
O
o
o
..SW3H10..
18
-------
6
-3
l--^x-i
-Q.
«* 00
. 0>
6
<
(O
o
10
oi
o
CD
KHl
B>
if)
cc -
Q.
of
-a.
of
-a.
(M
a:
-a.
8
6
Q.
en
O
0)
s_
^
en
(I/6T/) ,0,
19
-------
As noted, the output was saved and the MVP program was accessed to pro-
vide estimates of the verification status of the model. For chlorophyll,
some results of the verification analyses are shown in Figure 10 where re-
gression analyses and the relative error are shown across all years for the
two kinetic systems. Statistical comparison results such as these provide
the analyst with quantitative measures of model status. Figure 11 shows the
median relative error for 1967-1976 for each variable and across all vari-
ables. For the latter, the results show that the inclusion of the more com-
plicated LAKE1A kinetics did improve overall model status by about one-
third. The overall median error of 22% represents a useful measure when
discussing overall model credibility. Finally, MVP also produces verifica-
tion scores (using a Student's "t" test) and some results across all years
and variables are shown in Figure 12. At the given standard errors of the
mean, the verification score for LAKE1A kinetics was 70%. In other words,
for the months and variables where a comparison could be made between ob-
served and computed values, there was no difference between observed and
computed means 70% of the time at a 90% confidence interval. If the esti-
mated standard errors are in question, then a new score is easily computed
as shown in Figure 12 for values of 1/2 and 1-1/2 of-the given standard
errors. These results illustrate the utility of the WASP-MVP package - in
this case for a long term, multi-variable model analysis.
WASP has also demonstrated its flexibility as a modeling tool in in-
vestigating the spread and impact of hazardous materials in receiving
waters. A particular application was to determine the receiving water pH
response to a spill of a strong acid (9). The model's kinetic framework
was based upon the carbon dioxide-bicarbonate-carbonate equilibria (since
this buffer system is the predominant buffering system in natural waters),
and the exchange of carbon dioxide with the atmosphere under supersaturated
conditions. As a check on the model kinetics (carbonate-bicarbonate buffer-
ing), a simulation of the laboratory titration of a sample volume of river
water with a strong acid was performed. Figure 13 shows a good comparison
between lab data and model results. Figures 14 and 15 show receiving water
response to a 72 minute spill of a strong acid. Figure 14 also indicates
the type of segment grid which one may use to model spills of hazardous
materials, a fine mesh at the immediate point of discharge and an expanding
mesh as one is further removed from the point of discharge.
WASP was also employed in a recent study (10) for the State of New York,
Department of Environmental Conservation, concerning the need for remedial
action and the impact of such action on achieving polychlorinated biphenyl
(PCB) environmental objectives in the Hudson River ecosystem. A modeling
framework was developed which defined the major interactions that lead to
the PCB distribution in the biotic and abiotic sectors of a given aquatic
environment (Figure 16). A simplifying assumption was made which permitted
the decoupling of the biotic and abiotic sectors (i.e., if the food chain,
above the phytoplankton, is viewed as a whole, then its uptake and excretion
produce sinks and sources of PCB in the abiotic sector which are negligible
when compared to those within the abiotic sector alone). WASP was used then
to compute PCB distributions in the abiotic sector which were dependent upon
advective and dispersive transport, sedimentation and resuspension, sediment
release, evaporation, absorption and de-absorption kinetics and external
20
-------
10.
X 6
O.
O
X
o
o
Ul
oe
UJ
-------
CO
OJ
ns
r
S-
i
o>
fO
4->
CO
03
O)
>
Ol
CO
O)
o
o
E
CO
0)
co
S-
o
S-
s-
O)
Ol
n3
'QJ
O)
O)
s_
k|1VO-S80|
3AllV13d NVIQ3N
22
-------
< V)
I O
si
v>
_ o
KZ
o v
0 c
OJ
T3
Z3
4->
00
CXI
CD
S-
(03IJIW3A SH1NOW-3T8VIWVA JO %) 3HODS NOI1VDIJIM3A
23
-------
ID 15
ml 0.1026 N H2S04 ADDED
Zfl
2.5
Figure 13. Acid spill problem lab data versus model output.
24
-------
LATERAL DISPERSION
EQUALS LONGITUDINAL DISPERSION
AT END OF 72 MINUTE SPILL
ALKALINITY=I23 mg CoCOs/l
DEPTH = 5 FEET
FLOW =
72CFS
10 MINUTES AFTER SPILL CEASES
1
25 MINUTES AFTER SPILL CEASES
MOTC
CONTOURS IN pH UNITS
Figure 14. Spatial river pH response to a 72 minute spill of
96 percent sulfuric acid at 0.8 gprri.
25
-------
LATERAL DISPERSION
EQUALS LONGITUDINAL DISPERSION
DISPOSAL SEGMENT
0.8 GPM 96 %
Y/////7if7/77/77.
'//Y////A I
10
I
Q.
Q.
20 40 6O
TIME (MINUTES)
80
100
8O FEET DOWNRIVER OF SPILL LOCATION
I
20 40 60
TIME (MINUTES)
80
100
I6O FEET DOWNRIVER OF SPILL LOCATION
I
ZO 40 60
TIME (MINUTES)
80
100
Figure 15. Temporal river pH response to a 72 minute spill of
96 percent sulfuric acid at 0.8 gpm.
26
-------
r
c
to
_c
o
-o
o
o
O
_c
CD
2
o
rs
T3
CQ
Q-
QJ
3
P
(O
O)
O
<->
to
4-
O
O
01
_c
o
CD
27
-------
waste sources. Some model results for the upper Hudson river are presented
in Figure 17.
Although not initially thought of as a tool for modeling biological
treatment processes, WASP is capable of being programmed so as to provide
this capability. WASP has been applied to a number of research efforts
aimed at developing useful kinetic models for the anaerobic filter,
rotating biological contactor and and the pure oxygen treatment system to
name a few (11,12).
Even though WASP is essentially a stand-alone program, i.e., not exe-
cuted in conjunction with another program (that might produce hydrodynamics
or runoff waste loads), it can be user-programmed to interface with other
computer programs if desired. In the recent New York City 208 Study (13)
WASP was modified to accept time-variable non-point runoff flows and waste
loads from a landslide rainfall runoff model. Figure 18 shows the general
system logic and interfacing of the programs.
PROGRAM OVERVIEW
Before going into a detailed presentation of the WASP program and the
theories upon which it is founded and numerically implemented, a brief dis-
cussion will be presented of the basic philosophy and assumptions underlying
the program.
The key principle upon which the model equations of WASP are founded is
the principle of conservation of mass. This principle simply states that
the mass of each water quality constituent being investigated must be
accounted for in one way or another. WASP, conserving mass both in time
and space, accounts for and traces the water quality constituents from their
point of spatial and temporal input to their final points of export.
In order to perform the spatial and temporal mass balance computations
the user must supply WASP with input data defining the model segmentation,
advective and dispersive transport fields, boundary conditions, forcing
functions (waste loads), segment parameters, kinetic constants, time vari-
able kinetic functions, and initial conditions for the state variables
(water quality constituents). WASP utilizes this input data together with
the user supplied kinetic subroutine to construct the mass balance equa-
tions, which are then numerically integrated in time. At a user specified
time interval (print interval) WASP saves the current values of the state
variables, and other user selected variables of interest, and stores them on
auxiliary storage disk files for subsequent retrieval by the WASP graphics
subroutine and the MVP program. A simplified program flow chart is pre-
sented in Figure 19.
MVP may be executed at the user's discretion. If the user should chose
to perform MVP analyses he will need to supply field data for comparison to
the theoretical computations generated by WASP. MVP uses three statistical
tests for scoring the model verification. The scores are determined using a
Student's "t" test on a comparison of the means of the theoretical and ob-
28
-------
=8
°
I »
u »
5 5
X O
0»
v.
21-
o
a.
tu
_i
Si
O
M
o
a.
o
CO
o
to
o
CM
Oco
83d
(I/GUM sanos a3QN3dsns
0)
o
o
to
O)
Q.
to
"O
(O
CO
o
Q_
n3
o
ai
s-
O
to
o
-------
INPUT
OUTPUT
PERIOD
OF
INTEREST
^
r
SYNTHETIC RAINFALL
RECORD
(HOURLY)
t
QUALITY
MATRIX
LAND AREA
CHARACTERISTICS
DRAMAGE AREAS
LAND USE
RAIN GAGES
RUNOFF DISTRIBUTION
QUALITY
C$0 CAPTURE
DIURNAL EXCESS
STP CAPACITY
TIME
AVERAGING
INTERVAL
\
\
*
RAINFALL
RUNOFF
MODELING
PROGRAM
1
(RRMPI )
*.
RUNOFF LOADS
MOOCL
MBMCMT
TIMC
RUNOFF FLOWS
MOOCL
*E*MCNT
TIMC
FLOW ROUTE
MATRIX FILE
k
RUNOFF FLOWS
(FROM RRMPI )
BOUNDARY FLOWS
j BOUNDARY
TIMC
UNIT FLOWS
( BY MODEL
INTERFACE )
RAINFALL
RUNOFF
MODELING
PROGRAM
2
(RRMP2 )
TIME VARIABLE
FLOWS
MOO EL
INTERFACE
TIME
BASE FLOW
STEADY STATE
MODEL CHARACTERISTICS
GEOMETRY, INTERFACES,
DISPERSION, KINETICS,
CONTINUOUS LOADS, ETC
TIME VARIABLE
STORM LOADS
(FROM RRMPI)
TIME VARIABLE
FLOWS
( FROM RRMP2 )
TIME VARIABLE
BOUNDARY
CONDITIONS
WATER
QUALITY
MODEL
(WASP)
7V*f VAKIASLE
WATER QUALITY
MOOCL
MMCNT
TIME
Figure 18. Land side simulator - WASP systems diagram.
30
-------
-------
served data, a linear regression of theoretical and observed data, and a
comparison of differences between theoretical and observed data. It should
be noted, however, that even without extreme field data the MVP may prove
useful in aggregating computed output according to any spatial or temporal
averaging scheme.
Both WASP and MVP are written in FORTRAN IV, using a modular, subroutine
orientated approach. This has permitted both programs to be available for a
number of different computers with different core capacities. Table 2 con-
tains a list of the computers on which WASP and MVP have been implemented.
TABLE 2. WASP AVAILABILITY
Computer
DEC POP 11/45 (EPA Grosse He)
DEC POP 11/70 (Manhattan College)
DEC POP 11/70 (EPA-Athens)
IBM 370/168 (EPA-COMNET)1
CDC 6600 (NYU)1'2
Core Capacity
Words
32K
32K
32K
256K
TOOK
Program
Systems
16
12
16
19
23
Configuration
Segments
20
40
29
120
20
^Completely core-resident (all other versions require program overlays).
^Modified version of WASP - input structure not compatible with the IBM,
DSC, and DEC versions.
REFERENCES
1. Hydroscience. 1972. Mathematical Model of Phytoplankton Dynamics in
the Sacramento-San Joaquin Bay Delta, Preliminary report. Prepared for
the California Department of Water Resources.
2. Hydroscience. 1974. Western Delta and Suisun Bay Phytoplankton
Model - Verification and Projections. Prepared for the California
Department of Water Resources.
3. Thomann, R.V., .et ^1_. 1975. Mathematical Modeling of Phytoplankton in
Lake Ontario, 1. Model Development and Verification. EPA-660/3-75-005,
NERC, ORD, EPA, Corvallis, Oregon. 177 pp.
4. Thomann, R.V., .et _al_. 1975. Mathematical Modeling of Phytoplankton in
Lake Ontario, 2. Simulations Using LAKE 1 Model. EPA-660/3-75-005,
NERC, ORD, EPA, Corvallis, Oregon.
5. DiToro, D.M., et_ _§_! Lake Huron report in review.
32
-------
6. DiToro, D.M., et_ _al_. Lake Erie report in review.
7. Thomann, R.V., jjt _al_. 1978. Verification Analysis of Lake Ontario and
Rochester Embayment Three-Dimensional Eutrophication Models. In
review.
8. Thomann, R.V., et _al_. 1978. A Ten Year Modeling Analysis of Phyto-
plankton in Lake Ontario. In preparation.
9. Nusser, J.A., £t _al_. 1978. Mathematical Modeling for Impact and
Control of Hazardous Spills. Proceedings of the 1978 National Con-
ference on Control of Hazardous Material Spills, sponsored by the U.S.
EPA and the U.S. Coast Guard, Miami Beach, Florida, April 1978.
10. Hydroscience. 1978. Estimation of PCB Reduction by Remedial Action
on the Hudson River. Prepared for the State of New York, Department of
Environmental Conservation.
11. Mueller, J.A., et_ _al_. 1977. Application of Carbonate Equilibrium to
High Purity Oxygen and Anaerobic Filter Systems. Presented at the
173rd National Meeting of t^e American Chemical Society, New Orleans,
Louisiana.
12. Mueller, J.A., .e_t _al_. 1978. Nitrification in Rotating Biological Con-
tactors. Presented at the 51st Annual Conference of the Water Pollu-
tion Control Federation, Anaheim, California.
13. Hydroscience. 1978. Task Report 334, Time Variable Modeling, sub-
mitted to Hazen and Sawyer, Engineers for NYC 208.
33
-------
SECTION 4
WASP THEORY
MASS BALANCE EQUATIONS
The basic concept in writing a mass balance equation for a body of water
is to account for all of the material entering and leaving the water body
via direct addition of material (runoff and loads), via advective and dis-
persive transport mechanisms, and via physical, chemical, and biological
transformations. In order to formulate a mass balance equation for a speci-
fic water quality variable, consider a coordinate system as shown in Figure
20a, where the x- and y-coordinates are in the horizontal plane and the z-
coordinate is in the vertical plane. Several authors (1,2,3) have stated
the proper form of the mass balance equation around an infinitesimally small
fluid
.2-
A) COORDINATE SYSTEM B) INFINITESIMAL
VOLUME
Figure 20. Coordinate systems,
volume (Figure 20b) to be
3"t BX x 3x 3y y ay 8z z 9x 9x 5
- ^ Uyc - -2 Uzc + S (x.y.z.t) (2.1)
34
-------
where: c = concentration of the water quality variable [M/L ]
t = time [T]
S = vector of all other sources and sinks of the water quality
variable
2
E = diagonal matrix of dispersive coefficients [L /T]
U = vector of velocities [L/T]
The dispersion coefficient represents the overall phenomenon of mixing
due to the temporal variation of tidal velocity, the lateral and vertical
gradients in velocity and the density differences within the water body.
For the sake of clarity and brevity the derivation of finite-difference
form of the mass balance equation, to be presented in the following section,
will be for a one-dimensional rectangular estuary. Under the assumptions of
vertical and lateral homogenuity, and permitting variation, or "a gradient",
along only the length of the estuary, the mass balance equation for the one-
dimensional case is
3(E Ac)
3U Ac
FINITE DIFFERENCES
The fundamental method of solution employed in WASP is the use of
finite-difference approximations to the derivatives of Equations (2.1) and
(2.2). In order to keep this section from becoming overly long, the devel-
opment of the numerical procedures will of necessity be abbreviated, with
the occasional use of "it can be shown". However, if the reader is
interested in a mo^e detailed explanation of finite-difference methods, the
author recommends an excellent presentation by Smith (4).
Taylor's series expansion theorem states that for a function u, where
its derivatives are single-valued, finite, and continuous functions of x,
then
u(x+h) = u(x) + hu'(x)
h3u'
(x)
(2.3)
u(x-h) = u(x) - hu'(x) + hu"(x) +
h3u"'
(x) +
(2.4)
Assuming that terms containing the third and higher powers of h are negli-
gible in comparison to the lower powers of h, then Equations (2.3) and (2.4)
may be subtracted and added together and rearranged to produce Equations
(2.5) and (2.6) respectively,
35
-------
u'(x) =
(u(x+h) - u(x-h)}
(2.5)
u"(c) =
4 {u(x+h) - 2u(x)
\\
u(x-h)}
(2.6)
with an error terms of order h2. Referring to Figure 21, Equation (2.5)
can be seen as an approximation to the slope of the tangent centered at P
formed by chord AB, and is known as the central-difference approximatior
M(M+X)
x-h
x+h
Figure 21. Finite-difference approximations.
The slope of the tangent at P may also be approximated by the slope of the
chord PB giving the forward difference formula,
u'(x) - \ (u(x+h) - u(x)}
(2.7)
or the slope of the chord AP, giving the backward-differences formula,
u'(x) - \ {u(x) - u(x-h)}
(2.8)
Both Equations (2.7) and (2.8) may be obtained from (2.3) and (2.4) re-
spectively, by assuming the second and higher order powers of h are negli-
gible. The error term for both the forward-difference and the backward-
difference is in order h. To ensure positive stable solutions WASP employs
a backward-difference approximation in the spatial plane (5). For pro-
gramming simplicity WASP employs a forward-difference approximation for the
temporal plane.
36
-------
Formulate a grid, as shown in Figure 22, that subdivides the x-t plane
into sets of equal rectangles with sides equal to Ax and At, and let the
coordinates (x,t) of a representative of point P be
x = iAx and t = nAt,
where i and n are integers. Use the following notation to denote the value
of u at P
Up = u(iAx, nAt) = u.
Jk
P(ih.jk)
1,1
i+I.J
T
l.j-l
Ih
Figure 22. Finite-difference grid.
Then by Equation (2.6)
n
'
3X
Vl,n-2ui,n-Vl,n
Ax
(2.9)
Similarly, for the backward-difference approximation Equation (2.8)
ui,n~ui-1,n
Ax
(2.10)
and for a forward-difference approximation time equation
37
-------
Restating Equation (2.2)
1 9Qxc
O- _ _ * ' * 4. C fO ION
§1 'A 2 -- A is (2'12)
dX
where QY = U A, and
* x
A = the cross-sectional area
and using Equations (2.9) and (2.10) as finite-difference approximations for
the first and second terms of the right side of the equation respectively,
one may develop a finite difference approximation at the mesh point (i,n)
as follows:
dci i
ar = ' A - A* -
, (EAc).+1 n-2(EAc). +(EAc). ,
+ ' - l±ii!3 , LiD - !iiiH
+ s (2.13)
Letting V = AAx and re-arranging terms Equation (2.13) may be expressed
Vi 3T = ^c)i-l,n-^c)i,n
(EAc). , -2(EAc), +(EAc). ,
+ - HUD - LiH - llliJ] + S (2.14)
Dividing the water body into completely mixed finite segments as pictured in
Figure 23, and recalling that S represents the source-sink terms both for
the direct addition of the water quality constituent (point source dis-
charge, distributed runoff, etc.) and addition-removal through reactive pro-
cesses, one may see that Equation (2.14) represents a mass balance for seg-
ment i, and may be restated for a fixed point in time, as follows:
vi '
38
-------
1,1+1
Figure 23. Completely mixed finite segments,
39
-------
where o
c-j = segment concentration, M/L
3
V-j = segment volume, L
3
Qi-l,i = net advective flow, from segment i-1 to segment i, L /T
EI_I i = dispersion coefficient for the i-1, i interface, L /T
' ' 3 ' f\
Ai-l,i = cross-sectional area of the i-l,i interface, L
L-j-l,i = characteristic length or mixing length of segment i-1 and
i, L.
-
Li-l,i I
W-j = point, or distributed sources-sinks of the water quality
constituent, M/T
K-j = kinetic or reactive process rate, 1/T.
It can be shown (5) that Equation (2.15) can be extended to develop the
multi-dimensional, multi-constituent form of the equation as stated by Equa
tion (2.16)
dc m E A
m m
+ W. + K.V.c. + K..c, (2.16)
i ill i 11
where ,
K-J = the cross-coupling reaction term
c-j = the concentration of state variable 1 in segment i
It should be noted that for the one-dimensinal and multi-dimensional
cases, the advective and dispersive fields are assumed known. It should
also be noted that for irregularly shaped water bodies (such as embayments)
that V = AAx, where A is some representative cross-sectional area, as shown
in Figure 24a. Also note that for multi-dimensional water bodies a segment
may have different characteristic lengths for adjoining segments, see Figure
24b.
40
-------
V = AK
A) CROSS-SECTIONAL AREAS
Lj K
B) CHARACTERISTIC LENGTHS
Figure 24. Cross-sectional areas and characteristic lengths.
If we divide Equation (2.16) through by V-j and permit the new right-
hand side to equal c-jn, and then use the forward-difference approximation
for the time derivative, we can form the following expression
(2.17)
which states that the concentration at time n+1 is equal to the concentra-
tion at time n plus the derivative evaluated in time n times the time step,
At.
STABILITY AND NUMERICAL ERROR
Classical stability analysis (Smith) requires At/Ax^ <_ 1/2 to
guarantee stable solutions for (2.17). However, in the classical analysis
At and Ax are normalized variables and it is difficult to extend these vari-
ables to practicable applications. It has been shown (6) that for computa-
tional stability a necessary condition is that
At < Win (;
Vi
(2.18)
where
= exchange coefficients =
E . . A.
= 1J
(2.19)
However, since K-j is, for many applications, non-linear and time-dependent
itself, and therefore may be difficult to evaluate, the following criteria
may be used for choosing the integration step-size
41
-------
At < Min ( ; 1 (2.20)
providing that
Max (K.V.) « zQ.. + ZR...
j J j J
It should be recognized that the use of the completely mixed finite seg
ment approximation in conjunction with the backward difference spatial ap-
proximation, introduces a numerical error (sometimes referred to as numeri-
cal or pseudo-dispersion in the literature) into the model (7). The extent
of this effect is given by
Enum
where Enum is the numerical error expressed as a pseudo-dispersion coeffi-
cient. For some applications, especially intra-tidal models, where u may
be large, this error term may lead to highly distorted results (i.e., arti-
ficially spread spatial profiles). However, for many applications,
especially in well mixed estuaries, where the time scale of importance is on
the order of days to seasons rather than hours, and where u is the net ad-
vective freshwater velocity (usually small), the effect of the numerical
error, Enum, is generally not significant. It should be noted that one
cannot reduce or eliminate the numerical or pseudo-dispersion coefficient,
Enum, by adjusting E, the true dispersion coefficient. The only way to
reduce the effect of Enum in WASP is to reduce the segment length in the
direction of u. However, one must also be aware that arbitrary reduction in
the segment length effects the integration step size since
iQ + zR
or
Z +
and a reduction of Ax requires a like reduction of At. Choice of proper
spatial and temporal grid sizes is still somewhat of an art and considera
tion of computer core size and execution speed, the nature of the problem
being analyzed, and the degree of accuracy required, all influence the
choice.
42
-------
COMPUTATIONAL PROCEDURE
In order to integrate the aforementioned finite difference equations
WASP uses a second order Runge-Kutta or predictor-corrector method (8). The
operational procedure for the Runge-Kutta method may be best explained
schematically via Figure 25. First, using (2.17), cn+i/2 is computed for
the half-step inteval At/2 (predictor). Using cn+i/2, tn+]/2 (i.e.,
tn + At/2), the derivative cn+]/2 is evaluated, and this derivative is
used as the average derivative for proceeding over the whole interval At
(corrector). Experience has shown that for the general class of applica-
tions that WASP is used for, the second order Runge-Kutta method yields re-
sults comparable to the more commonly used fourth order Runge-Kutta method
at an execution time savings of fifty-percent.
WASP uses a slightly modified version of the Runge-Kutta method in that
WASP, in normal operations, prohibits any segment concentration from going
negative and causing either numerical instability or numerical oscillations
of the solution. In some applications, particularly in the nutrient limited
growth phase of eutrophication modeling, it is possible for a particular
segment derivative-timestep combination to cause a negative concentration.
This negative concentration, if permitted to occur, might degenerate the
true solution by causing either instability or oscillations. Rather than
permit this to happen, WASP, upon detection of a derivative-timestep com-
bination which would cause a negative segment concentration, maintains a
positive segment concentration by setting the segment concentration pro-
jected for timestep n+1 to half the concentration that was present at time-
step n. It should be noted that this procedure does not maintain a proper
mass-balance, i.e., mass is not conserved. However, experience has shown
this procedure to be acceptable within reason. WASP informs the user if
and when a half-concentration procedure was performed, and the user can
monitor the frequency of occurrence of the procedure and determine if the
simulation must be rerun at a smaller timestep.
REFERENCES
1. Pritchard, D.W. 1958. The Equations of Mass Continuity and Salt
Continuity in Estuaries. Journal of Marine Research, Vol. 17 (Nov.
1958), pp. 412-423.
2. Daily, J.W. and D.R.F. Harleman. 1966. Fluid Dynamics. Addison-
Wesley Publishing Company, Reading, Mass., Chapter 16.
3. O'Connor, D.J. and Thomann, R.V. 1971. Estuarine Modeling: An
Assessment, Chapter III EPA, 16070 DZV 02/71.
4. Smith, 6.D. 1965. Numerical Solution of Partial Differential Equa-
tions, Oxford University Press, Ely House, London, Great Britian.
5. Thomann, R.V. 1972. Systems Analysis and Water Quality Management,
Environmental Science Services Division, Environmental Research and
Applications, Inc., New York, New York.
43
-------
PREDICTOR STEP
»Cl4l/2 111* I
-------
6. Kent, R. 1960. Diffusion in a Sectionally Homogenous Estuary. Trans
ASCE, Jour. SED, Vol. 86, No SA2, Mar. 1960, pp. 15-47.
7. Bella, D.A. and W.E. Dobbins. 1968. Difference Modeling of Stream
Pollution. Proc. ASCE, 94, No. SA5, Oct. 1968, pp. 995-1016.
8. Carnahan, B., ejt _al_. 1969. Applied Numerical Methods, John Wiley and
Sons, Inc., New York, New York, pp. 361-366.
45
-------
SECTION 5
WASP PROGRAM LOGIC
INTRODUCTION
WASP was designed and written so as to permit its application to as wide
a variety of water quality modeling problems as possible without a user
being required to make extensive program revisions for each new application.
In order to achieve this goal four objectives were set down and met:
1. FORTRAN IV would be chosen as the programming language due
to its universality.
2. WASP would be written using a modular subroutine orientated
approach for both program clarity and the ability of most
computer operating systems to accomodate subroutine overlay
structures where core requirements are restrictive.
3. WASP would permit the user a great deal of flexibility in
structuring the physical setup of his model via a number
of user selectable input options.
4. WASP would, by splitting the kinetic portion of the mass
balance Equation (2.16) away from the remaining terms (the
advective and dispersive transport and the source/sink
terms), require the user to develop only a FORTRAN sub-
routine that would describe the kinetic interaction of the
state or water quality variables. This would remove from
the user the need to know how WASP handles the remaining
portion of the mass balance equation on a FORTRAN coding
level, although it is important to understand the mass
balance equation on a general level.
WASP is comprised of a mainline program and twenty-eight support sub-
routines (forty subroutines for the DEC POP 11/45 version, which includes
graphics output) and one user written kinetic subroutine. Before going into
a detailed explanation of the program logic and the procedure to follow to
develop a new kinetic subroutine some background explanations must be given.
46
-------
TIME VARIABLE FUNCTIONS
WASP permits the user to specify time variable input for any of the fol-
lowing parameters: exchange coefficients, advective flows, boundary condi-
tions, forcing functions (waste loads) and miscellaneous functions which
might be required by a user's kinetic framework. The user specifies the
time variable input data for any of the previously mentioned parameters as a
series of time and value combinations which WASP uses as a piecewise linear
function of time. Suppose Figure 26a presents the observed measurements of
daily solar radiation incident upon a body of water.
Figure 26b shows how the user might approximate the solar radiation with
a piecewise function of time. Table 3 presents the input data the user
would supply. Core considerations, especially for the small minicomputer,
require WASP to store the entire piecewise linear function on auxiliary disk
storage, and maintain core resident only the appropriate information needed
to evaluate the function during the current time step. How this is accom-
plished will be discussed later in the section entitled "Writing a WASPB
Subroutine".
TABLE 3. PIECEWISE LINEAR APPROXIMATION FOR FIGURE 26b
lav
(ly/day)
187.5
487.5
750.0
487.5
150.0
Time
(days)
0.
75.
165.
255.
345.
, Iav ,
(ly/day)
262.5
640.0
730.0
300.9
170.0
Time
(days)
15.
105.
195.
285.
365.
lav
(ly/day)
337.5
712.5
637.5
190.0
Time
(days)
45.
135.
225.
315.
KINETIC DATA - CONSTANTS, SEGMENT PARAMETERS AND PIECEWISE LINEAR FUNCTIONS
OF TIME
Kinetic data - constants, segment parameters, and miscellaneous func-
tions of time - are read as input data by WASP for use in the user supplied
kinetic subroutine, WASPB. The actual choice, of what the constants, seg-
ment parameters, and piecewise linear functions of time are to be, is deter-
mined by the modeler and the systems analyst as they develop a new kinetic
subroutine. The selection and application of kinetic data will be discussed
in greater detail in "Writing a WASPB Subroutine", but their basic concept
is presented here for background reference. The BOD oxidation rate at 20°C,
the saturated growth rate of phytoplankton at 20°C, ammonia to nitrate
nitrification rate at 20°C, temperature correction factor for reaeration,
and saturated light intensity may all be thought of as constants. Para-
meters are segment dependent and may include such factors as segment depths,
reaeration coefficients, and water temperatures. Kinetic piecewise linear
functions might include the daily incident solar radiation for a year, or
the fraction of daylight hours over a year. Due to core requirements it was
decided not to permit segment dependent piecewise linear functions. How-
47
-------
CO
c
O
+->
O
S-
rcs
O)
O)
co
i
tu
(J
Ol
CM
OJ
s-
48
-------
ever, a user may use segment parameters and piecewise linear functions to-
gether to achieve the same effect. For example, a modeler may aggregate ob-
served variations in water temperature over a year for different areas of
the water body into one normalized composite piecewise linear function. The
modeler would then use the segment parameters for input of the maximum
yearly water temperature for the individual segments and a piecewise linear
function for input of the normalized aggregate variation during the year.
UNITS
WASP has been programmed using the following internal units conventions.
The units of concentration of the state variables (or water quality vari-
ables) are assumed to be in mg/1 or parts per million parts (ppm). The seg-
ment volumes are read in as million cubic feet (MCF). Time is in days. Ad-
vective flows are nominally read with units of cubic feet per second (cfs)
and internal.ly converted to units of million cubic feet per day (MCF/day)
using the following conversion factors:
or
= 0.0864 x Q cfs
Therefore the term Q-jjC-j has units MCF-mg/1/day. The user may read the
exchange coefficients as exchange coefficients with units of MCF/day, or as
dispersion coefficients, cross-sectional areas, and characteristic lengths
and use Equation (2.20) to compute exchange coefficients. If the latter op
tion is chosen the nominal input units for the dispersion coefficients,
cross-sectional areas, and characteristic lengths are square miles per day,
square feet, and feet, respectively. WASP uses the following factors to
convert the dispersion-exchange coefficient to MFC/day:
x A [ft2] 5
or
.
. . ay -- - x ,,_ x -3-3
,, [MCF/day] = 27.8764 x
E..[mi2/day] x A [ft2]
,, . - - ,
J
49
-------
Similarly the term Rij (cj-Ci) also has units MCF-mg/1/day. Forcing
functions, or waste source loads, are nominally read in units of pounds/day
(lb/da.y), and internally converted using the following convention:
Wn. [MCF-mg/1/day] = MI tlb/day] x 453.59 ^
v in3 "US v Lft3 v i MCF
gm 28.32 liters ~~"~
or
W. [MCF-mg/1/day] = -^ x W1 [Ib/dayl.
dc.
Looking at the term V -TT-, we can readily see that it also has units MCF-
mg/l/day. In order to have a consistent set of units for Equation (2.15),
the kinetic term must also have units MCF-mg/1/day. Table 4 presents some
TABLE 4. KINETIC TERM UNITS
Reaction Order Kinetic Term Units
Zero Kn-Vi Ki [ = ] mg/l/day
First Kicivi Ki t = l Vday
C.
Michaelis-Menton K^ c ]_ K V Ki [=] mg/l/day
Coupled (derivative term aK^-C,-1 C,-mV.j K,-[=] mg of cVmg of C1
for state variable 1) a [=] mg of C'/mg of Cm
of the possible forms of the kinetic term which meet this requirement. It
should be noted that a user may use units other than those nominally ex-
pected by WASP as long as he is careful and consistent in their use. Sup-
pose for example a modeler is investigating coliform bacteria as a state
variable. Normally coliform bacteria are measured as most probable number
per 100 ml of sample or MPN/100 ml. If the user wishes WASP to interpret
the pseudo-concentration of bacteria as MPN/100 ml, he should enter his
initial conditions and boundary conditions as MPN/100 ml, and use the fol-
lowing convention to enter waste loads as pseudo "Ibs/day"
W. ["Ibs/day"] = Qv [cfs] x 0.0864 x 62.43 x Cco]i [MPN/100 ml]
50
-------
or
W. ["Ibs/day"] = Qy [cfs] x 5.394 x Cco1i [MPN/100 ml]
or
Wi ["Ibs/day"] = Qy [MGD] x 8.34 x C [MPN/100 ml]
where
Qv is the waste discharge volume flow rate
Cco]i is the coliform bacteria "concentration".
WASP MAINLINE AND SUBROUTINE OVERVIEW
The WASP mainline is really just a program control module. As such it
performs no computations but does assign the logical units (disk) for
temporary scratch files and the dump-save files generated in the user's
kinetic subroutine WASP and does control the calling sequence of the WASP
subroutine.
In the following subroutine descriptions the references, within paren-
theses, to Card Groups are with respect to the card by card description of
the WASP input data to be found in Section 6, WASP Input Structure. Also,
in the following subroutine description occasional reference to internal
WASP program variables will be made. A complete definition of these vari-
ables is presented in the section entitled, "WASP Common" immediately fol-
lowing this section.
WASP1
WASP1 reads the model identification and system by-pass options for the
user's model (Card Group A). WASP1 also performs some variable initializa-
tion. Two variables of interest set by WASP1 are NBCPSY and NWKPSY. NBCPSY
and NWKPSY are used to indicate the maximum number of boundary conditions
and forcing functions that WASP is dimensioned for, respectively.
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).
51
-------
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 cfs to MCF/day. WASP4 also reads the flow by-pass options
for each system.
WASPS
WASPS, 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
WASPS, 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).
WASP 7
WASP7 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 i^1 update.
slope =
intercept = f(t).+,
next update time = t.+.
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.
52
-------
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.
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 - WA12A
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 pur-
poses will be noted as K-jVjC.,-"1. WAS12 then goes through the following
steps:
a. Using the IQ and JQ vectors as drivers, WAS12 computes
(ViC.m) = (V.C."1) + W-.C.m - ZQ.kC.m
or (V.C."1) = K.V.C.1" + ZQj1Cjm - SQ^"1
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 vec-
tor).
b. Using the IR and JR vectors as drivers, WAS 12 computes
(V.C^) = (V.C."1) + ZR... (C.m - C.m)
or (V.C.m) = (K.V.C/" + EQ.-.C.j"1 - ^C^} + ^ (C - C.m)
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).
53
-------
c. Using the IWK vector as a driver, WAS12 computes
W
11
or (V.C) =
- C/V
d.
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).
WAS12 completes the evaluation of the derivative, C-j"1 by dividing
through by the volume and multiplying by the time scale factor
(SCALT)
C.m = SCALT x (V1C.m)/V1
Note: C-j now has units M/L^/T - nominally mg/l/day.
WAS 13
WAS13 is used to print the intermediate system-segment pairs during the
simulation.
WAS 14
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 his simulation at any time other
than time equal to zero. The following real case history demonstrates the
use of TINIT. For the analysis performed in the Western Delta - Suisun Bay
area of San Francisco Bay, time zero of the simulation runs was January 1 of
the particular year being studied. However, due to lowered water tempera-
tures, low incident solar radiation, and high extinction coefficients, the
phytoplankton growth rate was almost zero for the first two months of the
simulation, therefore there was little or no algal growth. In addition this
first two month period was also the period of high flows (snowmelt and
rains) which, due to the stability criteria presented in Equation (2.21),
required the smallest integration step-size. Since the initial conditions
did not dramatically change over the first sixty days, it was decided to
start the simulation at day sixty for model calibration runs saving some 20
to 25 percent of the running time. Of course, for final verification and
projection runs, the simulation was correctly started at. time equal to zero.
54
-------
WAS 15
WAS15 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, and performs the second order Runge-Kutta integra-
tion. A brief flowchart is presented in Figure 27.
MS 16
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.
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 POP 11/45 at Grosse He, 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 subrountines called by WAS19 include SIR, PLOT, BLKPLN,
and ENCOD (for IBM version only).
Miscellaneous Subroutines
WASP also includes a number of subroutines which perform rather trivial
operations which the user need not concern himself with. These subroutines
include SCALP, WMESS, WERR, SETIA, SETRA and FMTER.
DEC POP 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
55
-------
VARIABLE
INITIALIZATION
t« TSTRT
Figure 27. Simplified WAS15 - WASPB flow charts.
56
-------
common disk buffer, reducing the excessive core requirements required for
separate buffers, at little cost to execution time.
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 sub-
routine to subroutine within the program. The R, I, and * contained within
parentheses after the variable name indicate, respectively, whether the
variable is a REAL (floating point), or INTEGER (fixed point), and whether
the variable is read as input data.
Variable Name
OUT(I)
NOSYS(I*)
ISYS(I)
ISEG(I)
ISIM(I*)
LISTG(I*)
LISTC(I*)
INITB(I)
IPRNT(I)
IDUMP(I*)
IDISK(I)
Definition
Device number for reading input data.
Device number for printer output.
Number of systems or water quality constituents in
the user's model.
System currently having its derivatives evaluated.
Segment currently having its derivatives evaluated.
Simulation type - currently only time variable is
permitted.
User selected option to print exchange coefficient,
segment volume, advective flow and boundary condi-
tion input data.
User selected option to print forcing function
(waste load), kinetic constants, segment parameters,
and miscellaneous kinetic time functions, and ini-
tial condition input data.
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.
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
57
-------
IREC(I)
MXDMP(I)
IDFRC(I)
NBCPSY(I)
NWKPSY(I)
SYSBY(I*)
RBY(I*)
QBY(I*)
NE6SLN(I*)
been reached, permitting the user to write the cur-
rent state variables or segment concentrations to
auxiliary storage (disk). Normally I DISK 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 since 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 incorpo-
rated 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).
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.
58
-------
TIME(R)
DT(R)
TZERO(R*)
SCALT(R*)
TEND(R)
PRNT(R*)
OMEGA(R)
ITCHCK(R)
MXITER(R)
C(R)
CD(R)
CMAX(R*)
CMIN(R*)
PARAM(R*)
CONST(R*)
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.
Not used in current version of WASP (although user
must include in his input data check).
Segment parameters for use in the WASPB kinetic sub-
routine.
Constants for use in the WASPB kinetic subroutine.
59
-------
EVOL(R*) Segment volumes (nominal input units MCF, inominal
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 equa-
tion
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 ex-
change, flow, etc., and the "M" vector (array) will contain zero
slope.
IR(I*), JR(I*) Contain the segment numbers between which exchange
is to take place.
IQ(I*), JQ(I*) 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
60
-------
IBC(I*)
IWK(I*)
IVOPT(I*)
NOV(I*)
IROPT(I*)
NOR(I*)
IQOPT(I*)
NOQ(I*)
IBCOP(I*)
NOBC(I*)
IWKOP(I*)
NOWK(I*)
NOPAM(I*)
NCONS(I*)
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 NOSE6).
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.
Number of segment parameters (for use in the WASPB
kinetic subroutine) read.
Number of constants (for use in the WASPB kinetic
subroutine) read.
61
-------
NFUNC(I*)
NVOLT(R)
NRT(R)
NQT(R)
NBCT(R)
NWKT(R)
NFUNT(R)
ITIMV(I)
ITIMR(I)
ITIMQ(I)
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 ISYS, 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.
Used as a breakpoint counter for obtaining correct
slope and intercept vajues 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.
62
-------
ITIMP(I) Used as a breakpoint counter for obtaining correct
slope and intercept values for the time-variable
functions required by the WASPB kinetic subroutine.
ITIMB(I) Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable bound-
ary conditions.
ITIMW(I) Used as a breakpoint counter for obtaining correct
slope and intercept values for time-variable forcing
functions.
WRITING (OR REVISING) A WASPB KINETIC SUBROUTINE
It is the intent and purpose of this section to define a procedure for
the modeler/systems analyst to follow in order to develop, program, and de-
bug a new WASPB kinetic subroutine or to revise an already existent kinetic
subroutine. The three steps - model development, programming and de-
bugging - are of equal importance in developing a new kinetic subroutine.
It is the responsibility of the modeler to develop the modeling or kinetic
structure. He must understand in sufficient detail the physical, chemical,
and biological principles that form the model structure, so as to be able to
write the basic mass balance equation, including kinetics, for each state
variable or water quality variable in the model. In addition he must, with
the possible assistance of the system analyst, determine what the kinetic
constants, segment parameters, and kinetic piecewise linear functions of
time are to be. It is the systems analyst's (or programmer's) responsi-
bility to program and debug the WASPB kinetic subroutine. He takes the
modeler's mass-balance equations and programs them within the WASP frame-
work. Finally, it is the joint responsibility of the modeler and the
systems analyst to "debug" and checkout the newly written (or revised) sub-
routine. This can be done by setting the print interval equal to the inte-
gration interval (stepsize) and taking a few integration steps. Then after
getting the state variables printed, the user can check the WASP results by
hand computations for a few system-segment-timestep combinations. A kinetic
subroutine should never be used unless this checkout procedure is performed.
It is recommended that the following procedure be used in developing anc
coding a WASPB kinetic subroutine:
Subroutine Design
1. Review and understand the kinetic structure proposed by the
. modeler - check for consistent units.
2. Assign WASP system numbers to the state variables.
3. Assign the order that the kinetic constants, segment parameters,
and piecewise linear functions are to be read in the user's input
deck (the internal assignments are facilitated by use of the FORTRAN
63
-------
EQUIVALENCE statement with the CONST, PARAM, MFUNC and BFUNC
arrays).
Subroutine Coding
4. Determine if any variable and/or program initialization need be per-
formed upon first entry into WASPB. If so, code the variable and/or
program logic, and use the variable INITB to determine primary
entry.
5. If any piecewise linear functions are required in the user's kinetic
subroutine, develop the code necessary to update the slopes and
intercepts (via a CALL to WASPS) and to compute the function values
using the appropriate slopes and intercepts.
6. Code the segment loop, which evaluates the kinetic portion of the
derivative.
7. Code the logic necessary for writing the state variables and other
variables of interest, that are computed in WASPB, to auxiliary
storage files. Using the variable, IDISK, as an indicator that a
print interval has been reached, the user may code the WRITE state-
ment in one of the following two ways:
a. for large kinetic subroutines where core requirements may be
restrictive, it is recommended that the user code the WRITE
statements in the segment loop.
b. for smaller kinetic subroutines or where core requirements
are not of concern, the user should code the WRITE statements
outside the segment loop (after the end of the segment loop),
in order to cut down in I/O time.
EXAMPLE WASPB KINETIC SUBROUTINE
As an example of how this procedure might be followed consider the
modeling structure for a simple BOD-DO model. This model will assume first
order linear kinetics for BOD removal, DO utilization and reaeration. The
BOD removal coefficient, Kr, the deoxygenation coefficient, Kj, and the
atmospheric reaeration coefficient, Ka, all at 20°C, will be assigned as
kinetic constants. However, these coefficients together with the DO satura-
tion will be temperature corrected using a segment parameter for inputting
the segment dependent temperatures. Following the previously defined pro-
cedure:
Subroutine Design
1. Review and understand the kinetic structure. The rate equations
are given by Equations (3.1) and (3.2)
64
-------
dBOD.
-3Tls-Kr(T1)*BOD1 (3.1)
dBOD.
= - Kd(T.) * BODi + Ka(Tl) * (DOsat(Tl) - DO,) (3.2)
The formula for temperature correcting reaction coefficient is
K(T) = K(200C)6T"20 (3.3)
and the following empirical non-linear equation will be used to com-
pute DO saturation (1 )
DO . = 14.652 - 0.41022T + 0.007991T2
s at
- 0.000077774T3 (3.4)
Equation (3.1) states that BOD removal is a first order linear func-
tion of the removal coefficient and the segment concentration (Note:
there are no "kinetic" sources of BOD). Equation (3.2) contains
both a sink (deoxygenation) and source (reaeration) term for DO.
Kr need not equal K
-------
1C
r.
c
LJ
Z
»-
^>
K
^j
C,
rs
'
i
* *
*
* *
*
* *
*
o
* o *
o
1
* *
o
z
<
« X
LJ O
4 C O
Z 1
LJ
IB Z
>- LJ
X IS
o >-
X
-10
«
1 i ^t
M LJ
X >
LJ -1
X O
u to
o to
M M
ID 0
X
LJ
»- M ->>
»- < « « 1 1 1
M O O C 1 1 1
Z ^ N. X
S M « <-l
u
o is
LJ LJ
ts o
LJ IS
O o LJ o- or or
cu o o o o
0 >-»->-
(VI »- 0 U LJ U
« ex e a e
»- L. U. U.
>- « Z Z Z
Z 0 0 0
K LJ M M M
M M LJ LJ LJ LJ
u u. M e; ir c-
M u. LJ or or or
U. LJ M o O O
U 0 LL LJ LJ LJ
LJ LJ LL
O LJ LJ LJ LJ
u z o or or or
O 0 => => 3
« t- z « < c
> « o o- or or
O Z M LJ LJ LI
x LJ »- c. c. a
LJ 19 « X X X
e >- o: LJ LJ LJ
X LJ 1- t- t-
00 f
O LJ LJ K O C
ID o or x x x
1 1 1 1 1 1
e o o M rv «n
(v rv rv < « <
or O « t- >- »-
K >C 1C LJ LJ LJ
m Xix
£
4
^
IT M (VI 1C * f. \C
o
u
OUOUOUUUOL
UNITS
*
0
4
JJ|
4|
41
0
/jl
4t
^
^
^
o
AC
LJ
LJ
X
or
e.
1 U UU
u
ts
LJ
O
in
LJ
or
3
^
^
or
LJ
IV
LJ
»-
Z
LJ
X
IS
LJ
n
1
u
e.
LJ
*
U U
(1?0).DOSAT(120).
^
1C
o
(VI
*
0
»
»*
o
4 (VI
M
or
1C
o
(VI
c
«
D
(V
O
1C
» Wk
e o
(V (V
KM
1C
u
a
_l LJ
< »-
LJ
or
*
U UU
P T IS NOT STANDARD
N OF WASP - AND SHOULD
P VERSIONS
O 0 O
M a
u to
or u
> o
f or
0
»- e
Z K 0
LJ LJ r-
x x m
LJ
»- u r.
« tO 03
»- O M
to
z t-
Z M z
K >-
or. z M
Ik O
O LJ U
tS LJ C
z in 3 in
X I) «
O 1 Z
. _l MX
LJ _)Z
>- O « LJ
O L. or cc
y >-
LJ or »- >-
X O O
»- u. z o.
0
u
UU U U
ASP
lao SESMENT MODELS
FUNT
50».WWKU2.HO).HFUNC(20)
EGtlSI^.LISTGtLISTC
ECtMXOIPtlDFRCd?)
SLN
MIN(12I
X Z to O" IS LJ
O »(Vf*»M*-!LJ »
C) (VI M tO JC Z (VI
icvft-ca
OX SC O M c <£
o: LJ x is - v z x~
<[ (- Z 0 LJ tD LJ LJ IT>
^ M(vt«NO**» «»^
ztn (vi zco>-»-ero>-
C) ^ O ^ CD ^ LJ (V (/"
XM *-.tt>-I>-LJ>--C
C) LJ * * * V) 2 C/"l V >< CV CJ
U O 030 C C » X M »
O »(VM»»i/;ciiroo
ZM ^'^^^Q.LJLJLJ*14
*I OCT^OZX^^**
X z»«DOorx»-Lj e
«o m»o.z»^e(M
LJ u. a o tr or »>-c M »
C3 Z(VI*LJCOtO»<»X
LJ X * M ^ M ^ O. LJ t5 ^ ^
^O K'^D'^^LJTLJ'^or
M^ GO MMZ^GLJH
" Z 2
-l-l Or.
« _J LJZZZZZZ?
LJM ISOCOCOOC
:X _»jLiy*»»ar^^
««<-rr> I y *^ >
LJLJZ-OOOCOCC
OrOrLJLJLJLJLJULJ
U L) U U
Figure 28. Simplified BOD-DO WASPB kinetic subroutine (Version 1)
66
-------
O I
O LJ I-
c\. o in «
o. o v>
O « ~ Z _)
Z O JC IO 3 O O
r> » « o >
u - t- o
C CM K> LJ K O «
. ,1 I Z LJ
fj t- t- in « z o. ~
o X >« vi»cc - c
» x z o en »- LJ o
(I O >- O \B >- I- 13 »- O
CX CE J U O ZO I
OCMZtt t-LJ IJCT CL
e?»->~ »cn jr tr
a- r, - - z «= o on LJ LJ
'jtooocrz*: MM t- in
r.rv*~LJLJ»-orrG* U. * *-i
4*«fr4G>O «*«» U. *4
- o t- LJ ru »- e> LJ en c\ »-
ou(L7i-iui I n i o o a c. e
c\/ico=>>z »-**«-< of, r» c^
» M u L. » in-~«u~' « o o
cv»czcuc/i ? rv o: e: s Z j» o o
H^k-4»^(»M OCVtftf tf«<» O «-
c v o - a o cr <: >- e>
OITOZK-CCJC I >- LJ I Z ct iC >- LJ
cctvoozzc) c/> i - a. or * a e
CD Z(_1D""OC^ X Z >- X I «^- T LJ
o»zut-_«irec»e L.O LJ au w
3O 3 U «4 «CO I»- «-"
o«-iro«-iz » >- ui«
»-«M oz KO.
» o. jc
i~in»>o3i-izz t-tnt-c.ccL~*f t- cv w '
os-rw-zzt-io Z X > cr «: u>- cg»- o o*
a>- tr--«-OLJ <«oowiLJ"-«ieto: u or 0.0.0.0* >
scoo t->oo»-ruLje:t-ZJjt« ct i- KZE^O. o o
-Xi-iZCVLJlO WLjZtr LJ»OI»Q. KZ ULJLJ»r DO
05:»»«[E3OLJt^«Jt>- U X ~ O LJ O »->-»- LJ _J CCT
ITi»^k» Z^«IJ^«XK»«yvH«ff) (JIT) O^- CD*
rvucD.X't-^eciv>"f'>-K* J««o***
crino3>-«-c LJ LJ » Zii cc»« o. K^f-Z «a« f- »- U>u
cr>> >-iZzcr«>-o»-«a.TLjir 3 t- LJ i >->-- «^ i- LJI^
ca-'-r^iCLtf- v. »-ctsc K »- o >- LJ LJ LJ « r- Z is IT. -.
-r~^> >-» Ox oZLJV-orcucr «zooo ill >cr- LJ «->
r . . o CV or - UO.X »OZZ«KZC CC >- -LJOt-h-l-l'l-O LJ K LJ O
'-*^7*^7»« 7 C Clfr-LJOOfc-O LJ li.i/^LJfcV^LJ^O I O »>frLIV
?c« *-u u * f t- o.z u o m o o o <^. «o IB x z e»-2c>i
«>-£.*-> > LJ V Z ->t» t/JU. I Id LJZ*-ryf\,fV!- t- » IB LJ >- I -
. ^ r, c O _l » OOO. »> » LJ O CL LjaoO' CtO< HC - IB »L.-Ub-e >-
. - C- be O - _> O JC »«/) LJ K »-n«O« « < II
> 3 > >- _J LJ a LJ >- LJ or f, LJ O LJ II a Q. M LJ t£ » l/> «- >
-~«->Z>- O2X O.OI/1CXOOZO X LJ 7 2- LJ (2 X II II II Z 19 I LJ»*-_JOr »*
U. I- SZZ3ZUOZ «^t-O«t-LJLJ LJ LJ M V> V. -CLJCSCCC
ILJOILJ »«LJ I ^^Ot^d*-^^'^*11'^-!/) LJ^- LJ f-4 C\.>OLJLJLJ
LJ tr jc_io»:_iLjt-_i ii a t- »- a - uidta - s LJ 11 5 o er irci^
^^?zzz I o z o « LJ c > er « i- ce»z 11 LJ LJ LJ 11 zto z o «J c " *-
'-oooooc H-L.O -> »">«:o> bjoocezoe v> v \f> * ». LJ cc »- on H o o
sy53ss» » c.»--a«-^a tn*or>>-oma»-«->-a« *-ct»-t- "-1 ' co « o r>,
^syysas z 3 3 3 o z Zioz«-z o _i
o c ~. o o o o o o o c x LV OLUCTO<|MC c zo ocoo c c
OCJL. OCJCJCJ O LJ UJ LJ Z ** OKKX^*"O U »«o OILC> O C)
O O O U O U L> O U O U O O U OOO
c.ooocoo o o o e o o o o o o c - c o oo ooeo o oc
OOOC3COO O O O O O O OOOCCCO O OO OOOO C: OC
Figure 28 (Cont.)
-------
OJ
n
* tc
* J?
o
ff> VI tC.
ui o.
UI * O
»- UI U. Z
« »r «
1 Z «- 1C
i »- u tr> v>
n x o * »-
* z o o z
se _i « «J <
tn u X at »-
*»-». o v>
O C Z U Z
M 7 K U O
a u o in u
x Z cr ui
»- »- u it) r>
tc < a. >- o
fr- IT UI Z 0
O 3 kl UI UI
_i o o ozzz»- n
< H- BC C> **** Z ***
> O UI >-»->- U
t- Vl 3 3 3 O
IB IS r> in O O O Zn
LJ LJ o. L. a a: c uj i
r r> ininzKDO) a> am v>
«iC.LJ«nOOOO33 3 UJ < IU
_>TCK ZZUO V) W " V. »OX -
»-«O"-iUI»-« (919 K
Z > CJ t- »- w< »« C« U Itl 4J LJ W EC O «
1KO.3 « II LJ O < «O V) MCO trt LJ »- K
CklLJ B. lb«S»>K O O MO 1-1 t- O
c. »- x x. LJ LJ cc o » _< _i _j » z> r o.
Z M LJ U O V>I/>LJV)LJ U U «4U CWCL.K XK
» »- . x u - >- in cr LJ v x v r n LJ LJ
»- o. »- - » LJor z z uz o o i- »- I
s z f- »- tn «o o i- i u ui ui o LJ
n.KZ»«utuieur-3ui- LJ c x o o v> z m
3>-orLJionrcDtoinx«(r «i ui o «>-i
OLJC.XLJOLJ ««-«-xc »- a: *-K o K a MZ *
e >- _t o «n M»»ec>e:*-xLj e. LJ ui LJ UI~K< " e
>c«»-ffi>kjar>«ni-ie.« u. u. to u. m LJ ui
u u oxti. r> C »- «r t->-">o«jecOorLJL,3 3 Z 2 Z JC CD -
. - .
O I LJ O. LJ U&*'Cj»«lJ»"«-ib.«-i»-lL.O»-
C L> a L. LJ X *- U * * *- X ~ LjlfoUaXollKCt.0
o o < ec u. x ec « ui LJ viuK
a «_> i »- «n o t- ou»OLJV>»ouiir»ui e
se o KO « CM u e 3 -J e> .j » .1 j » - _j»- .- _i *> >i LJ n c c.
LJ " > > 3 Ij. > u.li.« _JL. _JUO _JU.O U.IO Z Bt X
XO4CGLJ LJeruo_/J»-i->_j>««»-iJ>.J<«»- fftno Z u«> Z c-
_
t- o or OCKZO o o « a «OK « o or z ouz LJ or - v/r> o -J
Ck.L>> * » LJ U MI-I u> Cl O » O O B LJ MKU »-LJV »-* (L 3
or o 1-1 z u
e o* OLJUU>i»-f-
- -
CLZtruie
S >
a. tc *
O«J OUCI U U OUU O Cl O Cl U O Cl U UUUUUUU V> «-> * UI «->
O O 3 O 3
m or or O u Ui V>
UI O « UI >
erXOarjcxU.
r>yy»-o
»~LJ'CUi««
'^i-tt-iio
!_ o w. o a. ui z
u o or LJ
Figure 28 (Cont.)
68
-------
Line number 0002: dimensions or allocates storage for the tempera-
ture corrected rate coefficients and DO satura-
tion values.
Line number 0023: conserves storage by equivalencing the dimen-
sional vectors to the unused portions of the seg-
ment parameter array (PARAM).
Lines 0025-0034: computes the temperature corrected reaction co-
efficients and DO saturation. This will be per-
formed only once (saving some computer time) by
checking the status of INITB (line no. 0025).
Note the user must reset INITB himself (line no.
0034).
5. Piecewise linear functions. Since there are no piecewise linear
functions this procedure was not coded.
6. Segment loop. Line numbers 0036 through 0042 contain the body of
the kinetic subroutine. Note the inclusion of the segment volumes
since at this point we are computing Vj dC-j/dt.
7. Saving the state-variables. Line numbers 0043 through 0048 provide
the instructions necessary to store the state variables on disk
files for the IBM 370 series computer, while the Comment statements
immediately following show the code to be utilized on the DEC POP
series computers. Note the user must set MXDMP, as shown on line
number 0024.
Figure 29 shows a second version of the subroutine which has temperature
as a piecewise linear function of time. Since temperature is time variable
the reaction coefficients and DO saturation will also be time variable, ne-
cessitating their evaluation at every time step (thereby eliminating the
need for step 4). Lines 0025 through 0027 show the FORTRAN code necessary
to update the piecewise linear functions of time, while line number 0028
shows the computation of the piecewise linear function at time = TIME, using
the MFUNC and BFUNC vectors. As in the previous example there is sufficient
storage available (unused portion of the PARAM array) to compute and store
the temperature corrected reaction coefficients and the DO saturation values
and output these numbers, together with the state variables, outside the
segment loop. However, for demonstration purposes the code was included in-
side the segment loop. Note: the variable MXDMP is now set equal to four
to reflect the fact that three additional variables per system are being
saved; and the WRITE for system 2 also contains a dummy variable, DUMMY, to
"pad out" the WRITE to MXDMP equal 4.
REFERENCES
1. Solubility of Atmospheric Oxygen in Water. 1960. Twenty-ninth progress
report of the Committee on San. Engr. Res. of San. Engr. Div., ASCE,
Jour. San. Engr. Div., Vol. 86, No. SA4, July 1960, pp. 41-53.
69
-------
*
*
*
*
*
*
«
*
*
^
,
41
CD
ft.
V
tf
X
41
41
LJ
2 *
**
^ 41
3
O
or
CD
3
in
CJ LJ (J
o
o
o
u
o
ffi
1
o
z
<
X
LJ
0
Z
LJ
19
*-
X
O
_l
«
u
»*
z
LJ
I
CJ
o
4
CD
Z
LJ
fr- V*
in
^
in
CJ U UU
*
*
41
41
*
*
9
*
41
^
41
O
o *
1 *
z
LJ
es
>-
X *
0
c\
LJ *
>
_J
0
in 41
en
** 4>
O
0
41
CM
0
CJ 0 LJ 1
in v >- v
»- « C < 1 1 1
-> C O O 1 1 1
Z X V -V
3 »4 11 rH
CJ
CJ C
LJ LJ
IB O
LJ CB
o c LJ o- or or
01 n o n o
ce - i- -
iv t- o LJ o u
c CM « « a
U- L. U. U
»- < Z Z Z
z o o o
t- LJ i- «- >-
LJ LJ Z U O U
i l-t LJ LJ LJ LJ
CJ U. i- o: cr or
i u. LJ ce cc cc
U LJ -> C O G
L- O U. CJ U CJ
LJ LJ U.
O LJ LJ LJ LJ
U Z O CC IX CC
O LJ 3 3 3
> < C or or or
O Z ** LJ LJ LJ
z LJ t- a a a
LJ U C Z Z Z
a: > cc LJ LJ u;
X LJ t- t- t-
o o <
O LJ LJ or o e
CD o or se je X
i i i i i i
o o o * CM m
CM CM CM < a <
CC O < »- (- t-
X JC X. LJ LJ LJ
in XII
Z
a
|w
in « CM K> « «-, MI
z
o
u
J IJ CJCJ U CJ LJ U CJ CJ
*
*
41 in cj
C3
Z LJ
30
en
LJ
or
3
»-
0
or.
LJ
O,
Z
LJ
K
41
1-
Z
LJ
Z
19
LJ
I/I
Z
3
Z
»4
X
r
41
1
0
U
ft.
in x
LJ t-
* t-
LJ
X
* K
O.
41
CJ CJ LJ LJLJ
O
0
M
or
LJ
a.
z
0
»*
-
<
_)
3
Z
en
LJ
X
-
K
LJ
0
CJ U
*
*
*
*
41
,
41
4>
4^
4|
,,
in
* z
o
to*
1-
LJ
Z
3
L.
cr
«
LJ
Z
f*
_J
LJ
in
^
LJ
LJ
ft.
0 U U t
RIATION OVER TIME
e
>
CJ
or.
3
i-
«
or
LJ
ft.
Z
LJ
t-
O
LJ
ISI
**
J
<
Z
or
o
z
i
t-
a
z
LJ
^4
j u c
*
41
41
41
41
*
*
41
41
4)
41
4)
» CJ CJ
0)
P Y IS NOT STANDARD
N OF HASP - AND SHOULD
P VERSIONS
CM O 0 C\
o
z
LJ » or
1- 0
i- <
« z - o
* LJ LJ r-
x y «
a LJ
X 1- LJ T
» < en cr
cc t- c *~
K m
z >-
0 X " Z
CM « «
C CC >-
vr »- _i z
or z »*
e o o
CM L. O
O O LJ
if 19 LJ O
z in =>
0 » 3 _J
CM I CJ
K >O 1 Z
JC _l »
LJ _) Z
1- O c U>
O L. or CD
« LJ CC H-
LJ TOO
or »- u. z
u u u u
CM
f^
O
ftSP
1?0 SEGMENT MODELS
x
o
or ^
o
L. **
in .
o y
er LJ
«: t-
LJ I/I
>-
Z in
0
Z CM
X ->
o
CJ O
H-
o
Z »H
T
0
z cr cr
LJ U
LJ O
LJ X
en »- o
z _i
4 ^t i
«t
x «
J J
< J
LJ ««
or 2
>-
ft.
O
CJ
LJ LJ LJ LJ
o o e o
FUNT
?0) .MWK(l?.i»0) .MFUNC(2ni
CG.ISIM.LISTG TSTC
? - v,
CM !-.
CM i in
t-c LJ >-
~ d >~ V
t- s. cr; i-
JC . O -
» ir
Z 0 LJ
» in ru i/~.
CM -> O
CM ^
«1 Q >-
z m tfi
*- * a. t-
LJ U-
tr o c
Z in CM s
» C\* t-« »
»-- »-
c a >- 3
Z 5 cr o
» IT, -
t- >- r
tr o 10 c
Z CM 0
*H >- »-<
*- 3 -
J -I O ^
0 C -
> >
z *
or
LJ ?
13 C.
J _J L, 5
C t H"- E"
LJ LJ S C
CE or » t>
0 O O O
fO ^ If \£
O CD O 0
O O O O
o c> o c
Figure 29. Simplified BOD-DO WASPB kinetic subroutine (Version 2)
70
-------
o D. CL c. a a
r- D c o c o
r> c. a. a a c.
z u e. u (_ o
CC LJ LJ LJ LJ U-
-i O C. O D C
o eo 2 z
c. c & a c-
CVi CT> LJ LJ
o H c «-
u. CM co ~ _j «j MZ o
o »4 c rti ? or ~ >- bc_i«r
- - O 2 .H O C VI >- - V) LJ X Z 3 or
je CM « * 2 c u i-i »- >- t- O
C. * SO l> O or 2 OO i-it-ili.
> «D cvi O<^« LJ|i.«« «LJ C/l
x *4«~****c)KLj ey xcQfO looo:
v o o n. ^ 1-1 LJ I o i «« r^ » >- i- LJ LJ LJ
2 2cv>ccor>»2 »- ^r * K-OK> ore »->-»-
tJ_J «-i»l-Ob. CO - « LJ t- t- KK C«2
LJ f, I r\. ec z cv m 2 fv o: L. »- a ? 2?3 LJ
ori*; o*H^->-i»«^i-i ocvf« o«LjZ3 *»_j oo OIAO 2
»"LJ o tr, - u »- o. csoii. «-i« K« «» ci K
2»- oir. o2»-nbt i t- LJ i o v> in t-
Sf>2C O>-l LJLJI/I 3
lfl>-cr»«CB 2U3iiOcn I2>-T O K U. 2 LJ LJ OOCI O
M c a. « »o»2u.>-jir IBO-CS MU2X r>K KII LJ c
OO»X~"5*-»2««DLJ 3 O 13 »< »- (E l-i I-i I U.Z 3 UILJIT CD
^D«ec>»o.sr»» f- c~inc«- u >-LV cc>-" O.LJ ccmo 3
>- 2 » m2 c<>-L_2Z: or o: x 2 o o » » _j y i D 10
o uj L>incv:>oaQ.xr oa I >- x y 3 t- z»i-i * »- « o ^- _i _i «
*«a:>- »r-L.ri«-ica«-> *i K- u *- J_i LJ
»**>oz(v»z2KZO 2z o «-> »MCCQ. «- i o »< co
«o>->-eto^-ct» ot»»«-i o o viooinLj o: LJ n o i ^ ce LJ ui x» o: o
e_)Ljcuio»icoo >- o «9 - r>2CVUn tO LJ Z Or LJ LJ2 2O 3 ZixLJLJ tnifl U O
»">-«->*--oo*'»»«-iccy OLJIO «se t- 2 O r> u. er u. o -!>-» eror u V.
r3i/*&OK<^Ljir. »^.»-«^2»-i UJIMX ^-*»Vv4 MV^-LvZLj 2 fr-m o.*«>^^ LJLJ o:*-i 2
C»>X»«>CUOCCLX*»- _i cr z «r c _i u f * >- i r> >- »- »- >- i LJ
»->-C- CITIOJI*-*" LJ LJ 2 * QT >4 I O U 2 LJ O,O2 33 UILJ O.
(oa^ocarzr^ccoo- CD f< in o »< « ~ LJ t- e v. z 3 x i-t-i- ir« O.CL I-LJ o
»-o.LJ«ju> «r» 2 2 n- « >- u»->-i-t ox U2u»-ar * K OOD ts »- _i o -c«i0..ocvor»- u> a » MOZZCV LJ _i . x- ^ LJ K «t-cc. LJO.«-U
^^»»-»-»«V»^WO'^-2**y»-l Z KO^LJCLe CJUILJV>C»« »"^* VI OO OO-yOtBLJU-
*'^Cl "^ZO- » »«IJ LJ . Z . LJLJZ ZLj*«-« - O»-0: I^O*-O:UJL.3
tr i^; <«: -in>-£i.»-> XLJV z ""LS o i «-" -i c>>ec>ji>t- _/ LJ o. LJ »- LJ r o LJ to i- o LJ «- 3 ooearL.jr«cDLj WLJ
»-z»-CLjo.«r, r«->-»-2i-« oix cvuciaur x »-cruo.ccnt-»-u ii.ujLJ3O«-«»-o 1.1 n n c
Lv t- 3Z2323 Z»<*-»3C.«*-LJZ o: L> < »- >- OLJ
iLjoiuo I 3 Lv H- inLj u XX o or or o 3 -i
___. h?£_ JCJOjcj n_j zt-«ern u oi^LjuriiLJOO-Q- r*._Ji-
Z2222Z22ZZZZZ XOZ O « u « »- «»-3inXu * LJ > > 3 t-UU. T---IU.
COCOCOCOCOCOO ^U.O *^> ^*> LJO.>3U.3 ZU.O2 XOCtfCL.>O LJOLJLJ^
Z>779FZFyrryV X* CL»*»CL»"*« KZLjCr2V_J»^Z2O. U^-ZlTKT'ZLjLJ^-UCrL.LJ^
zz>zzF>rszzzr * 3 3»- o . - j ix Tioo^uar _)
COCOCOOCOOCCO O O O« X L. « t- LJ u. >- C O Cr CO <
ULJLJUOLJL;LJLJLJLJLJL> u u LJ O X ««-i LJ -> >- » o. LJ u jt > "-1 LJ
«J O U ^ U LJ U 4>U U U O U U O «J DUO «J LJ LJ LJ LJ U LJ LJ
OOOOOOOOO O O O O O
OOCXXMOM^^MXX M CM CWCM CV «V (V CW CW
oocoooooooooo o c oo o D c o e
ooooooooooooo o e oo o o c o o
Figure 29 (Cont.)
71
-------
OQ.Q.Q. O O. 0. O.
r-OOD »- O O O
«o a. o. o_ m a. a. a.
x LJ u LJ x. LJ <_> LJ
a ui ui LJ m LJ LJ LJ
MOOO «-t O O O
Ul
^
f-
4
cr
Ul
CL
X
UI
V-
K
ID
LJ
in
»«
^
X
X
CL
FICTENTS
"i.
LJ
O
°
Z
o
^.
u
c
LJ
nr
^.
U
LJ
CC
007991*TEMPSG«TE*P'
o
ts
in
0.
u
^
4
CM
IV
Q. a a. o
X X X 2 -
ui ui ui o a-
»- K >- ^-
IS
cn
a.
F
Ul
l-
*
CL
X
LJ
t-
in
a
x
U'
41
^
|^
f«.
^.
r-
0 »-
o <->
C LJ
O.KR.KO.Kft
UTIME
O 0
CD CC
m
in
94
Ul -H
4 in
4 O LJ
_l O
Ul LJ LJ
k- *>v _l
_l - 2 >
O 10 E Ul U
> LJ » a
* ' O *
O CD J
0 « K «
C(l.ISES) .KR.KD
»
^ .
u
cr
ra
^
»
v*
94
Ul
>-
b«
er
o
o
o
1
4
in
O
^
4
^
JC
*
o
o
CD
T-
3
O
U
Ul
o
o
0
tf
in
O
,
o
o
#*
CM
94
cn
Ul Ul
_J f-
eo -
e er
UTIME
o
or
CD
CO CM
94
LJ
in u
O 0
_l LJ
LJ _l
V t-
2 U
LJ
Ct _l
C _l
o
u
Ul
o
o
o
9-
«
in
o
o
o
o
»
LJ
or
U.
O
»
Cy
+*
^
LJ
^H
94
or o or
n
94
10
Jf
o
2
Cf
3
9»
Ul
K
O
2
or
C
CL
O
2
o
I
10
in in « en CM *
o.
O _ _ _ _
o
O f- O LJ _IO»-»_ILJ~2 CL
X X
LJ LJ
t- I
O
2 O CD _l "-I O CD - « 94 O
o _ _. __ _.____.___
v-xcro**«fr-o*9*o cc *4 ui LJ x p* x ui x 4< m
to ZLJOfM»»««ouOLJ»oru O cr u LJ
2 19 * I « « C O O I « X > 3 >l>3 * CT f"
IS LJ IS 2 Ul »- 9- 9- 9- CM 9-t«.-eO « I CD In CDC 2
O m er 10 >-
is - IT, ui >- CL 9- 9- t- in - « o v. o o M « < LJ r-
l/^ 94 » O r- CC LJ O G O r» *4 9»CLJI3(/)»O**CCtS(/>*O** O C^ _J
LJ cc C- d ccC)*~-J>-i-< «->H-_l^-~^l^ 2 CC - X
on iicmiw o»- ouc^eoeccc-eoeoo^ujKer Occ K o
O O O C LJ LJ Or D « C C O U U U C U U LL O CUI2 LJ Or >- fir C I
OCCO> K H>lclca£ O O LJ * *»« (J i^ ^MNMLJ »«eiUI t-U'C' K«-< CL 3
K O -" 2 U
o e ouiLjLjtvf'
«i e o-»-uiX2«rt^
i b.2cruiO»-ui
3 ». a cr * 2 LJ
LJ U U UUOUULJLJCILJLJULJ M >->£UILJ
O O 3 O 3
V) cr cr O LJ LJ
ooocoooooo o o o ee OCOC*-LJ«UI««
oooceoocoe o o o o o oeooc2t-ce-i-if3
ui r> i^ O e LJ z
Lv O CT LJ
Figure 29 (Cont.)
72
-------
SECTION 6
WASP INPUT STRUCTURE
INTRODUCTION
This section will describe, in detail, the input required to run a
user's WASP model. No attempt has been made in this section to detail the
job control language (JCL) commands required to execute the WASP program
since JCL is not only machine dependent but often is site-specific even for
the same series of computer.
To arrange the input data into a logical format, the data cards required
are divided into fourteen card groups, A through N. The card groups are
briefly summarized in Table 5. For each card group, a brief description of
each card is given to define the variables which appear within the group,
and any options which 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.
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.
INPUT DATA
Card Group /\
Model Identification and System Bypass Options
The variables which appear on each card are as follows:
1. Model Identification Numbers
73
-------
TABLE 5. 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
6. Parameters
1. Number of Parameters
2. Scale Factors
3. Parameters
H. Constants
1. Number of Constants
2. Constants
74
-------
TABLE 5. (CONT.)
Card Group
I. Miscellaneous Time Functions
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
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
75
-------
Card Group A
broup ^\
md Syste
Model Identification and System Bypass Options
5 10 15 20 25 30 35
MODEL ISER IRUN NOSEG NOSYS LISTS LISTC
FORMAT (715)
MODEL = model designation.
ISER = series designation.
IRUN = run number.
NOSEG = number of model segments.
NOSYS = number of systems.
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.
2. Title
1
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.
3. Simulation Option
Presently WASP permits only time variable simulations, therefore include
the following card:
1
24
TIME VARIABLE SIMULATION
FORMAT (6A4)
76
-------
4. Systems Bypass Options
SYSBYCTrSYSBY(2T ... 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).
Exchange Coefficients
Exchange coefficients may be inputted in one of two formats, actual ex-
change coefficients or, they may be calculated from inputted dispersion co-
efficients and accompanying cross-sectional areas and characteristic lengths
as per Equation (2.19).
1. Data Input Option Number; Number of Exchange Coefficients
5 10
IROPT NOR
FORMAT(2I5)
Data input options:
IROPT = 1, constant exchange coefficients.
2, all exchange coefficients proportional to one piecewise
linear approximation.
3, each exchange coefficient represented by its own piece-
wise linear approximation.
= 4, constant exchange coefficients calculated from the dis-
persion coefficient, cross-sectional area, and characteris-
tic lengths specified for each interface.
= 5, all exchange coefficients proportional to one piecewise
linear approximation, calculated from a piecewise linear
dispersion coefficient approximation and respective cross-
sectional areas, and characteristic lengths.
= 6, each exchange coefficient proportional to its own piece-
wise linear approximation, calculated from a piecewise
linear approximation for the dispersion coefficients, cross-
77
-------
sectional area, and characteristic length 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.
2. Scale Factor for Exchange Coefficients
10.
SCALR
FORMAT (E10.3)
SCALR = scale factor for exchange coefficients. Exchange coeffi-
cients are normally expressed as million cubic feet per day
under options 1, 2, and 3. Options 4, 5, and 6 normally re-
quire the following units:
Dispersion Coefficient - Square miles per day
Area - Square feet
Length - Feet
The conversion of sq. mi. - feet/day to MCF/day for options
4, 5, and 6 is handled internally in WASP. If units other
than the normal units are necessitated by alteration of the
space and time scales, SCALR should be set such that the
product of SCALR and the exchange coefficient (or the
equivalent computed from the dispersion coefficient) yields
MCF/day.
3. Exchange Coefficients
The data input format is determined by the option selected.
Option 1
Each card in this package contains the exchange coefficient information
for four interfaces. The number of exchange coefficients read is equal to
NOR. The information on each card is described below:
10 15 20 30 35 40
BRTKJ IR(KT JR(K) BR(K+1) IR(K+1) JR(K+TT
50 55 60 70 75 80
BR(K+2) IR(K+2) JR(K+2) BR(K+3) IR(K+3j JR(K+3)
FORMAT (4(F10.0, 215)
BR(K) = exchange coefficient between segments IR(K) and JR(K) in
million cubic feet per day.
78
-------
IR(K),JR(K) = segments between which exchange takes place -- the order of
the segments is not important; if a segment exchanges witn
a boundary, the boundary is specified as zero.
Option 2
The card package consists ot two sub-packages. Sub-package I contains
the exchange coefficient data, while sub-package II contains a detailed
specification of the piecewise linear approximation to which all the ex-
change coefficients are proportional.
Sub-Package J_ - Exchange Coefficients
Each card in this sub-package contains the exchange coefficient informa-
tion for four interfaces. The number of exchange coefficients read is equal
to NOR. The information on each card is described below:
10 15 20 30 35 40
BR(K) IR(K) JR(K) BR(K+1) IR(K+1) JR(K+1)
50 55 60 70 75 80
BR(K+2) IR(K+2) JR(K+2) BR(K+3) IR(K+3) JR(K+3)
FORMAT (4(F.10.0, 215))
BR(K) = ratio of the exchange coefficient between segments IR(K) and
JR(K) to the piecewise linear approximation.
IR(K),JR(K) = segments between which exchange takes place. NOTE: the
order of the segments is not important; if a segment ex-
changes with a boundary, the boundary is specified as zero.
Sub-Package _II_ - Piecewise Linear Approximation
The number of breaks required to describe the piecewise linear approxi-
mation is followed by a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; an exchange value,
and a time (normal time scale is days). The input is as follows:
5
NOBRK
FORMAT (15)
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 70 80
_ _ _
RT(K) T(K) RT(K+1) T(K+1) RT(K+2) T(K+2) RT(K+3) T(K+3)
FORMAT (8F10.0)
RT(K) = value of the approximation at time T(K), in million cubic
feet per day.
79
-------
T(K) = time in days; if the length of the simulation exceeds
T(NOBRK), the piecewise linear approximation will repeat
itself, starting at time T(l); i.e., the approximation is
assumed to be periodic with period equal to T(NOBRK), this
holds true for all piecewise linear functions time.
Option 3
Each exchange coefficient is defined by a package of cards consisting of
two sub-packages. The first sub-package identifies the two segments between
which the exchange will take place, and the number of values comprising the
piecewise linear approximation. The second sub-package defines the piece-
wise linear approximation which describes the exchange coefficient. The
input is as follows:
Sub-Package l_
5 10 15
IR(K) JR(K) NOBRKlK)
FORMAT (315)
IR(K),JR(K) = segments between which exchange takes place. NOTE: for
exchange only, order of segments is not important. If a
segment exchanges with a boundary, the boundary is speci-
fied as zero.
NOBRK = number of values and times used to describe the piecewise
linear approximation. All exchanges must have the same
number of breaks, and all breaks must occur at the same
time relative to each other.
Sub -Pack age ll_ - Piecewise Linear Approximation
This consists of a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; an exchange value,
and a time (normal time units - days). The input is as follows:
10 20 30 40 50 60 70 80
_
RT(K) T(K) RT(K+1) T(K+1) RT(K+2) T(K+2) RT(K+3) T(K+3)
FORMAT (8F10.0)
RT(K) = value of the piecewise linear approximation at time T(K) in
million cubic feet per day.
T(K) = time in days. All break times must agree for all segments,
i.e., T(l) must be the same for all exchanges, T(2) must be
the same for all exchanges, etc.
Option 4
Each card in this package contains the information to calculate the
exchange coefficients for two interfaces. The number of dispersion co-
80
-------
efficients is equal to NOR.
below:
The information on each card is described
10
E(K)
50
E(K+1 )
FORMAT
20
A(K)
60
A(K+1 )
(2(2F10
25
IL(K)
65
IL(K+1)
.0, 2F5.0,
30
JL(K)
70
JL(K+1)
215))
35
IR(K)
75
IR(K+1 )
40
JR(K)
80
JR(K+1)
E(K)
A(K)
IL(K)
JL(K)
IR(K),JR(K) =
Option 5
dispersion coefficient for the interface between segment
IR(K) and JR(K) in square miles/day.
the interfacial cross-sectional area between segments IR(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.
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.
The card package consists of two sub-packages. Sub-package I contains
the information necessary to calculate the exchange coefficients, while sub-
package II contains a detailed specification of the piecewise linear ap-
proximation to which the dispersion coefficients contained in sub-package I
are proportional.
Sub-Package l_
Each card in this sub-package contains the information necessary to
calculate the exchange coefficient for two interfaces. The number of dis-
persion coefficients is equal to NOR. The information on each card is de-
scribed below.
10
E(K)
50
E(K+1)
FORMAT
20
A(K)
60
A(K+1 )
(2(2F10
25
IL(K)
65
IL(K+1 )
.0, 2F5.0,
30
JL(K)
70
JL(K+1)
215)
35
IR(K)
75
IR(K+1 )
40
JR(K)
80
JR(K+1)
81
-------
E(K) = the ratio of the dispersion coefficient between segment
IR(K) and JR(K) to the piecewise linear approximation.
A(K) = the interfacial cross-sectional area between segments
IR(K) and JR(K), in square feet.
IL(K) = the length of segment IR(K) in relation to the IR(K)-
JR(K), in square feet.
JL(K) = 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),JK(K) = segments between which exchange takes place. NOTE: for
exchange only, order is not important.
Sub-Package U^ - Piecewise Linear Approximation
The number of breaks required to describe the piecewise linear approxi-
mation is followed by a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; a dispersion co-
efficient and q time (normal units are days). The input is as follows:
10 20 30 40 50 60 70 80
_
RT(K) T(K) RT(K+1) T(K+1) RT(K+2) T(K+2) RT(K+3) T(K+3)
FORMAT (8F10.0)
RT(K) = value of the piecewise linear approximation at time T(K),
in square miles/day.
T(K) = time in days.
Option 6
Each exchange coefficient is defined by a package of cards consisting of
three sub-packages. The first sub-package identifies the two segments be-
tween which the exchange will take place, and defines the number of values
comprising the piecewise linear approximation. The second sub-package de-
fines the piecewise linear approximation which describes the dispersion co-
efficient. The third sub-package defines the interfacial cross-sectional
area, and the characteristic lengths of the two segments involved. The
input is as follows:
Sub -Pack age _!_
5 10 15
IR(K) JR(K) NOBRK(K)
FORMAT (315)
82
-------
IR(K),JR(K) = segments between which exchange takes place. NOTE: for
exchange only, order is not important.
NOBRK = number of values and times used to describe the piecewise
linear approximation. All exchanges must have the same
number of breaks, and all breaks must occur at the same
time relative to one another.
Sub-Package _!!_ - Piecewise Linear Approximation
This consists of a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts; a dispersion co-
efficient, and a time (consistent with the normal time scale of the model).
The input is as follows:
10 20 30 40 50 60 70 80
RT(K) T(K) RT(K+1) T(K+1) RT(K+2) T(K+2) RT(K+3) T(K+3)
FORMAT (8F10.0)
RT(K) = value of the piecewise linear approximation at time T(K),
in square miles/day.
T(K) = time in days; all break times must agree for all segments,
i.e., T(l) must be the same for all exchanges, T(2) must be
the same for all exchanges, etc.
Sub-Package III
This card defines the interfacial cross-sectional area and the
characteristic lengths of the segments involved.
10 20 30
A(K) IL(K) JUKj
FORMAT (3F10.0)
A(K) = the interfacial cross-sectional area between segment IR(K)
and JR(K) in square feet.
IL(K) = the length of segment IR(K) in relation to the IR(K)-JR(K)
interface, in feet.
L]L(K) = 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.
83
-------
4. Exchange Bypass Options
2 4
RBY(1) 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).
Card Group C
Volumes
1. Data Input Option Number; Number of Volumes
5 10
IVOPT NOV
FORMAT (215)
Data input options:
IVOPT = 1, constant volumes. Currently WASP only permits constant
volumes.
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.
2. Scale Factor for Volumes
TO
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.
3. Volumes
The data input format is determined by the option selected.
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.
84
-------
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.
Card Group J3
Flows
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,
NOQ = number of flows.
If no flows are to be inputted, set NOQ to zero, and go to Card Group E.
2. Scale Factor for Flows
10
FORMAT (E10.3)
SCALQ = scale factor for flows, flows are normally read in cubic
feet per second (cfs).
3. Flows
The data input format is determined by the option selected.
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.
85
-------
10 15 20 30 35 __ 40
BQ(K) IQ(K) JQ(K) BQ(K+1) IQ(K+1) JQ(K+1)
50 55 _ 60 _ 70 _ 75 __ 80
BQ(K+2) IQ(K+2) JQ(K+2) BQ(K+3) IQ(K+3) JQ(K+3)
FORMAT (4(F10.0, 215)
BQ(K) = flow between segment IQ(K) and JQ(K) in cfs. AESOP conven-
tion is: if the flow value is positive, then flow is from
segment IQ(K) to JQ(K).
IQ(K) = upstream segment.
JQ(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.
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
Each card in this sub-package contains the flow information for four
interfaces. The number of flow specifications is equal to NOQ. The infor-
mation on each card is described below:
10 15 20 30 35 40
_ _
BQ(K) IQ(K) JQ(K) BQ(K+1) IQ(K+1) JQ(K+1)
50 55 _ 60 70 75 _ 80
BQ(K+2) IQ(K+2) JQ(K+2) BQ(K+3) IQ(K+3) JQ(K+3T
FORMAT (4(F10.0, 215 j)
BQ(K) = ratio of the flow between segments IQ(K) and JQ(K) to the
piecewise linear flow approximation.
IQ(K) = upstream segment.
JQ(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.
86
-------
Sub-Package j_I_ - Piecewise Linear Flow
The number of breaks required to describe the piecewise linear approxi-
mation 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
NOW
FORMAT (15)
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 70 80
_
QJ(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 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 J_
5 10 15
IQ(K) JQ(K) NOBRK
FORMAT (315)
IQ(K) = upstream segment, flow from segment IQ(K) to JQ(K),
assuming positive flow.
JQ(K) = downstream segment flow from segment JQ(K), assuming posi-
tive 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.
87
-------
Sub-Package JJ_
Sub-package II is a time series describing the piecewise linear approxi-
mation. Each time series element consists of two parts: a flow and a time.
The 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 cfs.
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. 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(1) 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).
Card Group £
Boundary Condition
All input is read NOSYS times; once for each system of the model.
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.
88
-------
= 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.
2. Scale Factor for Boundary Conditions
TO
SCALE
FORMAT (E10.3)
SCALB = scale factor for boundary conditions. Boundary conditions
are normally expressed as milligrams per liter (mg/1), or
parts per million parts (ppm).
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) IBC(K+1) BBC(K+2)
45 55 ' 60 70 75
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 ^
Each card in this sub-package contains the boundary condition informa-
tion 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)
89
-------
45 55 60 70 75
IBC(K+2) BBC(K+3) IBC(K+3) BBC(K+4) IBC(K+4l
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_ - Piecewise Linear Boundary Condition Approximation
The number of breaks required to describe the piecewise linear boundary
condition approximation is followed by a time series describing the boundary
approximation. Each time series element consists of two parts; boundary
concentration, and a time. The input is as follows:
5
NOW
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 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
IBCIK) NOBRK(K)
FORMAT (215)
IBC(K) = boundary segment number.
90
-------
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_ - Piecewise Linear Boundary Condition Approximation
The segment number and the number of breaks required to describe the
broken line approximation is followed by a time series describing the broken
line approximation. Each time series element consists of two parts: a
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 approxi-
mations. The input is as follows:
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 boundary approximation at time T(K) in mg/1.
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 equation to T(NOBRK). All break
times must agree for all segment, i.e., T(l) must be the
same for all exchanges, T(2) must be the same for all ex-
changes, etc.
Card Group £_
Forcing Functions
All input is read NOSYS times, once for each system of the model.
1. Data Input Option Number; Number of Forcing Functions
5 ]0
lUKOP(ISYS) NOWKdSYSy
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 I SYS. NOTE:
forcing functions may also be considered as sources (loads)
Qr sinks of a water quality constituent. If no forcing
91
-------
functions are to be inputted, set NOWK(ISYS) to zero, and
continue with next system or go to next card group.
2. Scale Factor for Forcing Functions
_ TO
SCALW
FORMAT (E10.3)
SCALW = scale factor for forcing functions. Forcing functions are
normally read as pounds per day.
3. Forcing Functions
The data input format is determined by the option selected.
Option 1
10 15 25 30 40
BWK(K)
45
IWK(K+2)
IWK(K) BWK(K+1)
55
BWK(K+3)
60
IWK(K+3)
IWK(K+1) BWK(K+2)
70
BWK(K+4) IW
75
K(K+4)
FORMAT (5(F10.0, 15))
BWK(K) = forcing function of segment IWK(K), in pounds/day.
IWK(K) = segment number to which forcing function BWK(K) is to be
applied.
Option 2
The card package consists of two sub-packages. Sub-package I contains
the forcing function data, while sub-package II contains a detailed specifi-
cation of the piecewise linear approximation to which all the forcing func-
tions are proportional.
Sub -Pack age J_
Each card in this sub-package contains the forcing function information
for five segments. The number of specifications is equal to NOWK. The in-
formation on each card is described below:
10 15 25 30 40
BWK(K) IwKTK) BHKlK+1) IWK(K+1) BWK(K+2)
45 55 60 70 75
IWK(K+2) BWK(K+3) IwX(K+3) BWK(K+4) IWKlK+47
FORMAT (5(F10.0, 15))
BWK(K) = ratio of the forcing function for segment IWK(K) to the
piecewise linear approximation.
92
-------
IWK(K) = segment number to which forcing function BWK(K) is to be
applied.
Sub-Package n_ - Piecewise Linear Forcing Function Approximation
The number of breaks required to describe the piecewise linear forcing
function approximation is followed by a time series describing the forcing
function. Each time series element consists of two parts; a function value
and a time. The input is as follows:
_ 5
NOBRK
FORMAT (15)
NOBRK = number of values and times used to describe the piecewise
linear approximation.
10 20 30 40 50 60 70 80
_ _
wKT(K) T(K) wXT(K+l) T(K+1) MKT(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/day.
T(K) = time in days, if the length of the simulation exceeds
T(NOBRK), the forcing function approximation is repeated,
starting at T(l), i.e., the approximation is assumed to be
periodic, with period equal to T(NOBRK).
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:
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 ll_ - Piecewise Linear Forcing Function Approximation
The segment number and the number of breaks required to describe the
piecewise linear forcing function approximation is followed by a time
series, describing the forcing function. Each time series element consists
of two parts: a function value and a time. The input is as follows:
93
-------
10 20 30 40 50 60 70 80
MKT(K) T(K) WKT(K+1) T(K) HKT(K+2) T(K+2) WKT(K+3) T(K+3f
FORMAT (8F10.0)
WKT(K) = value of the forcing function at time T(K), in pounds/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 boundary
conditions, T(2) must be the same for all boundary condi-
tions, etc.
Card Group £
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.
1. Number of Parameters
NO PAN
FORMAT (15)
NOPAM = number of parameters required by the model. If no para-
meters are to be inputted, set NOPAM to zero and go to
card group H.
2. Scale Factors for Parameters
10
SCALP(l)
20
SCALP(2)
30
SCALP(3)
SCALP(NOPAM)
FORMAT (BE 10.3)
SCALP(K) = scale factor for parameter group K.
3. Segment Parameters
5 15
ANAME(K) PARAM(ISEG.K)
30 35
PARAM(ISEG,K+1) ANAME(K+2)
50 60 65
ANAME(K+3) PARAM(ISEG,K+3) ANAME(K+4)
20
ANAME(K+1)
45
PARAM(ISEG,K+2)
75
PARAM(ISEG,K+4)
FORMAT (5(A5, F10.0))
94
-------
ANAME(K) = an optional one to five alpha-numeric character descriptive
name for parameter PARAM(ISEG,K).
PARAM(ISEG,K) = the value of parameter ANAME(K) in segment ISE6.
Card Group H.
Constants
The definition of the constants will vary, depending upon the structure
and kinetics of the systems comprising each model.
1. Number of Constants
5
NCONS
FORMAT (15)
NCONS = number of constants required by the model.
If no constants are to be inputted, set NCONS equal to zero and con-
tinue with the Card Group I.
2. Constants
5 15 20 30 35
ANAME(K) CONST(K) ANAME(K+1) CONST(K+1) ANAME(K+2)
45 50 60 65 75
CONST(K+2) ANAME(K+3) CONST(K+3) ANAME(K+4) CONST(K+4)
FORMAT (5(A5, F10.0))
ANAME(K) = an optional one to five alpha-numeric character descriptive
name for constant CONST(K).
CONST(K) = the value of constant ANAME(K).
Card Group _!_
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.
1. Number of Time Functions
5
NFUNC
FORMAT (15)
NFUNC = number of time functions required by the model.
95
-------
If no time functions are to be inputted, set NFUNC equal to zero, and go
to Card Group K.
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.
ANAME(K)
NOBRK(K)
VALT(K)
TOO
Sub-Package I
10
ANAME(K) NOBRK(K)
FORMAT (A5, 15)
= an optional one to five alpha-numeric character descriptive
name for the time function K.
number of breaks used to describe the time function K.
10
20
Sub-Package _!!_
30 40
VALT(K) T(K) VALT(K+1) T(K+TT
50
60
70
80
VALT(K+2) T(K+2) VALT(K+3) T(K+3)
FORMAT (8F10.0)
value of the function at time 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.
Initial Conditions
The initial conditions are the segment concentrations for the state
variables at time zero (or the start of the simulation).
1. Initial Conditions
5
ANAME(K)
15
C(ISYS.K)
20
ANAME(K+1)
30
C(ISYS,K.+ 1)
96
-------
35 45 50 60
ANAME(K+2) C(ISYS,K+2) ANAME(K+3) C(ISYS,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 appro-
priate 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.
Card Group K
Stability and Accuracy Criteria
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 concen-
tration (normal units mg/1 or ppm) for system K which if
exceeded by any segments in system K indicates that the
numerical integration procedure has become unstable. If un-
stability occurs an appropriate message is printed and the
integration procedure is terminated and a call is made to
the display subroutines.
2. Accuracy Criteria
10 20 30 80
CMIN(K) CMIN(K+1) CMIN(K+2) CMIN(NOSYS)
FORMAT (8F10.0)
CMIN(K) = originally WASP read the accuracy criteria for system K;
i.e., for time variable simulations the minimum concentra-
tion (normal units mg/1 or ppm) that governs integration
step-size control, if the user chooses option 1 (described
later). Since it is not recommended that the user utilize
this option due to numerical difficulties, just set CMIN =
0.0 for each system.
97
-------
Card Group j.
Intermediate Print Control
1. Print Interval
K)
PRINT
FORMAT (F10.0)
PRINT = print interval in days for time-variable applications.
NOTE: The maximum number of print outs = total prototype
time/print interval + 1 (for time zero) must be equal to or
less than 41.
2. Compartments (System - Segment) to be Displayed
ISYS(l) ISEG(l) ISYST2) 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 miscel-
laneous calculations may be displayed at the end of the
simulation, see Card Group N.
Card Group M
Integration Control Information
1. Integration Option - Negative Solution Option
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.
= 1, user will permit the integration procedure to go nega-
tive - used for special applications (ex., DO deficit, pH -
alkalinity).
98
-------
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).
3. Number of Integration Step Sizes
NOSTEP
FORMAT (15)
NOSTEP = number of integration step sizes to be used in the simula-
tion.
4. Integration Step Size History
10
DT(K)
20
TIME(K)
30
DT(K+1)
40
TIME(K+1) ...
FORMAT (8F10.0]
DT(K) = integration step size (normal units-days).
TIME(K) = time until which step size DT(K) will be used, then
switching to DT(K+1) until TIME(K+1).
Card Group JY
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 is 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).
99
-------
1. Variable Names
8 16 80
ANAME(I) 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
subroutine. NOTE: As presently written WASP permits a
maximum of 10 variables per segment to be saved in each
system for the IBM 1130, DSC META4, IBM 370 versions and a
maximum of 8 variables for the DEC PDR versions.
2. Variable Number, Segment Numbers
3 6 9 27
VARNO SEG(K) SE6(K+1) .... SEG(K+7)
FORMAT (913)
VARNO = the position of the desired variables, to be displayed, in
the WRITE file statement in the kinetic subroutine (see
previous note).
SEG(K) = segment number to be displayed. NOTE: Order of display
unimportant, i.e., need not be sequential.
A blank card terminates display for system, ISYS. Then another Vairable
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.
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 since 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.
100
-------
3b. Plotting Scales
10 20
PMIN PMAX
FORMAT (2F10.0)
PMIN, PMAX = minimum and maximum values, respectively, to be used for
this plot.
3c. Segment to be Plotted
36 9
SEG(1) SEG(2) .... SEG(INSPLT)
FORMAT (513)
SEG(K) = segment numbers to be plotted (a maximum of 5 segments/-
plot allowed).
A blank card terminates the plotting for system, ISYS.
At this point, unless the user has written additional display routines
(such as the POP graphics), WASP input is finished, and the user should end
his deck with the appropriate end of data set indicator.
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 TO
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+7l
FORMAT (8(15, F5.0))
101
-------
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
MXTIM
IVAR
FORMAT (215, 2F5.0, 215, 40A1)
= number of time selections to be included on this plot (maximum
of 5).
= 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 preceeding 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 sug-
gest a maximum of three; YSTR and YSTP values should be com-
patible for overlaid plots.)
= title for plot; when overlaying plots, the first two titles
and the last title will be printed.
TITL1
4d. Time Selections and Characters for Predicted Variable Plot
5 10 15 20
TIM(l) TIM(2) TIM(3) TIM(4)
27 28 29
SYMTAB(2) SYMTAB(3) SYMTAB(4)
rif
25 26
4(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).
102
-------
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 (de-
fault = 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 preceeding
plot, if OVRLAY on the preceeding 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+1) VALUE(K+lT
50 60 70 80
VALUE(K+3)
RIVMIL(K+2) VALUE(K+2) RIVMIL(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.
103
-------
4g. Format Specification for Data on Auxiliary Input File "IUNIT"
J 80
FMT
FORMAT (20A4)
FMT = 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 for-
mats, such as (2F5.0), (16F5.0), or (F5.0/F5.0)
5. Pen Plot Display Cards
The following cards permit the user to obtain plots of WASP results on a
digital pen plotter, should he have access to one on his computer system.
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
FORMAT (315, A14)
SYS = system number of the desired variable to be plotted.
SEG = segment number of the desired variable to be plotted.
IVAR = 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 ot
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, has 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 infor-
104
-------
mation: a survey number, the time (relative to time zero of the WASP simu-
lation) at which the data was collected, the value (mean value if several
samples were collected and are to be aggregated together) 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) SURVYT(3) Y(3) SD(3)
FORMAT (15, 3(F5.0, 2F10.0))
SURVEY = survey number that the data was collected for. Must be a
number from 1 to 5. In actuality SURVEY selects the appro-
priate plotting symbol for plotting the field data.
SURVYT(I) = time, relative to time zero of the WASP simulation, that the
data was 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 to-
gether than Y(I) is the aggregated sample mean.
SD(I) = the standard deviation of the aggregated sample mean (if
applicable).
The following description, used together with Table 6, should demon-
strate to the user how to prepare data input for the plots. The user will
be plotting phytoplankton and s.ecchi depth results for segment 1 through 3
of his model. The kinetic subroutine was written such that the phytoplank-
ton were written 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 is 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, while
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.
105
-------
LO LO
o o
iO LO
o o
<*
o
If) CO
oo
01
00
LO
in
o
r^ oo o o
CM CM OO CO
LO
oo
CSJ
LO
*d- Is*
O O
LO «3-
O O
O O
CL.
X
o
_J
Q.
UJ
d.
CO
0)
.o
T3
to
o
011 LO cr> r^ LO
Or CO LO f^. CO LO O~> LO LO
.. ..... . ..
i CM OCMi OCM i «:J-OO
COCO COCOLOOCM CM LOO
^" co ^* co ^~ LO oo co ^r LO
CM CM CM CO r r CO OO
i CO LO LO CO
r^ i ^- oo cfi ^* oo oo
o o o' o o d o
I I I D£ CC O£.
"*s. > UJ UJ UJ
I tO 10 I I I
iy p-x LO ^* LO ^. co i' i 111 111
OCMCMOI tDOi cr>t^.Lor-.*d- CDCM ^ o SOCMI^ SO
i .... ^^ ....... ^3 . ^*^. . *^ ... ^»^ .
LO r CO O LOi CMOCMO' r O OO'TCO O
r- «d- «* **
O CMO CO i OCMOO COO
* *
OLOOi OLOOi OOO O O Oi O O
i OOi i OOi r CM LO LO CMCVJ
i CO i CO CM OO OO
iii i r r- r CM CM OO r CO CM CM CSJ CM CM
^VJ CO ^^ *-O ^O ^** 00 O^ ^^ ^^ ^^1 CO ^tf" LO LD ^** OO O^ ^^ ^^ C\J CO ^d~ LO ^D ^*> CO ^T^ c^3
i i i i i r i i i i CMCMCMCMCMCMCMCVJCMCMOO
106
-------
SECTION 7
MVP THEORY
INTRODUCTION
The task of a modeler in calibrating and verifying a model is compli-
cated by a number of factors. Often two of these factors are the size
(number of compartments) and the quantity of data (either lack of or wealth
of). The eutrophication analysis performed on the Western Delta-Suisun Bay
was performed using an 11 system-39 segment model. The 429 linear and non-
linear equations were numerically integrated for a 12 month period. Using a
print interval of 10 days, some 15,000 numbers were generated for each simu-
lation run. The output generated from these simulations was no small task
for the modeler to assimilate. The limitations caused by a lack of data in
a model calibration are obvious, but what a wealth of data could mean. For
example, consider that as a result of the International Field Year on the
Great Lakes (IFYGL) some 200,000 observations encompassing some 75 water
quality parameters were made on Lake Ontario alone. The modeler attempting
to understand the behavior of such large scale models would find an almost
impossible task. Even with the most sophisticated computer graphics package
the user can only view certain portions of the model (either in variable,
temporal or spatial planes). Furthermore, since portions of the model are
so interactive, "adjustments" to improve one portion of the model often ad-
versely effect other portions of the model.
An important consideration in a model calibration and verification
analysis is the degree to which a model computes a "reasonable" representa-
tion of the real world. It is at this point that a considerable degree of
confusion and difference exists both in the realm of the model builder and
in the mind of the decision maker. What is "reasonable"? Is it sufficient
for a model to generate profiles which "look" like what is being observed?
Is it sufficient, for example, for a phytoplankton model to simply generate
a spring pulse which has been observed or is it necessary to ensure that the
magnitude of the pulse is correct using some quantitative measure?
The Model Verification Program provides a statistical framework (1) in
which to answer some of these questions in a quantitative way and to enable
the user to determine a model's relative behavior and its "reasonableness"
in representing the real world. Further, MVP, as will be discussed in
Section 8, provides a means for determining the temporal and spatial scale
over which a model may be applicable. That is, does the modeling framework
chosen by the modeler apply equally well over all temporal scales from hour-
to-hour, to day-to-day, to year-to-year, and/or all spatial scales from
local near-shore or tidal flats up to open-lake or estuary scales?
107
-------
STATISTICAL THEORY
Assuming that the user has a sufficient water quality data base and can
generate mean and variance statistics over segments and over time (daily,
weekly, monthly, or yearly, depending upon the time scale of the problem)
for the variables of interest, MVP can perform three simple statistical
tests to compare model output to observed field data. These tests are:
1. Test of the difference of means, using a Student's "t" test.
2. A linear regression analysis between observed data and model
output.
3. An evaluation of the model using relative error to provide a
gross measure of model adequacy.
1. Student's "t"
Let jT-jjk = the observed mean for variable i, segment j, time span or
time average (day, week, etc.) k, and c"-jj|< = the comparable model com-
puted mean. Permit d = cf-x", the difference of the calculated and observed
means where the triple subscript has been dropped, to be distributed as a
Student's "t" probability density function. The test statistic is given by
(2)
t = _« (5.1)
where
6 is the difference between the model and the data,
" is the standard deviation given by the pooled variance or
2 x
s/ = -T?- (5.2)
for sx^ as the data variance for the specific variable, segment, and time
averaged period.
Under the null hypothesis: 6=0, there is a critical d which deline-
ates the region of rejection of the hypothesis and is given by
d~c = ± tsj (5.3)
which for a 95% confidence region (5% change of making a Type 1 error) is
given by (5.4)
108
-------
(5.4)
The distribution of d and the critical regions are shown in Figure 30. As
4- MODEL MEAN-OBSERVED MEAN
TISTICALLV SKMiriCAM
OIMC*CNC[ UTWCCN
MODEL ««0 OMCHVCO
C*NS IM% CONFIDENCE
LIMITS)
VERIFIED
V SCO RE = 0
OVER-
ESTIMATION
Figure 30. Determination of verification score V.
indicated if (-dc d
(5.5)
Vijk - dc - dijk
(5.6)
The V score is therefore a measure of the degree to which the model de-
viates from the observed data, given the spatial and temporal variability
within a segment and timespan. More precisely, V is the extent to which the
analysis has penetrated into the region of rejection of the null hypothesis.
Of course, all the caveats of the application of such statistics apply,
especially the change (unknown) of making Type 2 error (that null hypothesis
is not rejected when it should be).
Another simple measure that may be used is the number of segments in a
given time span that have a V score of zero. Therefore let
109
-------
K
ijk
vijk =
A score, defined as the S score, for variable i and time span k is given by
Sik = Kijk/n
where n is the total number of segments where a V score can be computed,
either for the entire water body or just for regions of the water body (as
for example tidal flats vs. deep water channels for an estuary or near shore
vs. open water for a lake). The score then simply represents the fraction
(or percent) of the segments that "passed" the verification test. Since a
number of variables may be scored using this verification analysis, an over-
all aggregated score may be of interest and can also be computed. However,
equal verification of all variable may not be of concern. For example, in
an eutrophication model one may be willing to accept a lack of verification
of ammonia nitrogen in the water body, but be particularly concerned about
say, total phosphorus and chlorophyll. Therefore, a series of weights,
w-j, can be assigned to each variable i, representing the relative impor-
tance of each variable within the analysis. The aggregated S score for a
time span k is then given by
r n
S. = Z Z w. K.
* i j ' n
r
n Z w.)
(5.7)
where r is the number of variables in the aggregated score. S|< therefore
represents the weighted fraction of the total number of segment variables
that "passed" a Student's "t" test (V=0) for the time span k. It should be
noted that not all segments and variables can be tested every time span, so
that r and n are functions of the data availability for time span (daily,
weekly, etc.) k.
As an example of the utility of this scoring framework consider Figure
12. A comparison of verification scores is made for two different eutrophi-
cation kinetic structures that were used in the Lake Ontario analyses (3).
The three pairs of bar-histograms reflect the verification scores as a func-
tion of the standard deviation or standard error of the mean. The smaller
the estimated standard error of the mean the more stringent are the require-
ments for "passing" the "t" test. It should be noted that the minimum im-
provement in verification scores, irregardless of the estimated standard
error, was 40%, due solely to a change in the k-inetic structure. For the
readers information the principal kinetic changes between the LAKE1 and
LAKE 1A models were as follows:
a. The addition of silica as a state variable and the splitting of
chlorophyll into diatoms and non-diatoms.
110
-------
b. Use of threshold nutrient limitation in contrast to product ex-
pressions. The growth rate is limited by
Mln (*'.. '"I
KsN + IN1> KsSi
rather then
[P°4-P] [N] [Si]
KsP + [P04 " PJ KSN + LNJ K + IS1J
where Ksp, Ksfj, Ks$i are the half-saturation constants for
phosphorus [P04-P], nitrogen [N], and silica [Si], respectively.
The rate of mineralization of unavailable to available forms of the
nutrients, [PCty-P], [N], [Si], for uptake by the phytoplankton
is dependent, through a Michaelis Menton recycle expression, upon
chlorophyll. Therefore, the general expression for conversion of
unavailable forms is
R
[Unavailable] > [Available]
for R - KB1'20
f°rR K6
where
K = mineralization rate @ 20%.
K = half -saturation constant for chlorophyll [chl-a].
2. Linear Regression Analyses
An alternate perspective on the adequacy of a model can be obtained by
regressing the calculated values with observed values. Therefore, let the
testing equation be
x" = ct + BC" + e (5.8)
where a and 3 are the true intercept and slope respectively between the cal-
culated and observed values and_e is the error in "x. The model Equation
(5.8) assumes, of course, that c, the calculated value from the user's model
is known with certainty which is not the actual case. With Equation (5.8),
standard linear regression statistics can be computed, including
a. The square of the correlation coefficient, r^, (the % variance
accounted for) between calculated and observed.
Ill
-------
b. Standard error of estimate, representing the residual error between
model and data.
c. Slope estimate, b of 8 and intercept estimate, a of a.
d. Tests of significance on the slope and intercept.
In this work, the null hypothesis on the slope and intercept is given by
6 = 1 and a = 0.
Therefore, the test statistics
j^J and a/s
sb a
are distributed as student's "t" with n-2 degrees of freedom (4). The vari-
ance of the slope and intercept, s^ and Sa^ are computed according to
standard formulae. A two-tailed "t" test is conducted on b and a, sepa-
rately, with a 5% probability in each tail, i.e., a critical value of t of
about 2 provides the rejection limit of the null hypothesis.
Regressing the calculated and observed values can result in several
situations. Figure 31 (b) and (c) shows that very good correlation may be
obtained but a constant fractional bias may exist (bl); also Figure
31 (d) indicates the case of good correlation but for a > 0 a constant bias
may exist. Evaluation of r?, b and a together with the residual standard
error of estimate can provide an additional level of insight into the com-
parison between model and data.
3. Relative Error
An additional simple statistical comparison is given by the relative
error defined as
e =
x - c
(5.9)
for each variable, segment or time span. Various aggregations of this error
across regions can also be calculated and the cumulative frequency of error
over time spans or segments can also be computed. The difficulties with
this statistic are its relatively poor behavior at low values of x~ and the
fact that it does not recognize the variability in the data. In addition,
the statistic is poor when x~ > c" since under that condition the maximum re-
lative error is 100%. As a result, the distribution of this error statistic
is most poorly behaved at the upper tail. Nevertheless, if the median error
is considered, this statistic is the most easily understood comparison and
provides a gross measure of model adequacy. It can also be especially use-
ful in comparisons between models.
112
-------
10
o
UJ
h- 8
<
0 6
^J
<
0 4
2
/»
CALCULATED (fl)
OBSERVED
g ~- * x^
- s * r
i
TIME (day*) '
/
X|'
' '
/ r* = 1
/ b = I
/ 0 = 0
/ 1 1 1 1
10
O
Ul
< «
3 6
<
<-> 4
2
y-v
r«. i (b)
11 ^ *
J / 1 0« 0 /
» / *\ ,/
o J \ -
_ \ /]|
' . *. ^ 1
" / '
/
- /
/
~ x ^^^^^^^
' J--<***"^"^
^-*r*"*i^ 1 I I
02468 10 02468 10
OBSERVED OBSERVED
10
0
UJ
1- 8
-J
3 6
O
g 4
2
0
rz, , (c)
-0=0 / j
^^
/ f
_ X
/.
/ s *
/ / 8 ,/N!
/' TIME(doyt)
r i I i I
10
a
K 8
_J
^ 6
0
g 4
2
O
rt=| (d)
b * i /
- a > o 0' '
- / 4
/ s
/
-/ / o
X § .yv
^ TIME(doyc)
/ \ \ 1
0 2 4 6 8 10
OBSERVED
4 6 8 10
OBSERVED
Figure 31. Possible cases in regression between
calculated and observed values.
113
-------
Figure 11 presents a comparison of the relative error analyses for the
primary state variables incorporated in the LAKE1 and LAKE1A modeling frame-
works. Again, this type of analyses demonstrates the improvement in model
verification using a revised kinetic structure.
REFERENCES
1. Thomann, R.V. and R.P. Winfield. 1976. On the Verification of a Three-
Dimensional Phytoplankton Model of Lake Ontario. Environmental Modeling
and Simulation, ORD and 0PM, U.S. EPA.
2. Wine, R.L. 1964. Statistics for Scientists and Engineers. Prentice-
Hill, Inc., Englewood Cliffs, N.J. 671 pp.
3. Thomann, R.V., et _al_. 1978. Verification Analysis of Lake Ontario and
Rochester Embayment Three Dimensional Eutrophication Model, in manu-
script.
4. Carnahan, B., ejt
-------
SECTION 8
MVP PROGRAM LOGIC
INTRODUCTION
MVP was designed and written to be a direct compliment to the WASP pro-
gram, so that a user would have a statistical framework for judging the ade-
quacy of a model developed using the WASP program. As such the same design
considerations that were incorporated into the WASP program, were included
in MVP. A very important design consideration, carried over from WASP, was
that the user not have to make any programming changes to MVP no matter how
he formulated his model within the overall WASP framework.
However, the most important design philosophy included in MVP was to
provide the user with the greatest flexibility in analyzing the model veri-
fication over different temporal and spatial scales. Assuming a sufficient
water quality data base is available for the analyses, a user may look at a
model's verification or ability to reproduce observed data from week to week
up to seasonal time scale or from a segment by segment near field scale up
to a whole lake, river, or estuary scale within the same MVP run. This pro-
vides both the model builder and the decision maker with useful information.
For the model builder it will quickly inform him of areas or time periods
where additional modeling effort and/or data collection is required. For
the decision maker, the information provided might enable him to formulate
waste load allocation strategies that reflect the spatial and temporal
scales verified by the model. That is, waste load allocations might be
formulated on a regional or seasonal basis, if the model could only be veri-
fied over large spatial scales or long term time scales.
The purpose of this section is to review the program logic of MVP. How-
ever, unlike the WASP program, which requires the model builder (if not the
user as well) to understand the program structure as well as the theory
upon which WASP is based, it is sufficient for the model builder and user to
understand only the statistical theory which MVP utilizes. In addition,
sufficient flexibility has been incorporated into the program so as to make
program modification virtually unnecessary. For these reasons the level to
which the program logic is detailed will be limited to a rather simple over-
view. 'Instead greater emphasis will be placed on providing an understanding
of the types of temporal, spatial, and variable aggregation schemes avail-
able through user selectable options within the program.
115
-------
MVP MAINLINE AND SUBROUTINE OVERVIEW
The MVP mainline is really just a program control module and as such it
performs no computations but just consists of a calling sequence to various
subroutines in the MVP package.
In the following subroutine descriptions the references, within paren-
theses, to Card Groups are with respect to the card by card description of
the MVP input data to be found in Section 9, MVP Input Structure.
MVP01
MVP01 performs various initialization procedures and reads the model
configuration and user selectable options for all of Card Group A.
MVP2A
MVP2A is the heart of the model verification analysis procedure. MVP2A
controls the calling sequence for the various subroutines which read the
field data and the model theory, as computed by WASP, and perform the
various statistical analyses on both. A brief flow chart is presented in
Figure 32.
MVP02
MVP02 reads the. variable name card, score aggregate weighting card(s)
and field data to be used in the model verification statistical procedures
(Card Groups BB through BD). As the field data is read it is transferred to
a workfile for later retrieval when needed by the statistical analyses sub-
routines MVP03, MVP04, and MVP05. When necessary field data are aggregated
over time, if two or more data entries would fall within the same time in-
terval .
MVP03-RESCOR
The linear regression analyses on field observations and the WASP
model's theoretical results, are performed by subroutines MVP03 and RSCOR.
MVP03 permits the user to select a number of ways to aggregate field obser-
vation and model theory over space and time before performing each linear
regression analysis. An option to provide graphical display of the re-
gression analysis and the confidence intervals around the regression line is
provided. All options for MVP03 are selected by the appropriate cards read
from Card Group BF.
MVP04
MVP04 performs the relative error analyses. The percent relative error
between the WASP model theory and the field observations is computed and if
desired the user may obtain distributive and cumulative relative error
histogram plots. The user may, via the input options available through Card
Group BG, aggregate field observations and model theory over space and time.
116
-------
ISYS =
ISYS-H
PRINT
AGGREGATE
SCORES
ISYS =
ISYS-H
CALL MVP02
TO READ
FIELD DATA
Figure 32. Simplified MVP2A flow chart.
117
-------
MVP05-TSCOR
MVP05 and its associated subroutine RSCOR compute the Student's "t" test
scores on the differences in means between field data and theory results.
As with the linear regression analyses and relative error analyses, the user
is permitted (via Card Group BE) to aggregate field data and WASP theory
over space and time before performing the Student's "t" test analyses.
LOADR
LOADR performs the task of loading the WASP theory results into core
memory from the WASP disk storage files. Theory for the current system
variable combination is loaded into main memory for all segments and for all
times (every print interval specified in the WASP input) so as to minimize
disk storage retrievals. LOADR is the one subroutine in the MVP package
which is the most sensitive to disk input/output formatting of different
computer software operating systems (as is the WASPB kinetic subroutine disk
output structure), and so the system analyst must ensure compatibility
between LOADR and WASPB.
Miscellaneous Subroutines
The remaining subroutines in the MVP package (MVPHR, COAD2, SETA, and
SETRA) perform various heading, printer plot, and variable array initiali-
zation procedures.
MVP COMMON
The following list defines the relevant variables of importance con-
tained in blank COMMON of the MVP package. Blank COMMON is used by MVP as
the vehicle to pass information from subroutine to subroutine within the
program. The R, I, and * contained within parentheses after the variable
name indicate, respectively, whether the variable is a REAL (floating
point), or INTEGER (fixed point), and whether the variable is read as input
data.
Variable Name Definition
SCORE(R)* Used to compute the aggregate "t" test scores across
systems. SCORE (J,l) contains the number of "t" tests
passed, and SCORE (J,2) contains the number of "t"
tests performed. The final aggregate score is defined
as FSCORE(J) = 100. * SCORE(J,1)/SCORE(J,2).
LSGLG(I*) Used to indicate which segment or segments are to be
included in which aggregate score. Rather than use an
entire word to indicate whether or not segment ISEG is
to be included in aggregate score J, a single bit (bit
J of LSGLG(ISEG) is set to 1 if segment ISEG is to be
included, and set to 0 if segment ISEG is not to be
included.
118
-------
INTMT(I*)
VARWT(R*)
TMTRX(R)
TIMEX(R)
IBEG(I*), IEND(I*)
TMULT(I*)
DT(R)
NOINT(I)
HEADR(R*)
NOSYS(I*)
NOSEG(I*)
MXTIM(R*)
MXDMP(I*)
NOAGG(I*)
Used to indicate if the aggregated "t" test score is
to be computed over all time intervals (INTMT = 0) or
if the "t" test score is to be computed for a specific
time interval.
VARWT(J,K) are the weights to be assigned each vari-
able (K) for each scoring aggregate (J). This permits
the user to give more weight (or importance) to the
"t" test score of particular water quality variable
than another, when computing an aggregate score across
state variables.
Used to hold in core the WASP model theory results for
a particular water quality variable for all segments
over all time intervals. TMTRX is loaded with the
model results by subroutine LOADR.
Contains the times at which the WASP model theory re-
sults were saved on the auxiliary disk files.
Contain the beginning and ending time intervals over
which the verification analysis is to be performed.
Specifies the number of WASP print intervals that are
to be aggregated together to form one WASP theory data
point. For example if the user had chosen a print in-
terval of 10 days for his WASP simulation and wishes
to perform the MVP analysis using a monthly averaging
interval TMULT would be read as 3.
Contains the print interval used in the WASP simula-
tion.
Contans the number of aggregated WASP time intervals
that the MVP analysis will be performed on.
Contains the user selected heading or title for the
MVP run.
The number of systems or state-variables in the user's
WASP model.
The number of segments in the user's WASP model.
Is the ending time of the WASP simulation.
Is the blocking factor or the maximum number of vari-
ables at each print interval. MXDMP is determined
from the user's WASPB kinetic subroutine.
The number of segment aggregations that the user
wishes to compute.
119
-------
NOWGT(I*)
NOVAR(I*)
OUT(I)
SOPT(I*)
IOPT(I*)
SC(I*)
ISYS(I)
LSTVR(I)
IDF(I)
SYSBY(I*)
IRCRD
IDF40(I),
IDF41(I)
LO(R), HI(R)
IPTR(I)
The number of segment weighting vectors specified by
the user for spatial aggregation.
The current variable number in system ISYS that the
user is performing an analysis on,.
Device number for reading input data.
Device number for printer output.
Contains the variance option for the WASP theory as
chosen by the user.
Contains the degree of freedom option for the
Student's T and Snecedor's F distribution parameters
as chosen by the user.
Used in conjunction with the SORT variance option.
System currently undergoing verification analysis.
Used as an internal flag to control the calling of
subroutine LOADR. If NOVAR is equal to LSTVR then it
is not necessary to call LOADR since the necessary
WASP theory is already in core. If NOVAR is not equal
to LSTVR then a call to LOADR is made, to obtain the
WASP results for variable NOVAR and then LSTVR is set
equal to NOVAR.
Used in the DEC-POP version of MVP as the record
address pointers for the direct access WASP dump
files. Not utilized in the IBM 370 version since
sequential files are used.
User selected system bypass indicators. If a user
wishes he may choose to bypass the verification
analysis for a particular system (or systems) ISYS, by
setting SYSBY(ISYS) appropriately.
Not currently used.
Internal record address pointers for the direct access
files utilized by MVP.
Used to store the confidence intervals of the linear
regression analysis.
Used as an internal record address pointer for lo-
cating field data for each segment (similar in concept
to an ISAM table).
120
-------
JPTR(I) Used as an interval record address pointer for re-
trieving WASP theory results.
WGT(R*) Stores the weighting factors to be used in aggregating
across state variables for the "t" test scores.
MVP ANALYSIS OVERVIEW AGGREGATION
As stated earlier, MVP permits the user to analyze his model's verifica-
tion under different temporal and spatial scales, as well as across state
variables. It is the purpose of this section to outline the various reasons
why a modeler may choose to perform these aggregations, and the various ways
that these aggregations may be formulated using MVP. This second point is
especially important since, due to the flexibility built into MVP, there are
a number of options and aggregation strategies which the user may choose
from and unless the user understands what is available to him he may be con-
fused by the MVP input structure.
Segment Weights
If the user wishes to formulate a spatial aggregation, i.e., compute a
mean for a state variable or water quality variable that is representative
of a number of segments, he must decide how to weight the individual segment
concentrations in forming the overall area average. To do this MVP permits
a maximum of three segment weighting vectors (Card Group AE). Normally a
volume weighted average is what is desired, but other schemes may be in-
cluded by the user. For example, in a multi-layer model the user may wish
to obtain the average light extinction coefficient. In this case the user
may choose to use a weighting scheme which is a function of segment depth,
rather than volume average.
State Variable Aggregating
A user may determine how well his model verified against all state
variables and parameters of interest using the Student's "t" test on the
difference between theory and field means as the scoring criteria. The
user, via Card Group AG, specifies the segment and time aggregations for
which he will compute "t" scores accross variables. Then utilizing Card
Group BC (read for each non-by-passed system), he determines what state-
variable and/or parameters are to be included in the aggregate score, and
what weight each is to be given in forming the overall score.
Spatial Aggregation
Perhaps the best way to demonstrate the flexibility in spatial aggrega-
tion is to refer back to Figure 1, the segmentation map for the Western
Delta-Suisun Bay eutrophication study. Segments 1,2,3,4,5,6,14,24,25,26,
27,28, and 36 comprise the main stem of the Sacramento River (channelized
for shipping), while 7,8,9, and 11 are the main stem segments of the San
Joaquin River. The remaining segments are either tidal flats (regions of
high algal productivity) or turbid, low flow sloughs. The user may perform
121
-------
his verification analysis on a near field or segment by segment basis for
those segments for which is available. Going up the spatial scale he might
analyze: the upper Sacramento River by aggregating segments 1 through 6,
the lower Sacramento River by aggregating segments 14,24,25,26,27,28, and
36, and finally the entire Sacramento River by including both of the afore-
mentioned groupings. A similar procedure could be followed for the tidal
flats. A comparison of verification scores for the main stems or deep
channel segments vs. tidal flats could be computed by appropriate segment
groupings. Finally, an overall score is computed by aggregating all seg-
ments together.
Control of such groupings is provided through Card Groups BE, BF, and BG
for the "t" test, linear regression, and relative error statistical analyses
respectively.
Temporal Aggregation
With the exception of the linear regression analysis which has a little
more flexibility, the user has a relatively small degree of control over
temporal aggregation, in that he may aggregate over sequential time periods
only. For example, assuming that a user has performed a WASP simulation for
a prototype time of one year and used a print interval of ten days, he could
perform the verification analysis on a monthly basis by specifying TMULT
(Card Group AC) to be equal to 3. A user however (with the exception of
linear regression analyses) cannot aggregate over non-sequential time
periods, eg., days 0 to 30, with days 330 to 360.
122
-------
SECTION 9
MVP INPUT STRUCTURE
INTRODUCTION
The Model Verification Program, MVP, was developed and written for use
with WASP. It provides for the modeler a statistical framework for judging
model verification improvements (or lack of) resulting from adjustments to
model parameters or changes in the overall kinetic structure. Since MVP is
to be used in conjunction with WASP, many of the concepts and variable de-
finitions used in this section are direct carry overs from the section dis-
cussions concerning WASP. Therefore it is expected that the reader have
some familiarity with these sections before attempting to use MVP.
The input data required by the MVP program is divided into two major
card groups with each major card group having seven minor card groups.
These card groups are briefly summarized in Table 7.
TABLE 7. SUMMARY OF MVP CARD GROUPS
A^MVP Model Configuration Cards
AA Title Card
AB System Identification Card
AC Print Interal Aggregate Card
AD System Bypass Option Card
AE Segment Weighting Vector(s) Card
AF Variance Option Card
AG Score Aggregate(s) Card
B. MVP System Information Cards
BA System Card
BB Variable Name Card
BC Score Aggregate Weighting Card(s)
BD Field Data Card(s)
BE T-Test Card(s)
BF Linear Regression Card(s)
BG Relative Error Card(s)
For each minor card group, a detailed card by card description is presented
so as to define the variable fields which appear within the card group and
123
-------
to inform the user of any options which may be available. Depending upon
options to be selected by the user, certain minor card groups or cards with-
in a minor card group may not need to be inputted to the MVP program. Where
this circumstance might arise, the data input manual informs the user how to
avoid inputting the unnecessary information.
MVP INPUT
Model Configuration Cards
Major Card Group A is comprised of input data necessary to inform MVP
about the configuration of the user's model and the various aggregating and
variable weighting schemes the user wishes to apply to the verification
analyses. These input data or input variable may be thought of as being
system or state variable independent. The following is a card by card
description of Major Card Group A.
AA. Title Card
J 80
LINEAR REGRESSION OF AUG. 1973, RIVER SURVEY, REACHES 1,2&3
FORMAT (20A4)
Card columns 1-80 contain any information the user feels would be help-
ful in describing the run and identifying the ouput for later reference.
AB. WASP Model Configuration
5 10 15 20 25 30
NOSYS NOSEG MXTIM MXDMP NOAGG NOWGT
FORMAT (615)
NOSYS = number of systems in the user's WASP model.
NOSEG = number of segments in the user's WASP model.
MXTIM = maximum number of timesteps (print Intervals) generated by
the user's WASP simulation run. This parameter includes
time zero if it is in the simulation. Example: If a user
performed a 360 day simulation with printout every 10 days,
then the MXTIM = (360/10) + 1 (for time equal zero) = 37.
MXDMP = maximum number of variables dumped by WASP from each system.
MXDMP is determined by the user when writing the WASPB kine-
tic subroutine. See the WASP manual for further details, if
necessary.
NOAGG = number of aggregate scroes, aggregated across segments and
WASP variables, to be computed, present maximum is 26.
124
-------
NOWGT = number of segment weighting vectors specified for spatial
aggregation; at least one must be specified and there is
presently a maximum of three permitted.
AC. Print Interval - Aggregation Card
10 20 25
TBEG TEND TMDUT
FORMAT (2F10.5, 15)
Since the dates and times of the field statistical data will normally
differ from the WASP print intervals, this card allows the user to specify
parameters to structure (aggregate) the WASP results to parallel and corre-
spond timewise to the field data.
TBEG = WASP print interval time which the user selects as the
starting time point for the MVP analyses.
TEND = WASP print interval time which the user selects as the
ending time point for the MVP execution.
TMULT = number of WASP timesteps (print intervals) to be aggre-
gated over the time axis into each MVP Time Interval. If
not specified the default is one.
Note: These three parameters must satisfy the following equa-
tion for proper MVP execution.
TEND-TBEG =K where K - 1,2,3 ...
TMULT *{nWASP PRINT INTERVAL")
Example: Given the WASP print intervals ind ays to be:
0, 10, 20, 30, 40 340, 350, 360.
The user should specify a Print Interval - Aggregate Card of:
a. 0.0 360.0 3
to aggregate MVP results into twelve monthly averages
MVP interval: 1 2 3 12
WASP time: 0+^3030+-6060+-90 330+-360
or
b. 150.0 230. 2
to aggregate MVP results into four semi-monthly summer
averages
MVP interval: 1234
WASP time: 150+-170 170+-190 190^210 210+-230
125
-------
AD. System Bypass Option Card
SYSBY(l) SYSBY(2) SYSBY(NOSYST
FORMAT (1912)
SYSBY(K) = 0, perform any and all statistical tests specified by the
user for any of the variables associated with system K.
= 1, bypass all statistical tests specified by the user for
any and all variables associated with system K.
This option permits the user to bypass statistical analyses of any
system or combination of systems without having to restructure his input
data deck, i.e., it is not necessary to remove the input data cards for a
system if it is not to be included in the analyses.
AE. Segment Weighting Vector(s) Card
The user must specify from one to three weighting vectors according to
the user's input for NOWGT. The number of elements in each weighting vector
is equal to NOSEG, the number of segments in the user's WASP model. The
user inputs the first weighting vector, normally the segment volumes, from 1
to NOSEG, eight weighting vector elements per card, using as many cards as
necessary. If more than one weighting vector is to be input, start each new
vector on a new input card.
10
20
80
WGT(1,IWT) WGT(2,IWT)
WGT(NOSEG,IWT)
FORMAT (8F10.5)
WGT(K,IWT) =
weighting vector element K, for weighting vector IWT.
runs from 1 to NOSEG. IWT runs from 1 to NOWGT.
Note: The weighting vector(s) information is utilized by MVP
when statistical tests are performed over segment aggre-
gates. Thus if: C-j = concentration of a state vari-
able in segment i, and V-j = the volume of segment i,
then for mass conservation
n
where n is the number of segments to be included in this
aggregation.
126
-------
AF. Variance Option Card
Since the WASP program does not calculate standard deviations for its
theoretically calculated means, this input card allows the user to specify
what the standard deviation of the WASP data will be. The user must also
supply an option regarding degree of freedom calculation.
8 10 15 23
SOPT SC TOPT
FORMAT (2A4, 12, 5X, 2A4)
SOPT = 'S TIMES', WASP standard deviations will be a multiple of
the field data standard deviations dependent upon the value
of SC.
= 'S EQUAL1, WASP standard deviations are assumed equal to the
field data standard deviations (i.e., SC = 1.0).
'S CONST1, WASP standard deviations will be assumed to be
equal to a constant as specified by SC.
SC = integer value specified in conjunction with SOPT with the
meaning as indicated above.
Note: According to the above information SOPT = 'S TIMES' and
SC = 1 is equivalent to SOPT = 'S EQUAL1.
TOPT = 'T CONST', Student's T and Snecedor's F distribution pa-
rameters are given an average value throughout MVP analyses
independent of the actual number of degrees of freedom.
= 'T VARYS1, Student's T and Snecedor's F distribution pa-
rameters are exact and are a function of the degrees of
freedom for small values of degrees of freedom.
Note: Specifying the average value distribution parameter
option (T CONST) allows for faster MVP execution, while
specifying the exact parameters option (T VARYS) pro-
duces more accurate output results. However, result dif-
ferences become small if the degrees of freedom involved
are greater than ten. Thus, it is recommended to use the
degrees of freedom varying option (T VARYS) only if the
degrees of freedom are consistently small.
AG. Score Aggregate(s) Card
The user is permitted to form a number of spatial and/or temporal
aggregations to see how different areas of the model (ex., surface layers
vs. bottom layer, summer vs. winter) score relative to one another utilizing
the Student's "t" test as the scoring criteria. If NOAGG was specified as
zero than this section is bypassed. If NOAGG is greater than zero, then
NOAGG input cards are specified
127
-------
12 16 20 24
TINT NUMSG NSV(l) NSV(2) NSV(NUMSG)
FORMAT (F12.6, 17I4/(20I4)
TINT = a time within a MVP time interval over which the scope
aggregation will performed. If TINT is inputted as zero,
then all available time intervals will be aggregated for
this score aggregation.
NUMSG = number of segments within this score aggregation.
NSV(K) = segments which are to be included in this score aggregate.
For a given score aggregation, the various segments should
be unique.
Major Card Group J3
MVP Field Data and Statistical Test Options
Major Card Group B consists of the statistics of the observed field
data and various statistical option cards for any or all of the variables
(state variables or otherwise) in the user's WASP model. These cards, which
are read for each system (from 1 through NOSYS) in the model include the
WASPB variable names, additional aggregate weighting vectors, observed field
data, and user determined options to perform the three available statistical
procedures. The following is a card description of Card Group B.
BA. System Card
For each system of the NOSYS systems incorporated in the WASP simula-
tion, a system card must be inputted. This card has the character string
'*SYSTEM' in columns 1-8, optionally followed by an integer denoting the
current system. WASP systems are inputted in increasing numerical order
from 1 to NOSYS. The system-card must be supplied whether or not the parti-
cular systems is bypassed. However, the system card is the only input card
which must always be supplied. All other input cards are optional under
certain circumstances.
4 10
*SYSTEM ISYS
FORMAT (2A4, 12)
BB. Variable Name Card
This input card must always be supplied unless the system bypass option
for this system is set on (i.e., SYSBY = 1)
8 16
NAMEUT NAME(2) NAME(MXDMP)
FORMAT (10A8)
128
-------
The user supplies MVP with names for the MXDMP variables of the current
system.
NAME(K) = name of the kth variable in the ISYSth system dump file.
The user should refer to the appropriate WRITE statement
in the WASPB kinetic subroutine.
BC. Score Aggregate Weighting Card(s)
For a discussion of score aggregation weighting vectors refer to the
MVP Theory section. NOAGG score aggregate weighting cards must be supplied,
here, in each non-bypassed system. If NOAGG has been inputted as zero, by-
pass this section. The first score aggregate weighting card is utilized in
computing score aggregate number one, the second card for score aggregate
two, etc., up to and including the NOAGGth card.
8 16
VARHTTTT VARHT(2) VARwT(MXDMP)
FORMAT (10F8.4)
The ith card, where lo^NOAGG, should contain MXDMP weights. If the
user does not want to include the jth variable of the current system in the
ith score aggregate, a zero should be specified for VARWT(J) on the ith
score aggregate weighting card. If the user desires equal weights for all
variables for the ith score aggregate, a 1. should be specified for all the
variables from the current system that the user wishes to include in the
overall score.
Note: Since any of the variables in any of the systems may be in-
cluded (weighted or unbiased) in the score aggregate computa-
tion, the score aggregate results are printed at the end of an
MVP run.
BD. Field Data Card(s)
The user must now supply the observed field station data to the MVP
program for the current system. For each field station, the user must
supply MVP with two types of cards: (1) a Variable-Segment Card, and (2)
one or more Data Cards.
1. Variable-Segment Card
4 8
NOVAR NUMSG
FORMAT (214)
NOVAR = The variable number within the current system for which
the following data was collected. NOVAR may range from 1
to MXDMP.
129
-------
NUMSG = The segment number within the WASP model for which the
following field data was collected. NUMSG may range from
1 to NOSEG.
Note: A blank Variable-Segment Card signals the end of field
data for the current system to the MVP programs.
2.
Data
T|
Card
10
;i) xi
20 30
;i) S(l)
40 50
N(l) T(2)
60
X(2)
70
S(2)
80
N(2)
FORMAT (2(3F10.5, 110)
T(l) = time data was collected for NOVAR variable, NVMSG segment,
for current system.
X(l) = mean of the collected data, in same units as WASP output
results.
S(l) = standard deviation of the sample data.
N(l) = number of observations in this sample.
There is only one reduction on the input of field data. Data must be
entered in increasing time order. Otherwise execution of the MVP analyses
run will be terminated. Within one system, a Variable-Segment may be re-
peated, but only the latest entry will be used by the MVP program, all pre-
vious entries with the same Variable-Segment pair within one system will be
overridden.
Note: A blank time entry (with blank associated data) signals the end
of field data for the current Variable-Segment Card. MVP will
then expect a new Variable-Segment Card immediately following
the end of the data entries for the last Variable-Segment Card.
Two blank cards in a row terminates the input of field data for
the current system.
Statistical Performance Option Cards
Minor Card Groups BE through BG initiate the performance of MVP statis-
tical tests on the inputted field data and WASP results for the system under
current analyses. Each of these card groups has a title card which MVP re-
cognizes in order to implement the appropriate statistical test. Within
each system, the user has the option of specifying any, all, or none of the
statistical tests. The only restriction is that the user may specify a
statistical test only once within each system. Otherwise the MVP run will
be aborted without further analyses being performed.
BE. T-Test Card(s)
To initiate the Student's T-Test, the user must input the character
string '*T-TEST' in columns 1-7 of the title card for this statistical
130
-------
group. Then the single format for the control cards in this group are as
follows:
4 8 12 16 20
NOVAR IWGT NUMSG NSV(l) NSV(2) NSV(NUMSG)'
FORMAT (2014)
NOVAR = variable within current system, on which T-Test is to be
performed.
IWGT = weighting vector #1, #2 or #3 that the user has previously
specified, that is to be used for any spatial aggregations.
IWGT must always be specified.
NUMSG = 1, for T-Test of one WASP segment results versus a parti-
cular field station data that was inputted in the data
section of this system.
> 1, for spatial or segmental aggregation, i.e., the number
of segments to be aggregated.
NSV(I) = segment number or numbers (segment aggregation) which will
be contained in this T-Test. The number of segments speci-
fied must be equal to NUMSG. If segmented aggregation is
specified, it is permissible for more than one (even all) of
the segments to have associated inputted field data.
To exit from the T-Test process, input a blank T-Test control
card.
BF. Linear Regression Card(s)
To initiate Linear Regression, the user must input the character string
'*REGRESS' in columns 1-8 of the title card for this statistical group. To
exit from the Linear Regression process, input a blank card.
A single primary linear regression control card has the following for-
mat and may or may not be followed by secondary control cards depending on
the options chosen by the user on the primary control card.
1. Primary Control Card
1 19 22 25 28 31 43 44 47 50 53
OPTION! OPT2 OPT3 NOVAR IWGT OPT4
FORMAT (A19, 2X, A4, 2X, A4, 11X, 12, 2X, II, 2X, A4)
OPTION! = 'TIME POINTS BY SEGM1, this option indicates that the
linear regression graphs under consideration will consist
of one or more WASP model segments. The points for the
linear regression will have their y and x axis values cal-
culated respectively from station field data and WASP
theoretical data where MVP time intervals define the uni-
que points on each graph.
131
-------
'SEGM POINTS BY TIME1, this option indicates that the
linear regression graph(s) will consist of one or more MVP
time intervals as defined on the Print-Interval-Aggregate
Card. The points for the linear regression will have then,
y and x axis values calculated respectively from station
field data and WASP theoretical data where WASP model seg-
ments define the unique points on each graph.
OPT2 = 'ALLS': associated with the points of the linear regres-
sion, either TIME or SEGM (segment) points. This option
indicates all available points will be plotted separately.
For example, if the WASP model contained 25 segments and
OPTION! = 'SEGM POINTS BY TIME', was chosen by user than
there would be the possibility of 25 points on the linear
regression graph.
= 'SPEC' points, either time or segments, will be inputted on
secondary control cards. Each time or segment will gene-
rate a separate point on the linear regression graph.
= 'COMB' points will be inputted on secondary control card.
In this case each plotted point may be composed of more
than one MVP time interval or WASP segment.
OPTS = 'ALLS' associated with the option of specifying more than
one linear regression graph per primary control card. For
example, if the user has specified 11 time intervals and
OPTION1 = 'SEGM POINTS BY TIME', then 11 graphs will be
produced, one for each time interval.
= 'ALLC' all possible time or segment graphs, depending on
OPTION1 will be combined into one graph.
'SPEC' the user will input a secondary control card to
specify which time or segment graphs will be generated.
'COMB' same as SPEC except that user may combine specified
graphs into one graph.
NOVAR = The variable within current system, for which the current
linear regression is to be performed. The character string
'VARIABLE1 may be input in columns 34-42 as a mnemonic de-
vice to assist the user in reading or editing his input
data deck.
IWGT = the number of the weighting vector to be utilized for seg-
mental aggregation for this linear regression control card.
OPT4 = 'PLOT', plot all graph(s) associated with this primary con-
trol card.
132
-------
2. Secondary Control Card(s)
If the 'SPEC1 or 'COMB1 options are inputted by the user for OPT2 and/-
or OPTS, it becomes necessary for the user to supply the MVP program with
further control information in the form of secondary control card(s). In
all cases, any secondary control cards required by OPT2 are inputted before
any secondary control cards required by OPTS.
The 'SPEC' option for either OPT2 or OPTS has one of the two following
formats, depending on whether 'SPEC' is associated with segments or time
intervals.
A. Segments
4 8 12
NUMSG NSV(l) NSV(2) NSV(NUMSG)
FORMAT (2014)
NUMSG = number of segments in the 'SPEC' options.
NSV(l) = the segments to be included within the 'SPEC' options.
B. Time Intervals
4 11 20 21 30
NUMSG NIV(1) NIV(2) .... NIV(NUMSG)
FORMAT (14, 6X, 7F.10.5)
NUMSG = the number of times in the 'SPEC' option.
NIV(l) = the time intervals to be included within the 'SPEC' option.
For each MVP time interval wanted in the 'SPEC' option, a
time within that interval should b.e specified in NIV(I).
The 'COMB' option for either OPT2 or OPTS has one of the following two
formats depending upon whether 'COMB1 is associated with segments or with
time intervals.
NUMSG = the number of times in the 'SPEC' option.
NIV(l) = the time intervals to be included within the 'SPEC' option.
For each MVP time interval wanted in the 'SPEC1 option, a
time within that interval should be specified in NIV(l).
The 'COMB1 option for either OPT2 or OPTS has one of the following two
formats depending upon whether 'COMB1 is associated with segments or with
time intervals.
A. Segments
Card 1:
133
-------
NOCOB
FORMAT (14)
NOCOB = number of combinations, i.e., number of graph(s) or
point(s). Each graph or point may be composed of one or
more WASP segments.
Cards 2 through NOCOB + 1:
4 8 12
NUMSG NSV(1) NSV(2) NSV(NUMSG)
FORMAT (2014)
NUMSG = number of segments within the ith combination.
NSV(l) = segments within the ith combination. There must be NVSM6
SEGMENTS in this combination.
B. Time Intervals
Card 1:
4
NOCOB
FORMAT (14)
NOCOB = number of combinations, i.e., number of graphs or points.
Each graph or point may be composed of one or more MVP time
intervals.
Cards 2 through NOCOB + 1.
4 11 20 31 30
NUSMG NIV(1) NIV(2) .... NIV(NUMSG)
FORMAT (14, 6X, 7F10.5)
NUMSG = number of MVP time intervals within the ith combination.
NIV(l) = unique time within a unique MVP time interval within
'COMB1 option. For each MVP time interval wanted in the
ith combination, the user should include a time within that
time interval.
NOTE: Within any primary linear regression control card, no
WASP segment or MVP time interval should be specified
more than once.
Example of time interval specification:
Assume that the Print Interval-Aggregate Card (Card AC) supplied was
as: 0.0 360.0 3 so that the year is broken into 12 monthly average. To
134
-------
include January and February in a 'SPEC1 option or in the ith combination of
a 'COMB' option, the user should input the secondary control card as
4 11 20 21 30
where: T]: 0. £ T] <_ 30.
J2' 30. £ T2 £ 60.
Any values for T] and T£ which satisfy these constraints are valid.
We recommend a simple approach by splitting the desired time interval in
half. Then
TI = 15. and T2 = 45.
BG: Relative Error Cardfs)
To initiate Relative Error calculations, the user must input the
character string '*REL ERR1 in columns 1-8 of the title card for this
statistical group. There are two different types of relative error control
card: primary and secondary. For each relative error calculation performed
there is one and only one primary control card. However, there is at least
one and maybe more than one secondary control card.
1. Primary Control Card
4 8 12 15 18
NOVAR IWGT NOCOB POPT
FORMAT (14, 14, 14, 2X, A4)
NOVAR = Variable within the current system over which the Relative
Error calculation is to be performed. Relative error cal-
culations will be performed over each MVP time interval for
which the user has specified field data.
IWGT = The number of the weighting vector that the user wishes to
use with any segmental aggregation occurring within the re-
lative error calculation.
NOCOB = number of relative error calculations that will be performed
in creating a relative error distribution, (also determines
the number of secondary control cards).
POPT = 'PLOT' the user should specify this option if he desires a
one page plot of the distributive and cumulative relative
error distribution for the NOVAR variable in current system
over all MVP time interval for which the user has input
field data.
135
-------
= 'NPLT1 user does not wish a plot of the relative error dis-
tributions.
1 ' same as 'NPLT1
NOTE: To exit from the Relative Error process, the user should input a
blank primary control card.
2. Secondary Control Card
There are NOCOB secondary control cards for each primary control
card.
4 8 12 80
NUMSG NSV(lj NSVT2) ....NSV(NUMSG)
FORMAT (2014)
NUMSG = number of segments to be included in the ith relative error
calculation for the NOVAR variable of the current system.
NSV(l) = WASP segments to be included in the ith relative error cal-
culation which l£ij
-------
APPENDIX A
WASP LISTINGS
Appendix A presents supporting documentation for the Lake Ontario
eutrophication model, LAKE!. The supporting material includes:
1. A FORTRAN IV compilation listing of the LAKE1 kinetic subroutine,
with the appropriate code for generating the dump and plot files
for the IBM 370 and DEC POP computer systems.
2. The WASP input data for running the LAKE1 model, with card group
identifiers for following the input structure.
3. The output results, including printer and pen plots, from the
WASP - LAKE1 model.
4. The Job Control Language (JCL) for executing WASP on the IBM 370
series computers (JCL as used on the EPA - COMNET system).
5. The task builder and overlay command files needed to taskbuild or
link-edit WASP on the DEC POP computer systems (for a DEC POP 11/45
running under the RSX-11D operating system).
No attempt has been made here to explain the structure of the LAKE!
model or its associated equations. For this the user is asked to refer to
the EPA Report, "Mathematical Modeling of Phytoplankton in Lake Ontario",
EPA-660/3-75-005, March 1975.
Finally it should be noted that although the LAKE1 model input data set
and kinetic subroutine specify the model to have 10 systems, in actuality
the model used only has 8 systems. The last 2 systems, representing the
upper trophic levels of zooplankton, were bypassed.
The listings are availabe through the National Technical Information
Service, Springfield, Virginia 22161.
137
-------
APPENDIX B
MVP LISTINGS
Appendix B provides supporting documentation for the verification
analysis of the LAKE"! model using MVP. The supporting material includes:
1. The MVP input data for performing a verification analysis on the
LAKE! model. The input data listing has card group identifiers
for following the input structure.
2. The output results of the verification analysis.
3. The JCL needed for executing MVP on the IBM 370 series computer
(JCL as used on the EPA - COMNET system).
4. The taskbuilding and overlay command files needed to taskbuild or
link-edit MVP on the DEC POP computer systems (for a DEC POP 11/45
running under the TSX-11D operating system).
The listings are available through the National Technical Information
Service, Springfield, Virginia 22161.
138
-------
APPENDIX C
MAJOR PROGRAM MODIFICATIONS
(DEC PDP-11 VERSION)
For those modelers and programmers who maybe confronted with problems
requiring system/segment combinations different than those made available by
their current version, the dimensions of the appropriate arrays and direct
access files must be changed accordingly. As a guide in accomplishing this,
Tables C-l and C-2 have been developed.
When changing the dimensions of any of the arrays which appear in common
the entire program must be recompiled since each common block must be
changed identically. All direct access file changes should be either in
subroutine FILEOC or in WASPMAIN.
The modeler/programmer should note that system/segment configurations
are not the only model parameters which require changes in code. Included
below is a list of common changes and the program modifications needed for
each.
1. A change in input data such as number of flows, exchanges, etc. may
exceed the current storage capacity of arrays and/or direct access
files listed in Table C-l.
2. A longer simulation time or a shorter print interval may cause the
amount of output to exceed the size of the direct access files
listed in Table C-2.
3. When running WASP with MVP the following files must be defined with
the same attributes in both programs, i.e., record length, number of
records, etc.,
F10.XXZ: (FILEOC in WASP, MVPMAIN IN MVP)
all other files which occur both in FILEOC and LOADR.
Future updates of this documentation will include more explicit detail
according to responses from the user community.
139
-------
OO
LU
1
i i
U-
SX
CU.
^g
oo
CO
LU
C_3
O
cC
i
O
LU
CC
I i
Q
1
OO
LU
-J
CQ
<=C
i i
py*
^t
>*
<
O
oo
Q_
OO
cC
3:
r
1
CJ
UJ
1
GO
^£
i
C-5
O
LU
1
* t i
Xj 1 1
>- JZ
X tf- 4-> LU
< i o coz
£Z i i
II . O> h-
o -j ra
oi 2: o
LU II T3 Q£
h- >- S- CQ
Z. O ZD
1 1 r. (J OO
O S- 0)
O- a) o; 2:
i tt
-if- OJ
ZD -r- r- LU
11 ,,_ j
* C U- i i
fvl QJ II L<_
T3 M
«i i CO
>- *» OO
> QJ 00 LU
i "O C_)
X -t- S_ O
U. O et
LU II U
_l X O) I
KH a; o
U- LU
CH
1 1
Q
1 1
+J 'J3 -
SZ r Q) S-
<£ i i v. -
40 i i
C LO
zi o o> OJ
O i Zi CT
U "i C
o uo ro re
<: LO > S_
i i -.
c:
o
.,
40
r-
f
r
4
Ol
Q
OJ
^
rc
r
S-
re
>
i
QJ
s^
4-
0
S-
JO
E
c
-a
ai
S- 1^.
- n
ZI X
Cr
0) CU
s-
a) q-
-C
4-> C
i
l/>
r- ^ »
>-
- x* '
r
+ 00
> -a
t~*) i-
z. o
~~^ u
LU
s:
oo
o
to
o
ro
CU CO
i- LU
CO
oo o
CIJ 2T
E ii
o
o ^
>
^~
i| ^3
O S-
OJ
s- c:
QJ i
r
00 If-
1
c
' X T
r^
+ OO
c£. "O
O S-
z. o
> 0
Q
^£
LU
QL
Q
c^
LU
^
1
r
14
M-
ai
o
o
QJ
CO
c:
ro
c~
0
X
jQ C
E OJ
Zi !-
2: u
Qi
o
z.
i- " »
QJ 1 CO
_Q T~ ^~
E X CO
zi a; o
C i ^S
tl ^_-*
a
i- O
r- U_ +
ZJ >-
cr co
OJ O-
S- 0 0
OO CQ
OJ II Z.
^f- ^^ ^ J
JJ-)
a; - -r- o re
CO JZ JZ
Q 00 00 U
2: "O
- s- >- -a
. O !-
- - o >, o
i QJ -fJ >
+ i. ,- re
C5 l~~
2T 4 *i O
- o .a -f-"
>-
CO
Q.
CJ
CQ
z.
t
.._
oo
Q.
CQ
~^f
t
OJ
-C 1
4-> 1 -r-
E 4-> 0
OO O S- CQ
r o ro o
o. 2:
CO ^*^ "
2: s- re o
ro z.
LO "d S-
sr o
O) ZJ if-
CO E
J3 OO QJ
4-> C 4-J
ZJ if- O oo
O O !- >,
i_ +J 00- -
in s- ro - *
Z5 O) S- S_ CO
co ja 4-> re >-
E C r OO
C ZJ O) ZJ i i
l < C O O
,
ZS !Z ID
CO -1-
-4J 0-
ZI CO
2: s_ is
i- * *
O) 1 OO
_a - >-
B X co
ZI O) O
If- '
QJ S--~>
S- O
r- U_ +
ZI >-
cr oo
QJ . a_
s_ o c_>
LO CQ .
QJ II Z.*~*
JZ X ' ' CO
-fj >-
QJ Ol OO
00 i JD X
If. -0 v_-
oo c: -^^-~-
>- T- o i
OO JZ +
O oo oo >-
2: TO co
i. >- Q-
o ^
^^ (j >»S
i QJ 4-> 2:
+ S- .-
o
2: if- - s-
O JO O
>.
oo
Q-
^g
z.
t
^_
CO
D_
^f
z.
\
00
sz
QJ O
JZ -i- .
4_> 40 ^
(J QJ
00 C 4J
i- ZI 00
f- >>
O oo
2: co
c: s-
CJ I - -
QJ S- 3 00
c: o o >-
r- If- -i- CO
4-> 4-> 1 1
Zi Cf_ J- * '
o o ro ~a/
S- 0.2
JO S- Q
zi QJ ro 2:
oo jn i
E S- 0
c zi o 2:
i i C If-
1
Zi C r~
ZI .
CTCM
0) 1^.
S- II
X
.QJ
JZ Ol
r
oo 4-
c
^~
+ 00
cr-o
o s-
z o
^ - 0
Q
LU
a:
Q
(S^
LU
cc:
-^J
ro
O)
S-
oo
Q
r
f[
4-
o
OJ
ft
^
Zi
z.
o-
;z.
\
E O
ZJ CO
C II
X
Ol
JZ OJ
r- OO
oo 4- >-
i- CO
S- X
oo o s:
>- 4- II
CO >-
f"^ "^3
2: QJ QJ
S- >
r- re
ZJ JZ
O"
>> QJ T3
i S_ i
C Zi
O oo O
I 1 S- 00
'-O o
O U i
i QJ re
S- S-
LO QJ
LO 4- C
i i o OJ
CO
C S-
O QJ c:
to
,
-t-
LO
1
1
o
5-
Q.
JZ
QJ 4->
JZ !-
+-> s
00 !Z
E =J
QJ S-
00 O)
00
o
4- 4J
o
oo
S- "-
0)
"E re
Z! i-
Z. CO
CO
>-
oo
o
z.
140
-------
CO
CJ
o
UJ
_J
hj" u!
>- t~
X <+-+-> UJ
i i O CDZ
C l l
II QJ h-
O _l ID
ry* Z O
UJ II "O Qi
1 >- S- CO
z: o ID
i i « 0 OO
o s- 01
ex 01 c£ z
r i I
-4- Ol
ID T- r- LU
4_> »r |
« C U_ l l
NJ Ol II U.
"I i CO
> * CO
' ' QJ CO UJ
r- -0 0
X -r- S- O
u- o <:
LU II (J
1 X Ol 1
i i rv C_)
U- UJ
o;
i i
Q
l l
4-> CO - -
C r Ol O)
3 > 3 CT
O > C
U un ro ro
0 un > s-
4->l 1
c un *
3 O QJ O)
O 3 cr
O ««i C
U un ro ro
< LO > s-
i j \^^ f
C
o
1
4->
i
r
M
O)
l/> 1
r- -i ^
U gy
, ,-yN
*^^ UJ
01 S- 0
^ o «s
r M- II
ro >-
> X3
O) QJ
-l-> S- >
LO ! ro
O) 3 J=
CD CT
S- QJ T3
ro s-
- - CO O
o J=
O S- CO
UJ O
Qi O 1
i i QJ ro
S- S^
0 14- c
O O>
QJ O1
13 S_
r 0) C
ro o i i
> E
3
i C - -
ro OX
C O) i <
r- -C II S
U_ -U X
UJ
;>~
(^
CO
^
fV
fY*t
^^
^^
1
T3
CO QJ
i >
fO ro
> CO
O) C
4-> O)
C QJ
r- J2
4-> LO
C ro
r- JZ.
i.
Q.4J
^
M- Q-
O 4->
^
i- 0
01
_Q -H*
E ro s_
3 -C 0
Z -(-> <4-
(_}
UJ
1 1
1
Ol
S-
to
T3
S-
o
o
Ol
i-
4-
O
i.
01
_Q
E oo
3 r^.
C II
X
O)
-C QJ
+O ^^
VI If-
S
^^ 0
i <»_
o -a
^ QJ
rD s-
U_ -r-
Z 3
- - CT
l ^
e^
UJ
0£
f"^
^£
LU
CC
CO
Q.
CO
Ol =£
r 3
ro c
'T *T
^
ro O)
> to
^
Ol
E S-
r- 0
-1 ^ q
4- to
O C
O
S- T-
Ol 4->
fi (j
E c
z 14-
C_3
2T
^
U_
x s
un
Q_
CO
e^
"^
Ol
c
r
+->
ITU
O
i-
JD
CO
^
f~^
~y
QJ
QJ
CO
un
un
^) *
S_ E
ro QJ
"D -4->
C LO
3 >>
O LO
J3
S-
O Q.
S- to
0) C
JO 0
E -
3 +->
C ra
S-
E -!->
3 C
P *"
X C
ro O
s: o
CO
Q-
C_j
CO
^ X
ID
Q.
CO
<=t
3
O)
c~
r-
-I-"
O
r-»
^
CO
^^f
f~i
^^
O)
Ol
CO
un
un
oo
CO
_
r*n c~
C 4->
r- 3
(J-
0 E
<4- OJ
^_)
4 LO
0 >,
LO
I
QJ S-
JO QJ
E Q.
C LO
C
E 0
3 ("
E -t-1
- u
X C
ro 3
s: «4-
co
Q.
sx
Z
co
>-
CO
(>
"^
QJ
Ol
CO
CO
1
un
r-
LO -4->
E S-
QJ ro
-l-> Q.
LO
>> s-
co O
M-
o c
^
S- S-
QJ C
ja 01 o
3 CO
c c s-
ro QJ
§o >
E -c t.
i CJ ro
X <- i
0 JZ 3
s: s 0
CO
CO
X
o
CM
o
CO
CO
c
QJ
§>
QJ
to
'1
o
s_
cu
1
c.
E
1
X
ro
SI
0
UJ
CO
X
s:
141
-------
(*
8
i
O
LU
CO
<
I
H
tvl
>- t-
X 4- -I-)
i i O CT)
II . QJ
0 _J
Qi Z
LU II -O
fe'-fe
i i - CJ
O S- QJ
0- QJ Ct
"4- QJ
c-
2
i i
LJ_
LU
z.
t i
1
rD
O
o;
03
ro
GO
Z
h 4
^>-i-i LU
4-> -i- _J
« C Ll_ i i
fvl QJ || u_
-a hg
"I i GO
>- -GO
CD tO LU
TJ O
X -I- S_ C_>
u_ o >=c
LU II 0
i x QJ ^
j-i a: o
U_ LU
o;
> 1 1
c
r
^^ ^
M
-C-
40
CT)
c .
QJ O
i r->
u
ox
s^
QJ T-
.4-
^fe
4-> 4-
(/) t3
r- QJ
S-
i ^S
+ a-
^: a>
C£ S-
QQ
O to
'Z. -O
o
CM S
i i i
-M ;JD - -J
C i QJ Qj]
3 3 (
O « £
(J LO "3 f
ULO ^ C
U J *^» >
=3;i i v_
-i-)i i
c: LO --
3 o QJ a
31 Q
§ 2
-J C£
LJ
0 r- ^ Oi 0
o -i c i. c£
C
o
+->
C
4!
Q)
Q
QJ
r~-
.a
to
r
%
o
>
O
4->
~O QJ
Q) CT)
to c
3 C
4- QJ
O QJ -r-
-d O
S_ -1- -^
.
to
to <
-^ S sz
TCJ QJ O
QJ (J -r-
S- O) -l->
J3 -i- U
a. c
U_ T
Number o-
describe
volume fi
ro
0.
GO
<
2
s=
rvi
-C
+->
CT)
C .
QJ CM
i r^.
u
a x
S-
O QJ
(J i
QJ !-
S- 4-
QJ S-
-c o
4-> 4-
tO "O
r- QJ
i.
**" ~>. *r~
T o-
^ QJ
8S ^
O tO
Z. T3
^ 0
CM S
Q
<
LU
C£.
Q
<
LU
C£
' breaks used to
the approximation
function.
O QJ 5
.a o
4-
' *
=T
CD-
GO
et
2
C
r-si
-C
M
CT)
£Z
QJ O
i 00
II
a x
S-
O QJ
O i
QJ -r-
S- 4-
OJ t-
-C O
-(-> 4-
00 T3
^ QJ
S-
T o-
^ QJ
o: s-
03
O 00
^ T3
0
CM 2
f""j
1^1 J
f^x
Q
<:
LU
rv
breaks used to
the approximation
. cone, function.
4- -0
O QJ C
J& 13
i. -i- O
OJ S- JD
X3 U
E to s_
Z3 QJ O
^ T3 4-
^
LO
CL.
GO
C
i
rxi
_^
M
^CT)
C
OJ O
p LO
II
-o x
<_
O QJ
O '
(!) -r-
S- 4-
CU i.
-s: o
-t-J 4
00 T3
r- QJ
C
2-(NOBRK+l)
words requii
Q
<
LU
a;
Q
=c
LU
C£.
breaks used to
the approximation
orcing functions.
4- 4-
O Q)
-Q QJ
S- T- _C
0) S- +->
JD U
E S-
3 QJ 0
Z -D 4-
^ *.
to
0.
GO
§
C
txl
_£Z
-I-J
CT)
£Z .
QJ PO
»"" r^^
H
T3 X
S-
O QJ
(_) r
QJ -r-
S- '1
OJ i.
-c o
-M 4
to T3
^ QJ
&-.
^ ^ -r-
i 13
+ o-
i^ QJ
Q; s-
03
O to
^ TD
o
CM S
Q
>=C
LU
a:
Q
<
LU
C£.
breaks used to
the approximation
inetic functions.
4- .^
O QJ
-Q QJ
S- <- SI
OJ S- 4->
-Q U
E to s_
13 QJ O
Z: T3 4-
--^
1^
D.
GO
<:
142
-------
r-
>
>
H
F
<
*
c
-*
*
1
23
8
o
LU
1
CQ
UJ
-H o cnz
C i l
1 . QJ 1
0 _1 ^
si z o
LJ ii -o o:
- >- S- CQ
z o ^
-H ts O CO
O i. C Ll_ i <
-J QJ n LL.
-O r-J
«M 1 CO
>- *CO
QJ tO LU
-a o
X "- S- 0
l_i_ O <=C
LU II U
_J X QJ I
H-. CU <->
U. LU
o;
» *
Q
1 1
4-> <^5 '"^
C r QJ QJ
13 i zs cn
o « c
o uo "3 "°
O LO > S-
S-
i i >. -
^-
0
r~
-t-3
C
<4_
Q.
C
a
i
.£
n
S.
n
>
M
JZ
+->
cnLO
c r>>
QJ II
i X
-a QJ
o -r-
0 M-
QJ
S- S-
o
QJ 4-
^
2-(2NOBRK+l) is tf
in words required
Q
<:
LU
Oi
Q
to .
to
c -
O -
r 1 1
i ^ ^ ^
>o \
t_
cn<+-
QJ O
f-i
c Jr,
r- QJ
M E
O 3
C
S-
QJ «
B a]
E QJ
13
ZZ *i
^ ^
J i
'
S
3 Q-
OO
«=C
3 3
4->
(/) -
^ C3
El^
£% h-
CO - Z
i*, S
s: CL.
+ l-< h- ^.
i z s;
02S::1±
i 1 |O- r
&AI ^
3 A|
0 >i
i. >i
<- +J
+j ro a>
^ >
i +-> n3
i JC
o
to to +J
S- l/>
QJ -a =
J3 QJ E
E cn
^ c o
c: « r
^:
cu o QJ
^1 QJ -i-
U_ J3 U_
r
CU
-a -a
O QJ
E 1o
^: s-
O QJ
r- C
^: QJ
3 en
£ *"
fO J2
i O
'C -I-J
> .
S- >>
C -M "3
- rs -o
Q.
QJ | ' g~
g 3 T~
i o -. *
J ^^ -»__-
1
Qi
D-
O
-l-> T3
QJ
cn cn
c cu c:
r- _Q ^O
T3 JC
S- O O
O +J
0 h-
O to ^3
fO i > o
C h-
n -^- 1
^ *^~ _ '
QJ O O
cn a.
C QJ
fO M- C
JC O T-
(-) -4->
s- ^
QJ QJ O
ja ja s-
E J3
T3 = ^
i C to
3
o -o -o >>
_C QJ C i
to s- "3 cn
r- C
sj- to -a T-
LO QJ QJ T3
T3 +J S-
QJ -4-> O
r- QJ O O
i- -C 0
LU -4-> DL "3
t
O)
E
r
+->
C
o
-l->
fO
3
E
i
to
(O
4_l
o
1
o
^
LU
1
1
<4-
'i
-o
to
QJ
!Z
4->
3
0
S-
JO
to
+->
values in differen
-i->
C:
C/)
i
-M
**
21
O
Q_
CO
c.
1 ^
o
c
to
'""
s/
^^
Qi " ^
CQ to
0 C
Z O
i
4- -t->
O O
(^
QJ 3
Z5 4-
n3 -l->
> C
o>
QJ i-
_C QJ
1 1 4-
*
1 -K
143
-------
CO
>-
<
Q;
@g
r>
00
<:
CM
CJ
LU
I
ca
<=c
A
CO
QJ
r^
re
re
=> o-
a
S- r
ro -r-
i U.
ro
o
oo e
c,
«<
CO
>>
re s-.
i- T-
S- Q
<
-a
QJ
4->
re
"QJ
o;
' ^
E1 '
3 VO
H^s
X «> i
re LO Q
S^LO
^r i
3 LO
r y s
.Ji2s-
X "I 1
re LO Q
^tifj
o
i
-M
r
C
4^
QJ
Q
>1
re
S-
i.
<
i
r^
\S
Z:
CJ
CO
Z
4_
O
C
o
CO
c.
QJ
E
"O ' ^
CO
ro oo
S- X
QJ s:
c-
QJ
rj> n
C CO
1 1 "r
<£>
V
LO
f\
>,
+-> «r
X SZ QJ
QJ O >
C -P- -i-
4-> +J
QJ CJ CJ
^: c QJ
-»-> 3 Q.
4 00
-C QJ
O QJ i.
i^^3
C CJ
-t-> -r- S- O
re o o
.*:
QJ re . i
E QJ CJ r
1 J3 JO 5
[
I^~
\^x
Z.
n
H-
O
CO
Z.
VI
>-
QQ
o;
VI
>-
CO
OO
CO
l[
0
c
o
to
c
QJ
E
o co
>-
r OO
re x
s- s:
QJ -.
C
QJ II
CD
C CO
h- 1 CD-
CO
^^
LO
c
O
c
V. O
QJ -r-
OD+->
c n.
ro o
^:
O CO
X 00
QJ re
a.
« >>
E ^O
QJ
-t-> 3
oo O
CO 4_
>-
viCQ
>- CD-
CO VI
OO >-
>- CO
OO fV
4-
o
II C
o
_l !-
O CO
> c
SI QJ
4 'P~
o -a
C i
o re .
>r__ * ^^-^
CO QJ CD
f f-~ I^J 1
QJ QJ OO
E D1X
f- S!
T3 C-
p i
CM
re
i. II
QJx^
c CD cr
QJ LU s:
O5OO
X V,
^Sg
* >.X N
O O
CO =*
^IT CD-
CD 21
^^ VI
Slcg
--^o
O CM
VO r
_i cS-
os:
^^ ft
s: a:
| >
o> oT-x
re 3 4-
s- p o
4-> O
to > 3
^C54?
3 QJ
i cr cn
re QJ c
> re
QJ _C
= C 0
E -r- X
= QJ
A
g
«%
|
X ^
IE
,
c
QJ
E
"O
II x-^
>-
o re co
CO s- 0-
S QJ i-
i CO
00 X
£^S
E >-
i- CO II
0 Q.
0 :*i
r- 03 3
re z s:
QJ «q_
C CO O
QJ >-
cnoo c
X 0
C SI -r-
ll CO
LO LO
VI V.
CO CO
"-"
LO LO
r- CO
VI VI
LO LO
CO
C
O
-4-* *i
.C *->
cn . o
r "O C
re C 3
S- 3 4-
-(-> 0
OO .Q CT>
c:
c s- >,
r- O CJ i
4- S- QJ
QJ O >
3 It- T-
r CT +J
re QJ vi u
> . QJ
QJ CJ Q.
= c c: co
E -i- O QJ
= r O i-
|
*
O
CD
^~
^
u
c_>
o:
u.
Q
I i
4-
O
C
o
r
CO
C
QJ
E=
i
-o
re
QJ - s
C CO
QJ >-
CD CO
X
cs
CO
LO
1 QJ
-C* -C
=3 +->
CO
QJ C)_
Q) O
r- .C
4- -t-J S-
QJ
s- .c: _o
re +j E .
r - 3-0
3 3: c s-
o o
r- -o QJ CJ
*-> QJ _c QJ
i- (-> -l-> S-
re re
Q.-1- - ^4->
0 -l-> C
re o Q. QJ
CO T- S-
S- CO s_ S-
o re cj 3
LL. oo CJ
o
a:
LL.
o
t-H
II
Q
CJ
CJ
H-
o
c
o
1
CO
c
QJ .
JtD
-o uj
CO
i X
re 2:
ai v.
c oo
QJ >-
cnoo
X
c s:
, .
o
CM
CO
C?
CO
VI
LO
I
QJ
"O
>>
OO i
-t-> QJ
r- >
-0 -M
re QJ
CL
Concentration
rivative, res
o
0
*.
CJ
144
-------
^^^
h
8
Cvj
o
) 1 1
1
CO
CO
d)
-0
rO
i
rd
> CO
0)
i- f
rd T-
i ti-
ro
(J
CO O
O
\ «^
CO
rO C
S- T-
S_ Q
^c
"D
o;
4-*
rcs
r^
(D
a:
E ' '
3 <^O
E r-
i i y
X «l 1
t3 LO a
3 LO
E 0
-- s:
X "> i
rt3 LO c^
21 LO
C
o
f
4-3
'f~
C
1
<4-
a>
Q
ro
S-
s-
r*
X
=C
21
14-
O
c
o
1
CO
c*
O)
E -"^
~o ^^
oo
i X
'O ?T
S^- * »
O)
C II
O)
cnz:
i i
c: 2:
i i <_)
, X
*-O
^
^ ^
un
CO
c
o
E -M
3 ro
E S-
C C
i- a>
E 0
c
o o
c o
ro
a>
E -Q
3 cC
E 5
i~ O
X i
ro .
X "Z.
"^ t-H
§5
11
y
^^
r^
Q_
tj^_
O
C
O
r
CO
d
0)
"5 ^
Q-
o
i-
OJ n
O) LU
cnoo
X
c 2:
o"
N
o
CM
to"
f,
o
to
Ll 1 40
21 C
<=c cu
=C an
^_
C|_ ro
o a.
QJ *O
3
i S-
rO O
=> M-
^c
cv
Q-
^
y*
"~i
U.
M-
O -^
Q
C <
O UJ
i~ ry^
CO
C
QJ - v
££
^ ^~
-a o
U_
ro
S-
O) II
c~
OJ C_J
cnz
ZD
C U_
i I CO
,. X
o
^^
o"
S-
o
CO
4-> CO
Q. C
CD 0
(J !-
1 4_>
O) (_)
4-> C
C 3
r- C|_
t\ O
CO -r-
O> 4->
Q. O)
o c
I/O -^
00
z z
u_ u_
s: QQ
II
1
CO
8
M-
0
d
0
r-
CO
c
a>
E
-a
i>
ro
^_
CD
CD
0>
c*
i i
, ^
tn
LT)
i
UJ
z.
4->
c~
ro
4_j
CO
C
o
CJ
q
O
CU
3
rrj
1
CO
z.
8
oo
GOVERNMENT PRINTING OFFICE 1983/605-095/1935
145
------- |