WATER POLLUTION CONTROL RESEARCH SERIES • 16110 EGQ 04/72
  EXTENSIONS  OF  MATHEMATICAL
   PROGRAMMING  FOR REGIONAL
   WATER QUALITY MANAGEMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY

-------
      WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollution
in our Nation's waters.  They provide a central source of
information on the research, development,  and demonstration
activities in the water research program of the Environmental
Protection Agency, through inhouse research and grants and
contracts with Federal, State, and local agencies, research
institutions, and industrial organizations.

Inquiries pertaining to Water Pollution Control Research
Reports should be directed to the Chief, Publications Branch
(Water), Research Information Division, R&M, Environmental
Protection Agency, Washington, B.C.  20460.

-------
             EXTENSIONS  OF MATHEMATICAL PROGRAMMING

             FOR REGIONAL WATER QUALITY MANAGEMENT
                                by
                    University of California
                  Graduate School, of  Management
                  Los Angeles, California  90024
                              for the

                OFFICE  OF RESEARCH  AND MONITORING

                 ENVIRONMENTAL PROTECTION AGENCY
                       Project #16110 EGQ
                           April  1972
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402 - Price $1.00

-------
                  EPA Review Notice
This report has been reviewed by the Environmental
Protection Agency and approved for publication.
Approval does not signify that the contents necessarily
reflect the views and policies of the Environmental
Protection Agency, nor does mention of trade names or
commercial products constitute endorsement or recommenda-
tion for use.
                            11

-------
                          ABSTRACT
This report extends the work contained in the earlier work,
"Mathematical Programming for Regional Water Quality Manage-
ment."  A mixed integer, continuous variable non-linear pro-
gramming model is developed which promises to be much more
realistic and effective in selecting an optimal configuration
of regional treatment plants and pipe juncture nodes.  A new
Bender's type decomposition for the non-linear model is also
presented as an appropriate and feasible solution technique.
In the area of the physical response of the environment, the
earlier work linked the dissolved oxygen profile and the con-
trol options through a constant transfer matrix obtained as a
steady state solution to a system of differential equations.
This work shows how the transfer coefficients can be recalcu-
lated periodically in the course of the system optimization.
Full scale computations are carried out demonstrating this
capability.  The theoretical framework for including two addi-
tional treatment methods, mechanical re-aeration and flow aug-
mentation via reservoirs is also developed.'  In the social and
economic area, a new proposal for paying for a regional system
is presented.  Illustrative computations using the new proposal
are included.  This report was submitted in fulfillment of
project 16110EGQ under the partial sponsorship of the Environ-
mental Protection Agency.
                              111

-------
                   TABLE OF CONTENTS
Section    I.

Section   II.

Section  III.


Section   IV.

Section    V.


Section   VI.

Section  VII.
List of Figures

List of Tables

Conclusions and Recommendations

Introduction

Selection of a Waste Treatment
   Plant Configuration

Transfer Functions

Estimation of the Depth and Velocity
   of a River for a Given Rate of Flow

Transient Effects in an Estuary

Pollution Taxes and Charges in
   Perspective
Section VIII.  References
Page

 vi

vii

  1

  3


  5


 15


 37


 61


 67


 81
                             V

-------
                       LIST OP FIGURES

FIGURE   TITLE                                              PAGE
   1      Schematic of General System                          6
   2      Optimization Process,, Assuming Fixed Transfer       20
            Matrix
   3      Iterative Procedure for Finding Stable Transfer     21
            Matrix
   4      Solution of By-Pass Piping Problem Using New        33
            Fixed Transfer Matrix
   5      Solution of By-Pass Piping Problem Using New        35
            Variable Transfer Matrix
   6      Schematic and Notations for Spot Inflows            44
   7      Schematic and Notations for Derivatives of Spot     45
            Inflows
   8      Schematic and Notations for Spot Outflows           47
   9      Critical, Super-Critical and Sub-Critical Flows as  50
            Depicted by the  Depth-Specific Energy Relation
   10     Schematic and Notations for Continuous Inflow       57
                               vz

-------
                      LIST OF TABLES

TABLE    TITLE                                           PAGE
  1      Net Flows into Section i                         23
  2      Eddy Exchange Coefficients Between Sections      23
            i-1 and i
  3      Volumes of Section i                             25
  4      Reaeration Rate in Section i                     25
  5      Transfer Matrix                                  26
  6      Thomann's Transfer Matrix                        29
  7      Transient Period of Section i, Corresponding     65
            Starting Conditions and Computational Error

-------
           I.  CONCLUSIONS AND RECOMMENDATIONS

The benefits to be derived from the use of regional optimiza-
tion models for water quality management would seem to be enor-
mous.  Regional optimization models are applicable to bays,
estuaries and rivers.  They are capable of dealing with any
quality goal for which an appropriate system of diffusion, diff-
erential equations can be obtained to represent the physical
response of the environment.  In particular, they are capable
of dealing with dissolved oxygen, salinity and thermal contamin-
ation.  The control options include treatment at the source,
by-pass piping, regional treatment plants, mechanical re-aera-
tion and flow maintenance augmentation.  The'output from the
model gives the least cost solutions to achieve the pre-selected
goals using the full range of control options.  The optimal
solutions would specify the treatment levels at the major sour-
ces of contamination, the complete regional pipe network, the
location and levels of treatment of an optimal number of region-
al plants, the location and level of mechanical re-aeration and
the location and level of flow augmentation.  The model could
be used in conjunction with regional growth of demand projec-
tions to plan the long term treatment capacity expansion.  It
could ultimately be used in a real time regional control pro-
cess utilizing physical and chemical treatment as well as the
traditional biological treatment to capitalize on seasonal and
daily flucuation to greatly reduce yearly and daily treatment
costs.  This report develops the theoretical possibility of in-
cluding two additional treatment methods, mechanical re-aeration
and flow maintenance or augmentation via reservoirs.  In order
to deal with the question of the effect on the re-aeration co-
efficient of radical shifts in flow, Section VI investigates
the effect of flow on the depth and velocity of a river.  Fin-
ally, Section VII takes a preliminary look at the transient
effects in an estuary as opposed to the long term steady state
conditions.

In the area of system optimization, Section IV  examines the
question of an optimal configuration of regional treatment
plants and pipe junction points.  The original model deals with
this question by forcing a solution using only the regional
plants.  This forces the plants to operate at sufficiently high
levels as to utilize their economies of scale.  When the other
treatment methods are then introduced, they are compelled to
compete against the efficiently operating plants.  Where econo-
mically justified they then force out and close the non-compe-
titive plants.  This is computationally wasteful as the plants
must first be opened then some of them closed.  It is also
mathematically doubtful as to the global optimality of this pro-
cedure.  The mathematical foundations of a potentially much
more efficient procedure using an integer non-linear programming
model are presented.  Computational effectiveness hinges upon
the power of a new decomposition technique and the capability
of rapidly solving integer programming problems.


                            - 1-

-------
the realism of the models would seem to be adequate at this
time although they may ultimately be considereably improved.
The computational feasibility for full scale problems has been
demonstrated.  It would seem entirely practical to proceed to
actual pilot studies at this time with the tools already at
hand.  Further development of the theoretical possibilities
should proceed concomitant with the pilot studies.  In parti-
cular computer implementation of the new mixed integer, con-
tinuous variable non-linear model should be undertaken to sub-
stantiate its potential value.  The possibility of including
the additional treatment methods of mechanical re-aeration and
flow augmentation should also be exploited.  On a longer hori-
zon the dynamic or transient behavior of the physical environ-
ment should be studied more closely.  It would be entirely possi-
ble to realize substantial gains from an optimal holding and
discharge policy based on the transient behavior of the environ-
ment.  Finally the mathematical programming model for regional
water quality management should be coupled with a regional de-
mand model.  Rational planning over a realistic amortization
period requires fairly good estimates of projected demand over
the period.  The capital budgeting problem associated with the
optimal timing and amount of treatment capacity expansion should
also be carefully studied.
                            -  2  -

-------
                    II.  INTRODUCTION


This report continues and extends  the work of the earlier re-
port, "Mathematical Programming for Regional Water Quality
Management."  The goal of this work is  to develop a general
practical planning tool to provide optimal solutions for the
complex choices  involved in balancing alternative methods for
attaining water  quality goals9  The work encompasses three
broad related areas.  The first area is the social and econo-
mic question of  what are desirable and  feasible water quality
goals and what is the best way to  achieve and pay for them.
The second area  is the complex technical engineering problem
of the physical  response of the environment to various control
strategies.  The third area is the mathematical problem of ob-
taining least cost solutions in the total regional water system
that will satisfy the desired goals.

In the area of physical response our earlier work presented a
general mathematical model relating the control options of by-
pass piping, treatment at regional treatment plants and treat-
ment at the source to the-quality  goal  specified as a dissolved
oxygen profile on an estuary.  The original model was limited
in several respects.  The coupling of the dissolved oxygen pro-
file goal and the control options was through a constant trans-
fer matrix obtained as a steady state solution to a system of
differential equations first given by Thomann [10].  This is
inadequate when  by-pass piping and regional treatment plants
substantially alter the flow.  Section  IV of this report re-
examines the question of the transfer functions.  This work
shows how the transfer coefficients can be re-calculated period-
ically in the course of the system optimization.  It further
develops the theoretical possibility of including two additional
treatment methods, mechanical re-aeration and flow augmentation
via reservoirs.  In order to deal with  the question of the ef-
fect on the re-aeration coefficient of  radical shifts in flow,
Section V investigates the effect of flow on the depth and velo-
city of a river.  Finally, Section VI takes a preliminary look
at the transient effects in an estuary  as opposed to the long
term steady state conditions.

In the area of system optimization. Section III examines the
question of an optimal configuration of regional treatment
plants and pipe  junction points.   The original model deals
with this question by forcing a solution using only the re-
gional plants.   This forces the plants  to operate at suffici-
ently high levels as to utilize their economies of scale.
When the other treatment methods are then introduced, they
are compelled to compete against the efficiently operating
plants.  Where economically justified,  they then force out
and close the non-competitive plants.   This is computationally
wasteful as the  plants must first be opened then some of them
closed.  It is also mathematically doubtful as to the global
optimality of this procedure.  The mathematical foundations


                            - 3 -

-------
of a potentially much more efficient procedure using an integer
non-linear programming model are presented.  Computational ef-
fectiveness hinges upon the power of a new decomposition tech-
nique and the capability of rapidly solving integer programming
problems.

In the social and economic area. Section VII discusses economic
policy for regional water quality control.   A new proposal for
paying for a regional system is presented.   Illustrative compu-
tations using the new proposal are included.
                           - 4 -

-------
  Til.  SELECTION OF A WASTE TREATMENT PLANT CONFIGURATION


In our earlier report, "Mathematical Programming  for Regional
Water Quality Management/1 the model presented dealt  with  the
difficult question of the best treatment plant configuration
in a very limited manner.  Flow of effluent was restricted  flow
direct to a discharge point or flow to a single plant and then
from the plant to a discharge point.  There was no  flow between
plants and hence no way  to achieve say primary treatment for
a large percentage of the effluent at huge, cheap primary plants
and subsequent treatment for a much smaller percentage at smaller
more expensive secondary or tertiary plants.  Since pipe juncture
nodes can be considered  plants without treatment, the flow  pattern
was restricted to single juncture points.  This resulted in prac-
tice in manual pipe consolidation and weakened results.

The solution procedure for the original model required forcing
a solution using only the regional plants.  This  forces the plants
to operate at sufficiently high levels as  to utilize their  econ-
omies of scale.  When the other treatment  methods are then  intro-
duced, they are compelled to compete against the  efficiently
operating plants.  Where economically justified they then force
out and close the non-competitive plants.  This is  computationally
wasteful as the plants must first be opened and driven down and
closed.  There are also  difficult problems with infinite deriva-
tives in the cost function with the initial very  low or zero
flows.  These problems can be handled much more effectively using
discrete integer variables in the formulation of  the model.  This
possibility was not exploited originally becuase  of the limita-
tions of the solution techniques for large mixed  intefer and con-
tinuous variable non-linear programming problems.   Recent research
has greatly expanded our capability to solve integer programming
problems.  This extended capability coupled with  new and more
powerful decomposition techniques change the picture radically-
It now seems attractive  to switch to the mixed integer, continu-
ous variable formulation.

Mixed Integer Regional Treatment Model

In this model we can use a much more extensive and  realistic
supporting flow network.  This is illustrated by  the schematic
shown in Figure 1 on the following page.

The nodes in the network are of three types.  A node is a pollu-
ter or source, a sink or a discharge point, or an intermediate
node in the network consisting of a pipe junction or treatment
plant.  The sources are  fixed at the present sites  of the iden-
tified major polluters.  The discharge nodes are  the presently
used discharge sections  of the estuary plus any additional  po-
tentially attractive discharge sections.   The intermediate  nodes
consist of any existing  treatment plants or pipe  junctures  plus
any additional potentially attractive sites.  Combining existing
                             - 5 -

-------
   (alternate
     s i te s)
                                                        (alternate
                                                        capacities)
                      FIGURE 1

            SCHEMATIC OF GENERAL SYSTEM
and now sites into an optimal total network is particularly
nice in this approach because the fixed or capital costs can be
separated from the variable or operating costs.  This permits
a true comparison of the marginal cost of a new plant to an
existing plant.  The arcs in the network are pipes.  They can
consist of existing pipes and potentially attractive additional
pipes.  Again the fixed and variable costs can be separated.
Also it is possible to have discrete alternate capacities and
enforce a choice of pipe size.  Similarly, we can enforce a
                             - 6 -

-------
choice between alternate plant capacities at  a  given  site  or
between alternate plant sites.

In the mathematical statement of  the model we use  the following:

    Notation

       Physical Quantities

         A   - matrix of transfer coefficients

               (mg/i per Ib./day)

         C   - vector of D. 0.  goals (mg//)

         r.   -%removalofB.O.D. at node i

         q. •  - flow from node i to node j (MGD)

         J.   - concentration leaving node i (Ib/MG)

         m.   - mass of B.O.D. leaving node i  (Ib/D)

       Costs

         f.   - capital cost of plant i of a given capacity ($)

         v.   - operating costs of plant i ($/MG)/% removal)

         f..  - capital cost of pipe between node i and node j

               of given capacity  ($)

         v..  - operating cost of pipe between node i and node  j

               of given capacity  ($/MG)

         p.   - cost of additional treatment at source  ($/% removal)

    Index Sets

         R = {j I j is a section or reach of the estuary]

         S = {i I node i is a source]

         D = {i | node i is a discharge point]

         I = {i | node i not a source or sink]

         N = S U D u I set of nodes; a c N x N set of arcs

         a (x)  = {y e N | (x,y) e a] nodes "after" x

         b (x)  = {y e N  (y^x) e a] nodes "before" x

                             — 7  —

-------
Integer Variables

           r0  plant i is closed
     Zi  ~ *• 1  plant i is open

         _ rO  pipe (ijj) is closed
      ij   *-1  pipe (ijj) is open


The mathematical model is:

   Subject to

     Conservation of Flow

        Z  q .  -   Z   qk. = 0   i e I
     jea(i) 1J    keb(i) kl

        Z  q--=q.   i e S (q current discharge MG/D)
     jea(i) 13     L

     Conservation of Concentration
     J. -
           keb(i)
          m.                   _
     J. = ^-  (1 - r.)   i e S  (m .  current discharge Ib)

          qi

     Capacity at Plant  (&. - lower , (j, .  - upper capacity)

     2i • J&i ^   Z   qki ^ zi^i
      1    1   keb(i) kl    1 1


     (Note that when the plant is  closed or not chosen for
      the network all flow is shut off and when it is opened
      it functions in a restricted range.  By breaking up
      the range non-convexity can be eliminated and more truly
      global optimum achieved.  This would require merely add-
      ing some:

     Alternative Constraints  (Plantsj

      Z  z  =  1
     These constraints would force a choice of only one of
     the plants in the set E..)
                         _ Q _

-------
          Capacity of Pipes (i, .  - lower, ^. . - upper capacity)
          (Again  the  pipes can be eliminated or forced to
           function in a restricted range.  By breaking up
           the  range  non-convexity can be eliminated.  Again
           we would employ:

         Alternative Constraints (Pipes)

             2    z.   = 1
           These  constraints would force a choice of one of the
           alternatives in the set E...)

          The  quality goals are enforced by introducing the auxili-
ary variables
        m  =  J  •  (   2    q  )    t e  D
          ^   *•   keb(t)  Kt

        m  =  present discharge in Ib's

and the following:

        Quality Constraints

         2 a., (m -in.)  <^ -c.    i e  R
        teD   lt:  t   C

The objective is to  minimize  the total  cost function

        minimize

        T(zfc,zts,q,r)  = PL(zfc,q)  + PI(zts,q)  + PO(r)


which is broken down into:

        Plant Costs

                        zfc -  (ft+vt(r)-  kj(t)qkt)
Zts ' (fts
        Pipe Casts


        Pl(zts^) =  ,. Z,  a
                     V t , 5 ) 6A

        Polluter Costs
        P0(r) =2  r   - p
                teS
                             - 9 -

-------
Solution Technique

The foregoing complex mixed integer continuous variable non-linear
programming model presents a truly formidable computational chal-
lenge.  Although a direct implicit enumeration of the integer
variables coupled with the bounding off effect of solving  the
remaining non-linear programming problems might be attempted,
the outcome of such a venture would be extremely doubtful.  Re-
cent advances in Bender's type  (see [2]) decomposition offer a
much more promising avenue of attack.  A Bender's type decomposi-
tion applicable to a general non- linear programming problem

     Subject to
               gk(y)
           min
is possible.  In this type of decomposition a subset of the
variables say the y~ vector can be used to parametrize a re-
laxed Sub-Problem

     Subject to

               g1(Y1-Y2) £ o
                m /   —  \
           mm g (Y1,y2)

and the parameters y_ need not be continuous variables.

Whenever a fixed choice of parametes y9 satisfying the remain-
ing constraints

               gk(72) *o

permits a feasible solution to the sub-problem, we obtain a
feasible solution to the complete problem.  _In order to avoid
a random search through plausible choices of y~ , we employ a
class of lower bound functions D.  In a stepwise procedure we
generate members of this class.  By ensuring that at least each
lower bound function already generated indicates a better possi-
ble solution than the best known incumbent feasible solution we
guide the choice of the parameter variables.  In point of fact
since it is a necessary condition that each lower bound func-
tion be less than the incumbent when there exists a better solu-
tion, the lower bound functions provide us with a convergence
criteria.
                             -  10  -

-------
In order to generate the class of  lower bound  functions  of
interest here we require the  following expansions  of  the
functions.

Expansion 1   (With respect  to y, )

(i)  gi(y,y) =gi(YJy)  +  vgi(yJy)  -  AY
Now in order to eliminate the gradients and utilize  informa-
tion furnished by the dual variables  in the sub-problems, we
must make:

Expansion 2  (With respect to y )


(2)  Vg1(y°,y2) - AY]_ = vg1(y°,y°)  •  AYX
where                2 i
                    5 g
