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.
------- |