TECH. REPT. NO. 8
   CONCEPTS AND EQUATIONS
  FOR MULTILAYERED, VARIABLE
DENSITY ESTUARINE  HYDRAULICS
  ENVIRONMENTAL PROTECTION AGENCY
            REGION III

-------
     ^eosrx,
  US EPA Region in
    1650 Arch St.
Philadelphia, PA 19103

-------
   Concepts  and  Equations


 For Multilayered,  Variable


Density Estuarine Hydraulics
              by




       Robert L.  Grim


          May 1971


    Technical Report  No.  8
                Region III->^
                Centr.r v^":'•'";•'  .

-------
Introduction:




     The versatility of the "channel/junction" hydrodynamic and




quality models is well established.  Their abilities as predictive




tools have been verified in many widely varying locations.




     However, there are two basic assumptions inherent in the models




that limit their uses.  The first and primary assumption both hy-




draulicly and for quality is that of complete mixing vertically in




the junctions.  The second assumption is that the junctions are




connected with channels that are horizontal.




     It is the aim of this paper to present the derivation of the




equations of motion and continuity for a multilayered system of




channels and junctions in which channels slopes and density differences




play an integral part.




     Only with such a system may the effects  of the vertical placement




of heated discharges or of subsurface currents on contaminant distri-




bution be predicted.

-------

-------
                       TABLE OF CONTENTS
INTRODUCTION                                          1

DERIVATION OF THE EQUATIONS                           2

     Horizontal forces                                2
       Inertia                                        2
       Pressure                                       3
       Gravity                                        3
       Friction                                       3
     Summation of forces                              5

     Vertical forces                                  6
       Gravity                                        8
       Pressure                                       8
     Summation of forces                              8

     Friction forces and shear stresses               9
       Functional form                                9
       Finite difference expression                  10

     Finite difference form of the equations
       of motion and continuity                      12

     Computational scheme                            13
       Channel variables                             13
       Junction variables                            13
       Relationships between variables               15
       Density computations                          16

SOLUTION OF THE EQUATIONS                            17

     Stability consideration                         18

     Computation Algorithm                           21

CONCLUSIONS                                          24

REFERENCES                                           25

-------

-------
Horizontal forces:
                              figure A.
      Consider  the  two  layered  slice dx long of width B.  (figure A)
      Assuming   3y.dx     to be very  small compared to y., the mass
of each layer is


             MI  =

             M2  =
                                                              (1)

                                                              (2)
      Each layer  is  acted upon by  the forces of gravity, friction,

 pressure  and  inertia.  These forces must sum to zero.  The forces
 acting  on  layer  (1) are:
 Inertia;             du1
      FT=
                                                             (3)

-------

-------
or

     p  = PiylBdxf gul +u
      I            9t

Pressure:
     F  =
      p
which reduces to

     F  =  -gBdx
      P
Gravity;
                                            3x
             2       9x                                      (6)
     F  = +plyiBSidxg                                        (7)
      G
Friction:

     F^ = Bdx(T + T,)                                        (8)
      F        o

     where   is the shear stress on the layer bottom and  o is

the stress on the layer top.
                                                           2
     In general  T. is assumed to be quadratic.  (T.  = p.fu.) where

f is assumed to be independent of the variation in u with respect

to time and u is the relative motion between layers.


                   ) =  -PlylBSidxg - gBdx  3(Piyi )- Bdx(T0+ Tl)

                                                             (11)

-------

-------
Divide by plyiBdx   to get


                                       i2l-  (V  T0           (12)
 9t        3x               9y     3x
Note  that if    1 is  expressed  in  terms  of  the  Chezy coefficient

and   To =  0

Tl =  p  _!_ u/u/                                              (13)
         c 2

and  plis assumed constant, equation  (12) reduces to .the equation

of Barre' St. Venant.


 Ul  + u   Ul  +g 	 =-gS  -g.  u/u/                         (14)
9t        BX      gx
where R,  is the hydraulic radius of the channel, approximately equal

to yj for a wide channel.


     The forces exerted on block (2) are arrived at in the same

manner.

     The pressure force is handled by defining pressures.

        = pigyi                                               (15)
      a

and
9y
     P  +  3Pa dx = (pl -K-A) g fyi + A dx)                (16)
          ~                      lY          }

-------

-------
     Thus the resultant pressure force on  (2) is
     p   +
                                 Xlx) + g (Pz ,^-dx)
                                                            9y
                                                dx
                                                            3*
                                                             2 dx
                                                              (17)
