Front Cover
 Title Page
 Table of Contents
 List of Figures
 Governing equations
 Numerical methods
 Boundary conditions
 Model verification
 Model validation
 Model application
 Conclusions and recommendation...
 Biographical sketch

Group Title: UFLCOEL-97005
Title: A numerical hydrodynamic model for inlet-river systems
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00091092/00001
 Material Information
Title: A numerical hydrodynamic model for inlet-river systems
Series Title: UFLCOEL-97005
Physical Description: x, 83 leaves : ill. ; 28 cm.
Language: English
Creator: Seidle, Peter N., 1972-
University of Florida -- Coastal and Oceanographic Engineering Dept
Publisher: Coastal & Oceanographic Engineering Dept., University of Florida
Place of Publication: Gainesville Fla
Publication Date: 1997
Subject: Inlets -- Mathematical models -- Florida   ( lcsh )
Shore protection -- Mathematical models -- Florida   ( lcsh )
Ponce de Leon Inlet (Fla.)   ( lcsh )
Genre: government publication (state, provincial, terriorial, dependent)   ( marcgt )
bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (M.S.)--University of Florida, 1997.
Bibliography: Includes bibliographical references (leaves 81-82).
Statement of Responsibility: by Peter N. Seidle.
 Record Information
Bibliographic ID: UF00091092
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 37856747

Table of Contents
    Front Cover
        Front Cover
    Title Page
        Page i
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
    List of Figures
        Page vi
        Page vii
        Page viii
        Page ix
        Page x
        Page 1
        Page 2
        Page 3
    Governing equations
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
    Numerical methods
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
    Boundary conditions
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
    Model verification
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
    Model validation
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
    Model application
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
    Conclusions and recommendations
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
    Biographical sketch
        Page 83
Full Text




Peter N. Seidle










I would like to thank my advisor, Dr. Wang, for his

guidance and time I have used in my research. I would also

like to thank Professors Dean and Thieke for serving on my

committee and reviewing this thesis.

I would like to thank the Coastal and Oceanographic

Engineering (COE) Department for providing me with a research

assistant position and funding for my studies. With the

opportunity and resources the COE Department has provided me,

I am able to fulfill my personal goals and can better

contribute to the field of coastal engineering. I would like

to thank each and every professor in the COE Department for

piqued my interest and curiosity in their areas of research.

I would like to especially thank Dr. Lihwa Lin for his

patience, help, and guidance in my studies and research. I

truly appreciate the help from Professor Sheng, Dr. Eduardo

Yasuda, and Justin Davis for helping me with numerical

modeling. Gratitude is owed to Wendy Smith for her help as

well. I also thank each and every fellow student who has

helped me with class work.

Many thanks to Becky, Helen, John, Lucy, Sandra, and

Sidney for the administrative and computer support.

And saving the best for last, I would like to thank my

parents for their love, support, and encouragement, as well as

Micheal O'Shea and Pete's Wicked.


ACKNOWLEDGMENTS ................

LIST OF FIGURES ................

.......................... ii

.......................... vi

ABSTRACT ................................................. ix

1. INTRODUCTION ................
1.1 Background ............
1.2 Objective .............


3. NUMERICAL METHODS ...........
3.1 Linearized System ....
3.2 Grid System ...........
3.3 Non-Linear Terms ......
3.3.1 Friction Term ..
3.3.2 Convective Term
3.3.3 Diffusion Term
3.4 Numerical Stability ...


..................... .... 3

................... ......10
... ..... ............ .... .10
o ........................ 15
............... ..........17
.......................... 17
.......................... 19
.......... ............... 24

4. BOUNDARY CONDITIONS ...................................27
4.1 Open Boundary Conditions ........................29
4.1.1 Tide .....................................29
4.1.2 Neumann ..................................30
4.1.3 River and Canal ..........................31
4.2 Solid Boudaries .................................35


VERIFICATION ....................................36
Wind ............................................36
Seiching ........................................38
Massachisetts Bay ...............................41

6. MODEL VALIDATION ......................................49
6.1 Field Data ......................................49
6.1.1 Bathymetry ...............................50
6.1.2 Tide .....................................51
6.1.3 Current ..................................54
6.1.4 Drogues ..................................57
6.1.5 Wind .....................................60
6.2 Spectral Analysis ...............................60
6.3 Grid ............................................62
6.4 Open Boundary Conditions ........................64
6.4.1 Tidal OBC ................................64
6.4.2 Neumann OBC ..............................65
6.4.3 Canal OBC ................................66
6.5 Results .........................................68


APPLICATION .....................................71
Single Jetty ....................................72
Pair of Jetties .................................73

8. CONCLUSIONS AND RECOMMENDATIONS .......................76
8.1 Conclusions .....................................76
8.2 Recommmendations ................................78

LIST OF REFERENCES .......................................81

BIOGRAPHICAL SKETCH ......................................83


Figure Daage

2.1 Relationship between drag coefficient (shear stress

coefficient) and wind speed. ......................9

3.1 Rectilinear cell. .................................16

3.2 NumericalCells applied to friction term. .........18

3.3 Numerical cells applied to convective term. ......20

3.4 Numerical cells applied to convective term. ......22

4.1 Example of grid showing solid and open boundaries.


4.2 Example tidal and Neumann open boundaries. .......33

4.3 Example of a grid for a river or canal open boundary

condition. ........................................35

5.1 Comparison of analytical and numerically modeled

set-up. ...........................................38

5.2 Numerically modeled seich with 1 node. ..........40

5.3 Numerically modeled seich with 2 nodes. ..........41

5.4 Massachusetts Bay. ................................42

5.5 Geometry of 2D analytical solution. .............. 44

5.6 Analytical surface elevations at high tide and

maximim ebb velocity vectors. ....................45

5.7 Maximum velocity vectors from numerically modeling

Massachusetts Bay. ................................47

5.8 Numerically modeled surface elevations at high tide

in Massachusetts Bay. ............................48

6.1 Bathymetric survey of Ponce de Leon Inlet. .......51

6.2 Location of stations where tide and current data

were collected. ...................................52

6.3 Continuously observed tidal record and discrete

current measurements. .............................53

6.4 Flood velocity data taken inside the canal. ......55

6.5 Ebb velocity data taken in the canal. ............56

6.6 Flood current patterns studied using drogues. ....58

6.7 Ebb flow patterns studied using drogues. .........59

6.8 One tidal ccyle extracted from tide data. ........61

6.9 Discrete spectrum of the tidaal cycle. ........... 62

6.10 Area of Ponce de Leon Inlet which the grid is

applied to. .......................................63

6.11 Simulated tidal OBC signal compared to the observed

tidal signal. .....................................66

6.12 Maximum flood velocity vectors for numerical model

validation. .......................................68

6.13 Maximum ebb velocity vetcors showing numerical model

validation. .......................................69

7.1 Flood velocity vecors of a single jetty applied to

Punta Gorda, FL. ..................................72

7.2 Ebb velocity vecors of a single jetty applied to

Punta Gorda, FL. ..................................73

7.3 Flood velocity vecors of a pair of jetties applied

to Punta Gorda, FL. ...............................74

7.4 Ebb velocity vecors of a pair of jetties applied to

Punta Gorda, FL. ..................................75

8.1 Bathymetric smoothing. ...........................68

8.2 An idealized channel/bay system. .................80

Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science



Peter N. Seidle

May, 1997

Chairperson: Dr. Hsiang Wang
Major Department: Coastal and Oceanographic Engineering

A numerical model is developed for predicting tidal

induced near shore hydrodynamics of an inlet-river system or

inlet canal network. The model couples a two-dimensional

model for the inlet region and a one-dimensional model

representing the river and/or canal network. The model is

based on depth averaged equations of motion including the

convective acceleration, friction, turbulent diffusion, and

wind stress terms. The numerical method used to solve the

equations is the factorization scheme. The model can be

driven using tides and river discharges. A canal system is

treated as a river with zero discharge.

The model is verified by comparing known analytical

solutions to the numerical results for selected cases. The

numerical model is then applied to Ponce de Leon Inlet which