is a matrix of mixed partial derivatives.  Combining the error
terms
(3)
                     AY2
                          1
yields the desired representation of each function


(4)  g1(y1,y2) = g1(Yj,y2) + vg1(y°Jy°)  - AyL + Rx(Ayr *y2)


Using this representation of our functions and the following
relation,
(5)  Z x± vgy^y^) = vgm(y°,y°)
which  is derived  as  equation (14)  on page 17 of our previous
report,  "Mathematical  Programming for Regional Water Quality
                          - 11 -

-------
Management",;yieIds our class of  lower bound  functions.
The argument employed in establishing this  class  of  lower
bound functions is similar  to the argument  used to establish
the "Weak Duality Theorem"  of linear programming.  Any e-
feasible solution of the sub-problems satisfies


(6)  gi(y1,y2) = gi(y°;Y2)   + vgNy^y^) ' &i + ^(AY^AY;,) ^ el
or
 (7)
Using  (5) and (7) we demonstrate that as a function of

For any value of y,


 (8)  gm(y1,y2) = gm(Y^Y2) + vgm(yj,y°) • ^ + Rm(Ay-L
                                                         AY
                                , AY2)
                 + 2 x. -e.



For convenience we restate  the conditions  for  a stationary
point of the non-linear algorithm presented  in our previous
report  [4] on page 19.

Local Linear Stationary Point

a a set of dual variables x  3
                           P


1)  2 (-x ) > B-.    (insufficient  resolution)
    P   P
                            12 -

-------
a a w and H  3

2)   2  (-x )  •  gp(y°)  >  -gw(y°)  - e   (Inconsistency)
    peH   p
a a y  e F  ^

3)   m~^
    .2   (-Xi)  - g1(y°) >  -  e    (Optimality)
Now choose £-,_ = e/B^  and  assume  that condition 1)  of a local
linear stationary point is  not violated then


    f Xi >- -Bl   and   I  Xi ' el ^ -£

and when  the error expression


     i              '

we have our lower bound function


(10) gm(y°,y7) + 2(-x  ) - g1(y°y0)  - e^gm(Yl,y_)
         1  2    i    i       12             12


(Note that in the case where  the functions  are separable in
 y., and y« as in our  model,  the  convexity of the sub-problems
 in the continuous variables  is  sufficient to ensure condition
  (9) and the validity of  the  decomposition.)

The decomposition procedure then consists of alternately sol-
ving the relaxed sub-problems and  the following

Master Problem

Find a feasible solution  to,

        i)  gk(Y2) £  o



                                            (incumbent)

        3)   g (y-,,Y9)   +   2  (-x )g  (y-,,y0) <1 -e
              "         peH  p      L  2
                           — 13 —

-------
where a constraint of type  2)  is added when  an  optimal  solu-
tion is obtained for the sub-problem and a constraint  of type
3) is added when the sub-problem is inconsistent.

Convergence to an e-optimal solution in a finite number  of
steps using this decomposition  technique is easily  demonstrated.
Assume we return from a sub-problem without achieving  a  bound
ed gain in the incumbent solution when we obtained  a feasible
solution then

,..,    m,o   »  .   m,c  cx
(ii)   g (y1Jy2)  2 g  (Y1^2)

and

(12)   Z(-.x,.)  . gi(y°,y2)  > -e


if this set of dual variables is a repeat we contradict