which reduces to
F  =  -Bgdx
            g    3x
                             3x
                                                              (18)
Summing all the horizontal forces on  (2) gives




                                  T,          2
     .
    + u
3t
        :3x
           =  -S
                 2g
                      g
                   3p2y2
                                         .
                                                              (19)
                                 3x
     Generalizing  equation  (19)  to  a  system of  N layers,  the equation



of motion  can be expressed  as
3u
  n + u  "u n
 9t
       n
        9x
                -sng-
g
9pnyn
2 3(ynPan-l ),9(ynPn)~
g 3x 9X
(Tn_1 + TJ
P y
n n
                        n  =  1,2,3,...N
                                                             (20)
     Note that at  the  first  layer  (i.e.,  at  the  air/water interface)



                                           represents  the variation
0  represents wind stress and     a
                                  9
                                  x
    of atmospheric pressure along  the  channel.





     Further discussion of the  term  S^ and  its  relationship  to  the



depth y^ is needed.

-------

-------
     Define h^ as the elevation of the layer interface at any  time



t and at any point x.



     The slope of the interface and the quantities y^ and h^ are



related by the expression.





     Sir .  =  _Q . 4- 9h-i
     3x           dx



     Successive layers are related by
                                                              (22)
     The shear stress terms  Ti are very important to the distribution



of flow and will be discussed at length in a later sec', ion.




Vertical forces:



     The vertical forces acting on the layers of the channel slice



fall into two categories according to their origin.



     The first force is the local vertical acceleration causing  the



stream lines along dx to bend.  At the free surface this ac cLeration



is maximum with a value of


     A.. = 32h                                                 (23)

          at2


     and the vertical velocity V.,.., is expressed by




     VMAX - &                                               (24>
            91

-------

-------
     The Eulerian equations of motion give the vertical acceleration


as


     \ =  'I 31  -g                                         (25)
            p 3Z


     where P  is pressure


     and  z   is the vertical coordinate J f the pressure distribution


is hydrostatic, the right hand side of equation (25) vanishes and no


vertical velocities may exist.



     At any rate, in a periodic wave, the net transport vertically


will be zero.  Thus vertical pressure variations from hydrostatic


will be ignored.


     Equations (23) and (24) do provide a logical and convenient


basis of computing the vertical dispersion characteristics of the


channel .


     The vertical velocity at any depth may be derived from (24)


by assuming a linear distribution as
           =
                                                             (26)
              y    3t


     and the vertical dispersion EV may be formulated as


     Ev = f(VMAX^)}                                         (27)
     and may be determined emperically from observed data.

     In the second category of vertical forces is the buoyant force;

either positive or negative.  Returning to block (1) in figure A,

-------

-------
    the buoyant force is the resultant of the forces of



Gravity:
                                                            (28)
and



Pressure:
    FP-
                      3P
                        al dx)
                       3t
Bdx
2 LO
+ Pa + andx
o 3x
Bdx
2 (29)
which reduces to



    FT, = Bdx |2P,, + 8Pal
                          dx - 2P,
                                                            (30)
                                       3x
The inertia force is expressed as



    FT = PlylBdx 	1
     1            dt
                                                            (31)
where Wi is the vertical velocity positive in the upward direction.




    Vertical shear has been neglected due to the small velocities


involved.



    Summing forces vertically:
    FT + F  + F_ = 0
     1         Or
    PlylBdx dwl  = -PI yiBdxg +

             dt
                                Bdx  2

                                 2  L
                                            3P.
                                                            (32)
                                                            9P
                                                                 dx
                                                             (33)
and dividing by  p^Bdx  gives:




                           3P_
                                            3P
                              l dx - 2P
                                       a  " -g-Q-
                                        0    °X
                                                  dxl
                                                             (34)

-------
dw.
     Generalizing for the nth layer of an N layer system,
                       8P
dt
  n = -g +
2p y
  n n
               an dx -
                              3P,
                                     dx
                                                       (35)
     and the P  terms are as defined in equations (15)  and  (16)
             3.
     Friction forces and shear stresses:
     Von Karman's modification to Prandtl's mixing length theory
of shear stress in turbulent flow gives shear  stress at  any point
Z in  the vertical as
              2   2
             ^2)
             dZ
                                                       (36)
                     71
            yi