serves a canal system for the City of Punta Gorda, FL. The

model is used to study the effects on the hydrodynamic

processes when jetties are added to Ponce de Leon Inlet.


Coastal regions are important to communities for their

leisure purposes and resources. Access to these resources may

be through coastal waterways, residential canal system,

harbor, etc. It is important to keep these resources

accessible not only for leisure purposes, but also for

commercial, fish and wildlife, and ecological reasons. There

have been numerous studies and projects to learn how keep

these resources open, accessible, and stable in all aspects.

One area of study is numerical modeling. This thesis will

develop a numerical model such that hydrodynamic studies of

inlet improvement measures such as jetties, rip-rap fill,

breakwaters, groins, etc., can be evaluated.

1.1 Background

It is important to understand the hydrodynamics of water

bodies, large and small, so that the benefits gained from the

coastal resources can grow, or at least remain stable.

Technology has given us the tools so that these hydrodynamic


studies can be conducted in a faster, more efficient manner.

Numerical models, confirmed by analytical models and augmented

by physical models, can be developed to learn more about the

micro- and macro-scale processes which govern the

hydrodynamics of coastal waterways, harbors, canals, and


As an aid to the design of inlet and waterway

improvements, the effects of the tidal circulation, a

hydrodynamic process, near the inlet improvement measures can

be studied using numerical modeling. Assuming that the tidal

flow is not highly stratified and that the inlet width is

large compared to its depth, 2D numerical modeling should be

adequate to predict the tidal flows and velocities.

As previously mentioned, researchers have conducted

innumerable studies and published many papers about

hydrodynamics of estuaries and numerically modeling tidal

areas. Arcilla (1989) developed a general algorithm for

generating numerical models. Estuaries and bays are a focus

of study in two-dimensional circulation: Cheng (1982) studied

Tampa Bay, FL for circulation; Schomaker (1983) developed a

circulation model of Monterey Bay, CA; Reidel and Wilkinson

(1976) studied Cockburn Sound, Western Australia; Reeve and

Hiley (1992) studied the southern North Sea; Bathen (1976)

developed a model to study Keehi Lagoon, HI; etc. But these


models were developed for larger basins than a navigational

inlet and are only applicable to the water body for which they

were designed, so a more generalized two-dimensional numerical

model is needed which can be applied to different types and

sizes of inlets.

1.2 Objective

The objective of this project is to develop an implicit

two-dimensional numerical model which can easily be applied to

a multitude of inlets with different geometrical

configurations. The model will also be designed to aid in the

placement of jetties, groins, rip-rap fill, or other such

inlet improvements by simulating the tidal flows. Analysis of

this tidal process will allow for the optimum placement of

inlet improvement measures to ensure maximum benefits from the


This report will (1) present the equations of continuity

and motion used; (2) derive the numerical method used; (3)

show the differencing techniques applied to the governing

equations; (3) verify the numerical models efficacy; (4)

validate the model using data obtained from Punta Gorda, FL;

and (5) show how the model is applied as an aid to inlet

improvement design.


The governing equations of motion used for this model

were derived from the Navier-Stokes equations of fluid motion

for a non-compressible, Newtonian fluid. The simplifying

approximations applied were as follows:

1)By order-of-magnitude analysis, the vertical

accelerations are negligible compared to gravity.

Hence, a hydrostatic approximation is applied.

2)Since the Rossby number, the ratio of Coriolis

force to inertial forces, is small for inlets this

model is intended for, the Coriolis force is


3)Assume the turbulent Reynolds stresses are

proportional to velocity gradients by an eddy

viscosity coefficient and that the viscous

stresses are negligible compared to the turbulent


stresses. Also, for simplicity, assume the eddy

viscosity coefficient is constant for the inlet.

4)Assume atmospheric pressure is constant and

uniform for the region being modeled.

5)Assuming the inlet is well mixed and that there

are no density driven currents, the baroclinic

pressure gradient can be neglected.

6)Boussinesq approximation is applied to the

barotropic pressure gradient. This assumes that

the density gradient is small enough such that the

density is assumed to be uniform.

The continuity and resulting equations of motion are as


a +---uD+vD-= (2.1)
at ax ay

au au au a a2u a2u nd (2
-+u-a+v-=-g-- (-+-)+ (2.2)
at ax ay ax 8D ax2 ay2 pD

9v 9v av arn f 2 ua2v cv)y-wnd "
-+u- +v- = -g-i f v+A(-+- )+- (2.3)
at dx ay ay 8D x2 ay pD



u=depth averaged velocity in the x-direction,

v=depth averaged velocity in the y-direction,

u=depth averaged velocity vector,

D=total depth,

ri=water surface elevation,


x,y=right handed Cartesian coordinates,

f=Darcy-Weisbach friction coefficient,

A,=horizontal diffusion coefficient,

Tx=x-direction wind shear,

ty=y-direction wind shear, and

p=density of sea water.

There are two empirical coefficients in the above

equations, the Darcy-Weisbach friction coefficient and the


eddy viscosity coefficient. In numerical modeling they are

often adjusted to tune the model based upon field data. There

are numerous studies concerning the friction coefficient for

open channel flows. In engineering applications, the most

common methods of developing a friction factor use the Manning

coefficient or Chezy coefficient. The friction coefficient is

a function of many types of environmental conditions, such as

type and amount of vegetation on both the bottom and banks,

and some of these conditions have been quantified in terms of

the Manning coefficient (Henderson, 1966). For a sandy bottom

along an open coast, a Darcy-Weisbach friction coefficient of

approximately 0.01 is commonly used.

The value, as well as the range of variability, of the

eddy viscosity coefficient is far less well established. This

is because eddy viscosity is a flow property which is related

to turbulent mixing. Since flow turbulence is a complicated

phenomenon and the scale of mixing varies greatly from case to

case, it is not a simple matter to model its mixing effects

using only one representative coefficient. In an oceanic

environment the eddy viscosity coefficient is known to vary

enormously. The order of magnitude for the coefficient

suggested by Pedolsky (1979) ranges from 10 to 10000. For

each scenario numerically modeled, a stability condition

(discussed in chapter 3.4) must be considered in selecting the


eddy viscosity coefficient. The stability condition relates

the grid size of the numerical model and time step of the

numerical method to the coefficient.

The wind shear, tnd is somewhat empirical as well and

may be necessary to validate the numerical model. If field

studies were conducted on a windy day, then the wind must be

accounted for in order to validate the model. Once the

numerical model has been validated, the wind speed and

direction may be changed or set to zero depending on the long

term predominant conditions.

The wind shear forcing term can be written as:

x-wind P. C
air'd- 2 WsinO (2.4)
pD pD

y-Wfod P. dW cosO (2.5)
pD pD


Cr=drag coefficient,

W1o=wind speed at 10m, and

O=wind bearing.

Although the drag coefficient is empirical, researchers have



~ 0.002


fooi A

O .


..5 0 15 20 25 .. 30

0 10

20 30

Windsioeed 17.o

40 50

Figure 2.1
Relationship between drag coefficient (shear
stress coefficient) and wind speed. Figure was
taken from Roll (1965).

found that it is on the order of 10-3. The drag coefficient

could be found graphically from figure 2.1. Also, Wilson

(1960) concluded that for light winds (typically less that 10

m/sec), Cd=1.49 x 10-3 and for stronger winds, Cd=2.37 x 10-3.

For validating this numerical model, Wilson's findings were


The governing equations (2.1, 2.1, and 2.3) are solved

numerically for the problems posed in this thesis. The

numerical scheme and the applicable boundary conditions are

discussed in chapters 3 and 4.

-, ,. .. ,.

I _


.. .


The numerical method adopted for the model is the so-

called factorizationn" scheme. This method of implicitly

solving a set of linearized equations was developed to study

tidal currents in the vicinity of the Mississippi Sound in the

Gulf of Mexico (Sheng, 1983).

3.1 Linearized System

The governing equations (2.1, 2.2, and 2.3) are

linearized in order to develop the factorizationn" numerical

method. The governing equations now become:

By au av
L+D- +D-=O (3.1)
at ax ay