,,_\    m,o   >  ,  „ ,   ,i,o   \  ,   m/c  c\
(is)   g (y1,y2)  + s(-xi)g  (y-^y^ £ g (y^y^  ~  €


which has been established from the master problem.  Since
there are  only a finite number  of dual solutions between bound-
ed gains in a finite number of  steps we exceed any  lower bound
or fail to solve the master.  When we fail to solve the  master
we have established  from one of the lower bound  functions that
we are within e of the optimum.

When we return from  the sub-problems with an  inconsistent
solution we have

(14)    E (-x ) - gp(y° y )  > -gw(y° y )  - e
      peH   p        *-  *         ^  z

Again if w and H are a repeat we contradict

(15)   gW(yJ,Y2)  +  Z (-X )gP(y° y )  £ -e
                  peH   F


which has been established from the master.   So  in  at  most a
finite number of steps we  resolve the inconsistency or demon-
strate that the complete problem is inconsistent.

This type of decomposition has  recently been  employed  in a
large warehouse location problem  (see [3]) and proved  quite
successful.  It would seem to offer great promise  in the area
of regional water quality  control.
                           -  14 -

-------
                 IV.  TRANSFER  FUNCTIONS


Part 1:  Theoretical Considerations

We begin by considering  the  two fundamental  sets  of mass bal-
ance differential equations.  The  most  general  form of these
equations is given by Thomann [11] or as  equations  (1)  and
 (2) by Pence, Jeglic and Thomann [8, pp.  383-83] .   These basic
D.E.'s characterize a one-dimensional estuary divided into n
sections.  The  first set gives  the transformation of polluter
BOD to in-stream BOD and can be written as

  V.dL.
   "alT = Qi-l,i[5i.1,iLi.1-Kl-gi.lji)LiI-ffli_lji(Li_1-Li)
where the sections  are  numbered  from  upstream  to  downstream.
LJ_ is the ultimate  carbonaceous  BOD in  section i.   QI_I i i-s
the net flow  from section  i-1  to i.   Vi  is  the volume  of sec-
tion i.  §i-i  i  is  the  advection coefficient from section i-1
to i.  E^_;L ±' is the eddy exchange coefficient from section
i-1 to i.     d is the decay rate of BOD in  all sections.  K^  are
all direct sources  of BOD  discharged  to section i.

Since we require dL^/dt =  0 in the steady state (SS) the above
system can be  written as follows
where
       and r. . = 0   if  | i-j
              •J
In matrix notation we have

        RL +  K =  0

where R is a  square matrix  of  order  n.   Assuming R is non-
singular we have

        L = -R"^

The second set of D.E.'s  gives  the transformation of stream

                           - 15 -

-------
BOD into D.O. and can be written as
     dC.
                          -dLi)+Pi
where GJ is the D.O. concentration in section i.  r^  is  the
reaeration rate of D.O. in section i.  Cs^ is the saturation
concentration of D.O. in section i.  P^ is  all other  sources
and sinks of D.O. in section i.  Since we require dC^/dt = 0
in the steady state the above system can be written as follows
0 =
where
         s .
        and s..  = 0    if


 Let V be  a diagonal  matrix  of  order  n  with  the  V.  as the diagonal
 elements
         V =
V,   0

0    V,
                            0
                           V
                             n
 and  let W  be  the  column  vector
        W  =
               rics
               rnCS
                           -  16  -

-------
In matrix notation,  the above  system  is

        SC = dVL - V W + P

where S is a square matrix  of  order n.  Assuming S  is  non-
singular we have

        C = -dS^V R-1K - S~1V W + S~LP

The essential idea in calculating the SS  transfer matrix is
to inject 1 unit of BOD into each section and calculate  the
corresponding D.O. changes  in  each section.   A symmetric in-
terpretation is to remove 1 unit of BOD from  each section.
This gives a sign reversal  in  the final result.   Replacing
K above, by ej  (unit vector) and noting that  the 2nd and 3rd
terms on the right can be dropped we  define

                     R~1e .
                          3

where (f) . . is the jth column of the transfer matrix.  The trans-
fer matrix $ is then

        5 = -dS'-'-V R"1! = -ds'V R"1  = -d(R V~1S)~1

The transfer matrix  corresponding to  a 1  unit removal  in each
section is

        $ = d(R V^S)"1

The mass balance D.E.'s hold for non- tidal streams  as  well as
estuaries.  It is an interesting exercise to  show that the
classical Streeter-Phelps equation is contained in  these D.E.'s,
The Streeter-Phelps equation can be derived as follows:

In the absence of reaeration assume that  the  rate of change
of BOD at time t is proportional to the BOD at time  t
                      or
where yt is the demand  for  oxygen  (BOD)  at  time  t and L is  the
ultimate BOD at time zero.   The  general  solution to this D.E.
is [1, Thm.l, p. 40]

         (t) = ce~klfc


since L = (j)(0) = c, we  have

                  = Le"1^!1^
Similarly in the absence of  demand  for  oxygen assume that the
rate of change of D.O. deficit  (DOD)  is proportional to the

                           - 17 -

-------
DOD at time  t.


        diy
where Dfc  is  the  D.O.  deficit at time t, i.e., Dfc - Cs-Ct-

as before
Combining  the  above D.E.'s we have the Streeter-Phelps D.E.






        -r-7— =  k  y -k0D,    or   D!+k0D. = k,y,
        at      1 t  2  t         tzt    J_  t




The general solution is [ 1,  Thm. 2, p. 41]



                 ^"[^ «-» 4-    \f rt \r           •—T^" «""* "I"
          / _i_ \      JV / >— /*    JX ,' ^.T    T   ,    JV / •—
          (t) = e  z  /   e * k,y dx + ce  ^
                             n

                      0
using y  = Le      we  have
                1
               ,    -

          (t)  =   \\    f   e2-l    (v -k.)dx

                  k2~kl   0             2



               klL r -kit    -k2t,  ^  ^ -
              =  Z   t®    ~ e      + 'de
from the initial value  problem we have D_ = () (0) =  c  and the

classic Streeter-Phelps equation is
In the Streeter-Phelps D.E.Q. let D, = C -C.  giving
                                     t    w   t


        d(Cs-Ct)


           dt       klYt~ k2 (CS~Ct)



   or





(16)     -XT- = k2CS~ k2Ct~ klYt
        vj.1—      ^O    ^L-   J_l_




In  the mass  balance D.E. giving  the transformation of stream
                            - 18 -

-------
BOD to D.O. let E,.= 0 and  £.  ,  .=1 which gives
                  i j         • i— J.} i

     dCi    Qi-l  iCi-l    Qi  i+ici
 H?^   -    = —    '       -   3- , 3-TX 3. ,    _      _    -,T
 (17)  dt        v±          V..     +  i  Si   i i "    i


Note that  in  (16) and  (17)  the  following quantities are  equi-
valent:   k2~r, k-j^-d, and  yt~L.   Comparing (16)  and  (17)  we can
see that  the  terms of  (16)  are  contained in (17).

In the D.E. giving the transformation  of stream BOD into D.O.
we introducted the  term  P^ representing all other sources
and sinks  of  D.O. in section i.  Now let P^ represent sources
of D.O.  and Q^ sinks of D.O.  in  section  i.   From previous re-
sults we have

        C  = $K -  S  V W - S"1?  + S-1Q

where P and Q are column  vectors.   The same basic idea used
to construct  a SS transfer matrix  for  BOD can be used to con-
struct a SS transfer matrix  for  D.O. sources.   We inject 1
unit of D.O.  into each section  and  calculate the corresponding
D.O. changes  in each section.  Define

               „-!
where ty-j is the  jth  column  of  the  transfer  matrix.   The  SS
transfer matrix Y  is  then

        Y = -S"1!  = -S"1

This SS transfer matrix can  be  used to  account for  D.O.  changes
due to the following:

         1)  D.O.  in  the polluter discharges
         2)  D.O.  in  the regional plant discharges
         3)  in-stream aeration
         4)  flow  augmentation.

For 1) and 2)  let  QJ-,  and  R,  be  the  average D.O.  concentration
of effluent from polluter  i  and regional plant i respectively.
Thus 1) and 2) are handled by  the following  modification to
the quality constraints

             kt..(l-r.)p  J       	     kt..(l-r.)p  J
     AFJ + A 2	=H=	~	 + C - AF J - A 2	
                            -  19  -

-------
In-stream aeration 3) requires the introduction of a variable
specifying the amount of D.O. added in each section due  to
aeration.  Flow augmentation requires a conservation of  flow
constraint for each reservoir and additional flow variables.

Part 2:   Numerical Calculation of a Transfer Matrix

A general nonlinear programming problem for water quality con-
trol in a tidal estuary has been given by Graves., Hatfield
and Whins ton [4].  The flow chart [4, page 4] of the type of
operations involved in the overall approach is given below:
                      Basic Data from
                      Continuous
                          Monitoring
                      Basic Data
                      (costs,, removals,
                       etc.)
Compute Transfer

    Matrix
Solve Nonlinear
Programming
       Problem
                       Adjust Goals,,
                       Change Plant
                       Locations,, etc.
                        Acceptable
                         Solution
                          FIGURE 2

   OPTIMIZATION PROCESS,, ASSUMING FIXED  TRANSFER MATRIX
                           - 20  -

-------
In [4] the SS  (Steady State) transfer coefficients were  taken
as a matrix of fixed constants and the path  labeled  1  -»  2  - 3
of Figure 2 was not investigated.  There was no  need to  re-
compute the SS transfer coefficients since they  were assumed
constant.  Unfortuantely the SS  transfer matrix  A, is  not  com-
pletely independent of the variables of the  nonlinear  program-
ming problem.  The calculation of the SS transfer matrix de-
pends, among other things, on the estuarial  intersection flows.
Activating by-pass pipes and switching plant discharge points
can change the estuarial flow pattern.  This suggests  some type
of iterative process of:   1)  solving the nonlinear  programming
problem;  2)   using the current  solution to  update the transfer
matrix; and  3)   re-solving  the  nonlinear programming  problem
using  the updated transfer matrix.  With these  ideas,  we can
modify the path  1-2  - 3  of Figure 2 as in  Figure 3.

The purpose  of this section  is  to develop  the  details  of the
calculations implied by Figure  3.  We begin  by considering the
calculation  of the  transfer  matrix.   Let   Qi-j^i be the net
flow  from  section i-1  to  i.  V±  is  the  volume  of section i.
 F-_1  •  is  the  advection coefficient  from  section i-1 to i.
Ei-l  i is  the  eddy exchange  coefficient from section i-1 to i.
d is*the  decay rate  of BOD in  all  sections.   ri is  the reaera-
 tion  rate  of D.O. in  section i.   For  computational purposes,
 a convenient system of units is  Qi-i  i  and Ei_]^i in km^/day,
V-  in km3 and  d  and r± in I/day.  Note  that  ?j_-i}i_  is dimen-
 sionless.                   L
                      Compute Transfer

                          Matrix
                      Solve Nonlinear
                      Programming
                              Problem
                          Update
                      Transfer Matrix
                       Adjust Goals,
                       Change Plant
                       Locations,  etc,
                         FIGURE 3

  ITERATIVE PROCEDURE FOR FINDING STABLE TRANSFER MATRIX

                            - 21 -

-------
Let R and S be square matrices  of  order  n (n = the number of
sections).  If  reR and seS and
then
         r. .  =  c. .  - V.d
          11     11     i
         s . .  =  c . .  - V .  r .
                -
and     rij = sij = °   if  'i -  J I 2 2


Let V be a diagonal matrix  of order n with the  V^ as the dia-
gonal elements.  Then  the SS transfer matrix corresponding to
a 1 Ib/day BOD  input into each  section  is  given by

 (20)    A = -(4.536 •  10"7) d  (R V~1S)~1

With the system of units chosen as above  the elements of R V  S
will be close to one and accuracy will  be  maintained in cal-
culating the inverse.   The  units of A are  in mg/L/lb/day.

In order to calculate  a SS  transfer matrix for  the Delaware
Estuary j we need numerical  values for the  Oi_i  i, 5i-i i.>
Ei_i<;j_, Vj_j ri and d.   These are in turn  functions of the
physical characteristics of the estuary.   The necessary data
on the physical properties  of the estuary  is given by Pence ,
Jeglic,, and Thomann [8, page 390] .  The values  of QI_I 1= Q(i)
in km3/day for a base  flow  of 3000 cfs  are given in Ta6le 1.
Values of E-J__J_  i in km^/day are calculated from  (see [8] equa-
tion 6)

 (21)   E      = 1.57873 ^
where KI j_+± is the eddy-diffusion coefficient from section i
to i+1. 'AI^I+I is the cross-sectional area between section i
and i+1;  and  i^_ is  the  length  of section i.   Values of
Ki,,i+lj Ai,,i+l and  "i are given  in [8,  page 390].  The EI_I 1=
E(i)  are tabulated in Table 2.   The advection  coefficient is
taken to be PI- 1,1= ^(i)  = 0.5 = i=5 , 6, . . . , 31.  In sections one
to four^ the advective forces are  stronger than the tidal action
                               22  -

-------

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
Q (i) =Q±_lfi
7.34055E-03
7.43108E-03
7.67577E-03
7.69779E-03
7.77119E-03
7.86907E-03
8.35844E-03
7.69779E-03
7.74672E-03
7.77119E-03
8.29727E-03
8.37067E-03
8.37067E-03
8.48323E-03
8.95057E-03
9.92932E-03
;
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.


1.04040E-02
1.07123E-02
1.08591E-02
1.09081E-02
1.09081E-02
1.13778E-02
1.13803E-02
1.13925E-02
1.13925E-02
1.14047E-02
1.14170E-02
1.14170E-02
1.15393E-02
1.16127E-02
1.16861E-02

                       TABLE 1
              NET FLOWS INTO SECTION  i

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
E (i) = Ei_1 ±
9.77309E-04
2.43355E-03
3.37848E-03
3.88367E-03
4.49938E-03
5.38346E-03
8.71458E-03
1.56610E-02
2.02077E-02
2.18654E-02
2.40361E-02
2.56543E-02
2.60490E-02
2.86934E-02
2.32336E-02
1.93394E-02

17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.


2.07010E-02
2.23784E-02
2.49241E-02
2.81803E-02
3.95734E-02
5.99522E-02
6.42937E-02
6.97798E-02
1.00399E-01
1.04101E-01
1.08521E-01
1.18523E-01
8.65670E-02
7.03126E-02
8.49277E-02

                        TABLE 2
EDDY EXCHANGE  COEFFICIENTS BETWEEN  SECTIONS i-1 AND i
                             - 23 -

-------
and we take F ( 1 ) = 0.86686    5 ( 2 ) =  0.76969    F ( 3 )  = 0.61619,
and  p(4 ) = 0.52767.  Values  of Vj_ in km3  are  taken from [8]
pate 390 and are recorded in Table 3.

The reaeration  rate r-j_ and the decay rate d are  functions of
water temperature.  This has the effect of  introducing a de-
gree of  freedom in the calculation of  the SS transfer matrix
since each matrix will depend  on the particular  value of water
temperature chosen.  Values of r-j_  in I/day  are  calculated from
(see [8] equations 8 and 9)


             24 (8.1 • IP"5 u±)°'5   .018  (t _  20)
         ri = -  Ti -  e
where u^ is the mean tidal velocity  in section  i  and  h^  is  the
mean depth in section i.  Values of  u^ and h^ are given  in  [8]
page 390.  A value of t=20° C=68° F  was chosen  because:
 (a)  it simplifies the calculation of  (22) ; and (b) it is
roughly the average water temperature from June to November.
The r-j_ are tabulated in Table 4.  Values of d in  I/day are
calculated from  (see [8] equation 10)

        d = 0.23 •  1.099(t "20)


With t=20° C we have d = 0.23/day-

The SS transfer matrix for the Delaware Estuary calculated
from (20) is given in Table 5.  The  SS transfer matrix supp-
lied to us by Prof. R.V. Thomann of  Manhattan College , and
used in our earlier work is given in Table 6 for  comparison.

Updating the Transfer Matrix:

Recall that there are five sets of variables  (the F,  P,  T,
J and R matrices) in the nonlinear programming  formulation.
It is relatively simple calculation  to update the intersec-
tional flows Q^ i+i* given the current P and T  matrices  and
the original polluter discharges.  If the base  flow in the
stream is large, changing the polluter discharge  points  effects,
primarily, the intersectional flows.  For  the Delaware Estuary
problem we are assuming that all other secondary  effects are
negligible.  Some of the other quantities  indirectly  effected
by changing the flow pattern are:  1)  the advection  coeffi-
cient li^i+i/  and 2)  the average depth h^ and  consequently
the reaeration rate r^.  This assumption is apparently justi-
fied for the Delaware Estuary; however, this is not always  the
case,  especially in situations where polluter discharges make
up a large percentage of the base flow.  A detailed analysis
of this particular aspect of the problem is given in  Section  V.
                            -  24  -

-------


1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
V.
i
6.85344E-03
1.03085E-02
1.30272E-02
1.50662E-02
1.80115E-02
2.14099E-02
1.28856E-02
1.42733E-02
1.50945E-02
1.64822E-02
1.78416E-02
1.85496E-02
1.96541E-02
2.27976E-02
5.26752E-02
5.74896E-02


17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.




6.18509E-02
6.78547E-02
7.62374E-02
8.30342E-02
4.28198E-02
4.45757E-02
4.80874E-02
5.07494E-02
5.23920E-02
5.44877E-02
5.81693E-02
6.36634E-02
1.38655E-01
1.59158E-01


           TABLE  3
    VOLUMES OF SECTION  i


1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
r .
i
1.00116E-01
1.19462E-01
9.23170E-02 .
1.00116E-01
1.15604E-01
1.25953E-01
1.06598E-01
9.87044E-02
1.15604E-01
8.55555E-02
8.94839E-02
5.71866E-02
4.56965E-02
8.39498E-02
8.94839E-02
1.12356E-01


17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.




1.30556E-01
1.12356E-01
1.12356E-01
1.12356E-01
1.21359E-01
9.93308E-02
1.21359E-01
1.21359E-01
1.41016E-01
1.95081E-01
1.78124E-01
1.78124E-01
1.29738E-01
1.20993E-01


           TABLE  4
REAERATION RATE IN SECTION i
               - 25 -

-------
I
to
                                                 TABLE 5
                                             TRANSFER MATRIX

1
2
3
4
5
6
7
8
9
10
li
12
13
14
15
16
17
16
19
20
21
22
23
24
25
26
27
28
29
30
1
-1. 0180-05
-1.778D-05
-2.380D-05
-2.594D-C5
-2.4C1D-C5
- .9340-05
- .7210-05
- .4940-05
- .3C2D-05
- .1360-05
-1.CC4C-05
-8.9C60-C6
-7.633C-06
-6.267C-06
-4.581C-06
-3.033C-06
-1.9450-06
-1.262C-06
-8.0370-07
-5.0710-07
-3.5970-07
-2.7680-07
-2.0600-07
-1.4860-07
-1.1310-07
-8.3920-08
-6.213D-08
-4.6210-08
-2.8770-08
-1.245C-08
2
-1.5880-06
-1.2440-05
-2.0900-05
-2.4890-05
-2.4120-05
-1.9980-05
-1.8C20-05
-1. 5770-05
-1.381D-05
-1.2110-05
-1.0750-05
-9.5610-06
-8.2150-06
-6.7590-06
-4.9540-06
-3.2890-06
-2. 1130-06
-1.3730-06
-8.7580-07
-5.5310-07
-3.9250-07
-3.0220-07
-2.2490-07
-1.6230-07
-1.2360-07
-9.1730-08
-6.7930-08
-5. 0530-08
-3.1470-08
-1.3620-08
3
-1.113C-07
-1.0960-06
-1.4370-05
-2.2060-05
-2.3690-05
-2.0750-05
-1.9190-05
-1.704D-05
-1.509C-05
-1.333C-05
-1.1910-05
-1.0650-05
-9.1940-06
-7.5910-06
-5.5910-06
-3.7260-06
-2.4020-06
-1.566C-06
-1.0010-06
-6.33CO-07
-4.4970-07
-3. 4640-07
-2.5800-07
-1.8630-07
-1.4190-07
-1.054C-07
-7. 8050-08
-5. 8070-08
-3.6190-08
-1. 5670-08
4
-3.7840-09
-4.0820-08
-6.9190-07
-1.5C6D-05
-2.1270-05
-2.0870-05
-2.0200-05
-1.8380-05
-1.6550-05
-1. 4820-05
-1.3380-05
-1.2060-05
-1.0480-05
-8.7010-06
-6.4540-06
-4.3280-06
-2.8C40-06
-1.8350-06
-1.1760-06
-7.4610-07
-5.3070-07
-4.0920-07
-3.0500-07
-2.2040-07
-1. 6790-07
-1.2480-07
-9.2490-08
-6.8850-08
-4.2930-08
-1.86CD-08
5
-2.6160-10
-2.9550-09
-5.5450-08
-1.5780-06
-1.5420-05
-1. 9510-05
-2. 0470-05
-1. 9380-05
-1.7900-05
- .6350-C5
- .4990-05
- .3660-05
- .1980-05
- .0020-05
-7.5CCD-06
-5.0710-06
-3.3C60-06
-2.1740-06
-1.4CCD-06
-8.9050-07
-6.344D-C7
-4.8980-07
-3.655D-C7
-2.6430-07
-2.0160-07
-1.499D-C7
-1. 1120-07
-8.2790-08
-5.1670-08
-2.240C-C8
6
-3.347C-11
-3.8S4C-IC
-7.686C-09
-2.4350-07
-3. 0530-06
-1.52CC-05
-1.9UC-05
-1.949C-05
-i.eecc-05
- .773C-05
- .6630-05
- .541C-C5
- .3690-05
- .157C-05
-8.776C-06
-5.999C-C6
-3.944C-06
-2.6120-06
-1.69CC-C6
-1.C8CC-06
-7.7C9C-07
-5.96CC-C7
-4.453C-07
-3.225C-07
-2.462C-07
-1.6320-07
-1.36CC-07
-1.C14C-07
-6.3320-08
-2.747C-C8
7
-1.0350-11
-1.216D-1C
-2.465C-C9
-8.153C-C8
-1.107D-C6
-6.618D-Ct
- .496C-C5
- .77CD-05
- .838C-C5
- .819C-C5
- .7630-C5
- .672C-C5
- .5110-05
-1.294C-C5
-9.976C-C6
-6.911D-C6
-4.59CC-C6
-3.0630-06
-1. 995C-C6
-1.28CC-06
-9.163C-C7
-7.096C-C7
-5.31CO-C7
-3.85CO-07
-2.942C-C7
-2.1920-C7
-1.628C-C7
-1.215C-07
-7.598C-C6
-3.299C-C6
8
-5.3750-12
-6.367C-U
-1.3C9C-09
-4.4350-08
-6.2670-07
-4.C470-06
- .C650-C5
- .5360-05
- .7240-05
- .7850-05
- .7780-C5
- .718C-C5
- .574C-05
- .3620-05
-1.C63C-05
-7.4350-06
-4.973D-C6
-3. 3380-06
-2.1830-06
-1.4050-06
-1.CC8C-06
-7.8HD-C7
-5.8510-07
-4.2470-C7
-3.2470-07
-2.42CO-07
-1.7990-C7
-1.3430-07
-e.4CRD-OR
-3.652C-C8
9
-3.14SC-12
-3.752C-11
-7.791C-1C
-2.684C-C8
-3.896C-C7
-2.635C-C6
-7.484C-C6
- .167C-C5
- .516C-C5
- ,68tC-C=
- .75CC-C5
- .735C-C5
- .618C-C5
- .419C-C5
- .124C-C5
-7.961C-C6
-5.371C-C6
-3.629C-C6
-2.385C-C6
-1.541C-Ot
-1.1C7C-C6
-g.595C-C7
-6.445C-C7
-4.683C-C7
-3.584C-C7
-2.673C-C7
-1.969C-C7
-1.485C-C7
-9.3C7C-CP
-4.046C-C8
1C
-1.872C-12
-2.241011
-4.693C-1C
-1.64CC-C8
-2.43CC-C7
-1.7CIC-C6
-5.CE6C-C6
-6.3C6C-C6
-1.16CC-C5
-1.492C-C5
-1.661C-C5
-1.714C-C5
-1.642C-C5
-1.467C-C5
-i.ieec-c?
-6.544C-Ct
-5.83CC-C6
-3.972C-C6
-2.627C-C6
-1.7C6C-C6
-1.228C-C6
-9.5E1OC7
-7.173C-C7
-5.219C-C7
-3.998C-C7
-2.985C-C7
-2.222C-C7
-1.661C-C7
-1.C42C-C7
-4.533C-C8

-------
  TABLE 5 Cont.
TRANSFER MATRIX













1
to
^]

1














1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
11
-1.1390-12
-1. 3690-11
-2.8880-10
-1.0200-08
-1.5380-07
-1.1C50-06
-3.4240-06
-5.7590-06
-8.3990-06
-1.1600-05
-1.4780-05
-1.6290-05
-1.6260-05
-1.4910-05
-1.2430-05
-9.1290-06
-6.3180-06
-4.3500-06
-2.8990-06
-1.8930-06
-1.3670-06
-1.C65D-06
-8.0120-07
-5.8380-07
-4. 4770-07
-3.346C-C7
-2.4940-07
-1.8660-07
-1.1720-07
-5.1C40-08
12
-7.0230-13
-8.4710-12
-1.7980-10
-6.4110-09
-9. 7900-08
-7.1740-07
-2.2820-06
-3.9180-06
-5.8740-06
-8.4670-06
-1.1540-05
-1.4470-05
-1.5450-05
-1.4770-05
-1.2830-05
-9.6870-06
-6. 8290-06
-4.7640-06
-3.2050-06
-2.1070-06
-1.5270-06
-1.1920-06
-8.9880-07
-6.5620-07
-5.0400-07
-3. 7710-07
-2.8140-07
-2.1C8D-07
-1.3260-07
-5.7800-08
13
-4.2450-13
-5.1360-12
-1.0960-10
-3.939C-09
-6.0840-08
-4.5310-07
-1.4710-06
-2.56CC-06
-3.9250-06
-5.8250-06
-8.2840-06
-1.113C-05
-1.3630-05
-1.3990-05
-1.2960-05
-1.0180-05
-7.3590-06
-5.2220-06
-3.5550-06
-2.3570-06
-1.7150-06
-1.3430-06
-1.0150-06
-7.4290-07
-5.7140-07
-4.2830-07
-3.2010-07
-2.40CC-07
-1. 5120-07
-6.6030-08
14
-2.5900-13
-3.1420-12
-6.7320-11
-2.4360-09
-3. 7980-08
-2.866D-C7
-9.4580-07
-1.6680-06
-2.5900-06
-3. 9220-06
-5.7360-06
-8.0310-06
-1.0550-05
-1.2380-05
-1.2680-05
-1.0520-05
-7.8490-06
-5.6870-06
-3.9270-06
-2.6290-06
-1.9220-06
-1.5100-06
-1.1450-06
-8.3970-07
-6.4710-07
-4.8580-07
-3.6360-07
-2.7300-07
-1.7240-07
-7.5370-08
15
-1.2950-13
-1.5770-12
-3.3980-11
-1.2400-09
-1.9560-08
-1.5COD-07
-5.0460-07
-9. 0210-07
-1.4240-06
-2.2040-06
-3.3170-06
-4.8340-06
-6.7510-06
-8.7280-06
-1.1570-05
-1.0710-05
-8.4430-06
-6.3280-06
-4.4670-C6
-3.0350-06
-2.2350-06
-1.7640-06
-1.343C-06
-9.8860-07
-7.6380-07
-5.7480-07
-4.3130-07
-3.2450-07
-2.0550-07
-9.0010-08
16
-5.2980-14
-6.468C-13
-1.4CCO-11
-5.14SO-1C
-8.197C-09
-6.3670-08
-2.172C-07
-3.923C-07
-6.2650-07
-9.849C-07
-1.511C-06
-2. 2590-06
-3.268C-C6
-4.442C-C6
-6.529C-06
-9.7590-06
-9.C68C-06
-7.3S9C-06
-5.49CC-06
-3.851C-06
-2.878C-06
-2.2940-06
-1.760C-06
-1.3C50-06
-1.014C-06
-7.666C-07
-5.777C-C7
-4.3630-07
-2.7780-07
-1.222C-07
17
-2.1420-14
-2.6200-13
-5.693C-12
-2.104D-1C
-3.372C-CS
-2.642C-08
-9.1050-C8
-1.6550-C7
-2.664D-C7
-4.2290-07
-6.568D-C7
-9.9650-C7
-1.471C-C6
-2.054C-C6
-3.173C-C6
-5.5390-C6
-8.592D-C6
-8.291C-C6
-6.676D-C6
-A.906C-C6
-3.743C-06
-3.022C-C6
-2.343C-C6
-1.755D-C6
-1.372D-C6
-1.0440-06
-7.90SO-C7
-6.002C-C7
-3.8460-07
-1.7COO-07
18
-8. 6760-15
-1. 0630-13
-2.316D-12
-8. 5920-11
-1.3840-09
-1. 0920-08
-3.7890-08
-6.S210-06
-1.1200-07
-1.791D-C7
-2.8040-07
-4. 2980-07
-6.4320-07
-9.1310-07
-I. 4520-06
-2.7440-06
-5.C26D-06
-8.1260-06
-7.6740-06
-6.C85D-C6
-4.7870-06
-3.9370-06
-3.1COC-06
-2.3520-06
-1.R550-06
-1.4220-06
-1.C850-06
-8. 2070-07
-5.3550-07
-2.3810-07
IS
-3.5S9C-15
-4.415C-14
-•3.639C-13
-3.586C-11
-5.8CCC-1C
-4.598C-C9
-1.6C4C-C8
-2.941C-C8
-4.779C-C8
-7.t7SC-Ce
-1.21CC-C7
-1.868C-C7
-2.621C-C7
-4.C5CC-C7
-6.5ttC-C7
-1.3CCC-C6
-2.586C-C6
-4.916C-C6
-7.596C-C6
-7.C15C-C6
-5.817C-C6
-4.93CC-C6
-3.S72C-C6
-3.073C-C6
-2.454C-C6
-1.9C3C-C6
-1.467C-C6
-1.129C-C6
-7.381C-C7
-3.3C8C-C7
2C
-1.554C-15
-1.9C8C-14
-4.172C-13
-1.5560-11
-2.524C-1C
-2.CC8C-C9
-7.C36C-C9
-1.293C-C8
-2.1C8C-C8
-3.359C-C8
-5.376C-C8
-8.344C-C8
-1.268C-C7
-1.836C-C7
-3.C15C-C7
-6.152C-C7
-1.283C-C6
-2.635C-C6
-4.744C-C6
-6.S6SC-C6
-6.4690-Ct
-5.7S3C-C6
-4.857C-C6
-3.876C-C6
-3.158C-C6
-2.49CC-C6
-1.948C-C6
-i.siec-ce
-1.CC8C-C6
-4.566C-C7

-------
to
00
                                              TABLE  5  Cont.

                                            TRANSFER MATRIX

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IS
16
17
18
19
20
21
22
23
24
25
26
27
26
29
30
21
-8.4490-16
-1.0380-14
-2.272D-13
-8.4820-12
-1.3780-10
-1.C98C-09
-3.857D-09
-7.C970-C9
-1.1590-08
-1.872D-08
-2. 9670-08
-4.6170-08
-7.0410-08
-1.C23C-07
-1.6SOO-C7
-3.4S8D-C7
-7.447C-07
-1.5780-06
-2. 9940-06
-4.9220-06
-6.CC5C-06
-5.974C-06
-5.3540-06
-4.4820-06
-3.757C-C6
-3.C31C-06
-2.418C-06
-1.9130-06
-1.2950-06
-5.945C-07
22
-5. 3980-16
-6.6350-15
-1.4530-13
-5.4280-12
-8.8270-11
-7.0460-10
-2.4770-09
-4.5630-09
-7.4570-09
-1.2C6D-08
-1.915D-08
-2.985D-08
-4.5620-08
-6.6440-08
-1. 1020-07
-2.3C2D-07
-4.9650-07
-1.0720-06
-2. 0950-06
-3.6370-06
-4.9300-06
-5.6350-06
-5.4300-06
-4.7600-06
-4.0940-06
-3.3690-06
-2.7310-06
-2.1870-06
-1.5C30-06
-6.9700-07

-3
-4
-9
-3
-5
-4
-I
-2
-4
-7
-1
-1
-2
-4
-7
-1
-3
-7
-1
-2
-3
-4
-5
-4
-4
-3
-3
-2
-1
-8
23
.4300-16
.2180-15
.24CC-14
.455C-12
.624C-11
.4940-10
.5820-09
.9160-09
.7700-09
.7240-09
.2280-08
.9170-08
.934C-08
.2840-08
.131C-08
.5CCC-07
.271C-07
.1720-07
.4340-06
.5880-06
.745C-06
.6000-06
.0990-06
.8230-06
.308C-06
.6460-06
.0200-06
.4580-06
.7220-06
.083C-07
24
-2.1620-16
-2.6590-15
-5.8280-14
-2.1810-12
-3.5520-11
-2.8420-10
-I. 0010-09
-1.8470-09
-3.0240-09
-4.9010-09
-7. 8000-09
-1.2190-08
-1.8700-08
-2.7340-08
-4.5650-08
-9.6690-08
-2.1280-07
-4.7230-07
-9.6160-07
-1.7880-06
-2. 7080-06
-3.4780-06
-4.1510-06
-4.5430-06
-4.3210-06
-3.8140-06
-3.2570-06
-2.7090-06
-1.9470-06
-9.2770-07
25
-1.5220-16
-1.8720-15
-4.1C5C-14
-1.5370-12
-2.5040-11
-2.CC50-10
-7.0710-10
-1.305C-09
-2.137C-09
-3.467D-C9
-5. 5210-09
-8.6370-C9
-1.3260-08
-1.9420-08
-3.2490-08
-6.91CD-08
-1.5300-C7
-3.4230-07
-7.0500-07
-1.335D-C6
-2.C75C-06
-2.7280-06
-3.3750-06
-3.9250-06
-4.0830-06
-3. 8020-06
-3.3650-06
-2.8670-06
-2. 1150-06
-1.0240-06
26
-I.C62C-16
-1.3C6C-15
-2.8t5C-14
-1.073C-12
-1.749C-11
-1.4C1C-10
-4.946C-1C
-9.131C-1C
-1.496C-09
-2.428C-09
-3.869C-09
-6.05EC-C9
-9.3C9C-C9
-1.365C-08
-2.2E7C-08
-4.ee3C-08
-1.C87C-07
-2.448C-07
-5.089C-07
-9.772C-C7
-1.55CU-06
-2.C73C-06
-2.63CC-06
-3.180C-06
-3.483C-C6
-3.584C-C6
-3.367C-06
-2.975C-06
-2.28CC-C6
-1.127D-06
27
C.O
-9. 0050-16
-1.975C-14
-7.399C-13
-1.2070-11
-9.675C-11
-3.416C-1C
-6.31CC-IC
-1.034C-C9
-1.679C-C9
-2.678C-C9
-4.1950-C9
-6.452C-CS
-9.468C-09
-1.5890-C8
-3.404C-C8
-7. 6070-08
-1.724C-C7
-3.612C-C7
-7.017D-C7
-1.131C-C6
-1.533C-C6
-1.9810-Ct
-2.463C-C6
-2.7880-C6
-3.041C-C6
-3.1840-C6
-2.984U-C6
-2.418C-C6
-1.23CD-C6
28
C.C
-6.2130-16
-1.3630-14
-5.1080-13
-6. 3350-12
-6.6840-11
-2.3610-10
-4.3630-1C
-7.1540-10
-1.1620-09
-1.854C-C9
-2.-5C60-09
-4. 4720-09
-6. 5680-09
-1.1C4D-08
-2.371C-C8
-5. 3170-08
-1.2100-07
-2.552C-07
-5.CC50-07
-8.165C-C7
-1.1180-06
-1.4650-C6
-1.8560-06
-2.149U-06
-2.4290-C6
-2.6S8D-06
-2.B110-06
-2.4800-06
-1.3120-06
29
C.C
-3.3K8C-16
-7.434C-15
-2.787C-13
-4.55CC-12
-3.651011
-1.291C-1C
-2.386C-1C
-3.914C-1C
-6.361C-1C
-1.C15C-CS
-1.593C-C9
-2.454C-C9
-3.6C8C-C9
-6.075C-C9
-1.31CC-C8
-2.952C-C8
-6.762C-C8
-1.43RC-C7
-2.856C-C7
-4.7??C-C7
-6.566C-C7
-8.745C-C7
-1.134C-Ct
-1.347C-C6
-1.581C-C6
-1.86CC-C6
-2.116C-C6
-2.372C-Cf
-1.373C-C6
3C
C.C
-1.234C-16
-2.7C9C-ie
-1.C16C-13
-1.659C-12
-1.332C-U
-4.712C-11
-8.712C-11
-1.43CC-1C
-2.325C-1C
-3.713C-1C
-5.829C-1C
-8.986C-1C
-1.322C-C9
-2.22SC-C9
-4.819C-C9
-i.cscc-ce
-2.51CC-C8
-5.37ic-C6
-1.C76C-C7
-1.8C3C-C7
-2.523C-C7
-3.398t-C7
-4.471C-C7
-5.3S3C-C7
-6.476C-C7
-7.859C-C7
-9.333C-C7
-1.148C-C4
-1.1C3C-C6

-------
                                               TABLE  6
                                      THOMANN'S TRANSFER MATRIX
to

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
2^
24
25
26
27
28
29
30
1
7.834E-06
1.3872-05
1.7712-05
1.806E-05
1.616E-05
1.159E-05
9.7612-06
8.196E-06
6.6332-06
5.3832-06
4.5712-06
3.o$3E-u6
3.3S42-06
2.779E-06
2.016E-06
1.382E-06
9.3522-07
6.2032-07
4.028E-07
2.577E-07
1.835E-0?
1.3852-07
9.9232-08
6.T71E-08
4.7502-08
3.218E-08
2.1872-08
1.5122-08
8.650E-09
3.7662-09
2
9.331E-07
1.0692-05
1.706S-05
1.8922-05
1.7682-05
1.3022-05
i.noE-05
9.3892-06
7.643E-06
6.2342-06
5.311E-06
4.6502-06
3.951E-06
3.2502-06
2.3622-06
1.6212-06
1.0982-06
7.2912-07
4.7362-07
3.0312-07
2.1592-07
1.6292-07
1.1673-07
7.9672-08
5.5892-08
3.786E-08
2.574E-08
1.7802-08
1.0192-08
4.4322-09
3
5.7672-08
8.3562-07
1.2402-05
1.7702-05
1.8342-05
1.428E-05
1.2482-05
1.0702-05
8.8082-06
7.2512-06
6.2172-06
5.4692-06
4.6642-06
3.848E-06
2.8062-06
1.9312-06
1.311E-06
8.7102-07
5.6632-07
3.6252-07
2.5832-07
1.9502-07
1.397E-07
9.537E-08
6.6912-08
4.533E-08
3.081E-68
2.131E-08
1.2202-08
5.308E-09
4
1.3432-09
2.1302-08
4.084E-07
1.2742-05
1.7502-05
1.528E-05
1.399E-05
1.2292-05
1.0302-06
8.6122-06
7.4592-06
6.611E-06
5.6722-06
4.701E-06
3.447E-o6
2.3812-06
1.62UE-06
1.079E-06
7.0232-07
4.4992-07
3.2072-07
2.4222-07
1.7362-07
1.1852-07
8.313E-08
5.633E-08
3.8292-08
2.6482-08
1.5172-08
6.599E-09
5
3.271E-11
5.439E-10
1.1572-08
4.770E-07
1.330E-05
1.5292-05
1.5232-05
1.393E-05
1.2032-05
1.0292-05
9.0472-06
8.1082-06
7.0152-06
5.850E-06
4.3212-06
3.0022-06
2.051E-06
1.3682-06
8.9212-07
5.722E-07
4.0812-07
3.082E-07
2.2092-07
1.5082-07
1.059E-07
7.1742-08
4.8782-08
3.3742-08
1.933E-08
8.4102-09
6
9.865E-1.3
1.685E-11
3.7782-10
1.7382-08
6.325E-07
1.2272-05
1.5192-05
1.510E-05
1.378E-06
1.2272-05
1.1052-05
1.0082-06
8.8332-06
7.4372-06
5.555E-06
3.8912-06
2.672E-06
1.7892-06
1.169E-06
7.509E-07
5.359E-07
4.049E-07
2.904E-07
1.983E-07
1.392E-07
9.436E-08
6.417E-08
4.439E-08
2.5442-08
1.107E-08
7
2.0572-13
3.557E-12
8.166E-11
3.921JS-09
1.548E-07
3.7292-06
1.199E-05
1.4352-05
1.4462-05
1.371E-05
1.2802-05
1.1962-05
1.0662-05
9.0912-06
6.8872-06
4.873E-06
3.368E-06
2.265E-06
1.484V -06
9.549E-07
6.8202-07
5.1562-07
3.6992-07
2.5272-07
1.7742-07
1.2032-07
8.1842-08
5.663E-08
3.246E-08
1.413E-08
8
8.947E-14
1.559E-12
3.6262-11
1.782E-09
7.3262-03
1.924E-06
7.4892-06
1.221E-05
1.3832-05
1.3962-05
1.3462-05
1.283E-05
1.161E-06
9.9942-06
7.6552-06
5.4582-06
3.7902-06
2.5572-06
1.679E-06
1.0822-06
7.731E-07
5.8462-07
4.196E-07
2.8672-07
2.014E-07
1.3662-07
9.2912-08
6.4302-08
3.6872-08
1.606E-08
9
3.855E-14
6. 758E-13
1.589E-n
7-953E-10
3.3692-08
9.3722-07
4.044E-06
7.4092-06
1.177E-05
1.3562-05
1.386E-05
1.369E-06
1.2672-05
1.108E-05
8.630E-06
6.224E-06
4.3532-06
2.9502-06
1.9432-06
1.254E-06
8.970E-07
6.788E-07
4.873E-07
3.3322-07
2.341E-07
1.5882-07
1.0802-07
7.4792-08
4.290E-08
1.8692-08
10
1.687E-ll»
2.973E-13
7.0572-12
3.5872-10
1.557E-08
4.524E-07
2.089E-06
4.08UE-06
7.418E-06
1.1912-05
1.356E-06
1.416E-05
1.3562-05
1.213E-05
9.6722-06
7.0832-06
5.0002-06
3.4082-06
2.2532-06
1.4582-06
1.044E»06
7.9052-07
5.679E-07
3.8852-07
2.730E-07
1.852E-07
1.2612-07
8.7322-08
5.011E-08
2.184E-08

-------
                                            TABLE 6 Cont.

                                      THOMANN'S TRANSFER MATRIX
CO
o

1
2
3
4
5
6
7
8
9
10
n
12
13
Ik
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
n
8.9162-15
1.5782-13
3.7712-12
1.938E-10
8.555E-09
2.55TE-OT>
1.229E-06
2.4862-06
4.8102-06
8.6342-06
1.213E-06
1.3822-05
1.3892-05
1.278E-05
1.049E-05
7.8212-06
5.5602-06
3.6292-06
2.5422-06
1.6492-06
1.1822-06
8.9602-07
6.4402-07
4.4082-07
3.0992-07
2.1032-07
1. If 322-07
9.9232-08
5.6972-08
2.484E-08
12
4.8462-15
8.6052-14
2.0692-12
1.073E-10
4.8002-09
1.4662-07
7.2552-07
1.5022-06
3.0232-06
5.7672-06
8.8722-06
1.2352-05
1.3562-05
1.3112-05
1.1242-05
8.6032-06
6.2312-06
4.315E-06
2.883JE-06
1.8772-06
1.3472-06
1.0222-06
7.3522-07
5.036E-07
3.5422-07
2.4062-07
1.6392-07
1.1362-07
6.5262-03
2.847E-08
13
2.5982-15
4.6272-14
1.1182-12
5.8432-11
2.6442-09
8.2202-08
4.163J3-07
8.767^-07
1.8142-06
3.6002-06
5.8342-06
8.9092-06
1.2042-05
1.2742-05
1.1752-05
9.361E-06
6.9282-06
4.8592-06
3.271E-06
2.1412-06
1.5402-06
1.1702-06
8.4252-07
5.7762-07
4.0662-07
2.7632-07
1.6042-07
1.3062-07
7.5122-08
3.2792-08
14
1.4232-15
2.5402-14
6.1612-13
3.2412-11
1.4802-09
.4.6662-08
2.4032-07
5.1282-07
1.0822-06
2.2052-06
3.6902-06
5.9302-06
8.7822-06
1.1302-05
1.1802-05
9.9652-06
7.5932-06
5.4152-06
3.683JE-06
2.4242-06
1.7492-06
1.3312-06
9.5972-07
6.5882-07
4.641E-07
3.1572-07
2.1542-07
1.4952-07
8.6062-08
3.7602-08
15
6.1222-16
1.0962-14
2.6742-13
1.4182-n
6.5502-10
2.10XE-08
1.1042-07
2.3902-07
5.1522-07
1.0802-06
1.8652-06
3.l4i2-o6
5.0062-06
7.2872-06
1.0962-05
1.0472-05
8.421E-06
6.1822-06
4.2742-06
2.8442-06
2.0612-06
1.5722-06
1.1362-06
7.8142-07
5.5122-07
3.75^2-07
2.5652-07
1.7822-07
1.028E-07
4.497E-08
16
2.0852-16
3-7442-15
9.1702-14
4.8942-12
2.2802-10
7.4062-09
3.9502-08
8.640E-08
1.8902-07
4.0322-07
7.1052-07
1.2292-06
2.0372-06
3.1392-06
5.3282-06
9.8172-06
9.43IE-06
7.5002-06
5.4132-06
3.6912-06
2.7032-06
2.0752-06
1.5072-06
1.041E-06
7.3682-07
5.0332-07
3.4492-07
2.4032-07
1.391E-07
6.106E-08
17
7.2022-17
1.2952-15
3.1822-14
1.7052-12
7.9922-11
2.6172-09
1.4092-08
S.ioiE-08
6.8422-08
1.4762-07
2.6302-07
4.6202-07
7.8152-07
1.2382-06
2.2142-06
4.8542-06
9.1082-06
8.643E-06
6.7472-06
4.7942-06
3.5682-06
2.7662-06
2.025E-06
1.4o82-06
i.
-------
                                            TABLE  6  Cont.
                                      THOMANN'S TRANSFER MATRIX
U)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
21
1.8292-18
3. 2982-17
8.1372-16
4.3902-14
2.0752-12
6.8812-11
3-7552-10
8.3402-10
1.8642-09
4.0792-09
7.3842-09
1.3232-08
2.2972-08
3.7602-08
7.1182-08
1.8202-07
4. 6602-07
1.1322-06
2.477E-06
4.6';5^-o6
6.0742-06
5.9082-06
5.0202-06
3.8992-06
2.9793-06
2.1YT2-06
1.5862-06
1.1642-06
7.2542007
3.360E-07
22
1. 1042-18
1.9912-17
4.9132-16
2.6512-14
1.2542-12
4.1602-11
2.2722-10
5.0462-10
1.129E-C9
2.4732-0?
4.4792-09
8.032E-0?
1.396E-OS
2.289E-0?
4.344E-Cc
1.117E-07
2.8882-07
7.1232-07
1.602E-C-5
3-1782-C5
4 . 6992 -C-S
5.51.42-cS
5.1192-0?
4.1962-3;
3.306E-C=
2.475E-C-S
I.84i2-Co
1.3732-Cc
8.7352-:?
4.103E-07
23
6.5692-19
1.1892-17
2.934E-16
1.58*12-14
7.4952-13
2.4872-11
1.3592-10
3.0212-10
6.7592-10
1.1) 812-09
2.6852-09
4.8182-09
8.3822-09
1.376E-08
2.6162-08
6.7582-08
1.7592-07
4.3902-07
1.0072-06
2.0742-06
3.2902-06
4.2042-06
4.7472-06
4.2832-06
3.5432-06
2.7502-06
2.1052-06
1.6052-06
1.0482-06
5.0092-07
24
3.871)2-19
6.9892-18
1.7252-16
9.3162-15
4.4092-13
1.46iiE-ll
8.0022-11
1.7802-10
3.9822-10
8.7322-10
1.5832-09
2.8432-09
4.9502-09
8.1332-09
1.5^92-08
4.0152-08
1.0512-07
2.6462-07
6.1592-07
1.302E-06
2.1622-06
2'. 9012 -06
3.5872-06
3.9982-06
3.606E-06
2.9622-06
2.3632-06
1.8562-06
1.2552-06
6.1282-07
25
2.5422-19
4.5862-18
1.1322-16
6.1142-15
2.8942-13
9.6132-12
5.2562-11
1.1692-10
2.6172-10
5.7392-10
1.0412-09
1.8702-09
3.2572-09
5.3542-09
1.0212-08
2.6522-08
6.9662-08
1.7632-07
4.1412-07
8.8902-07
1.5152-06
2.0842-06
2.6892-06
3.2472-06
3.3872-06
3.0142-06
2.5332-06
2.0592-06
1.4452-06
7.2122-07
26
1.6442-19
2.9662-18
7.3232-17
3-955E-15
1.8732-13
6.2202-12
3.4022-11
7.5672-11
1.6942-10
3.717E-10
6.743E-10
1.2122-09
2.1112-09
3.4722-09
6.6232-09
1.7242-08
4.5412-08
1.1542-07
2.7302-07
5.9312-07
1.0292-06
1.44 22-06
1.9132-06
2.425E-06
2.7222-06
2.8662-06
2.6362-06
2.2592-06
1.6722-06
8.5872-07
27
1.0482-19
1.8922-18
4.6712-17
2.5232-15
1.1952-13
3.9692-12
2.1712-11
4.8302-11
1.0822-10
2.3732-10
4.3072-10
7.7412-10
1.3492-09
2.2202-09
4.2372-09
1.1042-08
2.9162-08
7.4382-08
1.7692-07
3.8802-07
6.8332-07
9.7012-07
1.3142-06
1.7202-06
2.0202-06
2.3182-06
2.5632-06
2.3992-06
1.9172-06
1.0232-06
28
6.7342-20
1.2152-18
3.0012-17
1.6212-15
7.6762-14
2.5512-12
1.3952-11
3.1052-11
6.9532-11
1.5262-10
2.7692-10
4.9792-10
8.6602-10
1.4292-09
2.7282-09
7.1202-09
1.8832-08
4.8172-08
1.1512-07
2.5432-07
4.5282-07
6.4942-07
8.9242-07
1.1952-06
1.44 52-06
1.741)2-06
2.1062-06
2.3522-06
2.1232-06
1.193E-06
29
3.3162-20
5.9832-19
1.4782-17
7.9832-16
3.7812-14
1.2572-12
6.8752-12
1.5302-11
3.4272-11
7.5232-11
1.3362-10
2.4562-10
4.2832-10
7.0522-10
1.3462-09
3.5232-09
9.3422-09
2.3992-08
5.7672-08
1.2872-07
2.3262-07
3-3792-07
4.7312-07
6.5132-07
8.1432-07
1.0362-06
1.3562-06
1.7202-06
2.2472-06
1.4162-06
30
1.177E-20
2.1242-19
5.246E-16
2.8352-16
1.3432-14
4.4632-13
2.4422-12
5.4362-12
1.2182-11
2.6742-11
4.8542-11
8.7312-11
1.5232-10
2.5092-10
4.7972-10
1.2562-09
3.3362-09
8.5932-09
2.0752-08
4.6652-08
8.521)2-08
1.2502-07
1.7722-07
2.4852-07
3.1732-07
4.1652-07
5.6912-07
7.6552-07
1.1292-06
1.4152-06

-------
 After updating the intersection flows to obtain a new set of
 Q-i i+1 t^ie updated transfer matrix is obtained from  (20) .

 Steady State Stability of the Transfer Matrix;

 Definition:   Let a-j_j  be any element of A and let a^j be the
 corresponding element of A1  (updated transfer matrix) .  Sta-
 bility of A (or A1)  implies that for all i, j   [a-J^-a^j !<16/dmax
 and
 In the  above  definition f, is an arbitrary parameter chosen to
 reflect the accuracy of calculation of DO changes and dmax is
 some  arbitrary BOD load.   To motivate the definition let
 €ij = aij-aij.   Intuitively we would like stability to imply
 that  £j_j  is so small that changes in DO in section i are negli-
 gible when some maximum BOD load, dmax,  is applied to section j.
 Theref ore,, we  want


 (24)     ^ij

 Since dmax>0,  6^>0 ;  and ^ij  = aij~aij we  have

 (25)     JCij|  = !a!.-a..]   6/d      for all
 In practive we  generally take  dmax equal to the maximum pollu-
 ter  input  in  all  sections.

 Experimental  Results:

 The  ability to  calculate the  transfer  matrix when needed adds
 quite a bit of  power and flexibility to  the water quality con-
 trol algorithm.   The first problem we  investigated was to re-
 solve the by-pass piping problem  using the  new transfer matrix
 given in Table  5.  Recall that this case has linear constraints
 A solution of the by-pass piping  problem using the transfer ma-
 trix of Table 6 can be  found  in [5]  (see Figure 8 page 31) .
 The new solution, with  no updating of  the transfer matrix,  is
 plotted in Figure 4 in  a manner similar  to  that of Figure 8 in
 reference  [5]  .  Upon comparison of the two  solutions we note
 that:  a)  the  overall  flow patterns are similar; and  b)  the
 final costs after eliminating  duplicate  piping are almost the
 same ($2,089,148. as compared  to  $2,070,150.  see Figure 9 of
 [5]).  Note that as we  are now applying  the piping cost equa-
 tion the unadjusted solution of $2,447,995.  is more indicative
 of the cost with consolidated  piping.

 The next problem we investigated  was to  update the transfer
matrix and continue the  solution  carrying each by-pass piping
problem to the optimum.   The optimal solutions obtained were
                             - 32 -

-------
   ^ X/^wWPV^/ /
 *—-  x:—/s'1^  i '
L_ --IZII" jL-il*""_^//
                                             Optimal Solution  f2,447,99s.
                        FIGURE  4


       SOLUTION OF  BY-PASS  PIPING PROBLEM

          USING  NEW  FIXED TRANSFER MATRIX
                              - 33  -

-------
               transfer matrix           cost

                  (Table 5)            $2,447,995.

                1st update            $2,436,462.

                2nd update            $2,436,573.

                3rd update  (stable)

The total computer running time was 132 seconds with one re-
start at 60 seconds (some calculating time was lost since the
restart occurred while solving a by-pass piping problem).

The final solution is plotted in Figure 5.  We can conclude
that for the by-pass piping problem on the Delaware Estuary
updating the transfer matrix is of little or no importance.
Since pure by-pass piping is the most extreme case of re-
arranging the flow pattern the above conclusion holds for any
problem configuration.
                           -  34  -

-------
                                  section    D.o. goal
                                   1        0
                                   1        1
                                   1
                                   14
                                   15
                                   It
                                  optiMl Solution ¥2,436,573.
                FIGURE  5

SOLUTION OF BY PASS  PIPING  PROBLEM
USING NEW VARIABLE  TRANSFER MATRIX
                 -  35 -

-------
    V.  ESTIMATION OF THE DEPTH AND VELOCITY OF A RIVER

                 FOR A GIVEN  RATE  OF  FLOW

Main reference  for this  study has  been:   F.M.  Henderson,
Open Channel Flow  (see  [6]).   No attempt has been made to
make specific references and  for further clarification re-
garding the physics of  the  problem,  the  reader is referred
to the above text or to  its equivalent.   The motivation for
this study stems from the general  equation for the behavior
of dissolved oxygen in an estuary  as  suggested by Thomann
[10],  One way  to improve the DO level  is by controlling  the
flow in the estuary which in  turn  determines the depth and
velocity.  The  rate of flow appears explicitly in Thomann's
equation.  The  depths are implicit in the volumes of  the  vari-
ous sections.   If during the  solution process  of the  general
model,, large changes of  flow  are generated,  the matrix of trans-
fer coefficients has to  be  updated to correspond to the new
flow variables  which in  turn  requires estimation of the new
depths resulting from the new rate of flow.

Four cases to be considered in this chapter  are:

         A.  No inflows  or  outflows.
         B.  There are spot inflows  (pipes) .
         C.  There are spot outflows  (pipes).
         D.  There is a  continuous inflow (tributary).

The underlying  assuptions in  all these cases  are:

         1.  Water level does not  vary across  any cross^section.
         2.  Velocities  of  flow are low  (or  equivalently—slopes
             are small)  and the steady state velocity is  sub-
             critical.
         3.  Velocity of flow is uniformly distributed in any
             cross-section.
         4.  Rectangular cross-section.

Additional assumptions do not apply to all four cases and they
are listed in their appropriate sections.

A.  Regular flow - fixed rate of flow in the whole sections

        The Bernoulli equation:


                   H = y +  z  +-
applies to small slopes.

        H =  total head

        y = depth of water  in  the particular  cross^section




                           - 37 -

-------
         v = velocity of flow (average over the cross-section)
         z = height of river bed above reference line
         g = gravitational constant
         For two consecutive cross-sections:
         H.  = H. ,-,  + i.  . .,, 1=0,..., n  and i=0 at the beginn-
          1    1"T"J_     1^ 1 i X
 ing of the first section and i=n  at the mouth of the river.
         n = number of sections into which the river is divided.

         !L^ i+1 = losses between cross-section i and i+1.

         Applying the Bernoulli formula to H.  and H.    gives:
                   v2                  2    1
         Yi + Zi +  2? =
         Defining:   S^  = average  slope  of river bed at section
                          = length of (i+l)st section
 leads  to:
        Substituting gives:,,
        Eliminating v.^ and vi+1 by observing  that
        Q = vA = vby

whe re:
        A = cross-sectional area
        b = width of cross-section^
results in:
                           -  38 -

-------
 or:

                                                       +          - 0.
 The loss function:
                                              2
          Ai,i+l =  (Sf} average   ' ^i+1  + CL \2g) average
                       over                       over
                    section  i+1               section  i+1

 where :

          Sf = friction slope = (-^J
                                          "R

          n = Manning's resistance coefficient

          R = hydraulic mean  radius = A/P

          P = wetted perimeter

          A = Cross-sectional area

         c  = Eddy  coefficient  =  2— for angles between  90
          LJ                        C
              and 180°

         r  = radius of curvature of center-line
          c

The arithmatic average used to  calculate Sf is a crude approxi-
mation but this is  also the nature of the whole loss expression.
The estimate of CL  is biased  upwards severely in most cases,
but in the calculation here,  the  model does not take into account
bridges or islands.  In addition, the larger the river,  the lar-
ger the Eddy Coefficient for  the  same b/rc.

Substituting gives:
                           29
                          TT^          —   O
                {—

                R

            2  '2g    2g


Substituting R = ^ = .^  and v =
                            - 39

-------
            1   /fin  ^ AV       r  (b,+2y.)4/3
          — 	 •  I r*i   I   AX       I   J-    J_       T~
          — i   VT  /id'   **"• •  -i    i-        •*•	
            2    1.4y     !+l           in/-3
                                    (by).
          + 	TTT-TT- ]  +     1+1    \-—^ + 	Y~'
             (by)^{3         4g       {by)i   (by)i+l


To approximate c_ for angles smaller  than  90 :
         -      L    b +b. ,          b.+b.^n
        2b   _   2    . i  i+l) ._0_  =   i  i+l -.  ff
   °-'-    -ri+1(  ^      ,0°      ^  90o
where:  3 = angle of curvature  (in degrees).
            T  ( Qn  )2-AX       [   1    (^_ + -L)
            2            i+l            b.   y.
                                  )4/3]  -
      .
by i+l  b


                   (^)  +  UJ
                       i+1     i+1


                                            2
              4gri,i+l     9°    bY i     bY 1+1

Equations  (26) and  (27) define  a recursive relation between
y^ and y^+i-  The initial  condition required to solve this
recursive  set of equations is obtained by assuming that at
the mouth  of the river

           yn = CONST.

y  has to  be known  to begin the solution process.

Equations  (26) and  (27) have  to be solved by using a numeri-
cal method.  To update  the longitudinal profile after a vec-
tor of flow-changes has been  induced,  there are two possibi-
lites :
         a.  Identify the  new flows in every section and re-
             solve  equations  (26)  and  (27) .
         b.  Approximate  the changes  in the y  values as
                          - 40 -

-------
Looking at equations  (26)  and  (27)  as  components of a vector
function which implicitly  defines
    ':   >
and using the implicit function  theorem^leads  to:

   Equation (26)   F, [ (
   Equation  (27)  F2[(y±+l>  Q) > (yi^i, i+l) ]  =0
or:
         = 0 = F
so that:
                                                      BQ
Only the first row of the matrix  is needed  to  obtain  yi  ,
since:
                 dy
                   i+1
   dQ     Byi+1    dQ
   [F1
                                                      dQ
                           - 41 -

-------
       = A =
            -
          = B = -y
                 *
           2  1.
                                                 +
        Q-n
1.49b.y.
                     y.
      b.y.
       i i
                              90
   dF,
          = 1
The determinant of  this matrix:
   A = A-l - BC. = A + y.  C.
   y                   *
                       F '          -1
Only the first row of  [ y.  )..  .  , ]    is required to obtain  the

                Yi        1}  1}

first row of f ' (       ) .   This  row equals:
                  ~^1       1      1   Yi
                  ( - i.    ) .  ii  = r.i   _i i
                            - 42 -

-------
     is obtained from Ci by  interchanging  i  and  i+1  in the  y
and b values in C..
   BQ
       = G =
           30
                                  D
E

                     E+yjro

                     A+y.C.
                        •* i i
These expressions may be  used to calculate
^i = ^i  .dyi+l  +
dQ    ay, ,-, dQ
This is a recursive  relation which may generate a large round-
off error in addition  to  a  high truncation error which accumu-
lates as i decreases.

It may be both simpler and  more accurate  to use method j. (see
p. 40)  to update  the profile.   The same is true in the other
three cases examined,  as  will  be shown later.
                           - 43 -

-------
B.  The Longitudinal Profile with  Spot  Inflows  (Pipes):

Additional assumptions  for  this  case:

1.  There is a negligible region where  the  change  in heights
and velocities spreads  across  the  width of  the  river and then
the cross-sectional water profile  is  independent of the  lateral
position.

2.  The pipe is perpendicular  to the  direction  of  the main
flow.  This means that  the  incoming  flow does not  add momentum
in the direction of the main flow.

                               AX,^,
            AX .
                                     AX '
                             I  I
                     v
                          I
      v .
        i
                               i
                                                       Q+AO
                                                 Vi+l
                              Q > 0

                        FIGURE 6

        SCHEMATIC AND NOTATIONS FOR SPOT INFLOWS

The solving equations is obtained by postulating that  the  mo-
mentum is preserved in the direction of the main flow:

         Q'V .  , =  (Q + AQ) •v •
            1,1            1,2

y^ 2 may k"2 calculated recursively from the known situation in
cross-section i+1.
          Q + AQ
   v.
AQ is assumed to be fixed and positive.
                         Q
Substituting:  y.   = 	
                ljl   b. .v. ,
                       Q +  AQ
                              -  44 -

-------
gives:

   y.
              Q-Q

      b. , (Q+AQ)  (Q+AQ)
       •*-}•*-      K   \7
                         £

                 = CQ+AQ} " Yi,2
 i.l =
dQ
  !  2 2- AQ-VJ,2
I+^Q) [Q2(i+^Q)
   Q         Q

    r2-AQ    ^ _
                  Q(Q+AQ)
                                dy.
                              +    1  i =
                                 dQ  J
                                   dy.
                                     i.2.
                                dQ
Consider now the  case where  AQ changes.   The  vector of flow
changes will define  the  new  flow in every section.   The new
dpeths and their  first derivatives  may be estimated as follows:


5



A(
i
T V V
O 1 1 , 1
I
\
i \
\
\

— Q0+AQ
Yl,2 ^2 3

'3
dyo dYl dyia dy1<2 dy2 dy3
dQ dQ dQ dQ dQ dQ
Vo vl Vl,l Vl,2 V2 V3
                           FIGURE  7

   SCHEMATIC AND NOTATIONS FOR DERIVATIVES OF SPOT INFLOWS


The situation shown  above  describes  the  flow before the change.
Suppose a change of  A(A.Q)  occurs.   Then  for all the cross-sec-
tions downs tream,  the  new  flow is  Q0+AQ+A(AQ)  (assuming as be-
fore that A.Q+A ( AQ)>°) -   For all those downstream cross-sections,,

         ~ dy
The new dy/dQ may be estimated from the  basic recursive rela-
tion.  This process of  updating will proceed from the mouth of
the river up to  the cross  section next to the transition re-
gion  (yi2) -
                             - 45 -

-------
For the cross sections upstream  from  the change:
                        A(AQ)
                                    d(AQ)
Again, since the profile is obtained from a  recursive  relation
starting at the river-mouth, all the sections  downstream from
the change in flow do not distinguish between  the  possible
sources of the change in flow, so:
            cLQ
and:
     i,l
                         i,2
                                 i 9
 (Note that for computational procedure,  the expression  of the
 derivative dy^ i/d(AQ) is similar to  the expression  dy^ ^/dQ. )

The new value of dy^_ i/dQ may be calculated from  the  expression
for this derivative, ' substituting all  new values  which  are known,

The rest of the sections upstream may  now be updated  using y^_  ^,
dyj_ i/dQ and the basic recursive relation up to the place
where the next change occurs.

If there are more than one inflow  (or  outflow) into the same
section and they may not be lumped together, the  section may
be further broken into sub-sections.   A  convenient notation
procedure may be :  All odd second indices are upstream  and all
even ones are downstream from the location of discharge.
                            -  46 -

-------
C.  The Longitudinal Profile with Spot Outflows  (Pipes)
       v.
                    V
                            Axi+l
                                         i +19
                                                            Q-AQ
                            AQ > 0

                           FIGURE 8
           SCHEMATIC AND NOTATIONS FOR SPOT OUTFLOWS
In this case; the specific energy  (energy per unit of weight,,
as represented by the Bernoulli equation) is preserved provided
some condition is met, as shown later.
                                   t
Assuming that preservation of specific energy is possible:
v2
                         v2
           vi,l
       2gb

         1  ,
        2? <
                                       2
   2 ^s known from the recursive process
                           - 47 -

-------
It will be shown that there are  two  possible situations:

1   There are three real  solutions  to the cubic equation , one
    of them negative.  Since  the flow is assumed to be sub-cri-
    tical  the larger positive  solution is the appropriate one.
    The sign properties of  the  roots of the equation may be ob-
    tained directly from  the  form of the equation.  Yet, since
    it has to be solved anyhow,  the  nature of the solutions will
    be demonstrated by evaluating them.

2.  Critical conditions develop before the location of change.

To develop the two situations,  the starting point is the solu-
tion of Situation 1.

The Cardan method  [13]  for  solving the cubic equation is used.
Accordingly, the  following  substitution:


    Yi,l  = H  +i tyi,2 + 2?
                               —} "


 yields:

     3    1 . r    .,  ,fiiAfi ,2 2   +JM£)2_J_[     + i (glAfi ,2,^0
    u>  ~  T LY-  o  + -I  \-K,,   i  \    M>    9<-r *lv     ?7L-*i  9    9a   T-^
    r**    j.   1  ^    ——  jiD y .            ^ y  J^      /.. /   J-, ^i    £* ^3
                   2g     1} £
 or:




 To obtain the solution:
     a -        +
     A ' 21 ~ 2 +  \l(2   27
         c  _ D _
         27   2      2   27
 or:
       =^-^+^       1  -
         el _ D   ci     //n    27DX2
         27   2   27
                             - 48 _

-------
To get all three roots  to be  real,  it  is  required that A and
B be imaginary, or, the  term  inside  the square-root sign has
to be non-positive:
        2c
            2c
or:
    1   27D
    i.     r ;>  o
       2c
   Inspection  of  the  equation shows  that:

                                         D > 0

                                         c > 0

       27D
and so - -  > 0 is satisfied.
       2c

   2.
       2c


Substituting  the  original, values  for  D and c:
or,                                    2
Substituting,  the problem is  now transformed into:


   _4_  3    r -,  . Vi.2   , 3  /"   XX 2
   27 Yi,2  L-L  + -      I   > 1  U.;
                            - 49 _

-------
When this condition is met, Situation  1 describes  the  case.
Namely:  Ignoring losses in a  transition  area,  the  specific
energy downstream from the change may  also be maintained up-
stream.  The diagram shows the depth-specific energy relation
for different levels of specific flow.   (The specific  flow is
the flow per unit of width.)
                                                       Curve  for
                                                       specific
                                                       flow   up-
                                                       stream
                                                       from out-
                                                       flow
                                               specific energy
 Curve for specific
 flow  downstream from outflow
                        FIGURE  9

      CRITICAL,, SUPER-CRITICAL AND SUB-CRITICAL  FLOWS
    AS DEPICTED  BY THE DEPTH-SPECIFIC ENERGY  RELATION
If the change  in  flow  is  large enough and y^_  9  is  close  enough
to the critical value, Situation 2 will prevail.   In  this  case,
the specific energies  cannot be equal because the  specific flow
upstream from  the  location of the change is too high  to  be able
to ever be in  the  energy-level of the lower flow.   Since y^ -i
and y^ 2 are determined freely  (they are not  forced by any don-
trolling device) ,  y.j_ 2 will be determined by  the energy  required
to move the specified  flow and to overcome  losses  (namely, the
basic recursive relation) and y^ ^ will determine  the lowest
energy level possible  for the specific flow Q/b.   This lowest
energy is the  critical one.  So, in Situation 2, critical  con-
ditions will develop upstream from the change and  there  will
be energy dissipation  at  the transition area  so that  the condi-
tions defined  by  the basic recursive equation prevail after the
transition.  The  flow  above the critical cross-section will a-
gain be sub-critical.
                             - 50 -

-------
Solution  of  Situation 1:
       c3  ,,    27D,  .  c3   / n    27D.2
   A = -TJ  (1	3)  + yy   /  (1	3)   - 1

       27       2c     2/  V       2c
       c   n    27D>,     . c^   /~    M    27D.T
       — U       ;  -t- i      /  J.  -  u      3)

       Zl       2c       ll  4           2c
       c^  ,-.    27D,    . cj   /  ,    ,,    27D,2
   B = 27-  (1  -- 3) -  i 27   /  1  - (1  - - 3)

       2/       2cJ      /7  V           2c



       c3   rT!   27DN2    "    7!    27D,2    iff _
   A = -^=-   /  (1 -- o)  +  1 -  (1 -- o)   ' A

       27            J                 J
        27
Similarly:



    B  = — IicT
    B       SL
where:           /         27D
                 /I -  (1  -
                           2^L            ,      !
    a = arc tg [    ^  _  27D      ]  = arc tg  / —	J^2 ~ X

                        _3                  /  ' 1  ""    •*'

                        2c                  V      2c3




        i  ^ 27D
    for 1 ^ 	3

            2c
                                               cos -^
            c      a  ,
            —  cos -  +
           f cos f
           -r cos -  - v-r  c sin
       =  -  — c sin  (-| +  p)


                        c_


           B = arctan •  3    = arctan ^

                     "/Tc/3



                              - 51 -

-------
So:
         2     .   / a  ,  TT\
      = -- c sin  (-  +  ^)
Similarly:
        2
                     6'
To determine the behavior of  the  roots,,  a is explored:

                        L-27D\2
     = arctan
    27D
                         2c
                            3J
                   1 -
                       27D

                       2c3
For - r = 1^ the value of  CT  is  determined from the original
    2c
expression for A and B.  In  that  case:
   a
   A =
        . c3   c3
       X 27 = 27
   B = -i



So^ CT = 2


Since -1 <
           7   27
          for
,  3

27D

2c3
                           0
the above expression for a shows  that

   0 < a £ TT
Substituting gives:

   0 < f £  H-L < §  c
   0 > -  <
                 f >  0
There is also ordering  as  follows

   M-i ^ M-3 > M-2
                            - 52 -

-------
Substituting
             +f
yields :
Since the  flow is assumed to be subcritical ,  (y.   )   is always
                                                 i , ± 1

chosen.



   y±jl  =  f (2 cosf + 1) =
= I ^1,2
                     i,2
                             cos  [   arct*n
                                                            - 1 ]+l]
where:
   27D =  27.   _1_


   2c3 "  2    2^
                                 2g
 i.l = 1  rdc (2  cos £
dQ     3  ldQU  C°S 3
                                    _ ±     .   a

                                         3   dQJ
        i,2   2g
        dy.
         ^
   dQ "  dQ

   dc   dyi, 2 _  (Q-AQ)     r_-|  ,  Q-AQ   . dyi,2.
   -« .~ ^   -i xv              ^  I  ^~   - -        ^3 /^
         dQ
                             - 53 -

-------
   dQ   	i	 _1+1   "  dQ [/ n_27D,2
        (1 _  27D.2

             2c3
      /i   27D.2   1          '
   =  (i - —3)  •  j  • /	r
                            27D.2 -1

                            o  3'
                            2c
      /-,   27D. 2   1         J-         , 0 x      ,       / 27. d  ,D .
   =  (1 - 	o)   • ^- /	i	 ^~2^	±	   (-^~");T^ L3)
          2c3'    z  / 	^-^^T -1
                             _3 dD   -2 dS
                                   "
     27
                                     6
      2    /          2 7D  2            c

         J          2c3
     27          •»•           dD   3D   dC

     9^3* /In  27D.2   L  dQ ~ C    dQJ
     Z(-   /   -L    V-L-  0;
  - has been  found  before.
dO
^a c _

d°   gb2


Substituting  these expressions  gives dy^ j_/dQ.


As mentioned  in Case  B,  changes in AQ will be regarded as
changes in Q  in the sections  downstream from the location  of
change.  The  only variables affected by a change in AO are
Yi,l and dYi.,l/dQ •   yi,l maY be updated using a first order
approximation, for which dy-j_  j_/d(AQ)  is necessary:
                                        .in
                                                    dyi.2
               (Q-AQ)     rl  ,  Q-AQ . dy   -,
                       9  L -1- + „        J-^ J
               g(by   )^       Yi,2
                   1}z
                             - 54 -

-------
Substituting:
   d(AQ)     dQ


yields  for-r-    ar* expression very similar to that for
Appropriate  substitution yields:

    dc   = dc    2 (Q-AQ)

           dQ
                       1          3 dD   _    2  dc
    da   =  21 .                 f   d(AQ)   "^ d(AQ) ,
   d(AQ)    2   H    ,n 27D,2  [         _6          ]
                         2c3
         V         2c

   dc
  , .  -^ has been found before.
  d(AQ)
                           L   dD   -    _dC_
                            d(AQ)     C  d(AQ)
  d ( AQ)


So:
        = 0
dCT   _ 21          X   	   r-3D dc_
 / ,x->\ ~  Q   r        O Tn  T    L   /-• r!/AA\J
      2c   y  i -
                         2c3
                   dy.  7
Substituting into        gives a value for the first derivative
which may be used  to  update y^ 2 •

Situation 2 :

The critical depth is:

                             - 55  -

-------
   dY,
   d(AQ)      *

Transition between Situations  1  and  2:

If the expression:


(28)   (£[y1>2  + J, (g)2  ()]-<)J
                  zg        i, £

changes sign as  a result of  the  change  in  /\Q,  then a transition
between the two  situations occurs.

If expression  (28) changes from  negative to  positive,  then y^_ ^
and dy-j_ ±/dQ are  calculated  as in Situation  1  with the new tff~

If expression  (28) changes from  positive to  negative,,  then re-
gardless of yj_ 2> Yl l  is critical and  dy^_ ]_/d^O = 0.   If AQ
becomes negative_, the outflow  becomes an inflow.   Then the up-
dating has to be  done in two  steps:

         1.  To  the situation where  /\Q  = 0;

         2.  Using the  routine described in  Case  B.

The procedure is  similar in  the  case of changing an inflow
into an outflow.

D.  The Longitudinal Profile with Continuous Inflow (Tributaries)


The notational convention is:

   Capital letters - for main stream

   Lower case letters - for tributary

   Q^q = rates of flow

   V,v = velocities

   M = total momentum

   M  = component of total momentum  in  direction  of  main flow
    X.

   B,,b = widths of main flow and tributary

   ot = angle between main flow and tributary at  intersection.
                            - 56 -

-------
                         FIGURE 10

       SCHEMATIC AND NOTATIONS FOR CONTINUOUS INFLOW
Assumptions:

1.  Transition  regions  have  negligible  dimensions  and  effects.

2.  The depth of  the  tributary  equals  that  of  the  main flow
    and it  is approximated by y^  2•

3.  The component of  momentum at  the direction of  the  main
    flow is maintained.   (Actually  this may work only  for  low
    velocities.   In reality, preserving the momentum would have
    created a path of flow which  is forced  to  follow the path
    of the main-stream  river-bed.   In  that  case, there is  a
    bending phenomenon  which may  tend  to preserve  the  velocity.
    There is little  theoretical or  empirical knowledge describing
    the behavior  of the water in  the junction.  In the case of
    high velocities the depth of  the water  is  not  uniform
    across the  river.   In addition, there are  energy losses
    which cannot  be approximated.)

4.  Both flow velocities are sub-critical.
                            - 57 -

-------
    M  = QV + gv cos
        q X
    q = -r— sin
    v = Const.  = qO = qO
                 by   by.  _
                        i, z


    V = BY
          9    2
          2   q X
2QY
    dX
             §- Q2 S   ,  /V sin 2
    dM
      •*£ 	 71 / Q    Q \
    dx  ~ M o ~ bf;


whe re:
       A = BY = cross-sectional area of main  flow corresponding
                to point X

       S = river-bed slope

       S= frictional slope
           n
                   4/3     1.49-B-Y       BY
                  K
                   .  (JL+  1)4/3
                     V      ;
where n  is Manning's  resistance coefficient
    j  = - -T— sin
    dx      b
 Equating the two expressions  for   . x gives:
                                   dx

    B2y3So _  (_n2_)2 y  (I + 1)V3  =


                            -  58 -

-------
= - 2QY    Sin a - Q2
                                  2yi
                                     x j •
           2 v/l  . i}4/3    /_IOx2. B sin 2 »  . V2
               W   n'    "*"  \v«/->'    ••»..        * JL   —
   dX    V1.49'   iVY  '  B'     T vbQ'     2y7 .
                                         •i- • ^
         t, o      -,   2q   sin a
      -  
-------
           VI.  TRANSIENT EFFECTS IN AN ESTUARY


A.  Introduction

The transient behavior of the DO resulting  from a change  in  the
system is not considered by  the approach  to DO control which
utilizes the transfer coefficient matrix.  This study attempts
to justify the above approach by demonstrating that under  "real
world" conditions, the system representing  the Delaware Estuary
converges essentially monotonically without developing any
pathological behavior.  It is also demonstrated that the  length
of the transient period  (i.e.,, the time required to get within
+ 1% of the steady state solution and stay  there) in the  slower
sections is approximately one month.

B.  Representation of -the Transient

Thomann's model [10] consists essentially of  two systems  of
differential equations, one  describing the behavior of B.O.D.
and the other one describing behavior of  DO.  The second  sys-
tem couples with the first one through the B.O.D. levels.

Theoretically, an analytic solution to the  system of differen-
tial equations is available.  Its derivation  in this case, how-
ever, involves inversion of  a 30 x 30 matrix  (for the 30  sec-
tions of the Delaware) with  analytic expressions  (as opposed
to a matrix with numerical entries) .  This fact provides a con-
vincing argument for numerical integration.   In addition,  the
result of an analytic solution is again a 30  x 30 matrix  which
does not have any intuitive  meaning, necessitating numerical
evaluation of the transient  behavior under  assumed conditions
of loading.

Alternative approaches to the numerical integration are either
to investigate the stability of the system  (Pole-Zero analysis)
or to load it with various patterns of Dirac  Deltas  (spikes).
Both possibilities have been explored by  Thomann  ([10], p.27,
[12], p.354).  The main shortcoming  of these approaches  lies
in the fact that they do not provide a simple answer to the
question:  What would be the consequences of  a change in  the
system and would the transient create undesirable effects?

The approach adopted in this study tries  to provide the answer
by sampling numerically a few responses under conditions  which
are plausible in reality.  This approach  provides also an indi-
cation on the sensitivity of the system to  changes in certain
parameters.

C.  Assumptions

The transient reaction of a  linear system does not depend on
the forcing function  (the transient is the homogeneous solution)
                             - 61 -

-------
However, since a numerical time-pattern of the system's behavior
was sought, a simple forcing function which would give intuitive
meaning to the result, had to be assumed.

The most obvious one would have been to remove all sources and
sinks of oxygen and to evaluate, under proper initial conditions,
the mathematical homogeneous solution.  The disadvantage of  this
procedure lies in the fact that the physical phenomenon of sat-
uration with DO is ignored.  Since, as explained later  (Part D),
the system may be allowed to oscillate between the zero-level
of DO and the saturation level of DO, the following alternative
was preferred:

By removing all sources and sinks of oxygen except for the re-
aeration from the atmosphere and by assuming that the water  a-
bove the first section and below the last one is saturated with
oxygen,  (together with an additional simplifying assumption  ex-
plained in the sequel), the system is driven to a steady-state
solution at the DO-saturation level.

The additional assumption concerns the in-flowing tributaries:

The system of differential equations describes conservation  of
matter.  Accordingly, the amount of flow into a section has  to
equal the amount of water flowing out.  This is easily achieved
by equating these two quantities in the computer program, how-
ever, this implies that the conditions in the in-flowing tri-
butaries are the same as those prevailing in the next section
upstream.  This is a somewhat artificial assumption but it is
the simplest one that could be made if only the main stream  is
to be solved for.

It is easily seen that under these assumptions the DO satura-
tion level is the steady state solution for Thomann's equation.

D.  Computation Method

As stated before, no analytic solution was sought for the sys-
tem.  Instead, a numerical integration subroutine was used to
generate a time pattern of the behavior of DO.  The integrat-
ing subroutine utilized Hamming's modified predictor-corrector
method.  This subroutine checks the possible truncation error
at each step of the integration and automatically adjusts the
size of the interval of integration if the error exceeds a
specified value.

The integration period was deliberately chosen as 200 days,
which is much longer than the transient period  (as it turned
out to be), in order to discover beats  (dampening of the oscil-
lation and then a renewed increase in the amplitude) if any
would have occurred.

It should be stressed that Thomann's major assumption was the
                            - 62

-------
linearity of the main system.   This  assumption  may  not  be  valid
when the system approaches  its  boundary  conditions,  namely,  no
DO or saturation, the reason being  that  at  zero DO  a different
physical process takes over, while  the saturation phenomenon is
generated by a stochastic process of exchange of oxygen between
the atmosphere and  the DO.  These conditions are not explicit
in the system of equations, so  there is  nothing inherent  in  the
system of equations  to prevent  the  solution from oscillating
above the saturation value  or below zero if the numerical  struc-
ture of the equations is appropriate.  As a result,  one of the
basic problems was  the adequacy of  Thomann's model  for  describ-
ing the transient.   This is the reason why  it was not clear
apriori if a beat would develop, or  if the  system would try  to
oscillate beyond its boundary values.

The solution to the  second  problem  was to bound the  DO  values
and their derivatives in case of a  bad oscillation.   As shown
in Part E, the actual runs  were designed to investigate the
necessity of constraining the DO levels, with the result used
to indicate the adequacy of the model in estimating  the transi-
ent behavior.

E.  Computation and Results

    1.  The conditions assumed  correspond to those prevailing
        in the summer since this is  the  critical period.

    2.  The starting conditions  for  the integration  of the  DO
        levels were  the average values in the summer of 1964
         (Table 3, Schaumburg  [9]).   The  numbers are  given  in
        Table 7.

    3.  The saturation value was calculated as:

             c =14.652-0.410222T+0.0079910T2+0.00007777T3
              s
         (Source —  Pence, Jeglic, Thomann  [8]).

        Substituting T=25°C=?77°F

        yields

             c =10.605976 mg/A
              S
        So the value used for the integration was 10,6.

    4.  Since the numerical integration  is  a sequential process
        and at every step values are computed sequentially (by
        a DO-loop)  for all  sections of the  river, starting with
        1 and ending at 30, the errors  (truncation  and  round-
        off) tend to accumulate towards  the later points  in
        time and towards the  sections indexed by higher numbers.
                             _ 63 _

-------
    The integrating subroutine estimates the error for each
    component of the solution vector generated at each step
    of the integration.  Then the sum of those error  terms
    is compared with a predetermined value.  If the error-
    sum turns out to be too high, the integration interval
    is halved and the step is repeated.

    An option is offered by the subroutine to assign  weights
    to the individual error terms, thus stressing the accur-
    acy of certain components of the solution vector.

    Utilizing this option, an attempt was made to improve
    the results by re-weighting the error terms of the
    individual sections.  The purpose was to increase the
    sensitivity of the error-sum to those components  cor-
    responding to sections of the estuary which were indexed
    by higher numbers.

    The re-weighting was done proportionally to r? and to
    x/rf ,  where n is the section number.  In the case  of /rf ,
    the error term corresponding to section n was multiplied
    by /n/ 30 -/ri .  In the case of n^  the weight was  n^/ 30 n2
          2j _                                             2-i .,
          n=l                                            n=l

5.  Since the system turned out to be insensitive to  the
    variations tried, they will be listed and the results
    given for one of them.

    Sensitivities checked:

       i)  increasing the saturation value to 12.6 mg/ 1. and
           decreasing it to 8.2
      ii)  forcing as opposed to not forcing the DO  levels
           to be between 0 and the saturation value .

6.  As stated before, the results were insensitive  to  the
    above variations.  Re-weighting as n^ improved  a little
    the numerical behavior.   (For explanation, see  comments
    following Table 7.)  Forcing the DO to be at most equal
    to the saturation value and the derivative to be 0 or
    lower in case the section is saturated, introduced addi-
    tional error and the system converged to a somewhat
    lower value than the saturation level.  The runs made
    without these constraints demonstrated almost no oscilla
    tion above the saturation value, indicating that the
    linearity assumption is adequate in this sense.

7.  Lengths of time for the system to converge to + 1% of
    the saturation value are listed in Table 7.  In  this
    run the saturation value was 10.6 mg/£, the system was
    not forced to remain below this value and the error
    terms were re-weighted as n2 .
                        - 64 _

-------
Section
Number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
-Transient
Period
in Days
4
4
6
10
14
17
21
24
25
27
29
30
31
32
33
33
33
33
33
32
31
30
31
29
28
25
22
23
21
16
Starting Error
Conditions Range
mg/ & mg/ &
7.5 nil
7.0
6.2
5.1
5.1
5.6
5.4
5.3
3.0
2.2
2.3
1.5
1
1
1
1
1
1
1
2
2
3
4
4
5
5
6
7
7
8
.2
.1
.1
.0
.0
.2
.4
.2
.9
.5
.1
.8
.5
.8
.2
.2
.6
.0







11
197
759
2386
7201
13230
15816
13395
6799
743
56
11
n
"
"
11
11
11
x
x
X
X
X
X
X
X
X
X
X







10
10~q
10 ^
10 ^
10 b
10
10~n
10~^
"Is
10 R
10 b
                           TABLE 7

        TRANSIENT PERIOD OF SECTION i, CORRESPONDING
        STARTING CONDITIONS AND COMPUTATIONAL ERROR
The error-range  (last  column in Table 7)  is the maximal absolute
value of the difference  between the following two numbers:

        i.  The  saturation value (10.6 mg/jO .

       ii.  The  numbers  generated by the  computer for the per-
            iod  between  the end of the transient phenomenon
            (column 2  of Table  7)  and the end of the integra-
            tion process (200 days) .
                              -  65  -

-------
These deviations represent the effect of the cummulative error
generated by the numerical integration process.  Decreasing
these deviations was considered as "improving the numerical
behavior"  (see E.6) and this was the reason for re-weighting
the error terms as explained in E.5.

The maximal deviations occur in sections 25, 26 and 27 and they
constitute  an error of about 1.5%.  The deviations in the rest
of the sections are well below 1%.

These observations may serve as indications on the accuracy of
the numbers generated by the integrating subroutine and on the
validity of the results of this section.
                             - 66 -

-------
     VII.  POLLUTION TAXES AND CHARGES IN PERSPECTIVE


Previous work [4], presented a model to find a minimal treat-
ment cost solution to a water quality problem.  A solution to
such a problem gives an efficient solution in terms of cost
minimization.  However, such a solution may be costly to imple-
ment and enforce unless each polluter accepts the regional solu-
tion and his share of its costs.  We now wish to discuss whether
there are decentralized mechanisms so that decisions by indivi-
dual polluters correspond to minimization of regional costs.
The possibility of the decentralization of decision making is
important in economics since enforcement costs are minimized
when the social optimum can be achieved through individual de-
cision making.  In economics, taxes or charges of various types
are often suggested as a decentralized means of making indivi-
dual behavior correspond to a social optimum.  Several taxa-
tion and charge schemes will be discussed below.

Examination of the literature concerning the economics of
water quality reveals the assumption that society has a right
to a certain level of water quality, or there is a socially op-
timum level of water quality-  By implication, the polluters
should either

    1.  be forced to treat waste material at a high enough
        level to maintain the quality required by society

or  2.  pay the costs society incurs to maintain the required
        water quality.

In this view, two types of tax schemes have been discussed,
those relating to enforcing a quality standard and those re-
lating to raising revenue.  The first type can be called an
enforcement tax.  One example of such a tax would be to charge
a polluter an extremely high tax  (or fine) on wastes dumped
above a certain level.  Thus, above this level, the polluter
would rather cease polluting than pay the fine.   (To do this,
he could either curtail production or build a treatment faci-
lity.)  Below this level, the polluter is in effect given free
rights to pollute.  Such a tax would not be intended to raise
any revenue.

Another enforcement-type tax would be to have a uniform tax
per unit of waste and the tax would be computed so that a pollu-
ter would dump no more than a certain amount.  This is illus-
trated below
                               T
  price percent
  effluent
D
                            Eo
                             - 67 -

-------
DD represents a polluter's demand to dump effluent as a func-
tion of the tax.  If  E  is the amount of effluent to be allow-
ed under a quality standard, then the per unit tax which would
result in Eo is given by the slope of T.  If the polluter would
want to produce at a higher level than that corresponding to Eo,
then again he would have to build waste treatment facilities.
A related concept would be to view the service of waste disposal
via a river as a public resouce to be sold at a price.  The
price and amount sold is set to ration use of the service cor-
responding to the quality goals.  In both these cases, the tax
revenue collected is "gravy" since the tax collecting agency
does not construct treatment facilities; the taxes serve only
to force achievement of the quality standards.  Any treatment
undertaken under these enforcement tax systems would be at an
individual level and would be individually paid for.  Because
of the individual emphasis, such tax systems would not result
in minimizing regional treatment costs.

The second type of tax scheme would be for the purpose of
raising revenue for treatment purposes.  For instance, a pollu-
ter could be charged for waste disposal according to the social
cost of his pollution and the revenues collected could be used
either to alleviate the damages or pay for treatment.  Here
again, this may not lead to the least cost regional solution
since it is always easier and cheaper to treat wastes before
they enter a stream than to have to remove them from the stream
once they are there.  Alternatively, revenue could be raised
through charging the polluter for his use of a public or "joint"
waste treatment facility.  It would not be necessary for the
charge to be collected by a governmental agency in this case.
A contractual agreement could be formed among several polluters
to operate a joint facility such as a regional treatment plant
or reservoir.  The charge for use of the facility would be
based on cost.  The advantage of this system is that the least
cost regional treatment system could be obtained on a voluntary
basis.  However, in this case, there would first have to be some
sort of enforcement mechansim, either legal or through enforce-
ment taxes, to give incentives to polluters to have wastes
treated in the joint facility rather than dumping them.  If en-
forcement taxes were used as incentives, then the treatment
charges would have to be lower than the tax; otherwise polluters
would pay the tax rather than treat their wastes in the regional
facility.  The charges would also have to be low enough for the
polluters to choose the joint treatment alternative rather than
treating their own wastes individually.  The problem is compli-
cated by the need to cover full costs of building and operating
the joint facility from the charges.

It is possible that a proposed tax could fall into both of the
above two categories as it affects polluters.  For example,
suppose a uniform effluent charge of $.05 per mg/d is to be
levied on effluent per unit B.O.D. concentration to pay for
treatment costs.  Because of differing manufacturing processes,
                             -  68 -

-------
locations and other parameters the cost of removing waste mater-
ial may vary for polluters.  For this reason one polluter may
have an on-site B.O.D. removal cost of $.04 per mg/d and another
on-site B.O.D. removal cost of $.06 mg/d.  In this case, the
second polluter would pay his effluent charges since his margin-
al removal cost is more than the cost of dumping it in the river
and the first polluter would treat his own waste.  This is not
inconsistent with the least cost regional solution since the
first polluter can remove his wastes himself more cheaply than
in a regional plant; such a situation can arise when a polluter
is so far from a regional plant that piping costs overcome possi-
ble economies of scale.  Here we see that the tax has served
two purposes.  In the case of the first polluter he has been
forced to treat his own wastes.  In the case of the second pollu-
ter revenue has been raised to defer the cost of treatment under-
taken by the "government" or tax collection unit.  However, it
is clear that unless treatment plant sizes and charges are based
on the least cost regional solution, in general the result of
such charges may not cover costs of the facilities.  That is,
in the above example, if a treatment plant were designed to
treat wastes of both the first and second polluter and the charge
set accordingly at $.05 per mg/d, then the plant would in fact
be too large since the first polluter would do his own treatment
rather than use the plant at that charge.

The determination of charges for the purpose of quality enforce-
ment is generally based on two criteria:

    1.  Equity between polluters;

    2.  quality resulting from implementation of charges meets
        quality goals.

Johnson in [7] determines several types of effluent charges
which meet the second criterion above.  The methods used were
called least-cost, single effluent charge, and zoned effluent
charge.

The least-cost method proposed by Johnson seeks to find the set
of effluent charges which will result in the minimization of the
sum of the taxation costs incurred by the individual polluters
subject to the constraint that the quality goals are met.  This
approach could possibly result in a different effluent charge
for each polluter.

The second approach seeks to impose a single effluent charge in
terms of $/lbs. to all polluters, again subject to the constraint
that the quality goals are met.  Another approach, the zoned
method divides the polluters between regions and assigns a single
effluent charge to each region to minimize the total taxation
cost over the region.  Actuarlly, all of these methods are vari-
ants of the zoned method since the single charge method treats
the entire basin as one zone and the least-cost method in effect
                             - 69 -

-------
assigns each polluter to a zone.  In all cases the least-cost
solution is strived for given the constraints on the charges
and the quality requirements.  Johnson's charges are all exam-
ples of enforcement taxes.  His cost minimization criterion
is a way of determining how polluting rights are to be distri-
buted and is thus an equity criterion.  It would also be possi-
ble to devise other methods of assigning effluent charges.  For
example, instead of zones the polluters could be divided with
respect to type, each type being assigned an effluent charge.
Assigning the same charge for each type gives another equity
criterion.

As already mentioned, there is no reason to assume that the
revenue obtained from a uniform effluent charge would exactly
meet the cost of treatment.  Also, an enforcement type tax
would not in general result in minimizing basic treatment cost.
For these reasons, recently the revenue collection type of tax
has become more interesting since the thrust has been to reduce
total basin costs significantly through the use of treatment
alternatives such as by-pass piping, regional plants and flow
augmentation.  Each of these alternatives creates problems of
joint funding and operation of facilities.  In such a case the
revenue necessary for the operations of these joint facilities
must be obtained from the polluters concerned.  Charges of the
revenue collection variety are usually based on three criteria:

    1.  Equity between polluters;

    2.  quality goals maintained;

    3.  revenue sufficient to cover costs incurred by tax
        collecting unit.

Graves, Hatfield, and Whinston [5] construct a pricing scheme
in their paper on by-pass piping which is of this type.  Their
scheme is based on an equity criterion of charging polluters
according to damages.  The scheme is not self-enforcing and so
would have to be preceded by some enforcement mechanism to give
incentives to meet the quality goal through by-pass piping.
The minimum cost piping arrangements is first solved for sub-
ject to meeting the quality goals and conservation of flow con-
straints.  The dual variables associated with the flow conser-
vation constraints measure the marginal cost of extra effluent
flow from an individual polluter to the least cost basin solu-
tion.  Suppose each polluter were charged according to the dual
variable per unit effluent.  This would lead to a charge of
                             -  70 -

-------
   denotes the optimal effluent  flow  from polluter k  from  solu-
tion of the least cost programming problem,  TC  is the least cost
solution, and the partial derivative  is evaluated for the  least
cost solution.  Total revenue  from this charge  would  be


                 y  2£!C  £E  =   £
which in general would not equal TC.   The expression  in  (29)  i
only guaranteed equal to TC  for  all values of F|, k=l,n   if  the
least total cost function is  linear homogeneous  in the effluent
flows .

In the small 5 polluter example  in the Graves, Hatfield,  Whin-
s^ton paper [5], the least-cost solution is $55,577.   The  total
revenue collected using an effluent charge e.qual  to  the value
of the dual variables associated with the flow consistancy con-
straints is $179,499.  In terms  of the notation  above

               TC =  $55,577


               5
               Z TP,  = $179,499
               k=l   k

or
                    5
               TC <  Z TP,  .
                    k=l  K


In the Graves, Hatfield and Whins ton  paper the revenue is forced
to equal cost  by redefining  the  total payment of  polluter j.

                          n
               Tp* - (TP ./ 2  TP, ) TC
                3      : k=l  *

This automatically  gives
              Z  TP*  = TC  .
             k=l k

This cost allocation  is based  on  achieving  the quality  goals
at least cost, costs  are  covered  from  the revenues,  and there
is an equity criterion.   The equity  criterion  is based  on  pollu
ters causing the same marginal damage  paying the same per  unit
                             - 71 -

-------
effluent charge, where marginal damage is defined to be the
increase in the basin solution cost due to an extra unit of
effluent.

Although the above is one type of equity criterion, it is clear
that there are other types,, such as the same type of pollutant,
or the same charge for polluters in the same zone.  The notion
of what is meant by equity needs to be further examined.  Equity
is included as a criterion for both revenue and enforcement
taxes.  Its importance is due both to our notions of fairness
and to the difficulty of enforcing and collecting the tax or
charge when it is not felt to be equitable by those who will
be charged.  In general, the notion of equity is based on two
principles:

    1.  Polluters who cause more damage should pay more.

    2.  Polluters who cause equal damage should pay the same.

The problem which is evident from these requirements is that
of interpreting "damage".  The proposal to charge polluters
of the same type the same amount per unit weight makes the
implicit assumption that polluters who dump equal weights of
effluent of the same composition are causing equal damage.
However, this assumption may be false, for instance when the
flows  (hence concentrations) of their effluents differ.  The
location of a polluter also determines how much damage he will
cause.  A polluter located in an area of high water quality
(and/or high flow)  will do considerably less damage than the
same polluter located where the water quality is poor  (and
stream flow is low).  Damage also depends on what is located
downstream of a polluter.  Damage will be much less if a pulp
mill is located downstream than if a municipal water plant is
located downstream from a polluter.  Under closer examination
it is not at all clear that any measure of effluent composition
alone reflects pollution damage.  In an economic sens.e, the
only measure of damage which is really relevant is the amount
of additional cost necessary to maintain desired quality goals
due to the existence of polluter k.  If polluter k dumps a
certain effluent but causes no additional treatment cost to
the system  (i.e., quality constraints are not violated), then
his existence is irrelevant as far as any treatment cost scheme
is concerned.

The Graves, Hatfield, and Whinston [5] paper recognized cost
in this way as a reflection of damage through their use of the
dual variables of the programming model to generate equitable
shares of the total basin costs.  The problem associated with
this measurement of damage, as with any sort of marginal con-
cept,  is that this information is only "local" in nature.  Since
the relationship between the optimal basin solution, and the
flows from the individual polluters is non-linear, the marginal
cost at the optimum solution of an additional unit of flow from
                             - 72 -

-------
polluter i is no reflection of the added cost  to meet  the stan-
dards if effluent from a polluter were  to be increased by more
than one unit.  Thus, the equity criterion used in the Graves,
Hatfield, Whinston paper is not completely satisfactory.

Incremental Cost Allocation Scheme

In the previous sections of this paper we have indicated several
faults associated with the pricing and cost allocation systems
presented in the literature.  These  faults can be summarized as
follows:

    1.  Failure of revenue collected to meet total costs of
        treatment.

    2.  Failure of some pricing systems to be based on the most
        efficient  (least-cost) combination of  treatment methods.

    3.  Failure of pricing schemes to satisfy  acceptable equity
        criteria.

In this section we will present a cost allocation system which
we feel adequately deals with these problems.   This approach takes
as given the present composition of  polluters  in the river basin
and the given quality goals and attempts to allocate the costs
of a basin wide system.  Both the costs of individual  polluters
and the joint costs of using regional treatment plants are con-
sidered in the allocation scheme.

After an optimal solution to a basin polluter problem, such as
that described earlier, has been obtained, it  is not at all
clear what part of the cost of that  system of  dams, treatment
plants, and pipelines can be attributed to the fact that a par-
ticular polluter is in the basin.  That is, when several pollu-
ters are present, it is not clear who has the  responsibility
for violating the quality constraints and hence the exact damage
costs attributable to any one polluter  is not well defined.  If
we could measure the cost that an individual polluter  causes to
the system, then we would have some  basis to decide what his
contribution to the total cost of the abatement system should
be.  We will now describe a way of allocating  the total system
costs to polluters based on the idea of each polluter  paying
the incremental cost he causes the system.

Consider a basin with only two polluters designated 1  and 2.
Denote the least-cost solution to the basin problem  (29) by
C(l,2).  The incremental cost caused by polluter 1 being in
the basin, given that polluter 2 is  in  the basin, is defined
to be

              IC(1) = 0(1,2) - C(2)

where C(2) is the cost of polluter 2  treating  his own  waste.
                             - 73 -

-------
Both the values of C(lj2) and C(2) can be  computed using the
programming model in the previous  section.

A logical tax based on each polluter paying  his  social cost
would then appear to be to charge  polluter 1 1C (1)  and pollu-
ter 2 1C (2) , where

              1C (2) = C(2) .

We note that

              1C (1) + 1C (2) = C (1,2)

which means that the total cost of the basin project has been
covered by this tax scheme.  However, the  following question
now arises:  What if the polluters are considered  in reverse
order?  Consider the following definitions of 1C* (1)  and 1C* (2 )
              1C* (2) = C(l,2) - C(l)

If polluter 1 is charged IC*(1) and polluter 2  is  charged
1C* (2) the total cost of the basin system  is covered  as  before,
but the charge to each polluter differs  in  these  two  schemes.
As an example, consider a case where  if  either  of  the polluters
are alone in the basin the quality goals are met without treat-
ment and if they are both in the basin the  goals are  violated.
Assume that the optimal solution when both  are  present is to
build a reservoir for flow augmentation  at  a cost  of  C(l,2)  =
$100.  The values of C(l) and C(2) are both zero since the stan-
dards are met if either polluter is absent.  if the first cal-
culation of incremental cost is used  the charges for  polluters
1 and 2 are :

              1C (1) = $100

              1C (2) = $0

If the second calculation is used the charges are:

              IC*(1) = $0

              1C* (2) = $100

Polluter 1 clearly would prefer the second  allocation over the
first cost allocation and by the same token polluter  2 would
prefer the first allocation.  Since either  ordering of the
polluters is equally plausible, we may average all  the possible
incremental costs which each polluter could cause  to  obtain a
cost allocation of

              p(l) = & IC(1) + & IC*(1)


                   = 4 (100) + £  (0) = 50
                             - 74 -

-------
and

              p(2)  =  | 1C (2)  + | 1C* (2)


                   =  £ (0)  + i (100)  = 50


The sum of the payments from the  polluters  covers  the total
basin cost of 0(1,2)  = $100.   However,  the  equity  problem has
been taken care of by considering all possible  orders of the
polluters, in this case  (1 2)  and (2  1),  and  averaging over
possible definitions  of incremental  cost for  each  polluter.

This approach is easily extended  to  a case  of three  polluters
where there are six possible  orders  for calculating  incremental
costs for each polluter.   For example,  considering polluter 1,
the following orders  are possible.

     Order 1:   (123)]
                       V polluter 1  before  polluters 2 and 3
     Order 2:   (132) )

     Order 3:   (213)} polluter 1  after polluter 2

     Order 4:   (312) } polluter 1  after polluter 3

     Order 5:   (231)]
                        V polluter 1  after polluters  2 and 3
     Order 6:   (321) I

Corresponding to each of these orders is  a  way  of  defining in-
cremental cost for polluter 1 (the subscript  denote s the order
under consideration) :
                = C(l,2)  -  0(2)

                = C(l,3)  -  C(3)

                = C(l,2,3)  -  C(2,3)

                = C(l,2,3)  -  C(2,3)

For polluter 2, the  following  incremental  cost definitions are
possible :

        IC1(2)  = C(l,2)  -  0(1)

        IC2(2)  = 0(1,2,3)  -  0(1,3)
                             - 75 -

-------
         IC3(2)  = C(2)

         IC4(2)  = C(l,2,3)  - C(

         IC5(2)  = C(2)

         IC6(2)  = C(2,3)  -  C(3)

And for  polluter 3,

                = C(l,2,3)  - C(
        IC2(3)  = C(l,3)  -

        IC3(3)  = C(l,2,3)  -  C(l,2)

        IC4(3)  = C(3)

        IC5 (3)  = C(2,3)  - C(2)

        IC6(3)  = C(3)

Giving each of  the  possible  orders equal weight, the average
incremental cost for  each polluter is obtained:
         p(i) = j-  2  1C. (i)   ,   1=1,2,3
               6  j=l  ^

Note that

         3
         Z  p(i) = C(l,2,3)
         1=1

So that  total costs  would be covered by charging each polluter
the amount  p (i) .

Generalizing to the  case  of n polluters where there would  be
n! orderings, the cost allocation formula is given by


(30)    p(i)  =i  2*IC. (i)    ,    1=1, ...,n
              n-  j=i   3

where ICj (i) is the  incremental cost due to i in ordering  j.
Also,

          n
          2  p(i)  = C(l,2,  ...,n)
         1=1

This cost allocation scheme is discussed in [14].  Also  in [15]
                             - 76 -

-------
there is a theoretical development of  the cost  allocation me-
thod described.  It is demonstrated  that formula  (30)  is a
unique formula which satisfies  the following  axioms:

    1.  Charges for the use of  the basin treatment  system must
        cover costs.

    2.  A user's charges are based only on  the  incremental
        cost caused by the user.

    3.  The user's charge is independent of the ordering or
        labeling of the users.

    4.  If a user's incremental costs  increase  by \,  then his
        charges will increase by  \.

The equity properties of this scheme are evident from the above
axioms 2-4.

The following example illustrates cost allocation for the five-
polluter case.  In the report [5] the  solution  of a small scale
pure by-pass piping problem is  given on page  90.  The dual solu-
tion for this problem is

        XL =  $3285/MGD            xg = -$1, 929, 780 .MG/L

        X2 =  $11970/MGD           X10= °

        x3 =  $6445/MGD            x^ 0

        X4 =  $2865/MGD

        X  =  $2381/MGD


The extremal value can be calculated as follows

        Q± =  19.7      Q1x1 = 64,714    K-jXg = -147,107
Q2 =
Q3 =
Q4 =
Q5 =
7.0
3.0
6.1
4.0
Q2x2 = 83,790
Q3x3 = 19,335
Q4x4 = 17,476
Q5X5 = 9,524
                                                  194,839

                                                 -147,107

        Kx =    .07623         194,839             $  47,732


A pricing scheme suggested  in [5] was  to use  the values of  the
xiQi to give the proportion of  the  optimal  solution  to be paid
by each polluter.  This gives
                            - 77 -

-------
        Polluter        % Paid          Payment/Year

           1             33.2             $15,643

           2             43.0              20,261

           3              9.9               4,664

           4              8.9               4,146

           5              4.8               2,261

It would be interesting to  compare the costs from the global
scheme to those given  above.   This is economically feasible
since we only have  five polluters.  The coalition costs were
obtained from the optimal solutions of 31 L.P. problems and
are given as

        C(l,2,3,4,5) = $47,119    €(3,4,5)  =  0

        C(l,2,3,5)   = 30,763    C(l,2,5)  =$25,238

        C(l,2,4,5)   = 29,365    C(2,3,4)  =  16,690

        C(l,3,4,5)   = 0         C(l,3,4)  =  0

        C(2,3,4,5)   = 18,939    C(l,2,4)  =  27,115

        0(1,2,3,4)   = 37,903    C(l,2,3)  =  27,844

        C(l,3,5)     = 0         C(l,2)    =  22,988

        C(l,4,5)     = 0         C(l,3)    =  0

        C(2,3,5)     = 14,812    C(2,3)    =  12,563

        C(2,4,5)     = 14,084    C(l,4)    =  0

        C(2,4)        = 11,902    C(5)      =  0

        C(3,4)        = 0         C(4)      =  0

        C(l,5)        = 0         C(3)      =  0

        C(2,5)        = $11,210    C(2)      = $10,381

        C(3,5)        = 0         C(l)      =  0

        C(4,5)        = 0

The polluters not in the coalition were allowed to return their
wastes to the estuary  at no cost which in effect made  them non-
existent.  The  global  scheme gave the following cost allocation:


                             - 78 -

-------
        Polluter         Local Method     Global Method

           1               $15,643          $10,372

           2                20,261           24,474

           3                 4,664            5,159

           4                 4,146            4,633

           5                 2,261            2,478

Application of Pricing Scheme

In the last subdivision we presented an alternate cost alloca-
tion mechanism to be used for regional treatment systems.  In
this subdivision we apply the scheme to be the Delaware Estuary
Mode1.

We now wish to divide the total cost of the regional system,
in this case the cost of the plants and pipes, between the
polluters using the cost allocation method explained in the
last subdivision.  This necessitates the calculation of C(G),
G <^ N where G ranges over possible subsets of polluters.  The
interpretation of C(G) in terms of river basin model is not
straightforward.  For example, how do we compute C(j) and what
do we assume about polluters I^j.  Are they non-existent, are
they treating at some uniform level?  Are they dumping their
effluent with no treatment?  An infinite number of assumptions
could be made concerning the behavior of the polluters not in
the coalition to be examined.

The only restraint on the behavior of polluters outside of the
group being considered is that they are treating in a manner
which will allow the polluters to gain feasibility of the basin
problem  (i.e., the quality constraints can be satisfied for
some level of treatment of wastes).  For example, if the pollu-
ters outside of the coalition are assumed to be dumping all
wastes directly into the stream it is very likely that several
sections of the river will violate the quality constraints no
matter what the polluters do.  One procedure which can guaran-
tee the feasibility of the problem is to find the least cost
solution to the basin problem, C(N), and assume that all pollu-
ters outside of G are required to treat their effluent to the
same level as they did in the solution for C(N).

The term C(N)  can now be interpreted as the least-cost the
polluters in G must pay to maintain the quality goals given
that the other polluters in the basin are treating on site at
the same level required for the least cost base solution.  The
polluters in G can make use of regional plants and piping.   In
the special case where G contains only one polluter piping is
                             - 79 -

-------
the only alternative method available.  In some problems no
savings over on-site treatment may result for some coalitions
of polluters.  For example, if polluters in the coalition are
so far apart that piping costs overwhelm any economies of scale
gained from a regional plant,  the best possible alternative for
these polluters is to treat on site.  In such a case, the cost
to the coalition would just be the sum of the individual on-
site costs, or

              C(G) =  Z  C(i)
                     ieG
                            -  80 -

-------
                      VIII.  REFERENCES

 [1]  Coddington, E. A.,  "An  Introduction  to Ordinary Differential
     Equations/"  Prentice-Hall,  1961.

 [2]  Geoffrion, A. M.,  "Generalized Benders Decomposition," Publica-
     tion  forthcoming in Journal  of Optimization Theory and Applica-
     tion,  1972.

 [3]  	, Graves, G.  W., "Multicommodity Distribution System
     Design by Benders  Decomposition,"  Forthcoming in Working Paper
     No. 181, Western Management  Science  Institute, October 1971.

 [4]  Graves, G. W., Hatfield,  G.  B. and Whinston, A.,  "Mathematical
     Programming for  Regional  Water Quality Management," Water Pollu-
     tion  Control Research Series, U.S. Department of  the Interior,
     Federal Water Quality Administration, August 1970.

 [5]  	, Hatfield, G.  B., Whinston, A.,  "Water Pollution Con-
     rol Using By-Pass  Piping," Water Resources Research, Vol. 5, No.
     1, February 1969.

 [6]  Henderson, F. M.,  "Open,Channel Flow," Macmillan  Series in Ci-
     vil Engineering, 1966.

 [7]  Johnson, E. L.,  "A Study  in  the Economics of Water Quality Man-
     agement," Water  Resources Research,  Vol.  3, No. 2, 1967.

 [8]  Pence, G. D., Jeglic, J.  M., Thomann, R. V., "Time-Varying
     Dissolved-Oxygen Model, "  ^Journal of  the  Sanitary  Engineering
     Division Proceedings of the  American Society of Civil  Engineers.
     Vol.  94, No. 5915, April  1968.

 [9]  Schaumburg, G. W.,  "Water Pollution  Control in the Delaware
     Estuary,"  Submitted to the  Department of Applied Mathematics
     in partial fulfillment  of the requirements  for the degree of
     B.A.  with honors,  Harvard College, May 1967.

[10]  Thomann, R. V.,  "The Use  of  Systems  Analysis to Describe  the
     Time  Variation of  Dissolved  Oxygen in a  Tidal Stream,"  A disser-
     tation in the Department  of  Meteorology  and Oceanography  sub-
     mitted in partial  fulfillment of the requirements for  the degree
     of Ph.D. at New  York University, February 1963.

[11]  	, "Mathematical Model  for  Dissolved Oxygen,"  Journal
     of the Sanitary  Engineering  Division Proceedings  of the American
     Society of Civil Engineers.  No.  3680,  October  1963.

[12]  	, "Recent Results  from a Mathematical Model  of  Water
     Pollution Control  in the  Delaware Estuary," Water Resources
     Research, Vol.  3,  pp.  349-359,  1965.
                                - 81 -

-------
[13]   Uspensky,,  J.  V.,   Theory of Equations. Me Graw-Hill, 1948.

[14]   Loehmarij  E.,  Whinston.,  A.,,  "A New Theory of Pricing and
      Decision-Making  in Public Investment," Forthcoming in Bell
      Telephone Journal of Economics and Management, Fall 1971.

[15]   	,	,  "An Axiomatic Approach to Cost Alloca-
      tion for  Public  Investment^" Forthcoming.
                               _ 82 _

-------
1
Accession Number
w
5
2

Subject Field & tjroup
06 A
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Organization
          California University, Los Angeles
     Title
          Extensions of Mathematical Programming  for Regional Water
          Quality Management
10

Authors)
Glenn Graves
16

21
Project Designation
16110 EGQ
Note
 22
     Citation
          Water Pollution Control Research Series, EPA,  16110 04/72 EGQ
 23
Descriptors (Starred First)

     *Regions, *Computer  Models,  *Mathematical Models, Networks,
     Water Quality  Control,  Water Pollution Treatment, Optimum
     Development  Plans
 25
Identifiers (Starred First)
     *0ptimal Configuration  of  Regional Treatment Plans,
     Non-linear Programming
 27
    Abstract
          This report extends the work contained  in the earlier work,  "Mathematical
          Programming for Regional Water Quality  Management."  A mixed integer,  con-
          tinuous variable non-linear programming model is  developed which promises
          to be much more realistic and effective in  selecting an optimal configura-
          tion of regional treatment plants and pipe  juncture nodes.   A new  Bender's
          type decomposition for the non-linear model  is  also presented as an  appro-
          priate and feasible solution technique.  In  the area of the  physical response
          of the environment, the earlier work linked  the dissolved oxygen profile and
          the control options through a constant  transfer matrix obtained as a steady
          state solution to a system of differential  equations.  This  work shows how
          the transfer coefficients can be recalculated periodically in the  course of
          the system optimization.  Full scale computations are carried out  demonstrating
          this capability.  The theoretical framework  for including two additional treat-
          ment methods,  mechanical re-aeration and flow augmentation via reservoirs is
          also developed.
Abstract
                              Institution
 WR:102  (REV. JULY 1969)
 WRSIC
                        SEND WITH COPY OF DOCUMENT. TO: WATER RESOURCES SCIENTIFIC I N F ORM A T I ON G E NT ER
                                                 U.S. DEPARTMENT OF THE INTERIOR
                                                 WASHINGTON, D. C. 20240
                                                                               * OPO: 1970-389-930

-------