A7     	I
AZi     	j~
1          r
                     Y2
         ii_ri
                     Y3
                      U2
                      U3

                                            TI

                        figure B.

-------
    where K is a universal constant equal to 0.4




    Applied to the three layered system of figure B,
     n
        2   (_^a)
              dZ
    In finite difference form
                                                               10
                                                            (37)
^E  =  un "un-l

dZ       AZn_i
                                                            (38)
    and
        un ~un-l
                      - un "
     dZ
             AZ
               n-1
                      AZ
                            n
where
          AZn-l  +  AZn




          y _1 + y




          ~   ~
                                                            (39)
                                                            (40)
    and  Azn = y^ + yn+1


               2      2
                                                        (41)
Then




    T  = o K2 I"Un~U n-1 1  /
     n    n   I	  I /

              L AZn-l   J/
                                      -  un"u
                          AZ
                                n-1
                                          AZ
                                            n
                                 ^ ,   + AZ
                                 n-1      n
                                                            (42)
 where Tn acts to make velocity differences between layers smaller.

-------
                                                               11
    Equations (20) and  (42) can be combined to describe the

horizontal motion of the fluid in the elemental slice.  Equation

(42) gives an approximate picture of the interchange between layers

due to buoyant effects.

    The velocity gradient 3un  in equation (20) can be transformed

                          IhT