u--g (3.2)
at ax

--g (3.3)
at ay

The above equations 3.1, 3.2, and 3.3 can be combined as

a system of equations by letting w=[r, u, v]T The system of

equations can be written as:

Iwt+Aw,+Bw =O (3.4)

where I is an identity matrix and A and B are coefficient

matrices as follows:

0 D 0
A=g 0 0

0 0 D
B= 0 0 0
g0 0

The subscripts t, x, and y indicate the first partial

derivative respective to the subscript. Applying a forward

finite difference approximation to the time derivative yields:

n+1 n
w n+ w +! n 1 nI I
At +Aw. +BWy =0 (3.5)
At x y

Multiplying by At and factoring out w"'l results in:

w"n I(+AtAx+AtB6 )=w (3.6)

where 6 is the finite difference respective to the subscript

and the superscripts indicate the time step. The left hand

side (LHS) of the equation can be expanded to:

w"'(l +AtA6x+AtB8 +(At)2ABx y)-02=W (3.7)

where 02 denotes a 2nd order error term. This term will be

resolved and removed in later steps of this derivation. A

term of AtB6 is added and subtracted to the right hand side

(RHS) of the equation with no net change. The LHS of Equation

3.7 is then factored into two terms to give:

wn"(1 +AtA6x)(1+AtB6)-02=(1 -AtB8)w "+AtB6w (3.8)

An intermediate time step may be introduced, which is

also referred to as time splitting. Letting w* represent the
intermediate time step, n+ w*=-A w" and equation 3.8
1 +AtA86
can be rewritten as:

w"n'(1 +AtA86)(1 +AtB8y)- 2=(1 +AtA6)w *+AtB8yw "

Moving the error term to the RHS of the equation and re-

expanding it the above equation then becomes:

w"+(1 +AtAWx)( +AtB6y)=(l +tA8x)w +AtB86w "+(At)2ABx yw "


Factoring out AtB6y from the two last terms on the RHS yields:

wn'(1 +AtA6 )(1 +AtB8y)=(1 +AtA8)w +AtB (1 +AtA8x)w"


Lastly, removing the common term, 1+AtA6x from the equation

results in:

(1 +AtB86)wn+l = +AtB8yW n


In summary, the time splitting equations used to solved

the linearized system of governing equations (3.1, 3.2, and

3.3) are written as follows:



(1 +AtAx)w *=(1 -AtBS y)

(1 +AtB y)w"n+=w +AtB8 "

Expanding equation 3.13 gives the following:

Tq +AtD8xu* ="-AtD6v "



u +Atg6xr*=" (3.16)

v*=v "-A tgyl" (3.17)

and expanding equation 3.14 gives:

rl' n+AtD8 vn"'+=i*+AtD8 v" (3.18)

unI=u* (3.19)

v "ln+Atg6 n+1 =v* +Atg6T"n (3.20)

Note that equations 3.15 and 3.16 are a system of

simultaneous equations where *' and u* are the two

unknowns, as well as equations 3.18 and 3.20,

where 1"n' and v"nI are the unknowns. Since equation 3.17 is

independent of q* and u* it can be substituted directly

into equation 3.20. Likewise, equation 3.19 states that

u" =u* so equations 3.15, 3.16, 3.18, and 3.20 can be

adjusted accordingly. The four resulting equations are as


S* +AtD6xu ''n n-AtD6v (3.21)
x f"fDyn

Atgx6* +un"+ =u"


11"l+AtD8 v"nl=*+AtD8 v n (3.23)

Atg68 n+n+n=V" (3.24)

The above factorizationn" method is similar to the

Alternate Direction Implicit (ADI) method in that equations

3.21 and 3.25 are solved in the x-direction (x-sweep) and that

equations 3.23 and 3.24 are solved in the y-direction (y-

sweep). It is also important to note that equations 3.22 and

3.24 are essentially the momentum equations in the x- and y-

directions, respectively. This will allow the non-linear

terms to be added. Moreover, when equation 3.21 is

substituted into equation 3.23 for i* ,the original continuity

equation (Equation 3.1) is recovered.

3.2 Grid System

A rectilinear staggered grid system is used for this

numerical model, shown in Figure 3.1. In applying finite

differences to equations 3.21 through 3.24 the staggered grid

system is examined more closely. A forward difference is



Figure 3.1
Rectilinear cell

applied to the velocities for the continuity equations (2.21

and 3.23) and may be written as:



*A / i n+ .
Atg( -1_)+uJ =Uij

A backward difference is applied to the surface elevation for

the momentum equations (3.22 and 3.24) and are as follows:

n+il n++ n+l-t D vA,,- n n
4j +alvtvi,,l-v )=I Vi+j)(viji vi.

A -g(, nl -n+1r_ n+1 n
4gr i -r l ),i + =v



3.3 Non-Linear Terms

The non-linear terms in the governing x- and y- momentum

equations (2.2 and 2.3) which have been neglected thus far can

now be added to the numerical x- and y- momentum equations

(3.26 and 3.28), respectively. It is necessary to examine the

adjacent elements of the grid shown in figure 3.1 in order to

properly examine each non-linear term of the momentum

equations individually.

3.3.1 Friction Term:

The frictional forcing term in the x-direction momentum

equation is f-uu This term is non-linear in that it is
quadratic in terms of flow velocity vector. Here the velocity

vector term, lul is expressed numerically in terms of the

known velocity components in the previous time step, such as u"

and v" Furthermore, since the x-direction momentum
equation, eq. 3.26, is solving for u,, at the grid center

depicted in Figure 3.2, lul must be calculated at the proper

location. The numerical equation used to calculate the

Figure 3.2
Numerical cells applied to friction

magnitude of the velocity vector for the x-direction momentum

equation is as follows:

I|I= .ci2 +i 4 ) (3.29)

The numerical expression of the frictional forcing term which

is added to the LHS of x-direction momentum equation, eq.

3.26, is given as follows:

n n n n
f ( i-+(--- ij+ i-Iij )2 n+l
(ui+2 +
8D1 4

Similarly, the frictional forcing term added to the LHS

of the numerical y-direction momentum equation, eq. 3.28, is:

n n n n
S+ + +Ui+- 2 n+1
(( -) + vi
8D1 4

3.3.2 Convective Terms:

The convective terms have directional characteristics

which must be appropriately accounted for in the numerical

scheme. That is, the down stream behavior is affected by the

upstream convection phenomenon (Hirsch, 1994). A preferred

differencing technique is the upwind method, which is similar

to a backward difference written with respect to the direction

of the velocity vector. Because of this directional

phenomenon, the direction of the velocity must be

predetermined in the numerical expression of these convective


For the x-direction momentum equation, the convective
au au
acceleration terms are u-+v- taken from equation 2.2. In
ax ay
order to be included with the numerical equation of motion

(equation 3.26), finite differencing must be applied to these

convective acceleration terms such that they are congruous

Figure 3.3
Numerical cells applied to convective

with the upstream convection phenomenon. Both convective
au au
terms, u- and v- need to be examined separately.
ax ay

To analyze the term u- first the velocity direction
needs to be determined. An average of the velocities between

the two cells (figure 3.3), u_-1, and u,+12, are used to

identify the direction of the velocity so that the proper
upwind differencing method may be applied. If uin2j is

positive, then the upwind difference is u "-u"_ The

numerical equation used to evaluate the first convective term


S ui ui+

n n
au n uijUi-j~ (3.30)
x Ax

However, if ui+2, is negative, then the upwind deference

is u-i-ul and the equivalent numerical difference is as


n n
u -" ui1 -uii (3.31)
ax /2Ax

The second convective term, v- can now be examined
with the grid shown in figure 3.4 using the same steps in

analyzing the first convective term. If vi"n/2j is in the

positive direction, then the numerical convection term is as


n n
v vn _i ij-1 (3.32)
ay v'-1/ Ay

And, if vil"/2,+1 is negative, the numerical differencing is:

n n
-Vu- n i ij+l (3.33)
-V = -V i-1/2j +1 1
ay .Ay

At every time step, the direction of the velocity is

determined for each element so that either equations 3.30 or




Figure 3.4
Numerical cells applied to
the convective term.

3.31 and either equations 3.32 and 3.33 are added to the

numerical x-direction momentum equation, eq. 3.26.

The same process used to derive the numerical equations
au au
for u-+v- can be used to find the numerical equations for
ax ay
the convection terms in the y-direction momentum equation,
av av
u-+v- For a positive velocity direction, the two
ax ay
convective terms are as follows:

n n
Uav V U (3.34)
ax i-1/2J Ax


n n
a yn Vij v-1
V ay Vij-/2y


Likewise, if the velocity is negative then the terms are as:

n n
av n l ij
-u-a u
ax12 Ax/2h



n n
SV n Vij -Vij+l
- y- i+1/2
ay i+I Ay


3.3.3 Diffusion Term:

The diffusion term is linear provided that the diffusion

coefficient is not a function of position, time, velocity, or

water surface fluctuation. The numerical expression for the

x-direction diffusive term:

a2u 82u
ax2 ay2

is obtained by simply applying a center difference to yield:

Aun i-2u u' -2u+,
n n n n n n
A i(.lj ij-z --lj +u Ij+ 4 Iu-1)
Ax2 Ay2

which is included on the RHS of the numerical x-direction

momentum equation, eq. 3.26. The terms added to the RHS of

the y-direction momentum equation, eq. 3.28 are:

A nI-2V n n-2v" n
Ax2 Ay2

3.4 Numerical Stability

Numerical stability criteria for the non-linear system

presented here are not established at present. However, the

stability criteria of a linear system are used instead. There

are two criteria which need to be met for numerical stability:

one related to tide propagation and another to diffusion.

For propagation stability, the well-known Courant

criterion is invoked. The Courant Number used in the present

model is defined as:

Cr= t ( +lul+vl) (3.38)


At=time step,

Ax=x-direction grid spacing,

Ay=y-direction grid spacing,

h=maximum depth,


lul=maximum x-direction velocity superimposed on

the tidal velocities, and

Ivi=maximum y-direction velocity superimposed on the

tidal velocities.

For a chosen grid size the time step, At, is then constrained

by the above equation for a critical Courant number, Cr.

Generally, for implicit numerical techniques the critical

Courant number can be greater than unity, whereas for explicit

methods the Courant number must be less than unity to maintain

numerical stability. For the (ADI) method, the Courant number

is suggested to be less than 5 to flows that are not

consistent with observed nature (Weiyan, 1992). The advantage

of using the factorization method is that stability is

maintained for Courant numbers of about 10 (Sheng, 1983).

For diffusive stability the eddy viscosity coefficient is

examined. The following condition has been suggested by

Hirsch (1994):

A < 1 (3.39)
At 2

This criterion is used in the present model to guide the

selection of the diffusion coefficient depending on the

spatial grid sizes and the time step size.

In summary, the grid size, Ax and Ay, are first selected.

Next, the time step, At, and the eddy viscosity coefficient

are then selected using the Courant number relationship, eq.

3.38, and the aforementioned diffusion criterion by Hirsch

(1994), eq. 3.39.


This chapter addresses "geometrical" boundary conditions

(BCs) used to define, drive, and control the numerical model.

Since the rectilinear grid is used, the numerical model must

be able to differentiate land and water cells. The land and

water cells are what define the model's geometry, so that an

inlet, lagoon, or harbor can be numerically simulated using a

rectilinear grid.

The bottom and surface BCs are addressed by the governing

equations as friction and wind stress forcing terms,

respectively, which were defined in chapter 2. The following

sections in this chapter will address the details of applying

the land and water boundary conditions.

The "geometric" boundary conditions (BCs) used in the

numerical model can be broken into two types, open and solid

An open boundary (OB) is a cell which is located on the edge

of the grid that covers the modeling area. The input to the

model, such as tide, river discharge and length, or canal

length, are applied to the OB cells. Inside the modeling



1,4 .6,4

2,3 .6,3

2,2 .6,2



i=l i=2 i=I

Figure 4.1
Example of grid showing solid and open

area, some cells are designated either as land or water to

define the model's geometry. The land/water interface creates

a solid boundary. Figure 4.1 shows an arbitrary example.

Cells (6,2), (6,3), and (6,4) are defined as OBs, so input

must be applied to these cells to drive the model. Whereas


cells (1,4), (1,5), (2,2), and (1,5) are defined as land or

solid boundary cells.

4.1 Open Boundary Conditions

Four types of open boundary conditions (OBCs) are

incorporated in the model; tidal, Neumann, canal, and river.

Though a Neumann OBC does not drive the numerical model as a

tidal OBC would, it is used to control the hydrodynamic

stability of the numerical model in certain tidal OBCs.

There are some considerations to be made when selecting

the model's grid and designating the OB cells as either a

tide, river, or canal BC. In practice, one edge of a cell may

be designated as a tidal OB and a different edge may be

defined as a canal or river OB. These methods of application

will be discussed in the following sections.

4.1.1 Tide

The tidal OB is defined by water level changes such as a

tidal fluctuation of the form ==A cos(at+() where A is the

amplitude, a is the period, and q is the phase angle. Though

it is not necessary for all scenarios, the numerical model can

accept a real tide as input using up to four Fourier series


harmonics. The Fourier series representing a mixed tide would

be as follows:

l= C (A,cos(ot+(,)) (4.1)

When applying the tidal OB, the velocities remain as unknowns,

otherwise, the model would be over-constrained and could not

be validated. The edge of the tidal OB cells should be of

sufficient distance offshore from the region of interest.

4.1.2 Neumann

If the area that is going to be modeled requires placing

tide OB cells on the corners, then the a Neumann BC must be

incorporated. Figure 4.2 shows the grid which would be

applied to such a situation. Where the cells marked with a T

indicate where a tidal forcing condition is imposed. The

cells marked with an N are where a Neumann BC is imposed. The

Neumann BC requires that the surface elevation gradient equal

zero between the OB cells and the cells in the J-1 row or the

2nd row of cells.


J-1 T




2 T

1 N N N N T

Figure 4.2
Example of tidal and Neumann open boundary

4.1.3 River and Canal

In order to apply the river or canal OB, the velocities

must be perpendicular to the edge of the grid. If there is

circulation or lateral currents, with respect to the river or

canal, then the model area must be extended out to where the

velocity vectors become perpendicular to the OB. Also, in


applying a river OB, the head of the river must not be

influenced by the tide. The discharge at the head of the

river is used to drive the model and should be constant,

though the model can be easily modified to allow for a time-

dependent river discharge. The model becomes unstable or may

produce flows that are not indicative of nature if the tide

propagates to the head of the river. A canal is treated as a

river, except that the discharge at the head is zero.

Moreover, since the discharge at the head of a canal is zero,

the model will be stable and valid if the tide propagates to

the head, or end, of the canal.

The length, depth, and discharge at the head, as well as

grid spacing and friction coefficient, are needed as input at

the river/canal OB (figure 4.3). In addition, the user needs

to define a grid, or cell, length for the river or canal

boundaries. The cell can be significantly larger compared to

that of cells in the model area. It is recommended that the

river or canal be divided into 15 to 30 cells. If any more

cells are used, then the computing time required to run the

model is compromised. The governing equations used by the

numerical code in the river or canal cells are as follows:

S+-uD=0 (4.2)
at ax

at head

N' Cell

River or Canal OB

Cell Size

Figure 4.3
Example of a grid for a river or canal open boundary

du+ u = l- f
at ax ax 8D



u=depth averaged velocity in the direction of the

river or canal,


D=total depth,

l=water surface elevation,


x=direction of the river or canal, and

f=Darcy-Weisbach friction coefficient.

In the previous chapter (chapter 3.3), statements

regarding the governing equations dealing with the friction

coefficient and convective term also apply to the above

equations. The model uses the surface elevation from the

previous time step and the river/canal discharge as boundary

conditions to calculate the new velocity at the model's OB

cells. The computed velocities at the OB cells are then used

as a boundary condition to drive the model. The new surface

elevations calculated for the cells at the OB are then used as

input into the river/canal OB subroutine for the next time

step. Since the flow at the OB is assumed to be perpendicular

to the boundary, the transverse velocity at the OB is set to

zero. For example, since the velocity vector at the OB in

figure 4.3 should be in the x-direction, the v-velocity is set

to zero.

At all OB cells, the model becomes unstable if diffusive

or convective terms of the momentum equations are included.

This prevents eddies from occurring at the boundaries which


would cause instability of subsequent computations. If field

data shows that there is circulation at the edges of the

numerical model's grid, the grid area should be enlarged to an

area where there is no circulation found by the field studies.

4.2 Solid Boundaries

Finally, at land boundaries, the flow condition is

satisfied by setting the normal velocity component to zero.

There is also a non-slip condition being imposed so that the

velocity parallel to the boundary is set to zero as well.


The numerical code must be compared to analytical cases

to verify the efficacy. The model's efficacy was tested using

wind set-up and seiching in a closed basin, and by replicating

tidal currents in Massachusetts Bay. Once the model is

verified, it can be confidently applied to various inlets.

5.1 Wind

The analytical set-up due to wind stress is (Choi, 1992):

(x)= wnd(x- ) (5.1)
pgD 2


uwnd =Shear stress by wind,

p=density of water,

x=distance from middle of basin,


D=depth, and

L=length of basin.

A closed, flat, rectangular basin was used to develop the

analytical solution. The length of the test basin which the

wind was blowing across, was 450 m long, with a grid size of

50 m. The depth of the basin was constant at 2 M and the wind

speed was constant at 5 m/s. The wind shearing force,

calculated using the respective drag coefficient discussed

previously in the governing equations in chapter 2, equals

0.0448 Pa. Figure 5.1 shows the surface elevation of a cross

section from the numerical model compared to the analytical

wind set-up. Only one cross section is presented since, as

expected, every cross section has the same wind set-up


The reason there is a difference between the numerical

and analytical set-ups is that the model's output is rounded

to 0.1 mm, where the analytical line does not include any

round off. For example, at x=25 m, the numerical set-up is

-0.0004 m and the analytical set-up is -0.0004386 m. If the

numerical model's output included more significant digits for

this test, then a closer comparison could be made. Moreover,

if the analytical solution was rounded off to the 1/10 mm,

then the two lines would correspond exactly.

x10- Wind Set-up



-1 -- Analytical
1 --- Numerical

-2 -


-4 *

0 50 100 150 200 250 300 350 400 450
Distance, (m)

Figure 5.1
Comparison of analytical and numerically modeled

5.2 Seichina

The model was also tested for its ability to replicate a

standing wave in a closed basin. The basin length required

for seiching is:

1Ln (5.2)


l=basin length,

n=number of nodes, and

L=wavelength of long wave.

For a channel which is 2 m deep and has a tide period, T, of

12 hrs, the wavelength is 191,352 m. The basin lengths which

would form a seiche with 1 and 2 nodes are 95,676 m and

191,352 m, respectively.

The rectangular test basin used for both cases is 3 m

deep and 450 m wide. Also, for both cases (1 and 2 nodes),

the number of cells used for the width and length were 9 and

50, respectively. Which means that Ax for the two different

cases are 1913.52 m and 3827.04 m, and Ay was 50 m for both

cases as well. The Courant number for both cases was 2. The

model was driven by specifying te water surface elevation at

the right hand side of the basin as a function of

time, i=Acos 21 where A, the amplitude, was set to 0.1 m

and T, period, was set to 12 hrs.

Figures 5.2 and 5.3 present the surface elevation for a

cross section of the basin for 1 and 2 nodal standing wave at

times 1/4T, 1/2T, 3/4T, and T. Note that the numerical

results are reasonably good, however, they did exceed the

0 10 20 30 40 50 60 70 80 90 100
0.2 -- i -- i -- i -- i -- i -- i -- i --- i --- i -

0 10 20 30 40 50 60 70 80 90 100

0 10 20 30 40 50 60 70 80 90 100

00.1 .. ..-

0 10 20 30 40 50 60 70 80 90 100
Distance, (km)

Figure 5.2
Numerically modeled seiche with 1 node at time,
from top to bottom; (a) =1/4T, (b) = 1/2T, (c)
=3/4T, and (d) =T.

input amplitude because friction was set to 0.0001 in the

numerical model. The small friction coefficient did not allow

enough energy to be damped from the system. Though it is not

evident for times 1/4T and 3/4T, there were surface

fluctuations which were less than 0.05 mm. The initial

surface conditions were evaluated using a cosine function to

represent a seiche, so the model quickly reached equilibrium

and became accurate.

I I 0 I
10 20 30 40 50

60 70 80 90 100

0 10 20 30 40 50 60 70 80 90 100

-0.2 LI
0 10 20 30 40 50 60 70 80 90 100
Distance, (km)

Figure 5.3
Numerically modeled seiche with 2 nodes at time,
from top to bottom; (a) = 1/4T, (b) = 1/2T, (c) =
3/4T, and (d) = T.

5.3 Massachusetts Bay

The final verification of the model's efficacy was to

reproduce the results obtained from an analytical solution of

the 2D equations of continuity and momentum with geometry

representative of Massachusetts Bay (shown in figure 5.4).

The analytical solution for tide, derived by Briggs and Madsen

0.2 -


I I i I

:. .... t . .AL


Cape Cod

0.0 :.5, .. .10.0 .

Scae (Nautical milss)

V . . .. ... .. .

. 7Ioston
SLigh, h p

. : :

(+0.2) .


:.M ::. .I c

S I "..... . .

n .x .: ...", yi.. :. : .j

: .. .
: I, ..

-. r. -
Figure 5.4
Massachusetts Bay. Figure taken from Briggs and Madsen

E.: ;:


C. C
..t.r 4


W" 6 )D,:: ,


/ / / 7 / / / /


(1973), to the linearized continuity and momentum equations

(eqs. 3.1, 3.2, and 3.3) is:

i-=Aoe t[osmY _-E 2mosinmyo(sinkx2 sinkjx) coshmyosk] (5.3)
n=1 mk, (x2 -x) sinhmyo



o2 2
M = ---k.
n gh "

and Xo, x,, x2, and y, are the geometric dimensions of the bay

as defined in figure 5.5. For simplicity, the depth, h, was

kept constant. The tide at a single location anywhere in the

region must be known and is used to solve for Ao. Briggs and

Madsen used the tide at the Boston Light House as the driving

force in their analysis. However, Boston Light House is

inside the bay, not at an open boundary, so the tide

measurement there could not be used to drive this numerical

model. From examining the tide (figure 5.6) at the open

boundary, found analytically by Briggs and Madsen, a value of



Figure 5.5
Geometry of 2D analytical solution. Figure taken from
Briggs and Madsen (1973).

4.30 ft can be used to drive the numerical model. The depth

that Briggs and Madsen used was 120 ft. They also used a

sinusoidal tide as input at Boston Light House, so it will be

assumed that the tide at the entrance to the bay will be a

sinusoidal tide as well. Figure 5.6 also shows the

analytically derived maximum ebb currents for Massachusetts


U I-

o e
o 0


o "

analytical surface elevations at high
and Madsen (1973)


a In0



and Madsen (1973).

For the numerical model, the bay was divided into 30

increments in the x-direction and 10 increments in the y-

direction. The Courant number used for the model was 2.2, the

Darcy-Weisbach friction coefficient was 0.05, and the

diffusion coefficient was set to zero since the analytical

solutions neglected the diffusive terms. Note that the

analytical solutions do not include the friction or convective

terms. Briggs and Madsen used order-of-magnitude comparison

to neglect the friction term so that the governing equations

could be solved easier. Friction had to be include in the

numerical model so that energy could be dissipated and the

model would remain stable. A high friction coefficient was

used because of the deepness of the Massachusetts Bay. This

caused a discrepancy comparing the numerical results to the

analytical solutions.

Figure 5.7 shows the maximum ebb velocity vectors

developed by the numerical model. By visual comparison the

numerical model is in good agreement with the analytical

results. Figure 5.8 shows the numerically modeled surface

elevations at high tide in the bay. Again, the discrepancies

could be due to the friction and convective terms which the

numerical model includes, and the analytical solution



Velocity Vectors of Massachusetts Bay

0 10 20 30 40 50 6(
Nautical Miles

Figure 5.7
Maximum velocity vectors from numerically
modeling Massachusetts Bay.

, ttt it1 \\\.... .....
. l4 %t l \ \ \ \ \ \\ ... .
I t t I\ .- . .

- - - - - .. .
l l l ~\1 ...sec

* .. ...... =1 ftlsec

Surface Elavation (ft) of Massachusetts Bay

0 10 20 30 40 50
Nautical Miles

Figure 5.8
Nummerically modeled surface elevations at high
tide in Massachusetts Bay.


Since the model's efficacy was verified, the numerical

model can now be applied to various inlets, but must be

validated for each inlet. Validation can be an iterative

process; the model's coefficients (friction, diffusion, and

wind drag) need to be tuned until the numerical results

simulate the actual hydrodynamic process (i.e. circulation and

eddies) and the currents of the inlet being modeled. Once the

model is validated, any inlet improvement measure, such as

jetties, groins, or rip-rap fill, can be confidently


6.1 Field Data

The Coastal and Oceanographic Engineering Department of

the University of Florida was commissioned by the City of

Punta Gorda, FL to examine the potential for shoreline erosion

due to the navigational channel at Ponce de Leon Inlet in the

vicinity of Punta Gorda City Park (Wang, et al., 1996). Ponce

de Leon Inlet is located on the east side of Charlotte Harbor


at the mouth of the Peace River. The channel serves as the

public boat ramp at Punta Gorda City Park as well as a

residential canal system. The data which was collected and

presented in a report to the City of Punta Gorda will be used

to show the validation and application (chapter 7) of this

numerical model.

Data were collected from Ponce de Leon Inlet in November

and December 1995 and March 1996. The data necessary to

validate the numerical model are: a bathymetric survey; a tide

record; current measurements; tidal flow pattern studies; and

wind speed and direction measurements. However, other field

data may be necessary for coastal structure design purposes.

6.1.1 Bathymetry

The bathymetric data collected should be complete enough

such that bottom contours can be drawn accurately. If

available, navigational charts may be used in lieu of, or to

supplement the field data. No navigational charts were

available for this inlet so a bathymetric field survey was

conducted. The bottom contours obtained from this data are

shown in figure 6.1. The depth for each cell in the numerical

grid can be interpolated from these contours.

1 -.51 0) Contour in m (NGVD)

Figure 6.1
Bathymetric survey of Ponce de Leon Inlet, figure
taken from Wang, et al. (1996).

6.1.2 Tide

Tide gages were placed at channel markers 8 and 18

(figure 6.2) to measure tide height for the period of data

collection, these data are shown in figure 6.3. The time lag

of the tide between channel markers 8 and 18 is on the order

of 5 to 10 minutes. The maximum difference in amplitude

between the two tidal records is approximately 1.5 cm, though

most of the tidal peaks show no difference in amplitude.

Figure 6.2
Location of stations where tide and current data
were collected, figure taken from Wang, et al.




A '

SOddmeter data
0 I I I I

03 2

Figure 6.3
Continuously observed- Tide at Channel Marker 8
- Tide at Channel Marker 18
-0.4current measurements, Current at inlet entrance
al. Oddmeter data

December 1995

Figure 6.3
Continuously observed tidal record and discrete
current measurements, figure taken from Wang, et
al. (1996) .

6.1.3 Current

Two instruments were used to collect current data: an

electro-magnetic meter, and a vane-type meter. The electro-

magnetic meter manufactured by Marsh-McBurney, was installed

on a docking pile just inside the inlet (station 1 in figure

6.2). The data collected from this meter are also given in

figure 6.3. Comparing the current data to the tide data, it

may be concluded that the maximum currents occur when the tide

is near zero. It is also easily seen that there is slack

water when the tide changes either from flood to ebb or ebb to


Other current data were collected from an Ottmeter, a

hand-held device which uses vanes driven by the flow to

measure the current. One minute averages were taken at about

0.5 m below the surface and 0.5 m above the bottom. Figure

6.4 shows the data collected during the flood tide and figure

6.5 shows the ebb velocity data measured at station 1, 2, and

3 (see figure 6.2) in the inlet.

Figure 6.3 shows the Ottmeter data taken at the same

location as the electro-magnetic meter. The flood velocities

between the two different methods are in agreement, though

there is a large discrepancy at the ebb tide. This

discrepancy can be attributed to the improper installation of

miT e=16:43



0.s 9 0.28
-0. s cu1ren)95
Flood current(m/s)

Station 2


e 0.51

0 0.54

0039 4.36

Station 3
M 0.74 0 0.50 0.26
0- 0.54

Time=1 6:55


Horizontal Distance (m)
Figure 6.4
Flood velocity data taken inside the canal using
the Ottmeter, figure taken from Wang, et al.

Station I








S t I... ...... .. .

a 0.25

8 0.32

Time=1 4:07
0 0.31

^ 0022 a0.26 -'
-S 0 -0.27 -- 1/30/95
Ebb current(m/s)

Station 2 Time=14:34

0 0.27




0 0.30 a 0.30

o 0.23 1-

Time=1 5:02

Horizontal Distance (m)

Figure 6.5
Ebb velocity data taken in the canal using the
Ottmeter, figure taken from Wang, et al. (1996).

tantinn I


-2 -

Station 3




the electro-magnetic current meter. The meter should have

been installed in deeper water, instead the meter was

installed such that it was less than 0.2 m of water during low

tide therefore giving the low, inaccurate measurements. Also,

the current meter was about 0.75 m from the pile, there could

have been some adverse affects on the current caused by the


6.1.4 Drogues

The ebb and flood current patterns were studied using

drogues. As the drogues were placed in the water, their

motions were video taped. From a video tape the velocities

were estimated from using estimated distances and time elapsed

on the video tape. When extracting the distances the drogues

traveled from relation to channel markers of a known distance

a problem is introduced with the perspective of the video

camera, because the video camera was at fixed, stationary

point on land. This is a crude method to obtain the current

velocities; therefore, the values may vary greatly from those

obtained through more non-Lagrangian methods. As a result,

the quantitative velocity values should only be used as a

relative reference between individual drogue path lines. The

main purpose of the drogues is to qualitatively study the flow

0 50m

Ponce De Leo

n Inlet



Figure 6.6
Flood current patterns estimated using drogues,
Wang, et al. (1996).


Figure 6.7
Ebb flow patterns studied using drogues, Wang, et
al. (1996).


patterns. Figures 6.6 and 6.7 show the flow patterns from the

drogue study for the flood and ebb tides. This is one of the

more important field studies conducted for validation of the

numerical model.

6.1.5 Wind

The days that the tide and current measurements were

taken, the atmospheric conditions were ideal, there was a high

pressure system, a slight breeze, and no synoptic events had

occurred in the previous days or were predicted. Hence, wind

effects were negligible and not included for validating the

numerical model.

6.2 Spectral Analysis

Spectral analysis is performed on the real tide, so that

the predominate frequencies and amplitudes can be used to

simulate the real tide using Fourier series (eq. 4.1). The

tide record observed at Punta Gorda (figure 6.3) shows that

the tide is mixed diurnal and semi-diurnal. Moreover, if the

record was taken for 30 days, a spring and neap tide should be

evident. For the model validation, spectral analysis will be







Observed Tidal Cycle


Figure 6.8
One tidal cycle extracted from the tide data.

performed on only one period of the recorded tide in figure

6.8. The rational for selecting only one tidal period is:

(1) excluding the scattered current measurements at

the beginning of he current record (figure

6.3), the current data reflects only one tidal

cycle; and

(2) since the model is being used as only an aid to

inlet morphological improvement measures,


residual currents and tidal flushing are not

of interest.

The continuous tide data taken were digitized to 1 hour

increments for a cycle during the time that current

measurements were being recorded. Twenty four data points

were taken from the observed data record. Figure 6.9 shows

the resulting spectrum of the tidal cycle.

6.3 Grid

The grid chosen for the Ponce de Leon Inlet as shown in

figure 6.10. The grid size Ax and Ay, is 6 m x 6 m, which

results in 59 rows and 29 columns. This resolution will be

sufficient to allow jetties to be inserted and show the

currents close to the shore. A higher resolution grid could

be used, but the computation time required will increase. The

information gained by using a higher resolution grid does not

justify the excess time needed to run the numerical model.

Both the northern and southern boundaries were extended,

so that the eddies which were studied using the drogues

(figure 6.7) do not cause a numerical instability at the open

boundaries. The drogue studies also show that the tidal flow





0 50 100 150

Figure 6.10
Area of Ponce de Leon Inlet
which a grid will be applied

becomes one-dimensional in the canal; henceforth, the grid

chosen should allow the model to be validated against the data

collected in the field studies.

Using the stability criteria discussed in chapter 3.4,

the time step and eddy viscosity coefficient can be

determined. The Courant number used was 10, which resulted in


a time step of 18 sec and an eddy viscosity coefficient of

0.12. The Darcy-Weisbach friction coefficient was set to

0.01, and, again, wind was not used as a parameter for this

model, since field studies were conducted on a calm day.

6.4 Open Boundary Conditions

The open boundary conditions (OBCs) needed to drive the

model for the case of Punta Gorda are tidal, Neumann, and

canal open boundaries (OBs). The tidal OBs will be placed on

the western edge of the grid, the Neumann OBs will be placed

on the northern and southern edges of the grid which extend

into Charlotte Harbor, and canal OBs will be placed on the

cells on the eastern edge of the grid. The following sections

will discuss the parameters used for these OBCs.

6.4.1 Tidal OBC

The amplitudes, circular frequencies, and phase shifts

for the most dominant harmonics, determined by spectral

analysis, are shown in table 6.1. These are the parameters

which will be used to drive the numerical model at the open

boundary using equation 4.1. The simulated tidal signal is

n-th Amplitude, Ai Circular Phase Angle, 4
Harmonic (m) Frequency, a0 (rad)
1 0.0847 7.2722 x 10-5 0.8178

2 0.1654 1.4544 x 10-4 2.7241
3 0.0119 2.1817 x 10-4 0.9671
6 0.0079 4.3633 x 10-4 2.5830
Table 6.1
Parameters used for the tidal OBC.

compared against the measured signal from figure 6.8 in figure


6.4.2 Neumann OBC

Since the grid extends into Charlotte Harbor where the

tidal OBC is, the appropriate measures discussed in chapter

4.1 must be applied here. In comparing the area the numerical

model's grid will cover (figure 6.10) to the bathymetric data

(figure 6.1), it is shown that the boundaries were extended

roughly 60 m to the north and 60 m to the south. This adds a

buffer zone so that any circulation in the area of interest

does not cause an instability in the OBC and so that any

effects caused by the Neumann OBC do not adversely affect the

flow patterns in the area of interest. During the model

validation or its application, if the eddies occur near the

Simulated Tidal Cycle


Figure 6.11
Simulated tidal OBC signal compared to the
observed tidal signal in fig. 6.8.

Neumann OB cells, the model's grid must be extended farther to

the north and south.

6.4.3 Canal OBC

The parameters necessary for this OBC are depth at the

end of the canal, the Darcy-Weisbach friction coefficient, the

size of the increment (Axcana), and the number of increments

(Icanal) These parameters are somewhat empirical. The


friction coefficient may be greater in the canal than in the

inlet due to the different bottom sediments characteristics,

excess submerged vegetation (e.g. sea grass), or overgrown

foliage on the banks (e.g. shrubs and mangroves). In the

interest of computational time, as previously discussed, the

total number of increments should not exceed 15. The measured

or estimated length of the canal can be divided evenly into a

number of increments. The right side of figure 6.2 shows the

beginnings of a complex finger canal system. This may add

damping, or friction, to the long wave, so the friction

coefficient should be increased. The finger canal system also

makes it very difficult to obtain an "effective" length if the

canal were straight.

For the Punta Gorda Isles finger canal system, the

surface area of the system was estimated from maps. Since

there are five canal OBCs used on the numerical model and 10

increments were used for each OBC, the surface area was

divided by 50 to obtain the size of increment, Axcanai, which

was 1500 km. The friction factor used in the numerical model

was 0.05. The canal system was assumed have uniform depth,

hence the depth at the end of the canal was set to the depth

at the OBC, or beginning of the canal.

- - - -
y) UU - - - - - - - -
E 0 . '--- -- -- - - -. ' - .' ' ---
E 80

0 0.3 mlsec

20 . . . . . . .

-20 0 20 40 60 80 100 120 140 160 180

Figure 6.12
Maximum flood velocity vectors for numerical
model validation.

6.5 Results

The numerical model was ran for 5 tidal cycles, a

sufficient time for the model to reach steady state from the

initial conditions.

Figure 6.12 shows the maximum flood

velocities and figure 6.13 shows that maximum ebb velocities.

The maximum velocities, both flood and ebb, are about 0.3 m/s,

this value is a little low compared to the electro-magnetic

current meter readings of about 0.4 m/s and the estimated

100 - - - - - -


40- .


-20 0 20 40 60 80 100 120

Figure 6.13
Maximum ebb velocity vectors showing
model validation.

drogue velocities of 0.5 m/s.

140 160 180


The difference in velocities

may be as a result of the canal OBC. One could adjust the

parameters used in the canal OBC to produce better results,

though that could take many iterations.

Comparing figures 6.13 and 6.7, the numerical model does

not show the circulation patterns that the drogue studies

show. The eddy viscosity coefficient was increased to the

point that the stability criteria (eq 3.44) was slightly

violated, though no instabilities became apparent. A possible


causes for this difference are that the bathymetric data may

have errors. Another possible cause is that the bathymetric

relief near the channel may be too great for the model's

resolution. This will be addressed later in chapter 8.

Though a more likely reason that the numerical model was not

able to produce the circulation is because the circulation may

be initiated by superimposing the wave conditions on the ebb

tidal flow. The predominant waves were found to come from the

north which could cause a circulation pattern on the northern

side of the channel during the ebb cycle. Also, the Peace

River empties into Charlotte Harbor just north of Punta Gorda.

The current from the river, a north to south current, or an

ebb current throughout Charlotte Harbor can initiate a

circulation pattern which was observed.


Recall, the objective of this research is to develop a

numerical model to aid in the placement of coastal structures

as an inlet improvement. When applying the numerical model,

any differences or inconsistencies between the field data and

the model validation results must be kept in mind. The

effects of these discrepancies on the numerical output must be

realized when scrutinizing the inlet improvement measures,

such as where a jetty or groin is placed and its length.

Ponce de Leon inlet in Punta Gorda, FL will be used to

show how the numerical model can be used. Two cases will be

demonstrated where a single jetty is installed and a pair of

jetties are installed. The results presented in this section

are a "zoomed in" portion of the model's grid. This allows

the currents in the area of interest to be shown in better


-20 0 20 40 60 80 100 120 140 160 180
Figure 7.1
Flood velocity vectors of a single jetty applied
to Punta Gorda, FL.

7.1 Single Jetty

Figures 7.1 and 7.2 show the maximum flood and ebb

velocities, respectively, for a single jetty case.


comparing figures 7.1 and 6.12, the ebb currents become

greater along the southern bank, the bank with no jetty. The

jetty also causes return flow during the ebb current on the

southern shoreline, though the velocities may not be large

enough to initiate shoreline erosion.

The improvement

measures recommended to the City of Punta Gorda were to instal

100 -

E80 -




-20 0

Figure 7.2
Ebb velocity
Punta Gorda,

vectors of a single jetty applied to

the single jetty as shown and to add rip-rap fill to the

southern shoreline, Wang et al. (1996).

7.2 Pair of Jetties

Figures 7.3 and 7.4 show the maximum flood and ebb velocities,

respectively, for the case with a pair of jetties. The pair

of jetties virtually eliminates the erosional currents in the

proximity of the shoreline. Again, the discrepancies between

-20 0 20 40 60 80 100 120 140

Figure 7.3
Flood velocity vectors of a pair
applied to Punta Gorda, FL.

of jetties

the field data and the validation results, namely the drogue

study and circulations eddies, must be kept in mind.

It is

possible that the pair of jetties may cause significant eddies

outside of the channel near the shoreline which may result in

higher, erosional velocities.





. . . . . . .. 0.3 m/sec


160 180

ioo -- -- i-- i --i --- i --- '--- -- i i -- -
160. . . . . . . .

140. .

E 80- ...... .... . .

6 0 . . . . . . . . . .
0. . . . .
0- 0.3 m/sec

4o0 : - iii. -

-20 0 20 40 60 80 100 120 140 160 180
Figure 7.4
Ebb velocity vectors of a pair of jetties applied
to Punta Gorda, FL.


A two-dimensional numerical model was developed using a

factorization method for numerical solution. The model

developed accounts for constant friction, convective

acceleration, turbulent diffusion, and wind stress. The

numerical model can be adopted to a variety of inlet-river and

-canal systems to study the tidal currents. This model can be

applied to coastal engineering design projects involving the

installation of jetties, breakwaters, groins, or rip-rap fill

where it is necessary to study the feasibility of such inlet

improvement measures.

8.1 Conclusions

The numerical model's efficacy was verified reasonably

well. There is a discrepancy when comparing the numerical

model to the analytical solutions of Massachusetts Bay. A

possible cause of the difference could be the friction factor.

The analytical solutions did not account for friction, where

as if friction is not included in the numerical model, then


the model may become unstable or produce erroneous results.

The analytical solutions also neglected the inertial terms.

Perhaps, if the convective terms were set to zero in the

numerical model, then the results might change.

Overall, the numerical model that was developed can be

applied to a multitude of tidal inlets and areas. The model

was shown that it could produce reasonable results. The

numerical model developed was not able to produce ebb tidal

patterns as measured in the field when applied to Ponce de

Leon inlet at Punta Gorda, FL. Possible reasons for the

differences between the numerical model and the field data

could be bathymmetry and the canal OBC.

Most navigable inlets consist of a deep channel dredged

across a relatively shallow area. When numerically modeling

these inlets, the model may produce better results if the

channel is "numerically smoothed". Figure 8.1a shows the

channel at Punta Gorda and figure 8.1b how one might "smooth"

the relief of the bathymmetry. The sudden relief in

bathymmetry may cause a large gradient in the friction and

convective terms with profound affects on the numerical model.




40 -30 -20 -10 0 10 20 30 40




-40 -30 -20 -10 0 10 20 30 40
Distance (m)

Figures 8.1 a & b
(a) shows the original bathymmetry of channel and
(b) shows the "smoothed" bathymmetry of channel.

8.2 Recommendations

The following studies and improvements to the numerical

model are recommended:

1)The affects of "numerically smoothing" a channel

should be closely studied. From testing the

numerical model as it was being developed, it is

reasonable that a "smoothed" channel will produce


more realistic results. Information was not

presented in this thesis on the grounds that no

analytical solutions or data were available at

the time to support this hypothesis.

2)Develop anther OBC which allows for a lagoon

