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