by the use of the continuity equation;

    _9A  + 3 (AD)  =  0                                       (43)
    3t      3x

    where A is the cross section yB and B is the constant width.

    1A  + u M  + A3u                                       (44)
    9t      3t     3x
    JW  =  -1 PM  +   3AJ
    9x      A [jit    U 3xJ
    and since A=yB ; 3A =
                     3x    3x
(45)
    and

    9un =  - _1  ^n - ^n ^n                                (46)
              n  at

    where again
                 9,
                   n
                                                             (47)
    ox


    Substitution of equation (47) into equation (20) makes a more

tractable equation for numerical solution.

    8un =  _!    n + ^2.  ^n  _ Sng
           An   at   B  3x

                     2
                      111     I T-    ,  ^  1
                                                             (48)
                                  pnyn

-------

-------
                                                              12
The finite element forms of the equations;



    Equation (48) can be expressed as

                                                        2

^n=_^^An + %_yn-s g -   g   j" 2 A(ynpan-l)  + A(yn Pn

At    A^  At    BL    n    8~ d  L g      L          L
                              pn n






           -1 vi+ g
                                                          (49)

                ~ d
                n n
    where
      d  = layer thickness at the channel midpoint
    p~  = average density in the channel

     n




  Equation  (35) becomes




AW                r

—-£•  = ~g +   1    9P3  + AP   -  8P     - AP_   I         (50)

At          2p y  I    n      n       n~1      r~
              n n
                                                  .I
    The continuity equation is written as
    Ay,,

   _ °  = ^Inflows - IQutflows                             (51)
                 A
                  sn

-------

-------
                                                             13
    Figure C represents three junctions,  I,  J, K  connected by




channels M,  N.  There are three layers in this example.
(D
M
O
N
                      figure C.






    Symbols  retain the same meaning,  but will now be doubly subscripted.




(i.e., yj £,  refers to the thickness  of layer 1 let junction I, and




Uj^ £  refers to  the velocity in layer 1 in channel N.)

-------

-------
                                                               14
    Variables are separated into junction variables and channel




variables as follows.




Channel variables:




    UN £                      horizontal velocity
JN,
     N
     HT I
     N,
Junction variables:
                          horizontal  flow







                          shear  stress







                          layer  crossectional  area at  the




                          channel  midpoint




                          layer  thickness at the  channel




                          midpoint




                          channel  width  (constant for  all




                          layers)




                          channel  length (constant for all




                          layers)




                          slope  of the interface  between layers




                          in the channel.




                          average  channel density











                          vertical velocity




                          vertical flow




                          layer  head




                          density




                          bottom elevation above  some  assumed datum

-------

-------
                                                               15
    A             =           surface area


    y-r £          =           layer thickness


    Vj £          =           layer volume


    Pa  £         =           applied pressure  at  the  top  of  the layer.
      *- »

                              (P     represents the  atmospheric pressure)
                                 J-»0
    Certain basic relationships between  the variables  can be ex-


pressed as:




                    >                                        (52)



                                                             (53)
         zi + HI,£+I  + yi.a                                 (56)
                                                             (58)


                                                             (59)




                                                             (60)

-------

-------
                                                               16
Density computations:



    Normally when a fluid is said to be incompressible, the



assumption is made that the fluid density is constant.  However



incompressibility means that volumes do not change with pressure,



not that density is constant.



Conservation of density;



    Since, in the multilayered equations of motion, density is



an integral part, the time and spatial variations must be computed



simultaneously with flow.



    This is different from the normal method of computing the hydraulics



of a system and the quality separately.



    The driving force for the density currents in most estuaries is



salinity.  That is to say horizontal salinity gradients causes a



sloping of the isopleths of pressure with distance.



    In this treatment density will be considered to be a conservative



quantity with the same characteristics as TDS or chloride.  Thus the



rate of change in density of a layer is merely the sume of all the



rates of contribution or estraction by advection and dispersion.
                         out         in     HK   out
  At



          ~EEHA ^- ~ E]EvAs — V VJ,L                      (61)
                 L           y



    where ER and EV are horizontal and vertical dispersion coefficients.

-------

-------
                                                              17
    In figure C, equation (61) can be written for layer 2 of

junction J as:
 Ap
   J.2 =
 At
(QM,
                               - (Pj'T2 " ^ AN,2 EHN,2]
                                     LN
                                            "rVsjEvj,;
                                              /
                                                           (62)
      The dispersion coefficients can be determined emperically

from observed data.  Since the time step length is generally

quite small, dispersion is probably not a significant factor

in the density calculation.
Solution of the equations:

    Solution of the horizontal motion equation (49) , the vertical

motion equation (50), the continuity equation (51) and the density

equation (61) must be done by a numerical integration technique,

-------
                                                              18
subject to the following boundary and initial conditions.


    The variables of at least one junction must be known for all


time.  This junction is called the driving head and is normally


the downstream or ocean end of the estuary.  At all other junctions


only the flow rates, and the density of the inflow need be known.


    Initial values of the variables at each junction are required.


The accuracy of their determination will affect the rate of con-


vergence of the solution, but will not affect the final values.


    The stability criterion for a one dimensional solution of the


hydrodynamic equations is that
    /  gd    ^ -^-                                          (63)




    where D is the depth in the channel.



    This equation relates the celerity of a long wave in shallow


water to the channel length and travel time.


    To take into account the tidal current as well as wave


celerity equation (63) can be written as
     V + / gd
                       At                                   (64)
    In the stratified case V in equation (64) contains both


horizontal and vertical velocity, or:
                               --
        +w        s    'MAX     At

-------
                                                                19
    A complicating factor in the choice of time step and channel


length is the problem of numerical mixing or numerical dispersion.


This is a common problem in the solution of differential equations


by finite elements and criteria have been established which can


minimize the effect.



    Stated for this problem the criteria is
                                                             (66)
    Assume W and Z are sufficiently small compared to U and L


so that the equation reduces to



    * =   uA t  =  1                                        (67)
           L


    Examination of equations (67) and (64) shows that the two


differ by the wave celerity term
    * +  /gd  At     = 1                                    (68)
             L


    By defining two separate time steps (At^) and (AtD) for the


hydraulic and densimetric solutions respectively, equation (68)


can be rearranged to give the ratio


    At
    _p  _ 11 r L   _  j—- \ i
    At

-------
                                                                20
    If the layer depth  (d) is assumed  to be  fairly  constant  over



a tidal cycle the terms within the parentheses will be  constant.
    r  -  v'id" )  -  K                                       (70)
    n
    and the ratio (!„) will be a function of advective velocity
                    K.



and the constant K.
    TR = K/U                                                 (71)





    The term T^ can range from a very large number  (U-*0)  to  a value



of 1 during the tidal cycle.  In addition, in large  scale networks



TR may vary significantly from channel to channel and  layer  to layer.



    Referring to the subscripting procedure of  figure  C,  the array  of



T-n values is expressed as






                    £                                        (72)

-------
                                                                21
The Algorithm;
    Given initial values of the channel and junction variables,

1.  Compute   Un,&  for each channel layer at  time  t=0
             At

2.  Project ahead to predict new U_  «  at  time  t +  At.
                                  11 , X,               ——


3.  Compute Qn £ , Tn,jj, f°r time t + At
                                       2

4.  Compute AWj    for each junction layer.

            ~~At

5.  Project ahead to predict new WT  .  at  time  t+At
                                    '*-           2

6.  Compute qj £ for time t + At
                              2

7.  Compute Ay     for each junction layer based  on flows  from
             At

    3. and 6.

8.  Compute new HT 0 , S  „  and P_T 0  for time  t+At.
                 j_,x    n, x       «*!,*•                —r



9.  Compute new  U  £   and average with  the values of  step  1.
                  n«
                 At

10. Compute Un £ for time t+At.

11. Compute Q_ n , T_ n for time t+ At.
             11 y *s    11 y **>

12. Compute new ^Wj-^ and average with those of step 4.
                 At'

13. Compute qj £  for time t+ At.

-------
                                                               22
             A
14. Compute   yl,£   for each junction based on flows from 11 and

            At

    13.




15. Compute H  £ , Sn £ , and Paj £ for time t+ At.
16. Compute new values for A  £ , d  £  , V-r £





17. Transfer arrays to be initial conditions for next time step,
18. Compute values of  /,„    for each channel layer
and add to
    previous value to form  t + At
19. Examine     /^R  o  to see ^ it is greater than or equal to 1.0.





    a.)  if L "*"/TR  £  is less than 1.0, compute and store Qn £Pn  £





    using old junction layer density and subtract from Qn £Pn £



    computed using present junction layer density.  Store the



    difference.



    b.) if 2, /Tp  „  is greater than or equal to 1.0, compute
               Rn, *•


    Q  ,,p  „  and add all the differences computed in (a) .  Set



    the 1  VTR    term equal to 0.0.





20. At each junction layer compute the new junction density, and



    store for next cycle.



21. Return to step 1.

-------

-------
                                                               23
    Steps 18 through 20 indicate one way of using the TR  «



terms.  What is being done is to insure that the affect of a



new upstream junction layer density will not be felt in the



downstream junction layer until the flow carrying that new



density has had time to arrive at the downstream junction layer.



    An alternative method which may be suitable for small scale



networks is to compute the network average of 1/T^  « and bypass



the density computation for all junctions until the    VTR  0
                                                          «-n,x'


term is greater than or equal to 1.0.  The density computations



can then be .made using averaged flows over the total time period.



This second method is the most efficient from a memory and time



standpoint.

-------
                                                              24
Conclusions;







    The equations of motion and continuity have been derived and




applied to stratified estuary networks.  An algorithm has been




presented which can be used to solve the equations.  This algorithm




must not be considered complete or totally correct in its present




form.  Nor is the treatment of the various terms (i.e., friction,




mass transfer) necessarily in final form.




    As is commonly the case, the methods and techniques used by the




computer program will reflect the limitations in time and space




imposed by the machine and the numerical difficulties encountered




by the programmer.  In addition, a sensitivity analysis on each




parameter may indicate areas where simplifications can be made or




where elaboration is needed.




    It must be emphasized here that this method is to give an




hydraulic solution.  Although density is a primary factor in the




solution, no indication is given as to the type of source of the




densities.  Once the vertical and horizontal flow patterns are




known for the system, a quality model can predict the concentrations




of various constituents at each junction layer.

-------
                                                               25


                          References
(1)  Sverdrup,  H.  U.,  Martin W.  Johnson,  and Richard H.  Fleming.
     The Oceans,  Their Physics,  Chemistry and General Biology.
     Prentice-Hall, Inc., New York,  1946.

(2)  Grim, Robert L.
     Mathematical Models of Estuarine Hydraulics and Water Quality.
     Duke University,  Durham, N.C.,  1970.

(3)  LeMehaute',  Bernard
     An Introduction to Hydrodynamics and Water Waves.
     ESSA Technical Report ERL 118-POL 3-2,  1969.

(4)  Lamb, Sir  Horace.
     Hydrodynamics.
     Dover Publications, New York,  1945.

(5)  Harlow,  Francis H. and Anthony  A. Amsden.
     Fluid Dynamics.
     Los Alamos Scientific Laboratory
     Los Alamos,  New Mexico, 1970.

(6)  Welch, J.  Eddie  et al
     The MAC  Method, A Computing Technique for  Solving Viscous,
     Incompressible,  Transient Fluid-Flow Problems Involving Free
     Surfaces
     Los Alamos Scientific Laboratory
     Los Alamos,  New Mexico, 1966.

(7)  Streeter,  Victor  L.
     Fluid Mechanics.
     McGraw Hill  Book  Company, New York,  1971.

-------