A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS by
Peter N. Seidle Thesis
A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS
PETER N. SEIDLE
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
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.
TABLE OF CONTENTS
LIST OF FIGURES..................
... . .. . .. . 1
... . . . . . .vi
2. GOVERNING EQUATIONS...........
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..
.1.. .. . .. . .
.1.. .. . .. . .
... . . .. . . .3
... . . .. . . .4
... . . . . . .10
... . . . . . .10
... . . . . . .15
... . . . . . .17
... . . . . . .17
... . . . . . .19
... . . . . . .23
... . . . . . .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
5.1 5.2 5.3
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
Single Jetty....................................... 72
Pair of Jetties.................................... 73
8. CONCLUSIONS AND RECOMMENDATIONS......................... 76
8.1 Conclusions........................................ 76
8.2 Recomminendations................................... 78
LIST OF REFERENCES.......................................... 81
BIOGRAPHICAL SKETCH......................................... 83
LIST OF FIGURES
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 wiith 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 measuremets ................................ 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 OC 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
A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS
Peter N. Seidle
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.
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 inlets.
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.
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 changes.
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:
l)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 follows:
89 BuD 8vD
Bt ax By
Bu Bu Bu 89 f 82U B2U x-wnd
a+u-+v=-g- fuu+A,(-+2)+ (2.2)
8t 8x ay ax 8D 8x2 ay 2 pD
8V 8V 8V 84a f 8-%>2 B2V \ -wn
-a +u-a+v--g l f V+A-v a2v +-)+- (2.3)
at ax ay ay 8D ax2 y 2 pD
u=depth averaged velocity in the x-direction, v=depth averaged velocity in the y-direction,
u=depth averaged velocity vector,
p=water surface elevation,
x,y=right handed Cartesian coordinates, j=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, Twind 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:
Z .id P.J Cd
- i WjosinO (2.4)
pD p D
-wind Pair Cd (2.5)
___ -WjcoO 25
pD p D
W10=wind speed at 10m, and
Although the drag coefficient is empirical, researchers have
Relationship between drag coefficient (shear
stress coefficient) and wind speed. Figure was
taken from Roll (1965).
found that it is on the order of i0-. 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 I0' and for stronger winds, Cd=2.37 x i0-3. For validating this numerical model, Wilson's findings were used.
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.
J iiiiiiii ii i!i
The numerical method adopted for the model is the socalled "factorization" 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 "factorization" numerical method. The governing equations now become:
a_+Dau +DA=0 (3.1)
at ax ay
at .-a, (3.3)
The above equations 3.1, 3.2, and 3.3 can be combined as a system of equations by letting w=[ u, V]T The system of equations can be written as: Iwt+Awx+BWy=O (3.4)
where I is an identity matrix and A and B are coefficient matrices as follows:
B= 0 0 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:
Wnt -+Awn +Bw i=0 (3.5)
At x11 y
Multiplying by At and factoring out wn"1 results in:
wn+1(1+AtA8x+AtB8 )=w n (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:
wnIl(1 +AtA8x+AtB86 +(At)2AB8x8 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+I(I +AtA6x)(1+AtBy)-02=(I -AtB8)w n+AtB8 w (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+2 w*-n and equation 3.8 1+AtA8x
can be rewritten as:
w n+(1 +AtAx)(1 +AtB6y)-02=(1 +AtA6x)w *+AtB8 yw n
Moving the error term to the RHS of the equation and reexpanding it the above equation then becomes:
w n+(1 +AtA8)(1 +AtB6 )=(1 +AtA8x)w +AtB8 w n+(t)2AB88,w "
Factoring out AtB8y from the two last terms on the RHS yields:
wn+(1 +AtA6x)(1 +AtB8y)=(l +AtA8 )w +AtB6y(l +AtA8x)w n
Lastly, removing the common term, 1+AtAx from the equation results in:
(1 +AtB8 )wn.+l=w*+AtB8yWn
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 +AtA8x)w *=(1 -At )WB8 (I +AtBy)w"+1=w* +AtB yw n
Expanding equation 3.13 gives the following:
7 +AtD8xu ="-AtD6yv "
u *+Atg6x8' =un (3.16)
v* =v "-Atg6 q]" (3.17)
and expanding equation 3.14 gives:
"n +AtD8v""l=p*+AtD v n (3.18)
v n+Atg6 pnI =v +Agyqn (3.20)
Note that equations 3.15 and 3.16 are a system of simultaneous equations where T' and u* are the two unknowns, as well as equations 3.18 and 3.20, where ,l'n and vn"i are the unknowns. Since equation 3.17 is independent of r* and u* it can be substituted directly into equation 3.20. Likewise, equation 3.19 states that unl =+=u so equations 3.15, 3.16, 3.18, and 3.20 can be adjusted accordingly. The four resulting equations are as follows:
' +AtD8xun+=I"-AtD8 yv n (3.21)
Atg8x'l*+un I=u n
rg n+1+AtD6 Vn+l=r*+AtD8V n (3.23)
A t n+I+vn I=V (3.24)
The above "factorization" 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 (ysweep). It is also important to note that equations 3.22 and 3.24 are essentially the momentum equations in the x- and ydirections, respectively. This will allow the non-linear terms to be added. Moreover, when equation 3.21 is
substituted into equation 3.23 for r* ,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
applied to the velocities for the continuity equations (2.21 and 3.23) and may be written as:
1, +AtD(unj-u )=T-AtD(v,-vin)
A g (l -T 7 1 + U 1= u n
tg~ij-) l +uj =ij
A backward difference is applied to the surface elevation for the momentum equations (3.22 and 3.24) and are as follows:
Atg( n+l_ n+l+n+l-v n grij l) =
n+l n+ n+I n
rij +AtD( vjl-i )=il* +AtD( vij, -
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 -Lulu This term is non-linear in that it is
quadratic in terms of flow velocity vector. Here the velocity
vector term, Jul is expressed numerically in terms of the known velocity components in the previous time step, such as u and vn' Furthermore, since the x-direction momentum
equation, eq. 3.26, is solving for un'at the grid center depicted in Figure 3.2, lul must be calculated at the proper location. The numerical equation used to calculate the
Numerical cells applied to friction
magnitude of the velocity vector for the x-direction momentum equation is as follows:
n n n
u (u )2 VI- I +1 I I V-1 ij)2 (3.29)
lI l = uid)+ 4
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 ,n 2 V I-Ii+I +Vi+il +Vi-lj+vifi2 n+l
+v++ 4 ij
Similarly, the frictional forcing term added to the LHS of the numerical y-direction momentum equation, eq. 3.28, is:
n n n n
( Id+ +J d ij-1 i+lj-12+(f. 2. +
+ 1 4 ) 2 v 7
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 terms.
For the x-direction momentum equation, the convective du du
acceleration terms are u-+v- taken from equation 2.2. In
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
Numerical cells applied to convective
with the upstream convection phenomenon. Both convective
terms, u- and v- need to be examined separately.
To analyze the term u- first the velocity direction ax
needs to be determined. An average of the velocities between
the two cells (figure 3.3), ui-112j and u+112, are used to
identify the direction of the velocity so that the proper
upwind differencing method may be applied. If ui 1/2j is
positive, then the upwind difference is u'n-u The
numerical equation used to evaluate the first convective term
u=u. uij-ui-f (3.30)
-ax -1 /2 Ax
However, if ui+1/2j is negative, then the upwind deference
is uij-utl and the equivalent numerical difference is as follows:
aU -Un Uij-Ui+lj (3.31)
ax i12j Ax
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 vin/2j is in the positive direction, then the numerical convection term is as follows:
y pin Uli/ i_ ii-1 (3.32)
Vay i-1/2' Ay
And, if Vi1/2,+1 is negative, the numerical differencing is:
u nH Uig-Uij+l (3.33)
iV_ /2;i+1 ay Ay
At every time step, the direction of the velocity is determined for each element so that either equations 3.30 or
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
for u-+v- can be used to find the numerical equations for
the convection terms in the y-direction momentum equation,
u-+v- For a positive velocity direction, the two
convective terms are as follows:
av u- i vi-Ij (3.34)
ax -Ui-112 Ax
v n Vi Vi;V-ay=E Vi-1/2 y
Likewise, if the velocity is negative then the terms are as:
aV n Vij -ViIlj
-ua-- u 1il/2 A
8x 1/2; Ax
(9v n Vii'-Pid+1 ay fi-1/2 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: A 2u + 2U
is obtained by simply applying a center difference to yield:
n n N N U n+n
A(ui+lzj i+u-lj uid+-jl ij u i-1)
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 n( 'n -2vyn" + v nzl -2nin +V- l) Av.i+lj 4j Vi-ij ti- I S V IJ )
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= At ( +lu l+lvl) (3.38)
Ax=x-direction grid spacing, Ay=y-direction grid spacing,
lul=maximum x-direction velocity superimposed on
the tidal velocities, and
IvI=maximum y-direction velocity superimposed on the
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):
AV AX, +AY 2 < 1 (3.39)
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 (ECs) 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
j=J M I;:.. .. .. ...
. .....63 ....
=! i "i2 ii i iIi
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 OpTen 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.
The tidal OB is defined by water level changes such as a tidal fluctuation of the form r=Acos(ut+( ) where A is the amplitude, a is the period, and 0 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:
il=E (A~cos(at+4,)) (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.
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 .: i i i .:'~ i i:.. . .. .. .. .. .. .. ..... . . . .... .
11.: :::::::: ...... .....
X ... . . . ..:::.. ...............
J -1 :* ": . . .... .. ................ T
.. !..'.. :!.... ...: ............... T
2 ..:......:. T
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 timedependent 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: L + auD=o (4.2)
River or Canal OB
Discharge at head
Figure 4.3 Example of a grid for a river or canal open boundary condition.
au + au a,, f JuJu at ax ax 8D
u=depth averaged velocity in the direction of the
river or canal, t-time,
Tl=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.
The analytical set-up due to wind stress is (Choi, 1992): ll=Tind~x L) 51
Twind =Shear stress by wind,
p=density of water,
x=distance from middle of basin,
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 mn 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 rn/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 profile.
The reason there is a difference between the numerical and analytical set-ups is that the model's output is rounded to 0.1 nun, 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 mn. 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.
200 250 Distance, (m)
Figure 5.1 Comparison set-up.
of analytical and numerically modeled
The model was also tested for its ability to replicate a standing wave in a closed basin. The basin length required for seiching is:
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 elavation at the right hand side of the basin as a function of time, i=Acos-) 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
1 .... ....I
0 10 20 30 40 50 60 70 80 90 100
-0 .1 I I I I I I
0 10 20 30 40 50 60 70 80 90 100 .1 I I I I I I I
0 10 20 30 40 50 60 70 80 go 1iO0 Distance, (kmn)
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-I1 I I I I I I I
o 10 20 30 40 50 60 70 80 90 10
o 0 i1 20 30 40 50 60 70 80 90 C
0 10 20 30 40 50 60 70 80 90 Ic
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
Massachusetts Bay. Figure taken from Briggs and Madsen (1973)
(1973), to the linearized continuity and momentum equations (eqs. 3.1, 3.2, and 3.3) is:
= c 2m~sinmy(sinknX2 sinkjx1) coshmy
=Ae i [Cos mYnE mk (x2-xI) sinhmyocosknx] (5.3)
and x., x, X21, and y0 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 A.. 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
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 Bay.
J (D P- P-.
O Q PO Q
rt (D i'
- H- X (1
(-- I- s
n 0 HSH- ~0 r
0.0 1.0 2.0
L I I
Velocity Scale CFt./sec.)
I I I I
For the numerical model, the bay was divided into 30 increments in the x-direction and 10 increments in the ydirection. 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 neglects.
Velocity Vectors of Massachusetts Bay
0 10 20 30 40 50 6(
Maximum velocity vectors from numerically modeling Massachusetts Bay.
L II I t i t tit
I I I I t t I I i I % '
. . . .. . .. =1 f/sec
Surface Elavation (ft) of Massachusetts Bay
I I I I
10- 4.8 4.85
I I I I I
0 10 20 30 40 50
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 simulated.
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.
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.
Figure 6. 1
Bathymetric survey of Ponce de Leon Inlet, figure
taken from Wang, et al. (1996).
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.
Location of stations where tide and current data were collected, figure taken from Wang, et al. (1996).
* I p
- Tide at Channel Marker 8 S-lide at Channel Marker 18
-0.4- ... Current at inlet entrance
) Oddmeter data
0 2 3 4 5 6
Figure 6.3 Continuously observed tidal record and discrete current measurements, figure taken from Wang, et al. (1996).
Two instruments were used to collect current data: an electro-magnetic meter, and a vane-type meter. The electromagnetic 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 flood.
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
" 0-- 39 90.28
0 0.42 0.51 0.54
Horizontal Distance (m)
Flood velocity data taken inside the canal using the Ottmeter, figure taken from Wang, et al. (1996).
Station 2 Time=1 4:34
Horizontal Distance (m)
Ebb velocity data taken in the canal using the Ottmeter, figure taken from Wang, et al. (1996).
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 pile.
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
Ponce De Leo
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.
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
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
(2) since the model is being used as only an aid to
inlet morphological improvement measures,
residual currents and tidal flushing are not
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.
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 0pen 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, i
Harmonic (M) Frequency, a1 (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
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
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, Ax,,,,,,,, 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.
, I VI I
-~ ~ -0. -mlse----c ----20 0 20 40 60 80 100 120 140 160 180 meters
Maximum flood velocity vectors for numerical
The numerical model was ran for 5 tidal cycles, a sufficient time for the model to reach steady state from the
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
80 . . . . . .
60 40 20
0II I I
-20 0 20 40 60
Maximum ebb velocity vec
drogue velocities of 0.5 m/s.
80 100 120 140 160 180 meters
tors showing numerical
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 detail.
. ................ .. o
40 .1 I I
*~ ~ ~ .. .. . .
-20 0 20 40 60 80 100 120 140
Figure 7.1 Flood velocity vectors of a single jetty 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.
measures recommended to the City of Punta Gorda were to instal
-20 0 Figure 7.2 Ebb velocity Punta Gorda,
80 100 120 140 160 180 meters
a single jetty applied to
vectors of FL.
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
Flood velocity vectors of a pair applied to Punta Gorda, FL.
the field data and the validation results, namely the drogue
study and circulations eddies, must be kept in mind.
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 ., -.0.3m/sec
l U I I I I I I I
160 140 120
E 80 .
6 0 .. ..
20 - 0.3 m/sec
-20 0 20 40 60 80 100 120 140 160 180 meters
Ebb velocity vectors of a pair of jetties applied to Punta Gorda, FL.
CONCLUSIONS AND RECOMMENDATIONS
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.
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 bathynunetry 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
bathyinmetry 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 bathymetry of channel and
(b) shows the "smoothed" bathymmetry of channel.
The following studies and improvements to the numerical model are recommended:
l)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:
d2 'B F AB dB d'IB gAc gA(
+ + p
dt 2 2Lc Ac di dt LCAB LCAB
'IB=surface elevation of bay,
F=frictional losses, Lc=length of channel,
AB=surface area of bay,
Ac=cross sectional area of channel,
To=surface elevation at entrance.
aAREA- A- .. BAY FRESH
PMM. c'--TH.h WATER
OCEAN OLUME A'
ELEVATION Lc SURFACE AREA AO
ELE TIN 70 ELEVATION
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.
LIST OF REFERENCES
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
Institue 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.,
Pedlosky, Joseph. Geophysical Fluid Dynamics. SpringerVerlag, 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.
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-832, 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,
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.