connected by a small canal as shown in figure

8.2. The governing equation presented by Bruun

(1978), but developed by earlier researchers, for

this case are as follows:

d'lB F AB dfB dflB gA gA
t2 LcAC dt dt LA B (8.1)


lB=surface elevation of bay,


F=frictional losses,

Lc=length of channel,

AB=surface area of bay,

Ac=cross sectional area of channel,

g=gravity, and

io=surface elevation at entrance.


Figure 8.2
An idealized channel/bay system. Figure was
taken from Bruun (1978).

This lagoon/canal system can be applied to the

numerical model as an OBC the same as the canal

and river are discussed in chapter 4.1, except

that the ocean elevation shown in figure 8.2 can

be the surface elevation can be used to drive the

model as at an OB.

3)Morphological changes and sediment transport

should be added to the numerical model. This

will allow the long term sediment transport,

shoaling, and bank erosion/accretion to be

evaluated for a given inlet improvement measure.


Arcilla, A.S. "An Integrated Numerical Approach for Coastal
Engineering Problems." Journal of Coastal Research.
Vol. 5, 1989: 603-616.

Bathan, Karl H. "A Tidal and Wind Driven Numerical Model of
a Lagoon Before and After a Major Modification." J.K.K.
Look Laboratory of Oceanographic Engineering,
University of Hawaii, Vol. 6:2, 1976.

