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
>,
4—1
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
OrOr—LJLJLJLJLJULJ

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-i—io
                                                                                              !_• 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




    i—i—i—    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 — t—t
-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-
l—l 	 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

-------