- Email: [email protected]
Implicit solution for estuarine hydrodynamic modelling
Mathematical Mode|ling, Vol. 3, pp. 543-558, 1982 Printed in the USA. All rights reserved.
0270-02551821060543-16503.00[0
Copyright © 1983 Pergamon Press Ltd.
IMPLICIT S O L U T I O N FOR E S T U A R I N E HYDRODYNAMIC MODELLING MYRON
H.
YOUNG
Center for Wetland Resources Louisiana State University Baton Rouge, Louisiana 70803, U S A Communicated by J. C. J. Nihoul
Abstract--This is a two-dimensional, finite difference, vertically averaged hydrodynamic model for shallow water estuaries. It utilizes three time levels and solves for current components in alternating direction exclusively by an implicit scheme. The purpose of using the implicit approach is to minimize the instability problem by which the usual implicit-explicit alternating scheme is inherently plagued when grid size is small or when water depths are deep. This is due to the explicit scheme limitation, whose time step is governed by the Courant-Friedrich-Lewy criterion. By the implicit-implicit scheme, a large time step would greatly reduce the CPU time. An effort was also made to model the effects of canals in the estuarine system. A different grid system than Leendertse's was used to lessen the lateral smoothing in Chezy coefficients caused by the large gradient of water bottom features. This grid system has both water elevation and depth measured at the center of the grid and is especially applicable for uneven grid-size models. The computer program is written to ensure easy-to-use features so that this model can be readily applied to different estuaries. Louisiana coastal features are dominated b y shallow estuarine systems. Their hydrodynamic and ecological characteristics are balanced a m o n g Gulf of Mexico tidal influences, wind tides, and the upland watershed. T h e s e delicate balances have also been upset by h u m a n activities in canals and pipeline proliferation. The resulting saltwater intrusion has b e e n a m a j o r factor in wetland destruction. One of the tools used to understand wetland interacting p r o c e s s e s is mathematical modelling. M a n y h y d r o d y n a m i c models [1] have been developed in recent years as attempts to simulate some of the processes. Both finite difference and finite element [2] models have b e e n successful. Some are easier to use than others. Most of these finite difference models [3, 4, 5] e m p l o y an alternating direction, implicit-explicit scheme in solving the differential equations of m o m e n t u m and continuity. The time steps chosen are critical to their computational stability. Instability b e c o m e s quite a problem when grid sizes are narrow or w h e n w a t e r depths are deep. The working time step, governed by the C o u r a n t - F r i e d r i c h - L e w y (CFL) criterion, in some cases m a y have to be reduced to just a few seconds. The resulting small time step, in turn, costs both in c o m p u t a t i o n time and cumulative error propagation. The finite difference model we d e v e l o p e d is similar to the L e e n d e r t s e - G r i t t o n [6] a p p r o a c h , which is a two-dimensional, vertical averaged, three-time leveled, implicitimplicit c o m p u t a t i o n in alternating directions. Since the C F L criterion applies only to the explicit step, the implicit-implicit scheme allows greater flexibility in the choice of time step. H o w e v e r , a grid s y s t e m other than L e e n d e r t s e ' s was used in an effort to better model the effects of canals in the estuarine system. In this model, both the water
Presented at the Third International Conference on Mathematical Modelling, July 1981, University of Southern California, Los Angeles. 543
5~
MYRON H.
YOUNG
elevation and the water depth are measured at the center of the grid square while in L e e n d e r t s e ' s they are, respectively, measured at the four corners and at the center of the grid. The use of a single point for both the elevation and the depth is more convenient when defining a canal, and the Chezy coefficients tend to be less smoothed out by the large gradient of the water bottom features. For a homogeneous system, the continuity and the m o m e n t u m equations, after being vertically integrated, are given, respectively, as follows: a~ + O H U + a H V ax
= o
au au aU _ f v + a?, u VUr-7a--i- + u - - f f + v g + g HC 2
aV
aV
a--t- + U ~ where
(1)
pH = o
a V + I u + g a~ v'x/-U-r-+-~ ~~ = o , -~ + g HC 2 - pH
+ V ay
(2) (3)
= water elevation with respect to a reference water level, H=h+~,
h U V f
= water depth from reference water level to the bottom, = vertically averaged velocity in x-direction, = vertically averaged velocity in y-direction, = 2II sin ~b = Coriolis parameter, = Earth's angular velocity in rad/sec, ~b = center latitude of model area in radians, g = gravitational constant, p = water density, ~s = wind stress at water surface (x and y represent components), C = Chezy function for bottom friction, = ( 1 . 4 8 6 / m ) H ~16(Manning's formula with H in ft), m = Manning's roughness coefficient. Several approximations were made for the above equations. Because of the low velocity gradients of the estuarine system, the diffusive transport of momentum, i.e., the second-order derivative of velocity terms, can be neglected. T h e r e f o r e , both the laminar and eddy viscosity terms are dropped. The atmospheric pressure gradients are negligible within the modelled area, and rainfall and evaporation are not considered and are dropped. For simplicity, the water density is considered to be homogeneous and unchanging, which gives rise to a first-order approximation in the salinity distribution study. In this model, the finite difference approach to solving these partial differential equations divides the time variable into three levels. T h e y are designated as n - 1, n, and n + 1 levels. In the first half time step, the known variables are ~,-I, U,-1, ~, and V". The continuity equation (1) and the U m o m e n t u m equation (2) are solved implicity for ~+1 and U "÷~, one row at a time until the entire modelled area is covered. Advancing one time level, at the second half time step, the known variables are ~", V", ~"÷t, and U "÷~. The continuity equation (1) and the V m o m e n t u m equation (3) are solved implicity again for ~,+2 and V "÷2, one column at a time. The procedure of alternating directions for e v e r y advancing half time step in U and V is repeated for all subsequent calculations until the allotted time is reached. At the first half time step, the difference formula for the continuity equation and the U m o m e n t u m equation are given, respectively, as follows:
Estuarine hydrodynamic modelling ~+1-- - A t
fin
~ij
+
V ~ - I -~- V ~
~T.n.+l __,.[
n-1 --
n-I
--n
n
= 0
n-, __ u n - I
(4)
C;+I
U,+,,, Ui-,,j+ V*" Ui,~+, i'~-'-fV*" +~x-x [ 2Ax 2Ay
2At :'"+'
545
u n + l __ l~n ]rTn+l --n n ui+Lj i+l,JAx " ~ u o " ' i J q Hvij+~Vi'j+t-HvliVqAy
"-']
_ si-Li + ~i-~,~j +g 2
U~ + ~ + U ~ -1
2
+
,.2
-.
C~-I
2
-.~
.
X/(U~-~) ~ ( V i i ) I(Hu~Cu,)- ~'~/(PH u~) = 0, (5)
where the overbar denotes a two-point spacial average and the superscript * denotes a four-point spacial average. Since the grid system of this model is different from Leendertse's, the definition of these spacial averages are also different. A comparison of these averaged terms are tabulated in Tables 1 and 2. The major difference lies in the definitions of H and C, through which the averages of this model lessen the lateral smoothing effect for better canal modelling. To solve Eqs. (4) and (5) implicity, they are rewritten by grouping the coefficients as follows: on+l
--[- rv !"/" n+!
n+t
(6)
--a~i ]i -- ~ i j ~ i j + a~ij = Aij /3...~ +' ¢~+] a ,~.+' IJ~l] + "4- I"~i+l,jU i+t,j --~" Bib where
(7)
At
a =gAx a 0
( V i..i ) /2( H -u.j C v-.2 )
.-1 = 1 + AAxt (Ui+,,jU i. L- 1j ) + g A t V ( U ~ ) 2 +
At
~'J =
a---x FI ~,~
-
~n 2 -- n --n 2 n-1 A 0 = [1 - gAt~/(Uijn--I ) 2 + (Vo)/(H%Cv,j)]Ui~ At-.t
+g_A_~(~,Tt,j_~_t)+EAt.x/(pffI.vo)+[2fAt At
.
Bij = ~,j - ~
-.
.
( H v i + V ,j+ ~ - / 4
_ ~At (U,.j+I"-~- U,"j2~O] V*"
~,,jV ,])
Similarly, at the second half time step, the difference equations for the continuity and the V momentum equations are, respectively, q ,[~. Ui+IjLd l + l J
At V~+2 -
Vii",, + ~[T*n+lij
2A t yn+2
-- .t~LOii" LYq
q"
.
V i + 1'1
_
]
--si'J-=+~i'J-lJ+g 2
ax
Vq v i i
,, Vi-l'J -~" --q vn+2
.
Vi'j+l -
. Vi'j-1 _4- z~f[T*n+l
2A y
"
V~+2+ V~ .x/trT,,+l.t2 + gl[l.n~2l(l~n+|~,(n+l)2~
"~'-'0
2
= 0
Ay
2Ax n
Vi,l+l~'i,j+l
Ax
J
(8)
o r r ~[ +"~q2 + r . "aq z
+ ~y s
--
n+l'~
~--ip,~"v,j "~v,j J-'ry/(pHvoj-O.(9)
Rearranging them by grouping the coefficients: •" l t r n+2 -I-
- - L¢ S i,j--I -
I vn+2.4_
~ij--O
t Trn+2
OL ij V ij
n+2
"~" a t ~ ; ' j +2 = A i j
/~t
i7n+2
t
t
- ;~ij :4" ~,ij+l--ij+~ = Bij,
(10) (1.1)
546
MYRON
H. Y O U N G
Table 1. Leendertse-Gritton's grid system (3 time levels: n -½, n, n + ~) Ax j+l
r
Vi.j+(l/2)
(+) (0) (-) (I)
hi+(ll2).i+(ll2)
,~ Water level h Depth UVelocity V Velocity
I Ay
0
~i~
Ui +(I/2).j
0
0
i-1 i-I
i+1
The 4-point average and the 2-point averages for the point (i + ~, j): ,
l
Vi+¢1/2).i = ~ (Vi,j+(I/2) + Vi.i 0/2) + Vi+l.i+(I/2) + Vi+l,i o/2)) AU
I [rrn+(i/2)
ITn_(I/2) X
A ~ - = TA ~ L/i+(l/2).j -- ~ i+(I/2).j] AU
1
2Ax
irrn
(1/2)
rrn-(l/2) x
2Ax ~ ' ; i+~3/2)j- .-~ i-~l/2~,jy
AU
1
2Ay
{¥Tn-(II2 )
rr a-(l/2)
X
2Ay ~LJ i+(l/2)'j+l- Ll i+(l/2)'j-l]
1 ¢lrn+(l/2) ± rrn_(i/2 ) x ~Jt ~ ~ g w i+(I/2).j T W i+(I/2).j] ~
l
--
__ - -
Ax
~ yn-l/2
/bi+l'j/Yn+(I/2)2 hi+l.]
~+('/2)
2 ~ ~j--(l/2)) -,,
Ax
/~.
1
--n
l
°
ui+o/:l.i = ~ (hi+tl/2).j+(IF2) + hi+ol:),j-o/2) + ~i+l./+ ~ ) n
where C~
1.486 F 1
q--
where
| - ~ ( h , (i/2)d+(i/2)+ h,+(i/2).j+(i/2,h,+(l/2).j (1,2)+ h,-(l/2,.j (1/2))+
m
L~
"l
¢,~/116 3
At a ' = g -A-~ At a~ i = 1 + -~y
[J -
.
(Vi,j+
1-
V~,j_I) --1-g A t V ( U * " + I )
2 + ~,-ijltX"n~21trC-tn+lT"(n+lVli )2xt~,"' ~-~Vii ]
At fi.+l Ay v,j
A~j = [1 - g A t ~ v / ( U ~ " + l ) 2 + r~ ,x1.~:/r r ~ . + l r ~ ( . + l ) ~ l ~r. v i i J I~,~tx Vii 1~-'Vij ] J V i j At
B'~j = W ' -
n
s
A-Zt ',q"+' "~"+' AX
~"
Ui+l,j Ll i+l,j -
-n+l
~ " + ' U ~ ÷') Uii
[
At
- .
_ Vn
] l[..+ .
547
Estuarine hydrodynamic modelling Table 2. Grid system used in this model (3 time levels: n - 1, n, n + 1)
1¢
Ax
)
(+) (0) (1) (-)
Water level h Depth U Velocity V Velocity
Ay
j+l
Vi. j + i
~ii
U~j
U i+ i.i
,
V~j i
i+1
The 4-point average and the 2-point averages for the point (i, j): ,
1
Vi, =~(Vi I.,+l+ Vi.,+,+ Vi,+ Vi tj) AU 2At
1
.-t
.+,
2At(U~'
AU_
-U~j
1 "U" I
2Ax
2~X
AU 2Ay
1 2Ay
~
) .-I
i+ld-
Ui-ld)
n-I
n-I
u,-=~-I(uU'+u~') n+l
Ax
n
Ax
n+l
n
2
and --n
1
n
n
Cu,+,.~ = ~ (C~j + C~+,.~) _ 1.486
~),/6
-~m [(hi,+
r"
-,1/61
+(hi+i.i+si+,.jJ
J
The boundary condition treatment for each of the alternating direction steps is the same. For illustrative purposes, only one set of boundary conditions is discussed here. There are four possible cases involved in solving the implicit equations. 1. O p e n - o p e n ends
- 1 denotes an open end, 0 denotes a closed end, 1 denotes a wet grid element.
548
MYRON H. YOUNG
According to the grid system used by this model, the U velocity is at the left boundary of the grid element. The unknown velocities located at the × marks are, respectively, U2, U3, and U4. They are to be solved implicitly along with the water elevations ~2 and if3. The driving forces are the tidal inputs at grids 1 and 4. The tidal inputs are designated as ~] and ~4. The implicit equations for this 4-grid system are derived from Eqs. (6) and (7), dropping jth row subscript: I'a2U] + a~] = A2+ aff~ i = 2//32U2+ ~ + / 3 3 U 2= B2 2
2 2 2 2 2 ~- / 3 3 U 3 Jr- if3 -}-/34U4 = B 3
(a) (b)
i = 3 [ - a ~ 2 + a3U3+ a~3 = A3
(c) (d)
i = 4 { - a t 2 + c~,U~ = A n - OC2.
(e)
The unknowns are with superscript 2, i.e., time level n = 2, and are solved by the following tridiagonal matrix system:
J':~/
I a j [A2a1 1
/33
/33
B2
~ /34 /~,~/ -,, ,,4 Lt,4j
B; ,,,;-,,c~J
2. Open-closed ends
,1' 1
'
I /A
2
3
4
Grid 4 is closed. U4 and Eq. (e) are dropped from the above equations, yielding the following matrix system:
~,
-a
]lIA2+a1 /~2/
A~
3. Closed-open ends j
r//////A ~/_/~
1
1
-1
1
2
3
4
Grid 1 is closed. U2 and Eq. (a) are dropped in this case, yielding the following matrix system:
-a
or3
a
JU~|=
A3
Estuarine hydrodynamic modelling
549
4. Closed-closed ends
1
2
3
4
Both U2 and U4 and Eqs. (a) and (d) are dropped, yielding the following matrix system:
o33
/33
L #~J
[B2 IA3 B3
Whenever tidal input is not specified, a flow input should be given at the open ends. The flow velocities at the boundary are set prior to the solving of the matrix system. The procedure of the implicit-implicit scheme gives staggered solutions for U and V. With the known values at time level n - 1 and n, the modelled area is solved row by row for velocity U and water elevation ~v at the first time level n + 1. They are solved implicitly according to the boundary conditions at each row. At the second half time step, i.e., the next time level, the known values are at time levels n and n + 1. The modelled area is solved again, this time column by column, for velocity V and water elevation ~v at time level n + 2 according to the boundary conditions at each column. The procedure repeats itself every two, subsequent, half time steps. Graphically, the alternating scheme can be represented as follows:
Time level n- 1
C'uU ' ~
n n+ 1
/j~,, v'
~ ~tiU ~ _ , ~
n +2 n+3
~ " ~
n +6
~vV
~ ' v U ' ~ ~~
n+4 n+5
(known)
~'vV' ~vU
~
~vV
The implicit-explicit scheme, however, works with only two time levels. With the known values at time level n, the modelled area is first solved implicitly row by row at time level n + 1 for velocity U and water elevation ~. It is then solved explicitly column by column for velocity V at the same time level. At the second half time step, the modelled area is solved again implicitly column by column for velocity V at time level n + 2. It is then solved explicitly row by row for U at the same time level. Graphically, the alternating scheme can be represented as follows:
550
MYRON H. YOUNG
Time Level
n
&
V.
U.
n+l
~
U~V
n+2
G
V,~U,
n+3
~
U~V
n+4
G
V,~U,
(known) ( $ ) implicit ("*) explicit
Since the time step size for the explicit scheme is restricted by the CFL criterion, the step size is usually four- to eighffold smaller than that of the implicit-implicit scheme so that the computational instability may be avoided. Therefore, the implicit-implicit scheme is in general about four to eight times faster than the implicit-explicit scheme solely because of the savings resulting from the large step size. Both computer programs for implicit-implicit and implicit-explicit schemes are verified against a one-dimensional theoretical case. A rectangular prismatic channel of length L and depth h is used in the test. The standing wave equations are OZU
02U
w - - gh ~ - r a2~ _
a2~
~ F - gh ~ and the standing wave solution is U(x, t) =
,~(x, t) -
h cosac (to L ) sin[to ~L( - ~x- 1 ) ] a L cos(to C ) c o s [ t o ~ ( L - 1 ) ]
cos(tot)
sin(tot),
where a is the amplitude of the input sinusoidal tide with a period T, C = X / ~ and to = 2rrlT. A 30,000-meter-long channel with one end closed is divided into 13 grids. A sinusoidal tidal input of one meter amplitude is placed at grid 1. The theoretical velocity U2 is plotted against the output of both schemes on Fig. 1. The outputs of the implicit-implicit and the implicit-explicit schemes are identical and are plotted with the same symbol. However, both models have a built-in bottom friction function whose Manning's roughness coefficient is set at 0.026. The result is therefore shifted slightly to the right of the theoretical curve, as expected. A one-dimensional implicit model with friction was also written. The result is plotted alongside the two-dimensional models. The result of a one-dimensional model, written by Harleman and Lee [7], is also plotted for a general comparison. Figure 2 shows a water elevation comparison at the closed end, grid 13, between the theoretical prediction and results of the implicit-implicit model. This model has been applied to several geographic locations in order to check out the computer program. One of these locations is the Atchafalaya-Vermilion Bay area
Estuarine hydrodynamic modelling Tide
1
~.p,,~-,--q
,:;
2
551
3
13
F I I I I I I I I I I 1 I
!~
30,000 M
'"
o
4
~o~
i-]x d~ (HR.)
^
\ x
5-
~
Z
/
I ow
o o
Z~ o~
+
o
2~D imp-imp (w£th friction) l~e imp (with frtetlon) 1 D llarleman and Lee.
\
/ ~
/ ~
o
/
Fig. 1. Velocity at x = 2307 m vs time.
O.
C3
Tide~ I I I i I I I i I I l i~'l
//
!-~
3o.o0o.
4
~(:3 0
0 C)
(-DO_ ~CO
TIM
_j 01 00
2[ 00
II o X o
q. 00 I
tHR. )
6
0 I
~-
8.00 I
I
"
OLO ~o
(D I~od. EEl
-7~_
1",.4
CD I
--
Theoretical
d
2-D lmp-lmp.
6 C3
c~ o I
Fig. 2. Water elevation at last grid of channel vs time.
552
MYRON H. YOUNG
Fig. 3. Atchafalaya-Vermilion Bay.
(Fig. 3). After steady state is reached, the mass conservation is checked at the beginning and ending of a 24-h tidal cycle. The result is found to be less than a 0.002% error. The C F L criterion dictates a At <<.As['X/2ghmax for a stable time step in an explicit calculation. For the Atchafalaya-Vermilion Bay area, the implicit-explicit model had to use a maximum half time step of 112.5 sec to avoid running into instability, whereas the implicit-implicit model may have an equivalent time step of 450 sec. This is a fourfold gain in computer CPU time. Some of the results are shown in Figs. 4 through 9. An attempt was made to add a salinity dispersion feature to this model without rewriting the computer program. A first-order estimate of the result was obtained by assuming the water density to be constant and neglecting all sources, sinks, rainfall, evaporation, seepage, and slip of the dispersed substance. The results for the Atchafalaya-Vermilion Bay area under the no-wind condition are plotted and shown in Figs. 10 through 13. In Fig. 10, salinities are initially set uniformly to 10 ppt throughout the area. Figures 11 through 13 show the salinity distribution after steady state has been reached. Inputting data is relatively simple for this computer program, and that simplicity is one of the advantages of this model. For example, changing the geometry of the plot from one estuary location to another is accomplished essentially by changing one input card. Detailed explanations for input cards will be provided in a User's Instruction Manual (to be published).
Estuarine hydrodynamic modelling
553
/~)~))-~i71
ATCHA.
~,6, s 6 ~
BAr
If ~ ~ ~ ~ ~ ~,~,I,/ 5
6
?
~b~e.~f
8. ~
I
f
.1
I
-1~ ~
Z
5~
B
g
~
te
1
~.IFZ'•
3
"~
~lt
5
6
7
9 ~
lO II ~ ~
12 12 13 l¼ I~ 1~ L~ ~ , ~1' ~ ' ~ " ~1' . V ~V I v ~ /m ~v ~ ~ ~d~ ~d~ ~ ~ [6 16 l~ J6 IG IG I~ 1~, 16 17 l e
I
I
l
1¼
I
~
l
l
111
8
8
lo
12
It
/~'7"
f / ' l l
3
1111
VERMIL.
CURRENTS (CM/SEC)
~
__~,
7
6
~
-
ESTUARY
1
I 1-
1-
I
1-
Z7
1-
1-
3
O
•
r
9
R,
•
/~
9
7'
8
5
9
§
0
-4~ / "
5
9 fl
~
,~'
B
f
f
t / l
3 3 /'
187
0
7`
F
[
[
f
1-
7`
F
~
,~ l
p
7"
/!
7'
7"
1" 1" f /" 1" [ ~ 2 3 I 5 3 3
611116
$
/~
9 6 / ~a TM
~ j ~ / l ~ / ~ 7 ` / ~ / . f f @ [
I
Fig. 4. Current distribution by implicit-implicit scheme.
~-~,.,
,, ~ ,~
If ~_ ~ ' . ' . , ~ , - , / 0
1
|
1
I
9
ATCHA.
~,,a~'~
BAY
~ 7
(f,f,#,py,~,o,
CURRENTS
_
"~
~/
/ . . . . . . . . .
!2..~ 7~ T"~"-""Z~
•
f T - - . ,~ .~ T /" T T "f "I
~ ~ f_~'..'_[_~.,"~
-
T
/f_ ~ 8 Y ! i P f ' t " f' f ' V
VERMIL,
ESTUAAT (CM/SEC)
NO W I N D DAY HOUR
q
MIN.
0
q
L
t3 3
q
I
I
I
l
I
I
I
~
"it
~
I
17
l
17
I
I
I
15
18
119 18
l
I~
I
2
I
I
"
I
1
27
~
1
O
8
9
1
8
6
~
f,
fl
2;
1
1
1
I
1
_~, j
~
,q
17
17
8 t6
Fig. 5. Implicit-explicit results for comparison with Fig. 4.
lS
13
554
MYRON H. YOUNG
ATCHA. - VEAMIL. BAY ESTUARY
m ~ gs ='6 ~"s~s ~s ~ '¢s
(CM/SEC)
CURRENTS
+, 40 ¢, ¢,,'0,',"°%~
NO NIND
+,1~,~,2¢,,~,Y,T,,.', : i. i~s+,d,~¢,Y,~. ?i-'~ 4" ~ Zq I
DAY HOUR MIN.
¢. ¢. ,ts ,~ ,t~ ,t~ '('t
~, +,,~22~,,¢12"111
q
12 0
/I ~'
o
V V
4, ~+ ~, ~s ~, ~+ 1, 4o",~',,v2~'
,, ~/2,v2,~,'~,+21,, 4s",(+,,q~2,.',,~.,+
I~o¢,+¢,212,v, ¢,,'qi ,~o G+,~+ +~ "~,,~,+¢,,v,ov,++o
p p f2 It 12 P ~+ ,~ )" ;+ •~ ~_+.3.j.f.z ~ , / , ' q ~ 2 ~ , ° + . 4 , ~ . ~ , . 4 . 4 , 4 o 4 , ~ , . ' q , " , / , # " ,~2.~ . ,~. . .+. . '+. ++. "+ + r , , . . ,o ,., . .,~.,. . . ,~. . ,~. ..+. .-. ,,-22,.22. ..~,,+ ... . . .~';''; /r+ +,+ ..,,, ~, .~ #'~ +,~,q ~, ~, ,~ + .',,.',.+,. +°+..'+i,+,r,:",~,:.,ll , , > . . + .,,, +T< <~ + ,3, , + + + ,
r~S )13 )i I
. I~. . . i./ . . I~./
i#
.
.
" .Y"'P
"st
.
+p,,++,.,
.
),' ~ ,+ ~ ,+ ,+ ~, .., ~. ~, L~, , + ,11, +o ,e,',+.~,o ,, ~,,,-',, ,+ ~ ~ ~' ,+ ]+ ~ ~. ~, ~, ~, ~ __,.+ ,~ ,q ¢0 ,,212.',°+,°,~,.,,.'+,22,,,',,.',,,.'++,, 8
9
7
;"
/"
++" ; '
/~ ;"
)"
r)'
;+ ]"
7
6
6
~e, "~,
-!
f:!
"4
+++~ ; e . ~ . . . + . . . + . . . ~ . + e _ . ~ 9
9
e
o
8
Fig. 6. Current distribution at 4th day, 12th hour, with no wind.
ATCHA.
-
BAY g
I
X
T
~sq
9
7
g
~ ~/
|
CURRENTS(CM/SEC)
9 ,~....
~
11
15• 2 2 ~
l
i~'._~g" 1 r ~ _ ~ _ - - 8 - - 6 " ~ ' 3
.'1
~
~
20KI20D OAT HOUR
~"
14]ND q
12 O
MIN. 16 26 25 20 2 ~ 2 ~ 2 ~ 3 ~ 1 ~ 3 ~
~ "ft.. ~lt._zll._lS~.l_'i,.. ,,~J~s~9 t~" - j
~
~
~
4-- ~
~
i-.. i . . t - - 4 - o 4 . _ is i t ~29 9--,
# ,o
i
~ 9
.,._9 ~z ~2- ~ a
,i.-. ~ . ?
•~ //
~,,e ~s
/VV
~jl ,t~s II
A 2 9 ~ 2 8 1 2 7 ~ 2 8 1 2 9 1 2 9 = 2 8 ~ 2 7 = 2 6 = 2 5 = 2 3 - 2 3 5 ~ "~ ~ ~O~29ajN]~30~31h30h30h29.2e=27
29 2 q ~
~30~30~30~3Z~3Z~3ZL3Z~31L30--28
28
~30~32--32--33~33--33--$2--31--29 ~3~33~33~3~32~31~$0.29 32 3q
33 3]
31
a7 ~5
26.2~,23
29 27 25 2q 23 21
~31~3q~$3~31~2B=27=25=2u~_23._2~._20 ~32~35~$2--2g~27~2B=2q=23=~2 2
33
30
27
26 2q
VERMIL.
ESTUARY
23 22 21
~
4r~ ~
a ++, ~ " ~
1~"
"~h/ 4~"
~--1~'17
2~-2(~"
I~--
16
1~-19
16
23 2~f'-18
I6P"I61"-I~ 17
L7 18
L? 18
21
lq
~" 15 1~I ~
19
17
16
16
19 ]B
I?
1~
16 15
I
18
17
17
16
16 16
1617
18101919
17
17
17
17
16
17
101819
I
15
17
I~
17
17
I
I
2(I
19
19
L9 18
iOI7
15 17 17
18 l e
~7
i8
16
21r" 4-~ 4"+
16t~" 't~ ~/" 1~-2~1~
4~
1~2~
le
18
19 19
19 2 ~ 2 ~ ' -
lO
I
)9
19
1
~
19 19
i~I - ~
19 20" 2 P "
19
19 19
19
L§ 19
19
ig
le
16
eo ~ , , f , , p , ~ , , ~ . , ~ , , ~ + , , ¢ , , e . , d % , o . l o . . , % , % , ~ , ~ 2 o . ~ o . , % , , . , , ~ , ~ , ~ , % , ~ , % , k , % , ,
tqJ
2~t~
]~r 2 r
t9
20
18
Iq
1916
20
19 3e
26 21i~
1~
18
15=2~ 20
dC
dr"
~"
20 20~Or
19 20
~
~
4"
|
19
4r-
a(" 4"+
,
4~
~/ I("
2~~" ~
19
1~'-2~--2~-
15
19 18
16 L 9 1 8
17 IG
•. .,,. . . .10, k , , . .
Fig. 7. Constant wind results for comparison with Fig. 6.
10
\
Estuarine hydrodynamic modelling
555
z .," ~'.,,-'--, F , ,.~,,~:,,,,,,,.,},.~, p .dq•Cq.,,j~'4v •~q~4~_a
is's ¢" g' ~ )
FITCHFI VERMIL BFIY E S T U R R Y -
~
"
"
HouRCURRENTS (CM/SECI.o
~.4¢,#,,-',3",.". ~. ~ ~.2'e e e e e 't' t' ? _!J 5
g
5
9
7
.3 ~k 7 q~--8~'k9 ~',,-7 ~ .
.,-V_
"~o'~,'%o',.',.',.',o° h ~"
~q 2 3 ~
21
l~r & 7 1 6
/~,/,#',/. . . . . . .
15 1~ II 14 ~ t 9 13 ~ ~F ,,~" ~r" # ~F
~f
F5
15-~1-17
~L'.°._
'
,
~t",~ ,°¢3~J" )o ~
,,
¢'
al~ 23 22i~21
20
20
21 23 ,2~25 27 aS 2 6
...¢'...E.~. . . . . . . . . . . . K" ~"
30~30
30
~" ~" ~," ~-" ~
30
30
31
30
30
30
31
2q 2q 23 a~J~2a 24
23 2q 26 2;;)7 ~3 18 2 2 7 23
25
,°
.E.,.E.# ................... s" ~/ ~,- ~ ' . - " ~. ,,. ~" w"
31
30
*-" ~" ~," ~," ~
3t
30
30
29
I~
19 1 5 - 1 6 2 6 - 2 3 lgl? 1514 14 IG 10 16 16 1 2 . 1 1 9~1._ ,( ~" "" ~ ~'~ '~'- '~'- '~'- '~" * f ~" , ~ ..~ a l ' v 'm~
F././././././,E,#,#,#',#,°¢..~.,~..z.j.36.z.o~./,..'.r'.,.',.M ~..~.o¢.o%y.#.r.. . . . . 2
28
27
26
28
27
26
25
24
24
20
17
I
~ ~ 1o
~
',%(,~_,#,. 16
14
21~1~
~F~
~'.¢ ....... E ~ / ...... t.. ~ K
2~
24
23
~" ~
23
~3
~"
2~
24
63
~/
22
d.
tl~
~3
,
~11
20 ~S
\ Fig. 8. Current distribution at 4th d a y , 20th h o u r , with s i n u s o i d a l tidal input.
f! ~,W'~,.I,~[ q 4 66 ~ f~'~.~.~.~.~,,/~
/I~ ~' '~ #' ~' #" ~" ~' ~./ ff°f°t°J'°~}~,, 5
6
6
5
1
11 ' ' ~
7 o ~ ~ 71
''~
~a ~
RTCHFI. BRY
~
.
.
'J ~ 3, # # ? I~ I~ # ~ # '.l
NO NIND
I
.
"~ '~ ' ~ ' ~ r ,
,
,
,
.
.
.
.
.
.
.
.
8
.
¼
5
,
,
.
.
.
HOUR
r,
.
, 20~
qs~of,,t,s,pt, s#s ¢~ 6}
"~
.
ESTUflRY
CURRENTS (CM/SECJ
#( .
-'VERMIL.
J I /~
~
.
7
#+,} 8
)
9
7
U
5
0
)
0
8
0
19
20
6
3
7
II
8
g
7
,UA
I
0
I
3
!
8
~
~5
6
6
~
Fig. 9. Real tidal input for c o m p a r i s o n w i t h Fig. 8.
~
2
lq
I°
556
MYRON H. YOUNG
,o ,o ,o i . lO
10
10
~ ~
~I0
....
10
10
o ....
I0
/"
10
lO
~ 10 t o
tO
10
tO
It)
tO
I0
tO
10
|0
10
tO
10
I0
,o,o
,o
,o,o
10 10
10
10
t0
l0
10
t0
10
t0
l0
10
t0
10 10
10
10 I0
tO
10
ESTUARY
5ALINITT
,oo1,: :Oo :: ,,:; ~o zo ~o ~o to to to to tot/ to to ~o ~o ~o to ~o to toI/ tO
[0
tO
lO
tO
Io
1o
IO
IO
to
tO
tO
tO
tO
IO
10
IO
tO
tO
10
10
tO
tO
HOUR HIN.
l0
10
10
10
10 I0
I0
t0
10
Io
t0
tO
10
l0
10
l0
1()
10
10
10
v
I0-10
10
10
1o
l0
IO
tO
10
10
10
10
10
1O
t0
i0
10
IO
t0
10
1O 10
10
10
10
10
l0
10
l0
10
tO
1O 10
t0
10
10
l0
10
10
10
1O
tO
I0
10
10
10
I0
IO
I0
IO
10
10
l0
l0
10
l0
l0
10
l0
10
tO
I0
tO
tO
10
tO
10
10
10
//f'~
#
tO
10 l 0
,,
10
. . . . . . . . . . . . . . . . . . . . . . 10
(0/001
0 0/.~
ii ii !i iiiii ii iii/ t0
VERMIL.
NO HIND DAY l
/
!i. . .!i. . .i/. . . . . . . . . . . . . . . . . . . . . . . 10
-
BAY
~
I0
o o 10
ATCHA.
10
t0
t0
i0
10
10
l0
t0
10
10
10 10
10 t0
10 tO
10
10
l0
10
l0
10
l0
10
l0
tO
I0
tO I0
tO tO
tO
I0
tO
tO
10
l0
l0
10
10
10
10
l0
10
l0
l0
10
V, , , V ,
::: ::
10
I
t
It)
10
10
10
l0
10
1O
t0
10
l0
10
10
10
l0
10
10
l0
10
10
l0
10
l0
1O
t0
10
10
l0
l0
l0
tO
11) tO
10
10
tO
tO
ID
10
10
10
10
l0
10
10
IO
10
l0
to
10
1o
l0
10
l0
l0
l0
IO
l0
10
10
l0
l0
10
10 10
10
l0
l0
I0
10
l0
10
10
It)
i0
to
10
1()
10 10
l0
10
t0
10
t0
10
10
l0
10
l0
10
10
10
10
to
tO
1o
10
10
t0
10
10
l0
l0
i0
10
10
10
10
l0
J0
10
10
10
10
10
10
10
It)
tO It)
IO
It)
to
It)
10
It)
to
Io
1o
It)
It)
1o
It)
IO
It)
to
to
to
to
1010
It)
to
tO to
1o
10
It)
to
to
to
It)
tO
It)
Io
It)
to
It)
10
10
10
It)
It)
It)
It)
10
10
It)
It)
It)
10
10
It)
I0
Io
1o
It)
It)
ID
10
10
10
10 It)
It)
10
10
It)
It)
l0
It)
19
10 ~ k l 0
lt)
l0
It)
It)
It)
t0
It)
lid 10
it)
10
10
It)
It)
It)
10
It)
10
1O 10
t0
t0
tO
10
10
tO
l0
t0
tO 10
10
10
It)
10
It)
It)
tO " ~ l O
l0
it)
10
It)
10
It)
10
It)
10
1O
l0
It)
10
It)
10
t0
It)
It)
10
t0
10
tO
t0
|0
It!
10
It)
It)
It)
10
10
10
it)
l0
It)
It)
10
It)
t,,, It) it) It) tt)/
10
IO
It)
It)
it)
It)
IO
it)
I0
It)
It)
It)
it)
IO
It)
It)
1o
It)
it)
It)
It)
it)
Io
lO
In
It)
tO
to
1o
It)
it)
tO
it)
It)
tO
tO
~
10
I0
It)
It)
It)
10
10
It)
It)
It)
it)
It)
10
10
I0
It)
It)
It)
It)
It)
10
it)
it)
It)
I0
It)
10
It)
It)
IO
JO
I0
10
10
lO
10
I0
10
10
t0
It)
lt)
l0
t0
10
It)
It)
it)
It)
It)
10
It)
It)
l0
It)
It)
It)
10
It)
It)
It)
i0
It)
10
It)
It)
It)
Io
10
1o
It)
10
It)
10
It)
It)
It)
It)
10
tO
10
lO
I0
I0
It)
it)
it)
I0
It)
to
It)
It)
It)
It)
10
It)
it)
IO
It)
10 lO
lO
10
I0
It)
It)
It)
10
tO
It)
It)
Fig. 10. Initial uniform salinity setup.
II
10
10
10
10
10
l0
10
i0
//~1~
II
10
10
l0
10
10
tO
10
II
II
tl
l0
10
tO
10
10
I1
II
It
il
11
10
tO
tO
30
tO
10
10
10
l0
IQ
10
l0
10f
I1
II
II
II
l[
|l
|0
IO
Jl /
10
JO 10
10
10
I0
lO
IO
lO
IO
l
11
11
II
it
I1
11
11
Jl
it
lO
10
IO
10
I0
tO
10
tO
tO
IO(
II
10
]1
lJ
]I
J]
l[
tO
IO
lO
IO
lO
It)
IO
IO
lO
10 -10
IO
10
I1
11
II
II
I1
tO
10
10
IO
50
IO
I0
10
10
~Pl l ~ , I12 I1 k
% 12/
,, ,1 11 1~
'~2 ~111~1I I
It
/
12 13
l
)
l0
,o 1o 1o,o ,o , o , o , o I0
10
IO
I0
11 t)
iS
15
15
t6
17
17
I
15
15
16
[7
17
18
18
/
17
17
17
18
18
17
18
18
Ie
t8
19
19191919
18
19
191920202020
18
19
18
19
IS 18
19
tO
10
IG
IO
I0
10
tO
IO
II
I0
tO
It)
IO
IO
I0
I0
[0
II
10
I0
JO
Io
tO
tO
I0
10
II
it
10
l0
lfl
I(1
10
l0
10
It
it
II
11
tO
l0
1O
lq
Iq
|3
13
~,22222
13
IO
I09[~
8
7
6
6
t~
2 2 2 2 2
10
10
9
9
t)
7
7
6
5
3
~1
12
tl
to
to
9I
O
8
7
6
5
~
8
8
7
6
6
t!
I0
9
8
8
It
I09
8
8
8
7
17
17
16
15
lq
Iq
lq
le
18
[8
17
16
15
16
Ill
tq
13
ltl
I
19
lib
18
17
16
I~
13
131212
13131212
212120202019191817171616 21212120201919
1817171616
222221212020191918
17
171716
[5
iq
I
15
t5
IrA Iq
1312121l
IS
iS
I~
lt~
131312
Iq
lq
1~1 13
lq
Iq tq
l
i5
IG
201919181817171616
18
18
t7
8
6
7
7
9
6
7
7
3 3 !
I0
131312
II
II
I0109
9
8
8
7
13131215
I]
II
I0109
9
8
8
7
;
II
I1
I0
I09
9
8
8
7
7
II
II
10
I0
10
9
8
8
8
7
1¼13~212
II
Io
10
I0109
9
8
8
8
7
6 6 5 5
1¼131211
I0
I(]
1
9
8
8
7
7
6 5 6 6 6
t5
IrA 1 3 1 2 1 2
[5
t5
lq
t6
9 9
3 3 2
It
15
21202020191918181717161616
I111109
h
II
2221212020191916181717161615 222121202019191818171716
II
15
3 3 2 2 ~ 3 2 2
l!
is
515
O
12
t7
45 56 IS Iq
q 0
11
18
19191818171716
2 212120
13
VERMIL.
~
10
12
ESTUARY
HOUR HIN.
j
10 L-/
13131~
~0191919191817161
212020202020202Q 2121212i
15
-
lo 1o ,o 1o io
13121211 18
9
I0
2
-
BAY
SALINITY (0/00) NO HIND DAY q
tO
ltl
ATCHA.
13J12
Iz
I
9
9
Fig. I 1. Salinity distribution at 4th day, 4th hour.
]
I
I
2
|
t
2
2
I
~
Estuarine hydrodynamic modelling
557
VEAHIL. ESTUflRY
ATCHA.
// ,: ',', ,',' ,',' ',i ',', ','," '/kZ~'
I
BA~ lO 1O 1O 10 t0
I0
l0 )0
I0 t0
[0
[0
10
IO~
to
to
~o 1o zol
:'o,: :: :: ,o :: :: ::
,,
I
''
,I
11
tl
I$
11
11
10
I,
[I"II
[0
i
I
H ~
tl
tt
.
H
H
H
11
H
to 1o Eo Jo ~o
k
.... )1
II
to
10
I0
~
111
10 ] 0 1 0 1 0
l0
IO I 0 1 0 10 ~ ~
~0
10
i
!
!1
~
l0
l010
lO
tO l 0 1 0 1 0
l01010101010 Ii1 0
lO
12 12
15
15
15
IS
16
15
1717
II
15
16
lfi
17
17
16
!7
12
17
I!
18
IB
1tl
18181817
171718
12
16
HIN.
0
12
l01010
l0 t 0 1
11 I0 IO IO lO I0 10 11 tO I0 10 l0
12 1 2 1 1 1 o
L0
IO 1 0 1 0
12
11
I0
10 IO
12
12 I I
tO I0
0 3
0
I
2
2
2
2
3
3
~
2
1
3
3
2
I
II 10
1212
20 20 20 20
lg
[6171616
19191
212121212020201919
12 II II
18 I 1 1 1 7 1 6 1 6
I
I1118
I1117171616
l0
11
1212
I110
1212
I11010
I
2
3
0
0
2
I
I
0
0
0
0
I
I
I
3
2
2
2
I
I
I
3
2
2
2
I
I
2
~2
21
@t
2D
~a
20
19
19
18
18
17
17
16
16
12 IZ
II
tO
10
2
7~
~11 23
22
22
21
21
20
20
19
19
18
18
17
17
18
16
12 12 I I
10
lO
2
2
2.
2
2
25
2 o, 2
22
21
21
20
20
19
19
18
18
17
17
16
16
12 I I
11
I0
IO 9
9
?
7
6
6
5
2S
2 I& 23
2
22
21
20
20
19
19
16
18
17
17
16
16
12 I I
II
18
I0
11
9
6
7
7
6
5
25
211 24
23
22
21
20
20
19
19
19
L8
16
17
16
16
12
I0
I0
I0
9
9
8
7
7
6
5
S
25
21~ 2 o. 2tA 23
21
20
20
19
19
19
18
111 17
17
16
12 I1
I0
IO
I0 9
9
6
7
7
6
5
5
12 I I
I0
|
0
ii! i!ilii iiiii iiIi!ii lil i iilil}
1919191919191919191818 2120 ~3
IZ
q
I0
~i ~ ~i 0 1 ] ' iii 15
IS
IS 1 5 1 5 1 5 1 5 1 6 1
DflY HOUR
IO 1 0 1 0 1 0
010101010 I0
No
.o
............................ ~ 1 1 1 1 II 10 t 0 1 0 1 0 1 0 l01o to tO 1o 1o
~
-
S
Fig. 12. Salinity distribution at 4th day, 12th hour.
ATCHA, - VERM!L. BAT ESTUART S A L I N I T Y (O/OO) NO HIND
DRY
q
HOUR
20
HIN.
0
IO 10 15
Io
15 15
[7
17
18
15
15 16
16
16
17 L7
18
I~1 18
|~
17
17
18
18 18
18
Io
12 18
1~
11
I0
I0
;o
IO
Io ] )
12
11
It
I0
10
9
I
12
12
11
I0
IO 9
9
I
2
12
l]
I0
9
I
12 I I
9
I
9
18 18 18
15
17
16
IB
17
IG
16
16
12
18
17
17
16
16
2
18
18
17
17
16
18
18
17
17
16
16
I11 17
17
18
18
17
18
[B
18
IS
V
14
3
3
3
2
2
]
0
0
0
V
S
q
3
3
?
Z
1
0
0
0
~
5 S
~
3
2
I
I
0
0
6
0
12
II
I0
I
!
I
6
16
]2
II
II
I0 !
I
,~
6
~
I 3
3
~
2
2
16
12 I I
Io !
I
?
6
5
I,I
~
3
2
2
16
16
IZ
II
I1
I0 S
S
e
?
6
5
q
Z
2
16
16
IZ
I]
11 to I
9
8
7
2
2
17
15
16
12
12 I I
tO IO 9
9
8
7
7
6
fi
6
17
17
16
12
H
tO z0 l0 g
9
g
?
?
6
6
6
18
18
17
17
16
18
18
17
17
16
2
H
2
12
II I0
to
to e
g
9
8
7
7
6
5
6
S
12
II
io
IO 9
|
9
8
7
7
IS
6
5
5
g
g
e
6
6
5
I0
tO
Fig. 13. Salinity distribution at 4th day, 20th hour.
5
0
0
I
L
1 j
I
L
2
2
558
MYRON H. YOUNG CONCLUSIONS
T h e i m p l i c i t - i m p l i c i t s c h e m e is d e s i g n e d to r e d u c e c o m p u t e r C P U t i m e b y t a k i n g a l a r g e r t i m e s t e p as c o m p a r e d to t h e i m p l i c i t - e x p l i c i t s c h e m e . C o m p a r i s o n s w e r e m a d e w h e n t h e p r o g r a m w a s a p p l i e d to s e v e r a l d i f f e r e n t g e o g r a p h i c l o c a t i o n s . T i m e s a v i n g s v a r i e d f r o m f o u r - to e i g h t f o l d . T h e g r i d s y s t e m c h o s e n f o r this m o d e l will b e t t e r m o d e l c h a n n e l s a n d c a n a l s . H o w e v e r , t h e grid s i z e is still t o o l a r g e to t r u l y m o d e l a c h a n n e l . T h e n e x t l o g i c a l s t e p is to i n c o r p o r a t e v a r i a b l e grid s i z e s into t h e c o m p u t e r p r o g r a m . P r o v i s i o n f o r v a r i a b l e boundaries caused by flooding should also be made. T h e s a l i n i t y s t u d y o f this m o d e l s e e m s to g i v e a r e a s o n a b l y g o o d f i r s t - o r d e r e s t i m a t i o n . T h e s i m p l i c i t y o f u s i n g this p r o g r a m m a k e s t h e m o d e l a t t r a c t i v e to u s e r s . REFERENCES 1. Tracor, Inc., Estuarine modeling: An assessment. 16070 DZV, EPA Water Pollution Control Series, Washington, DC (1971). 2. J. D. Wang and J. J. Connor, Mathematical modeling of near coastal circulation. Report No. 200, MIT R75-19, Massachusetts Institute of Technology, Cambridge, MA (1975). 3. J. J. Leendertse, Aspects of a computational model for long-period water-wave propagation. RM-5294-pr, The Rand Corporation (1967). 4. S. Hacker, R. Pike, and B. Wilkins, Transport phenomena in estuaries. LSU SGP-RFL-73-1, Louisiana State University, Baton Rouge, LA (1973). 5. G. F. McHugh, Development of a two-dimensional hydrodynamic numerical model for use in a shallow, well-mixed estuary. LSU-T-76-008, Louisiana State University, Baton Rouge, LA (1976). 6. J. J. Leendertse and E. C. Gritton, A water-quality simulation model for well mixed estuaries and coastal seas: Volume II, Computation procedures. R-708-NYC, The Rand Corporation (1971). 7. D. R. F. Harleman and C. H. Lee, The computation of tides and currents in estuaries and canals. MIT Technical Bulletin No. 16, Massachusetts Institute of Technology, Cambridge, MA (1969).
0270-02551821060543-16503.00[0
Copyright © 1983 Pergamon Press Ltd.
IMPLICIT S O L U T I O N FOR E S T U A R I N E HYDRODYNAMIC MODELLING MYRON
H.
YOUNG
Center for Wetland Resources Louisiana State University Baton Rouge, Louisiana 70803, U S A Communicated by J. C. J. Nihoul
Abstract--This is a two-dimensional, finite difference, vertically averaged hydrodynamic model for shallow water estuaries. It utilizes three time levels and solves for current components in alternating direction exclusively by an implicit scheme. The purpose of using the implicit approach is to minimize the instability problem by which the usual implicit-explicit alternating scheme is inherently plagued when grid size is small or when water depths are deep. This is due to the explicit scheme limitation, whose time step is governed by the Courant-Friedrich-Lewy criterion. By the implicit-implicit scheme, a large time step would greatly reduce the CPU time. An effort was also made to model the effects of canals in the estuarine system. A different grid system than Leendertse's was used to lessen the lateral smoothing in Chezy coefficients caused by the large gradient of water bottom features. This grid system has both water elevation and depth measured at the center of the grid and is especially applicable for uneven grid-size models. The computer program is written to ensure easy-to-use features so that this model can be readily applied to different estuaries. Louisiana coastal features are dominated b y shallow estuarine systems. Their hydrodynamic and ecological characteristics are balanced a m o n g Gulf of Mexico tidal influences, wind tides, and the upland watershed. T h e s e delicate balances have also been upset by h u m a n activities in canals and pipeline proliferation. The resulting saltwater intrusion has b e e n a m a j o r factor in wetland destruction. One of the tools used to understand wetland interacting p r o c e s s e s is mathematical modelling. M a n y h y d r o d y n a m i c models [1] have been developed in recent years as attempts to simulate some of the processes. Both finite difference and finite element [2] models have b e e n successful. Some are easier to use than others. Most of these finite difference models [3, 4, 5] e m p l o y an alternating direction, implicit-explicit scheme in solving the differential equations of m o m e n t u m and continuity. The time steps chosen are critical to their computational stability. Instability b e c o m e s quite a problem when grid sizes are narrow or w h e n w a t e r depths are deep. The working time step, governed by the C o u r a n t - F r i e d r i c h - L e w y (CFL) criterion, in some cases m a y have to be reduced to just a few seconds. The resulting small time step, in turn, costs both in c o m p u t a t i o n time and cumulative error propagation. The finite difference model we d e v e l o p e d is similar to the L e e n d e r t s e - G r i t t o n [6] a p p r o a c h , which is a two-dimensional, vertical averaged, three-time leveled, implicitimplicit c o m p u t a t i o n in alternating directions. Since the C F L criterion applies only to the explicit step, the implicit-implicit scheme allows greater flexibility in the choice of time step. H o w e v e r , a grid s y s t e m other than L e e n d e r t s e ' s was used in an effort to better model the effects of canals in the estuarine system. In this model, both the water
Presented at the Third International Conference on Mathematical Modelling, July 1981, University of Southern California, Los Angeles. 543
5~
MYRON H.
YOUNG
elevation and the water depth are measured at the center of the grid square while in L e e n d e r t s e ' s they are, respectively, measured at the four corners and at the center of the grid. The use of a single point for both the elevation and the depth is more convenient when defining a canal, and the Chezy coefficients tend to be less smoothed out by the large gradient of the water bottom features. For a homogeneous system, the continuity and the m o m e n t u m equations, after being vertically integrated, are given, respectively, as follows: a~ + O H U + a H V ax
= o
au au aU _ f v + a?, u VUr-7a--i- + u - - f f + v g + g HC 2
aV
aV
a--t- + U ~ where
(1)
pH = o
a V + I u + g a~ v'x/-U-r-+-~ ~~ = o , -~ + g HC 2 - pH
+ V ay
(2) (3)
= water elevation with respect to a reference water level, H=h+~,
h U V f
= water depth from reference water level to the bottom, = vertically averaged velocity in x-direction, = vertically averaged velocity in y-direction, = 2II sin ~b = Coriolis parameter, = Earth's angular velocity in rad/sec, ~b = center latitude of model area in radians, g = gravitational constant, p = water density, ~s = wind stress at water surface (x and y represent components), C = Chezy function for bottom friction, = ( 1 . 4 8 6 / m ) H ~16(Manning's formula with H in ft), m = Manning's roughness coefficient. Several approximations were made for the above equations. Because of the low velocity gradients of the estuarine system, the diffusive transport of momentum, i.e., the second-order derivative of velocity terms, can be neglected. T h e r e f o r e , both the laminar and eddy viscosity terms are dropped. The atmospheric pressure gradients are negligible within the modelled area, and rainfall and evaporation are not considered and are dropped. For simplicity, the water density is considered to be homogeneous and unchanging, which gives rise to a first-order approximation in the salinity distribution study. In this model, the finite difference approach to solving these partial differential equations divides the time variable into three levels. T h e y are designated as n - 1, n, and n + 1 levels. In the first half time step, the known variables are ~,-I, U,-1, ~, and V". The continuity equation (1) and the U m o m e n t u m equation (2) are solved implicity for ~+1 and U "÷~, one row at a time until the entire modelled area is covered. Advancing one time level, at the second half time step, the known variables are ~", V", ~"÷t, and U "÷~. The continuity equation (1) and the V m o m e n t u m equation (3) are solved implicity again for ~,+2 and V "÷2, one column at a time. The procedure of alternating directions for e v e r y advancing half time step in U and V is repeated for all subsequent calculations until the allotted time is reached. At the first half time step, the difference formula for the continuity equation and the U m o m e n t u m equation are given, respectively, as follows:
Estuarine hydrodynamic modelling ~+1-- - A t
fin
~ij
+
V ~ - I -~- V ~
~T.n.+l __,.[
n-1 --
n-I
--n
n
= 0
n-, __ u n - I
(4)
C;+I
U,+,,, Ui-,,j+ V*" Ui,~+, i'~-'-fV*" +~x-x [ 2Ax 2Ay
2At :'"+'
545
u n + l __ l~n ]rTn+l --n n ui+Lj i+l,JAx " ~ u o " ' i J q Hvij+~Vi'j+t-HvliVqAy
"-']
_ si-Li + ~i-~,~j +g 2
U~ + ~ + U ~ -1
2
+
,.2
-.
C~-I
2
-.~
.
X/(U~-~) ~ ( V i i ) I(Hu~Cu,)- ~'~/(PH u~) = 0, (5)
where the overbar denotes a two-point spacial average and the superscript * denotes a four-point spacial average. Since the grid system of this model is different from Leendertse's, the definition of these spacial averages are also different. A comparison of these averaged terms are tabulated in Tables 1 and 2. The major difference lies in the definitions of H and C, through which the averages of this model lessen the lateral smoothing effect for better canal modelling. To solve Eqs. (4) and (5) implicity, they are rewritten by grouping the coefficients as follows: on+l
--[- rv !"/" n+!
n+t
(6)
--a~i ]i -- ~ i j ~ i j + a~ij = Aij /3...~ +' ¢~+] a ,~.+' IJ~l] + "4- I"~i+l,jU i+t,j --~" Bib where
(7)
At
a =gAx a 0
( V i..i ) /2( H -u.j C v-.2 )
.-1 = 1 + AAxt (Ui+,,jU i. L- 1j ) + g A t V ( U ~ ) 2 +
At
~'J =
a---x FI ~,~
-
~n 2 -- n --n 2 n-1 A 0 = [1 - gAt~/(Uijn--I ) 2 + (Vo)/(H%Cv,j)]Ui~ At-.t
+g_A_~(~,Tt,j_~_t)+EAt.x/(pffI.vo)+[2fAt At
.
Bij = ~,j - ~
-.
.
( H v i + V ,j+ ~ - / 4
_ ~At (U,.j+I"-~- U,"j2~O] V*"
~,,jV ,])
Similarly, at the second half time step, the difference equations for the continuity and the V momentum equations are, respectively, q ,[~. Ui+IjLd l + l J
At V~+2 -
Vii",, + ~[T*n+lij
2A t yn+2
-- .t~LOii" LYq
q"
.
V i + 1'1
_
]
--si'J-=+~i'J-lJ+g 2
ax
Vq v i i
,, Vi-l'J -~" --q vn+2
.
Vi'j+l -
. Vi'j-1 _4- z~f[T*n+l
2A y
"
V~+2+ V~ .x/trT,,+l.t2 + gl[l.n~2l(l~n+|~,(n+l)2~
"~'-'0
2
= 0
Ay
2Ax n
Vi,l+l~'i,j+l
Ax
J
(8)
o r r ~[ +"~q2 + r . "aq z
+ ~y s
--
n+l'~
~--ip,~"v,j "~v,j J-'ry/(pHvoj-O.(9)
Rearranging them by grouping the coefficients: •" l t r n+2 -I-
- - L¢ S i,j--I -
I vn+2.4_
~ij--O
t Trn+2
OL ij V ij
n+2
"~" a t ~ ; ' j +2 = A i j
/~t
i7n+2
t
t
- ;~ij :4" ~,ij+l--ij+~ = Bij,
(10) (1.1)
546
MYRON
H. Y O U N G
Table 1. Leendertse-Gritton's grid system (3 time levels: n -½, n, n + ~) Ax j+l
r
Vi.j+(l/2)
(+) (0) (-) (I)
hi+(ll2).i+(ll2)
,~ Water level h Depth UVelocity V Velocity
I Ay
0
~i~
Ui +(I/2).j
0
0
i-1 i-I
i+1
The 4-point average and the 2-point averages for the point (i + ~, j): ,
l
Vi+¢1/2).i = ~ (Vi,j+(I/2) + Vi.i 0/2) + Vi+l.i+(I/2) + Vi+l,i o/2)) AU
I [rrn+(i/2)
ITn_(I/2) X
A ~ - = TA ~ L/i+(l/2).j -- ~ i+(I/2).j] AU
1
2Ax
irrn
(1/2)
rrn-(l/2) x
2Ax ~ ' ; i+~3/2)j- .-~ i-~l/2~,jy
AU
1
2Ay
{¥Tn-(II2 )
rr a-(l/2)
X
2Ay ~LJ i+(l/2)'j+l- Ll i+(l/2)'j-l]
1 ¢lrn+(l/2) ± rrn_(i/2 ) x ~Jt ~ ~ g w i+(I/2).j T W i+(I/2).j] ~
l
--
__ - -
Ax
~ yn-l/2
/bi+l'j/Yn+(I/2)2 hi+l.]
~+('/2)
2 ~ ~j--(l/2)) -,,
Ax
/~.
1
--n
l
°
ui+o/:l.i = ~ (hi+tl/2).j+(IF2) + hi+ol:),j-o/2) + ~i+l./+ ~ ) n
where C~
1.486 F 1
q--
where
| - ~ ( h , (i/2)d+(i/2)+ h,+(i/2).j+(i/2,h,+(l/2).j (1,2)+ h,-(l/2,.j (1/2))+
m
L~
"l
¢,~/116 3
At a ' = g -A-~ At a~ i = 1 + -~y
[J -
.
(Vi,j+
1-
V~,j_I) --1-g A t V ( U * " + I )
2 + ~,-ijltX"n~21trC-tn+lT"(n+lVli )2xt~,"' ~-~Vii ]
At fi.+l Ay v,j
A~j = [1 - g A t ~ v / ( U ~ " + l ) 2 + r~ ,x1.~:/r r ~ . + l r ~ ( . + l ) ~ l ~r. v i i J I~,~tx Vii 1~-'Vij ] J V i j At
B'~j = W ' -
n
s
A-Zt ',q"+' "~"+' AX
~"
Ui+l,j Ll i+l,j -
-n+l
~ " + ' U ~ ÷') Uii
[
At
- .
_ Vn
] l[..+ .
547
Estuarine hydrodynamic modelling Table 2. Grid system used in this model (3 time levels: n - 1, n, n + 1)
1¢
Ax
)
(+) (0) (1) (-)
Water level h Depth U Velocity V Velocity
Ay
j+l
Vi. j + i
~ii
U~j
U i+ i.i
,
V~j i
i+1
The 4-point average and the 2-point averages for the point (i, j): ,
1
Vi, =~(Vi I.,+l+ Vi.,+,+ Vi,+ Vi tj) AU 2At
1
.-t
.+,
2At(U~'
AU_
-U~j
1 "U" I
2Ax
2~X
AU 2Ay
1 2Ay
~
) .-I
i+ld-
Ui-ld)
n-I
n-I
u,-=~-I(uU'+u~') n+l
Ax
n
Ax
n+l
n
2
and --n
1
n
n
Cu,+,.~ = ~ (C~j + C~+,.~) _ 1.486
~),/6
-~m [(hi,+
r"
-,1/61
+(hi+i.i+si+,.jJ
J
The boundary condition treatment for each of the alternating direction steps is the same. For illustrative purposes, only one set of boundary conditions is discussed here. There are four possible cases involved in solving the implicit equations. 1. O p e n - o p e n ends
- 1 denotes an open end, 0 denotes a closed end, 1 denotes a wet grid element.
548
MYRON H. YOUNG
According to the grid system used by this model, the U velocity is at the left boundary of the grid element. The unknown velocities located at the × marks are, respectively, U2, U3, and U4. They are to be solved implicitly along with the water elevations ~2 and if3. The driving forces are the tidal inputs at grids 1 and 4. The tidal inputs are designated as ~] and ~4. The implicit equations for this 4-grid system are derived from Eqs. (6) and (7), dropping jth row subscript: I'a2U] + a~] = A2+ aff~ i = 2//32U2+ ~ + / 3 3 U 2= B2 2
2 2 2 2 2 ~- / 3 3 U 3 Jr- if3 -}-/34U4 = B 3
(a) (b)
i = 3 [ - a ~ 2 + a3U3+ a~3 = A3
(c) (d)
i = 4 { - a t 2 + c~,U~ = A n - OC2.
(e)
The unknowns are with superscript 2, i.e., time level n = 2, and are solved by the following tridiagonal matrix system:
J':~/
I a j [A2a1 1
/33
/33
B2
~ /34 /~,~/ -,, ,,4 Lt,4j
B; ,,,;-,,c~J
2. Open-closed ends
,1' 1
'
I /A
2
3
4
Grid 4 is closed. U4 and Eq. (e) are dropped from the above equations, yielding the following matrix system:
~,
-a
]lIA2+a1 /~2/
A~
3. Closed-open ends j
r//////A ~/_/~
1
1
-1
1
2
3
4
Grid 1 is closed. U2 and Eq. (a) are dropped in this case, yielding the following matrix system:
-a
or3
a
JU~|=
A3
Estuarine hydrodynamic modelling
549
4. Closed-closed ends
1
2
3
4
Both U2 and U4 and Eqs. (a) and (d) are dropped, yielding the following matrix system:
o33
/33
L #~J
[B2 IA3 B3
Whenever tidal input is not specified, a flow input should be given at the open ends. The flow velocities at the boundary are set prior to the solving of the matrix system. The procedure of the implicit-implicit scheme gives staggered solutions for U and V. With the known values at time level n - 1 and n, the modelled area is solved row by row for velocity U and water elevation ~v at the first time level n + 1. They are solved implicitly according to the boundary conditions at each row. At the second half time step, i.e., the next time level, the known values are at time levels n and n + 1. The modelled area is solved again, this time column by column, for velocity V and water elevation ~v at time level n + 2 according to the boundary conditions at each column. The procedure repeats itself every two, subsequent, half time steps. Graphically, the alternating scheme can be represented as follows:
Time level n- 1
C'uU ' ~
n n+ 1
/j~,, v'
~ ~tiU ~ _ , ~
n +2 n+3
~ " ~
n +6
~vV
~ ' v U ' ~ ~~
n+4 n+5
(known)
~'vV' ~vU
~
~vV
The implicit-explicit scheme, however, works with only two time levels. With the known values at time level n, the modelled area is first solved implicitly row by row at time level n + 1 for velocity U and water elevation ~. It is then solved explicitly column by column for velocity V at the same time level. At the second half time step, the modelled area is solved again implicitly column by column for velocity V at time level n + 2. It is then solved explicitly row by row for U at the same time level. Graphically, the alternating scheme can be represented as follows:
550
MYRON H. YOUNG
Time Level
n
&
V.
U.
n+l
~
U~V
n+2
G
V,~U,
n+3
~
U~V
n+4
G
V,~U,
(known) ( $ ) implicit ("*) explicit
Since the time step size for the explicit scheme is restricted by the CFL criterion, the step size is usually four- to eighffold smaller than that of the implicit-implicit scheme so that the computational instability may be avoided. Therefore, the implicit-implicit scheme is in general about four to eight times faster than the implicit-explicit scheme solely because of the savings resulting from the large step size. Both computer programs for implicit-implicit and implicit-explicit schemes are verified against a one-dimensional theoretical case. A rectangular prismatic channel of length L and depth h is used in the test. The standing wave equations are OZU
02U
w - - gh ~ - r a2~ _
a2~
~ F - gh ~ and the standing wave solution is U(x, t) =
,~(x, t) -
h cosac (to L ) sin[to ~L( - ~x- 1 ) ] a L cos(to C ) c o s [ t o ~ ( L - 1 ) ]
cos(tot)
sin(tot),
where a is the amplitude of the input sinusoidal tide with a period T, C = X / ~ and to = 2rrlT. A 30,000-meter-long channel with one end closed is divided into 13 grids. A sinusoidal tidal input of one meter amplitude is placed at grid 1. The theoretical velocity U2 is plotted against the output of both schemes on Fig. 1. The outputs of the implicit-implicit and the implicit-explicit schemes are identical and are plotted with the same symbol. However, both models have a built-in bottom friction function whose Manning's roughness coefficient is set at 0.026. The result is therefore shifted slightly to the right of the theoretical curve, as expected. A one-dimensional implicit model with friction was also written. The result is plotted alongside the two-dimensional models. The result of a one-dimensional model, written by Harleman and Lee [7], is also plotted for a general comparison. Figure 2 shows a water elevation comparison at the closed end, grid 13, between the theoretical prediction and results of the implicit-implicit model. This model has been applied to several geographic locations in order to check out the computer program. One of these locations is the Atchafalaya-Vermilion Bay area
Estuarine hydrodynamic modelling Tide
1
~.p,,~-,--q
,:;
2
551
3
13
F I I I I I I I I I I 1 I
!~
30,000 M
'"
o
4
~o~
i-]x d~ (HR.)
^
\ x
5-
~
Z
/
I ow
o o
Z~ o~
+
o
2~D imp-imp (w£th friction) l~e imp (with frtetlon) 1 D llarleman and Lee.
\
/ ~
/ ~
o
/
Fig. 1. Velocity at x = 2307 m vs time.
O.
C3
Tide~ I I I i I I I i I I l i~'l
//
!-~
3o.o0o.
4
~(:3 0
0 C)
(-DO_ ~CO
TIM
_j 01 00
2[ 00
II o X o
q. 00 I
tHR. )
6
0 I
~-
8.00 I
I
"
OLO ~o
(D I~od. EEl
-7~_
1",.4
CD I
--
Theoretical
d
2-D lmp-lmp.
6 C3
c~ o I
Fig. 2. Water elevation at last grid of channel vs time.
552
MYRON H. YOUNG
Fig. 3. Atchafalaya-Vermilion Bay.
(Fig. 3). After steady state is reached, the mass conservation is checked at the beginning and ending of a 24-h tidal cycle. The result is found to be less than a 0.002% error. The C F L criterion dictates a At <<.As['X/2ghmax for a stable time step in an explicit calculation. For the Atchafalaya-Vermilion Bay area, the implicit-explicit model had to use a maximum half time step of 112.5 sec to avoid running into instability, whereas the implicit-implicit model may have an equivalent time step of 450 sec. This is a fourfold gain in computer CPU time. Some of the results are shown in Figs. 4 through 9. An attempt was made to add a salinity dispersion feature to this model without rewriting the computer program. A first-order estimate of the result was obtained by assuming the water density to be constant and neglecting all sources, sinks, rainfall, evaporation, seepage, and slip of the dispersed substance. The results for the Atchafalaya-Vermilion Bay area under the no-wind condition are plotted and shown in Figs. 10 through 13. In Fig. 10, salinities are initially set uniformly to 10 ppt throughout the area. Figures 11 through 13 show the salinity distribution after steady state has been reached. Inputting data is relatively simple for this computer program, and that simplicity is one of the advantages of this model. For example, changing the geometry of the plot from one estuary location to another is accomplished essentially by changing one input card. Detailed explanations for input cards will be provided in a User's Instruction Manual (to be published).
Estuarine hydrodynamic modelling
553
/~)~))-~i71
ATCHA.
~,6, s 6 ~
BAr
If ~ ~ ~ ~ ~ ~,~,I,/ 5
6
?
~b~e.~f
8. ~
I
f
.1
I
-1~ ~
Z
5~
B
g
~
te
1
~.IFZ'•
3
"~
~lt
5
6
7
9 ~
lO II ~ ~
12 12 13 l¼ I~ 1~ L~ ~ , ~1' ~ ' ~ " ~1' . V ~V I v ~ /m ~v ~ ~ ~d~ ~d~ ~ ~ [6 16 l~ J6 IG IG I~ 1~, 16 17 l e
I
I
l
1¼
I
~
l
l
111
8
8
lo
12
It
/~'7"
f / ' l l
3
1111
VERMIL.
CURRENTS (CM/SEC)
~
__~,
7
6
~
-
ESTUARY
1
I 1-
1-
I
1-
Z7
1-
1-
3
O
•
r
9
R,
•
/~
9
7'
8
5
9
§
0
-4~ / "
5
9 fl
~
,~'
B
f
f
t / l
3 3 /'
187
0
7`
F
[
[
f
1-
7`
F
~
,~ l
p
7"
/!
7'
7"
1" 1" f /" 1" [ ~ 2 3 I 5 3 3
611116
$
/~
9 6 / ~a TM
~ j ~ / l ~ / ~ 7 ` / ~ / . f f @ [
I
Fig. 4. Current distribution by implicit-implicit scheme.
~-~,.,
,, ~ ,~
If ~_ ~ ' . ' . , ~ , - , / 0
1
|
1
I
9
ATCHA.
~,,a~'~
BAY
~ 7
(f,f,#,py,~,o,
CURRENTS
_
"~
~/
/ . . . . . . . . .
!2..~ 7~ T"~"-""Z~
•
f T - - . ,~ .~ T /" T T "f "I
~ ~ f_~'..'_[_~.,"~
-
T
/f_ ~ 8 Y ! i P f ' t " f' f ' V
VERMIL,
ESTUAAT (CM/SEC)
NO W I N D DAY HOUR
q
MIN.
0
q
L
t3 3
q
I
I
I
l
I
I
I
~
"it
~
I
17
l
17
I
I
I
15
18
119 18
l
I~
I
2
I
I
"
I
1
27
~
1
O
8
9
1
8
6
~
f,
fl
2;
1
1
1
I
1
_~, j
~
,q
17
17
8 t6
Fig. 5. Implicit-explicit results for comparison with Fig. 4.
lS
13
554
MYRON H. YOUNG
ATCHA. - VEAMIL. BAY ESTUARY
m ~ gs ='6 ~"s~s ~s ~ '¢s
(CM/SEC)
CURRENTS
+, 40 ¢, ¢,,'0,',"°%~
NO NIND
+,1~,~,2¢,,~,Y,T,,.', : i. i~s+,d,~¢,Y,~. ?i-'~ 4" ~ Zq I
DAY HOUR MIN.
¢. ¢. ,ts ,~ ,t~ ,t~ '('t
~, +,,~22~,,¢12"111
q
12 0
/I ~'
o
V V
4, ~+ ~, ~s ~, ~+ 1, 4o",~',,v2~'
,, ~/2,v2,~,'~,+21,, 4s",(+,,q~2,.',,~.,+
I~o¢,+¢,212,v, ¢,,'qi ,~o G+,~+ +~ "~,,~,+¢,,v,ov,++o
p p f2 It 12 P ~+ ,~ )" ;+ •~ ~_+.3.j.f.z ~ , / , ' q ~ 2 ~ , ° + . 4 , ~ . ~ , . 4 . 4 , 4 o 4 , ~ , . ' q , " , / , # " ,~2.~ . ,~. . .+. . '+. ++. "+ + r , , . . ,o ,., . .,~.,. . . ,~. . ,~. ..+. .-. ,,-22,.22. ..~,,+ ... . . .~';''; /r+ +,+ ..,,, ~, .~ #'~ +,~,q ~, ~, ,~ + .',,.',.+,. +°+..'+i,+,r,:",~,:.,ll , , > . . + .,,, +T< <~ + ,3, , + + + ,
r~S )13 )i I
. I~. . . i./ . . I~./
i#
.
.
" .Y"'P
"st
.
+p,,++,.,
.
),' ~ ,+ ~ ,+ ,+ ~, .., ~. ~, L~, , + ,11, +o ,e,',+.~,o ,, ~,,,-',, ,+ ~ ~ ~' ,+ ]+ ~ ~. ~, ~, ~, ~ __,.+ ,~ ,q ¢0 ,,212.',°+,°,~,.,,.'+,22,,,',,.',,,.'++,, 8
9
7
;"
/"
++" ; '
/~ ;"
)"
r)'
;+ ]"
7
6
6
~e, "~,
-!
f:!
"4
+++~ ; e . ~ . . . + . . . + . . . ~ . + e _ . ~ 9
9
e
o
8
Fig. 6. Current distribution at 4th day, 12th hour, with no wind.
ATCHA.
-
BAY g
I
X
T
~sq
9
7
g
~ ~/
|
CURRENTS(CM/SEC)
9 ,~....
~
11
15• 2 2 ~
l
i~'._~g" 1 r ~ _ ~ _ - - 8 - - 6 " ~ ' 3
.'1
~
~
20KI20D OAT HOUR
~"
14]ND q
12 O
MIN. 16 26 25 20 2 ~ 2 ~ 2 ~ 3 ~ 1 ~ 3 ~
~ "ft.. ~lt._zll._lS~.l_'i,.. ,,~J~s~9 t~" - j
~
~
~
4-- ~
~
i-.. i . . t - - 4 - o 4 . _ is i t ~29 9--,
# ,o
i
~ 9
.,._9 ~z ~2- ~ a
,i.-. ~ . ?
•~ //
~,,e ~s
/VV
~jl ,t~s II
A 2 9 ~ 2 8 1 2 7 ~ 2 8 1 2 9 1 2 9 = 2 8 ~ 2 7 = 2 6 = 2 5 = 2 3 - 2 3 5 ~ "~ ~ ~O~29ajN]~30~31h30h30h29.2e=27
29 2 q ~
~30~30~30~3Z~3Z~3ZL3Z~31L30--28
28
~30~32--32--33~33--33--$2--31--29 ~3~33~33~3~32~31~$0.29 32 3q
33 3]
31
a7 ~5
26.2~,23
29 27 25 2q 23 21
~31~3q~$3~31~2B=27=25=2u~_23._2~._20 ~32~35~$2--2g~27~2B=2q=23=~2 2
33
30
27
26 2q
VERMIL.
ESTUARY
23 22 21
~
4r~ ~
a ++, ~ " ~
1~"
"~h/ 4~"
~--1~'17
2~-2(~"
I~--
16
1~-19
16
23 2~f'-18
I6P"I61"-I~ 17
L7 18
L? 18
21
lq
~" 15 1~I ~
19
17
16
16
19 ]B
I?
1~
16 15
I
18
17
17
16
16 16
1617
18101919
17
17
17
17
16
17
101819
I
15
17
I~
17
17
I
I
2(I
19
19
L9 18
iOI7
15 17 17
18 l e
~7
i8
16
21r" 4-~ 4"+
16t~" 't~ ~/" 1~-2~1~
4~
1~2~
le
18
19 19
19 2 ~ 2 ~ ' -
lO
I
)9
19
1
~
19 19
i~I - ~
19 20" 2 P "
19
19 19
19
L§ 19
19
ig
le
16
eo ~ , , f , , p , ~ , , ~ . , ~ , , ~ + , , ¢ , , e . , d % , o . l o . . , % , % , ~ , ~ 2 o . ~ o . , % , , . , , ~ , ~ , ~ , % , ~ , % , k , % , ,
tqJ
2~t~
]~r 2 r
t9
20
18
Iq
1916
20
19 3e
26 21i~
1~
18
15=2~ 20
dC
dr"
~"
20 20~Or
19 20
~
~
4"
|
19
4r-
a(" 4"+
,
4~
~/ I("
2~~" ~
19
1~'-2~--2~-
15
19 18
16 L 9 1 8
17 IG
•. .,,. . . .10, k , , . .
Fig. 7. Constant wind results for comparison with Fig. 6.
10
\
Estuarine hydrodynamic modelling
555
z .," ~'.,,-'--, F , ,.~,,~:,,,,,,,.,},.~, p .dq•Cq.,,j~'4v •~q~4~_a
is's ¢" g' ~ )
FITCHFI VERMIL BFIY E S T U R R Y -
~
"
"
HouRCURRENTS (CM/SECI.o
~.4¢,#,,-',3",.". ~. ~ ~.2'e e e e e 't' t' ? _!J 5
g
5
9
7
.3 ~k 7 q~--8~'k9 ~',,-7 ~ .
.,-V_
"~o'~,'%o',.',.',.',o° h ~"
~q 2 3 ~
21
l~r & 7 1 6
/~,/,#',/. . . . . . .
15 1~ II 14 ~ t 9 13 ~ ~F ,,~" ~r" # ~F
~f
F5
15-~1-17
~L'.°._
'
,
~t",~ ,°¢3~J" )o ~
,,
¢'
al~ 23 22i~21
20
20
21 23 ,2~25 27 aS 2 6
...¢'...E.~. . . . . . . . . . . . K" ~"
30~30
30
~" ~" ~," ~-" ~
30
30
31
30
30
30
31
2q 2q 23 a~J~2a 24
23 2q 26 2;;)7 ~3 18 2 2 7 23
25
,°
.E.,.E.# ................... s" ~/ ~,- ~ ' . - " ~. ,,. ~" w"
31
30
*-" ~" ~," ~," ~
3t
30
30
29
I~
19 1 5 - 1 6 2 6 - 2 3 lgl? 1514 14 IG 10 16 16 1 2 . 1 1 9~1._ ,( ~" "" ~ ~'~ '~'- '~'- '~'- '~" * f ~" , ~ ..~ a l ' v 'm~
F././././././,E,#,#,#',#,°¢..~.,~..z.j.36.z.o~./,..'.r'.,.',.M ~..~.o¢.o%y.#.r.. . . . . 2
28
27
26
28
27
26
25
24
24
20
17
I
~ ~ 1o
~
',%(,~_,#,. 16
14
21~1~
~F~
~'.¢ ....... E ~ / ...... t.. ~ K
2~
24
23
~" ~
23
~3
~"
2~
24
63
~/
22
d.
tl~
~3
,
~11
20 ~S
\ Fig. 8. Current distribution at 4th d a y , 20th h o u r , with s i n u s o i d a l tidal input.
f! ~,W'~,.I,~[ q 4 66 ~ f~'~.~.~.~.~,,/~
/I~ ~' '~ #' ~' #" ~" ~' ~./ ff°f°t°J'°~}~,, 5
6
6
5
1
11 ' ' ~
7 o ~ ~ 71
''~
~a ~
RTCHFI. BRY
~
.
.
'J ~ 3, # # ? I~ I~ # ~ # '.l
NO NIND
I
.
"~ '~ ' ~ ' ~ r ,
,
,
,
.
.
.
.
.
.
.
.
8
.
¼
5
,
,
.
.
.
HOUR
r,
.
, 20~
qs~of,,t,s,pt, s#s ¢~ 6}
"~
.
ESTUflRY
CURRENTS (CM/SECJ
#( .
-'VERMIL.
J I /~
~
.
7
#+,} 8
)
9
7
U
5
0
)
0
8
0
19
20
6
3
7
II
8
g
7
,UA
I
0
I
3
!
8
~
~5
6
6
~
Fig. 9. Real tidal input for c o m p a r i s o n w i t h Fig. 8.
~
2
lq
I°
556
MYRON H. YOUNG
,o ,o ,o i . lO
10
10
~ ~
~I0
....
10
10
o ....
I0
/"
10
lO
~ 10 t o
tO
10
tO
It)
tO
I0
tO
10
|0
10
tO
10
I0
,o,o
,o
,o,o
10 10
10
10
t0
l0
10
t0
10
t0
l0
10
t0
10 10
10
10 I0
tO
10
ESTUARY
5ALINITT
,oo1,: :Oo :: ,,:; ~o zo ~o ~o to to to to tot/ to to ~o ~o ~o to ~o to toI/ tO
[0
tO
lO
tO
Io
1o
IO
IO
to
tO
tO
tO
tO
IO
10
IO
tO
tO
10
10
tO
tO
HOUR HIN.
l0
10
10
10
10 I0
I0
t0
10
Io
t0
tO
10
l0
10
l0
1()
10
10
10
v
I0-10
10
10
1o
l0
IO
tO
10
10
10
10
10
1O
t0
i0
10
IO
t0
10
1O 10
10
10
10
10
l0
10
l0
10
tO
1O 10
t0
10
10
l0
10
10
10
1O
tO
I0
10
10
10
I0
IO
I0
IO
10
10
l0
l0
10
l0
l0
10
l0
10
tO
I0
tO
tO
10
tO
10
10
10
//f'~
#
tO
10 l 0
,,
10
. . . . . . . . . . . . . . . . . . . . . . 10
(0/001
0 0/.~
ii ii !i iiiii ii iii/ t0
VERMIL.
NO HIND DAY l
/
!i. . .!i. . .i/. . . . . . . . . . . . . . . . . . . . . . . 10
-
BAY
~
I0
o o 10
ATCHA.
10
t0
t0
i0
10
10
l0
t0
10
10
10 10
10 t0
10 tO
10
10
l0
10
l0
10
l0
10
l0
tO
I0
tO I0
tO tO
tO
I0
tO
tO
10
l0
l0
10
10
10
10
l0
10
l0
l0
10
V, , , V ,
::: ::
10
I
t
It)
10
10
10
l0
10
1O
t0
10
l0
10
10
10
l0
10
10
l0
10
10
l0
10
l0
1O
t0
10
10
l0
l0
l0
tO
11) tO
10
10
tO
tO
ID
10
10
10
10
l0
10
10
IO
10
l0
to
10
1o
l0
10
l0
l0
l0
IO
l0
10
10
l0
l0
10
10 10
10
l0
l0
I0
10
l0
10
10
It)
i0
to
10
1()
10 10
l0
10
t0
10
t0
10
10
l0
10
l0
10
10
10
10
to
tO
1o
10
10
t0
10
10
l0
l0
i0
10
10
10
10
l0
J0
10
10
10
10
10
10
10
It)
tO It)
IO
It)
to
It)
10
It)
to
Io
1o
It)
It)
1o
It)
IO
It)
to
to
to
to
1010
It)
to
tO to
1o
10
It)
to
to
to
It)
tO
It)
Io
It)
to
It)
10
10
10
It)
It)
It)
It)
10
10
It)
It)
It)
10
10
It)
I0
Io
1o
It)
It)
ID
10
10
10
10 It)
It)
10
10
It)
It)
l0
It)
19
10 ~ k l 0
lt)
l0
It)
It)
It)
t0
It)
lid 10
it)
10
10
It)
It)
It)
10
It)
10
1O 10
t0
t0
tO
10
10
tO
l0
t0
tO 10
10
10
It)
10
It)
It)
tO " ~ l O
l0
it)
10
It)
10
It)
10
It)
10
1O
l0
It)
10
It)
10
t0
It)
It)
10
t0
10
tO
t0
|0
It!
10
It)
It)
It)
10
10
10
it)
l0
It)
It)
10
It)
t,,, It) it) It) tt)/
10
IO
It)
It)
it)
It)
IO
it)
I0
It)
It)
It)
it)
IO
It)
It)
1o
It)
it)
It)
It)
it)
Io
lO
In
It)
tO
to
1o
It)
it)
tO
it)
It)
tO
tO
~
10
I0
It)
It)
It)
10
10
It)
It)
It)
it)
It)
10
10
I0
It)
It)
It)
It)
It)
10
it)
it)
It)
I0
It)
10
It)
It)
IO
JO
I0
10
10
lO
10
I0
10
10
t0
It)
lt)
l0
t0
10
It)
It)
it)
It)
It)
10
It)
It)
l0
It)
It)
It)
10
It)
It)
It)
i0
It)
10
It)
It)
It)
Io
10
1o
It)
10
It)
10
It)
It)
It)
It)
10
tO
10
lO
I0
I0
It)
it)
it)
I0
It)
to
It)
It)
It)
It)
10
It)
it)
IO
It)
10 lO
lO
10
I0
It)
It)
It)
10
tO
It)
It)
Fig. 10. Initial uniform salinity setup.
II
10
10
10
10
10
l0
10
i0
//~1~
II
10
10
l0
10
10
tO
10
II
II
tl
l0
10
tO
10
10
I1
II
It
il
11
10
tO
tO
30
tO
10
10
10
l0
IQ
10
l0
10f
I1
II
II
II
l[
|l
|0
IO
Jl /
10
JO 10
10
10
I0
lO
IO
lO
IO
l
11
11
II
it
I1
11
11
Jl
it
lO
10
IO
10
I0
tO
10
tO
tO
IO(
II
10
]1
lJ
]I
J]
l[
tO
IO
lO
IO
lO
It)
IO
IO
lO
10 -10
IO
10
I1
11
II
II
I1
tO
10
10
IO
50
IO
I0
10
10
~Pl l ~ , I12 I1 k
% 12/
,, ,1 11 1~
'~2 ~111~1I I
It
/
12 13
l
)
l0
,o 1o 1o,o ,o , o , o , o I0
10
IO
I0
11 t)
iS
15
15
t6
17
17
I
15
15
16
[7
17
18
18
/
17
17
17
18
18
17
18
18
Ie
t8
19
19191919
18
19
191920202020
18
19
18
19
IS 18
19
tO
10
IG
IO
I0
10
tO
IO
II
I0
tO
It)
IO
IO
I0
I0
[0
II
10
I0
JO
Io
tO
tO
I0
10
II
it
10
l0
lfl
I(1
10
l0
10
It
it
II
11
tO
l0
1O
lq
Iq
|3
13
~,22222
13
IO
I09[~
8
7
6
6
t~
2 2 2 2 2
10
10
9
9
t)
7
7
6
5
3
~1
12
tl
to
to
9I
O
8
7
6
5
~
8
8
7
6
6
t!
I0
9
8
8
It
I09
8
8
8
7
17
17
16
15
lq
Iq
lq
le
18
[8
17
16
15
16
Ill
tq
13
ltl
I
19
lib
18
17
16
I~
13
131212
13131212
212120202019191817171616 21212120201919
1817171616
222221212020191918
17
171716
[5
iq
I
15
t5
IrA Iq
1312121l
IS
iS
I~
lt~
131312
Iq
lq
1~1 13
lq
Iq tq
l
i5
IG
201919181817171616
18
18
t7
8
6
7
7
9
6
7
7
3 3 !
I0
131312
II
II
I0109
9
8
8
7
13131215
I]
II
I0109
9
8
8
7
;
II
I1
I0
I09
9
8
8
7
7
II
II
10
I0
10
9
8
8
8
7
1¼13~212
II
Io
10
I0109
9
8
8
8
7
6 6 5 5
1¼131211
I0
I(]
1
9
8
8
7
7
6 5 6 6 6
t5
IrA 1 3 1 2 1 2
[5
t5
lq
t6
9 9
3 3 2
It
15
21202020191918181717161616
I111109
h
II
2221212020191916181717161615 222121202019191818171716
II
15
3 3 2 2 ~ 3 2 2
l!
is
515
O
12
t7
45 56 IS Iq
q 0
11
18
19191818171716
2 212120
13
VERMIL.
~
10
12
ESTUARY
HOUR HIN.
j
10 L-/
13131~
~0191919191817161
212020202020202Q 2121212i
15
-
lo 1o ,o 1o io
13121211 18
9
I0
2
-
BAY
SALINITY (0/00) NO HIND DAY q
tO
ltl
ATCHA.
13J12
Iz
I
9
9
Fig. I 1. Salinity distribution at 4th day, 4th hour.
]
I
I
2
|
t
2
2
I
~
Estuarine hydrodynamic modelling
557
VEAHIL. ESTUflRY
ATCHA.
// ,: ',', ,',' ,',' ',i ',', ','," '/kZ~'
I
BA~ lO 1O 1O 10 t0
I0
l0 )0
I0 t0
[0
[0
10
IO~
to
to
~o 1o zol
:'o,: :: :: ,o :: :: ::
,,
I
''
,I
11
tl
I$
11
11
10
I,
[I"II
[0
i
I
H ~
tl
tt
.
H
H
H
11
H
to 1o Eo Jo ~o
k
.... )1
II
to
10
I0
~
111
10 ] 0 1 0 1 0
l0
IO I 0 1 0 10 ~ ~
~0
10
i
!
!1
~
l0
l010
lO
tO l 0 1 0 1 0
l01010101010 Ii1 0
lO
12 12
15
15
15
IS
16
15
1717
II
15
16
lfi
17
17
16
!7
12
17
I!
18
IB
1tl
18181817
171718
12
16
HIN.
0
12
l01010
l0 t 0 1
11 I0 IO IO lO I0 10 11 tO I0 10 l0
12 1 2 1 1 1 o
L0
IO 1 0 1 0
12
11
I0
10 IO
12
12 I I
tO I0
0 3
0
I
2
2
2
2
3
3
~
2
1
3
3
2
I
II 10
1212
20 20 20 20
lg
[6171616
19191
212121212020201919
12 II II
18 I 1 1 1 7 1 6 1 6
I
I1118
I1117171616
l0
11
1212
I110
1212
I11010
I
2
3
0
0
2
I
I
0
0
0
0
I
I
I
3
2
2
2
I
I
I
3
2
2
2
I
I
2
~2
21
@t
2D
~a
20
19
19
18
18
17
17
16
16
12 IZ
II
tO
10
2
7~
~11 23
22
22
21
21
20
20
19
19
18
18
17
17
18
16
12 12 I I
10
lO
2
2
2.
2
2
25
2 o, 2
22
21
21
20
20
19
19
18
18
17
17
16
16
12 I I
11
I0
IO 9
9
?
7
6
6
5
2S
2 I& 23
2
22
21
20
20
19
19
16
18
17
17
16
16
12 I I
II
18
I0
11
9
6
7
7
6
5
25
211 24
23
22
21
20
20
19
19
19
L8
16
17
16
16
12
I0
I0
I0
9
9
8
7
7
6
5
S
25
21~ 2 o. 2tA 23
21
20
20
19
19
19
18
111 17
17
16
12 I1
I0
IO
I0 9
9
6
7
7
6
5
5
12 I I
I0
|
0
ii! i!ilii iiiii iiIi!ii lil i iilil}
1919191919191919191818 2120 ~3
IZ
q
I0
~i ~ ~i 0 1 ] ' iii 15
IS
IS 1 5 1 5 1 5 1 5 1 6 1
DflY HOUR
IO 1 0 1 0 1 0
010101010 I0
No
.o
............................ ~ 1 1 1 1 II 10 t 0 1 0 1 0 1 0 l01o to tO 1o 1o
~
-
S
Fig. 12. Salinity distribution at 4th day, 12th hour.
ATCHA, - VERM!L. BAT ESTUART S A L I N I T Y (O/OO) NO HIND
DRY
q
HOUR
20
HIN.
0
IO 10 15
Io
15 15
[7
17
18
15
15 16
16
16
17 L7
18
I~1 18
|~
17
17
18
18 18
18
Io
12 18
1~
11
I0
I0
;o
IO
Io ] )
12
11
It
I0
10
9
I
12
12
11
I0
IO 9
9
I
2
12
l]
I0
9
I
12 I I
9
I
9
18 18 18
15
17
16
IB
17
IG
16
16
12
18
17
17
16
16
2
18
18
17
17
16
18
18
17
17
16
16
I11 17
17
18
18
17
18
[B
18
IS
V
14
3
3
3
2
2
]
0
0
0
V
S
q
3
3
?
Z
1
0
0
0
~
5 S
~
3
2
I
I
0
0
6
0
12
II
I0
I
!
I
6
16
]2
II
II
I0 !
I
,~
6
~
I 3
3
~
2
2
16
12 I I
Io !
I
?
6
5
I,I
~
3
2
2
16
16
IZ
II
I1
I0 S
S
e
?
6
5
q
Z
2
16
16
IZ
I]
11 to I
9
8
7
2
2
17
15
16
12
12 I I
tO IO 9
9
8
7
7
6
fi
6
17
17
16
12
H
tO z0 l0 g
9
g
?
?
6
6
6
18
18
17
17
16
18
18
17
17
16
2
H
2
12
II I0
to
to e
g
9
8
7
7
6
5
6
S
12
II
io
IO 9
|
9
8
7
7
IS
6
5
5
g
g
e
6
6
5
I0
tO
Fig. 13. Salinity distribution at 4th day, 20th hour.
5
0
0
I
L
1 j
I
L
2
2
558
MYRON H. YOUNG CONCLUSIONS
T h e i m p l i c i t - i m p l i c i t s c h e m e is d e s i g n e d to r e d u c e c o m p u t e r C P U t i m e b y t a k i n g a l a r g e r t i m e s t e p as c o m p a r e d to t h e i m p l i c i t - e x p l i c i t s c h e m e . C o m p a r i s o n s w e r e m a d e w h e n t h e p r o g r a m w a s a p p l i e d to s e v e r a l d i f f e r e n t g e o g r a p h i c l o c a t i o n s . T i m e s a v i n g s v a r i e d f r o m f o u r - to e i g h t f o l d . T h e g r i d s y s t e m c h o s e n f o r this m o d e l will b e t t e r m o d e l c h a n n e l s a n d c a n a l s . H o w e v e r , t h e grid s i z e is still t o o l a r g e to t r u l y m o d e l a c h a n n e l . T h e n e x t l o g i c a l s t e p is to i n c o r p o r a t e v a r i a b l e grid s i z e s into t h e c o m p u t e r p r o g r a m . P r o v i s i o n f o r v a r i a b l e boundaries caused by flooding should also be made. T h e s a l i n i t y s t u d y o f this m o d e l s e e m s to g i v e a r e a s o n a b l y g o o d f i r s t - o r d e r e s t i m a t i o n . T h e s i m p l i c i t y o f u s i n g this p r o g r a m m a k e s t h e m o d e l a t t r a c t i v e to u s e r s . REFERENCES 1. Tracor, Inc., Estuarine modeling: An assessment. 16070 DZV, EPA Water Pollution Control Series, Washington, DC (1971). 2. J. D. Wang and J. J. Connor, Mathematical modeling of near coastal circulation. Report No. 200, MIT R75-19, Massachusetts Institute of Technology, Cambridge, MA (1975). 3. J. J. Leendertse, Aspects of a computational model for long-period water-wave propagation. RM-5294-pr, The Rand Corporation (1967). 4. S. Hacker, R. Pike, and B. Wilkins, Transport phenomena in estuaries. LSU SGP-RFL-73-1, Louisiana State University, Baton Rouge, LA (1973). 5. G. F. McHugh, Development of a two-dimensional hydrodynamic numerical model for use in a shallow, well-mixed estuary. LSU-T-76-008, Louisiana State University, Baton Rouge, LA (1976). 6. J. J. Leendertse and E. C. Gritton, A water-quality simulation model for well mixed estuaries and coastal seas: Volume II, Computation procedures. R-708-NYC, The Rand Corporation (1971). 7. D. R. F. Harleman and C. H. Lee, The computation of tides and currents in estuaries and canals. MIT Technical Bulletin No. 16, Massachusetts Institute of Technology, Cambridge, MA (1969).