Briggs, Douglas A. and Madsen, Ole S. "Analytical Models for
One- and Two- Layer Systems in Rectangular Basins."
Report No. 172, School of Engineering, Massachusetts
Institute of Technology, Cambridge, MA, 1973.

Bruun, Per. Stability of Tidal Inlets. Elsevier Scientific
Publishing Complany, Amsterdam, 1978.

Cheng, Ralph T. "Modeling of Tidal Residual Circulation in
San Francisco Bay, California." Two-Dimensional Flow
Modeling. U.S. Army Corps of Engineers, The
Hydrographic Engineering Center, Davis, California,

Choi, Jei-Kook. "Three-Dimensional Curvilinear-Grid Modeling
of Baroclinic Circulation and Mixing in a Partially
Mixed Estuary." Dissertation, University of Florida,
Gainesville, FL, 1992.

Henderson, F. M. Open Channel Flow. Macmillan Publishing
Co., Inc., New York, 1966.

Hirsch, Charles. Numerical Computation of Internal and
External Flows. Volume 1, John Wiley & Sons Ltd.,
Chichester, 1994.

Pedlosky, Joseph. Geophysical Fluid Dynamics. Springer-
Verlag, New York, 1979.

Reeve, D. E. and Hiley, R. A. "Numerical Prediction of Tidal
Flow in Shallow Water." Computer Modelling of Seas and
Coastal Regions. Elsevier Science Publishers Ltd., New
York, 1992: 197-208.

Reidel, H. P. and Wilkinson, F. L. "Numerical Modeling An
Aid to Assessing Field Data." Coastal Engineering.
1976: 3243-3262.

Roll, H. U. Physics of the Marine Atmosphere. Academic
Press, New York, 1965.

Schomaker, Christine W. "A Model for Tidal Circulation
Adapted to Monterey Bay, California." Thesis, Naval
Postgraduate School, Monterey, California, 1983.

Sheng, Y. P. Mathematical Modeling of Three-Dimensional
Coastal Currents and Sediment Dispersion: Model
Development and Application. Technical Report CERC-83-
2, Aeronautical Research Associates of Princeton,
Princeton, New Jersey, 1983.

Wang, Hsiang; Lin, Lihwa; and Wang Xu. "Ponce de Leon
Channel Improvement Studies, Punta Gorda, Florida."
UFL/COEL-96/004, Coastal and Oceanographic Engineering
Department, University of Florida, Gainesville,
Florida, 1996.

Weiyan, Tan. Shallow Water Hydrodynamics. Elsevier Science
Publishing Company, Inc., New York, 1992.

Wilson, Basil W. "Note on Surface Wind Stress over Water at
Low and High Wind Speeds." Journal of Geophysical
Research. Vol. 65, No. 10, 1960: 3377-3382.


Peter N. Seidle was born in Binghamton, NY on 28

January, 1972. His family moved to Maryland in suburban

Washington DC where he started kindergarten and graduated

from Thomas S. Wootton High School, Rockville, MD, in 1990.

In 1994 he graduated from Widener University, Chester, PA,

with a Bachelor of Science in Mechanical Engineering.

Shortly after graduating from college, he became employed by

the ARIST Corporation to work as a program analyst in the

Ballistic Missile Defense Organization, which is part of the

Department of Defense. In August 1995 he began his graduate

studies at the University of Florida in the Coastal and

Oceanographic Engineering Department. In May 1997 he

graduated with a Master of Science in the field of study he

has and will always enjoy and be fascinated with.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs