NU14ERICAL CALCULATION OF THE RESPONSE OF COASTAL ',ATERS
TO STOPI SYSTEMS WITH APPLICATION TO HURRICANE CA,1IlLLE
OF AUGUST 1722, 1969
By
Bryan Rowell Pearce
A DISSERTATION PRESENTED TO THE GRADUATE
COUNCIL OF THE UNIVERSITY OF FLORIDA lIN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PH ILOSOPHY
UNIVERSITY OF FLORIDA
1972
DEDICATION
This dissertation is dedicated to my son,
Nicholas James, born August 14, 1972.
ACKNOWLEDGEMENTS
The author wishes to thank Dr. R. G. Dean, who initiated
and supervised this investigation, and Dr. B. A. Christensen.
The hours of discussion and field trips with Dr. Dean have provided
a unique educational experience and contributed to a better under
standing of the subject. The helpful criticism and conferences
with Dr. Christensen were invaluable and provided the extrospective
view that is necessary to any extended work.
Also to be thanked are Susan Phillips, Elizabeth Maltz and
Marilyn Morrison who typed the manuscript and Denise Frank who
drafted the figures.
The study was financed by the National Science Foundation.
All computations were carried out using the University of Florida
IBM 360/65 computer.
TABLE OF CONTENTS
Page
AC IOWLEDGEIMENTS iii
LIST OF FiGURES vi
LIST OF SYMBOLS x
ABSTRACT xv
CHAPTER
1. INTRODUCTION 1
2. STOPRM SURGE CALCULATIONS THEORY 7
a. Volume Transport Equations 7
b. Assumptions and Boundary Conditions 13
c. The Hurricane Model 20
d. Analytical Treatment Steady State,
Without Bottom Friction 29
i. Problem Definition 29
ii. Uniform Wind Stress Acting Over a
Shelf of Uniform Depth 30
iii. Uniform Wind Stress Acting Over a
Shelf of Uniform Slope 34
e. Analytical Treatment Moving Wind Stress
and Pressure Field 38
i. Problem Fonnulation 38
ii. Solution 40
Page
3. THE NUMERICAL MODEL 56
a. OneDimensional Numerical Model 56
b. Development in Two Dimensions 64
c. Effects of Discretization and Checks of
the Numerical Model 73
4. RESULTS OF NUMERICAL CALCULATIONS 87
a. Introduction 87
b. Numerical Calculations and Comparison
with Available Data 94
5. FIELD DATA AND INSTRUMENTATION 118
6. SUMMARY AND CONCLUSIONS 132
APPENDIX
A. DERIVATION OF GOVERNING EQUATIONS FOR SURGE
HYDRODYNAMICS 134
B. TRANSFOPSATION OF GOVERNI;IG DIFFERENTIAL
EQUATIONS INTO FINITE DIFFERENCE FORM 140
1. Transformation Without Convective Terms 140
2. Effect of Convective terms 146
REFERENCES 146
BIOGRAPHICAL SKETCH 153
LIST OF FIGURES
Figure Page
21 PERSPECTIVE OF NEARSHORE ZONE SHOWING COORDINATE
SYSTEM AND JNOMENICLATURE 8
22 COMPARISON OF VAN DOPN 'S k WITH OCEANIC WIND STRESS
MEASUREMENTS 15
23 TRANSIENT WATER VELOCITY DISTRIBUTION IN A CLOSED
CHANNEL DUE TO IMPULSIVE WIND 18
24 DEFINITION SKETCH FOR STORM TRAVELING AT CONSTANT
VELOCITY 22
25 STORM TRACK FOR HURRICANE CAMILLE, AUGUST 1618,
1969 23
26a PRESSURE CHARACTERISTICS OF HUPRICANE CAMILLE AS
A FUNCTION OF TIME 24
26b RADIUS TO MAXIMUM WINDS AND FORWARD SPEED OF
HURRICANE CAMILLE AS FUNCTIONS OF TIME 24
26c MAXI.MUIM WIND SPEED IN HURRICANE CAMILLE AS A
FUNCTION OF TIME 25
27a COMPARISON OF 30 FOOT SURFACE ISOTACHS 8/17/1200 26
27b COMPARISON OF 30 FOOT SURFACE ISOTACHS 8/17/1800 27
27c COMPARISON OF 30 FOOT SURFACE ISOTACHS 8/18/0000 28
28 NOTATION FOR ONEDIMENSIONAL PROBLEM 30
29 ONEDIMENSIONAL SHELF WITH UNIFORM DEPTH 31
210 SURGE HEIGHT VERSUS FETCH/DEPTH FOR UNIFORM DEPTH 33
211 ONEDIMIENSIONAL SHELF WITH UNIFORMLY SLOPING BOTTOM 34
212 SOLUTIONS TO UNIFORM SLOPE CASE AT THE COASTLINE 36
213 IDEALIZED SHELF BATHYM.ETRY 40
214 IMPULSIVE FORCING FUNCTION MOVING AT CONSTANT SPEED 43
Figure Page
215a ARBITRARY FORCING FUNCTION 44
215b TRIANGULAR FORCING FUNCTION 44
216 DEFINITION SKETCH FOR SHEAR STRESS DISTRIBUTION 47
217 RESPONSE TO TRIANGULAR SHEAR STRESS DISTRIBUTION 51
218 DIMENSIONLESS SURGE HEIGHT MAXIMA AS A FUNCTION
OF DIMEISIONLESS STORM WIDTH AND VELOCITY 52
219a DEFINITION SKETCH FOR NUMERICAL CALCULATIONS 54
219b FINITE DIFFERENCE MODEL OF TRIANGULAR SHEAR
STRESS DISTRIBUTION 54
220 COMPARISON OF NUMERICAL AND ANALYTICAL SOLUTIONS 55
31 DEFINITION SKETCH FOR ONEDIMENSIONAL NUMERICAL
CALCULATIONS, SECTION NORMAL TO THE COAST 58
32 MOVEMENT OF STORM 61
33 GRID ELEMENT DESCRIPTION FOR TWODIMENSIONAL
SCHEME 66
34a COMPONENTS OF q (i,j,n) 68
34b COMPONENTS OF qx(i,j,n) 68
35 SCHEMA OF GRID SYSTEM AND BOUNDARY CONDITIONS
EMPLOYED IN TWODIMENSIONAL CALCULATIONS 69
36 FILTERS USED TO CALCULATE WATER VELOCITIES 71
37a COMPARISON OF MAXIMUM SURGE HEIGHTS CALCULATED
BY THE ONE AND TWODIMENSIONAL MODELS 74
37b PROFILES OF STORM CHARACTERISTICS USED IN
COMPARE I SON 74
37c YCOMPONENT OF WIND VELOCITY 75
38 COMPARISON OF THREE ONEDIMENSIONAL MODELS 77
39 IDEALIZED SHELF RESPONSE PROBLEM WITH UNIFORM
DEPTH 79
310a ISOLINES OF DIM'ENSIONLESS STORM TIDE, ANALYTICAL
SOLUTION 81
Figure Page
310b ISOLINES OF DIMEN SIONLESS TRANSPORT IN
xDIRECTIO;, ANALYTICAL SOLUTION 82
310c ISOLINES OF D1;;ENSIOILESS TRANSPORT IN
yDIRECTIOil, ANALYTICAL SOLUTION 83
311 DEFINITION SKETCH FO, ;NU:MIERICAL MODEL 84
312 EFFECT OF LATERAL BOUNDARY PROXIMITY AND
RELATIVE GRID SIZE ON MODEL RESPONSE 86
41 COASTAL REGION MODELED 88
42a LARGE GRID MODEL 89
42b DEPTHS USED IN LARGE GRID MODEL 90
42c SMALL GRID MODEL 92
42d DEPTHS USED IN SMALL GRID MODEL 93
43a CALCULATED STORM SURGE ELEVATIONS 8/17/0920 95
43b CALCULATED STORM SURGE ELEVATIONS 8/17/1640 96
43c CALCULATED STORM SURGE ELEVATIONS 8/17/2000 97
43d CALCULATED STORM SURGE ELEVATIONS 8/17/2200 98
43e CALCULATED STORM SURGE ELEVATIONS 8/17/2340 99
44 COMPARISON OF OBSERVED AND CALCULATED SURGE
HEIGHT MAXIMA 100
45 SENSITIVITY OF SURGE CALCULATIONS TO VARIATIONS
IN f AND Xs 101
46 SCHEME USED TO EXTRAPOLATE SURGE HEIGHTS TO COAST 102
47 COMPARISON OF EXTRAPOLATED AND NONEXTRAPOLATED
SURGE HEIGHTS AT BILOXI BAY, MISSISSIPPI 104
48a COMPARISON OF OBSERVED AND PREDICTED SURGE AT
DAUPHIN ISLAND, ALABAMA LARGE GRID MODEL 105
48b COMPARISON OF OBSERVED AND PREDICTED SURGE AT
DAUPHIN ISLAND, ALABA'A SMALL GRID MODEL 107
49a WATER VELOCITY AND WIND VELOCITY RECORDS FOR
HURRICANE CAMILLE AT EGLI;I AFB, FLORIDA 108
vii i
Figure Page
49b COASTAL TOPOGRAPHY AT SITE OF WATER VELOCITY
MEASUREMENTS 109
410 COMPARISON OF CALCULATED AND MEASURED WATER
VELOCITIES 110
411 LONGSHORE CURRENT AND BEACH PROFILE MEASURED
AT MUNICIPAL PIER, PANAMA CITY, FLORIDA 112
412 COMPARISON OF MEASURED AND CALCULATED WIND
SPEEDS AT SITE OF WATER VELOCITY MEASUREMENTS 113
413 COMPARISON OF PREDICTED AND ACTUAL FLOODING
ON MISSISSIPPI COAST 115
414 COMPARISON OF PREDICTED AND MEASURED FLOODED
AREAS PER UNIT LENGTH ALONG COAST 117
51 SUMMARY OF CRITICAL WARNINGS 120
52 HINDCAST WAVE FIELDS FOR HURRICANE CAMILLE 122
53 NUMBER OF TIMES DESTRUCTION WAS CAUSED BY
TROPICAL STORMS, 19011955 124
54 PROPOSED FIXED CONTINUOUSLY MAINTAINED
INSTRUMENTATION SITES 125
55 SCHEMATIC OF RECORDING UNDERWATER TIDE GAUGE 126
56 RECORDING UNDERWATER TIDE GAUGE 127
57 TIDE MEASUREMENTUNDERWATER GAUGE 128
58a CURRENT METER DEPLOYED 130
58b PHOTOGRAPHIC RECORDING MECHANISM OF CURRENT
METER 130
LIST OF SYMBOLS
an Fourier coefficients
Bn Fourier coefficients for forcing function
c wave celerity, Jgh
c1,c2 constant coefficients used in moving wind stress development
C storm or forcing function velocity
CDT central daylight time
c* C/c
D total depth
f Darcyleisbach friction factor
fL linear friction factor
F,F' forcing function
Fx,F resistance coefficients for x and y equations of motion
g gravitational acceleration
G;.;T Greenwich ;.lean Time
h water depth, referenced to Mean Sea Level (';.SL)
h water depth at shelf break
hi water depth at coastline
H11,12 definitions used to simplify moving wind stress development
HF hindcast wave height
i indicial functional parameter
k wind stress coefficient
z shelf width
m bottom slope
MSL mean sea level
n indicial functional parameter
p pressure
pO surface pressure at center of hurricane
p surface pressure far from storm
pn surface pressure
q water transport
q water transport component normal to boundary
qx,qy water transport components in x and y directions
r distance from center of hurricane
R radius to maximum winds
t time variable
t* dimensionless time variable, t/T
t' dummy variable
t initial time
tl,t2,t3 time that zeros of idealized forcing function cross shelf break
T function of time; 4z/C
TF hindcast wave periods
u component of water velocity in x direction
u time average of u
u' instantaneous fluctuation of u
U wind velocity at 30 feet above the surface
Uy x and y components of U
Ucr U critical, in Van Dorn's expression for .,ind stress coefficient
Uc cyclostrophic wind
U geostrophic wind
UG gradient wind
v component of water velocity in y direction
v time average of v
v' instantaneous fluctuation of v
V forward velocity of storm
V' V sine
w component of water velocity in z direction; wind. field breadth
w time average of w
w' instantaneous fluctuation of w
w* lateral dimension of shelf incorporated into numerical model
x,xl horizontal coordinate, normal to coast
x*,x' dimensionless horizontal coordinates, x/Q
xc coordinate of coastline
X function of x; dimensionless depth, x/xB
Xs starting point of storm as distance from landfall
y,X2 horizontal coordinate, parallel to shore
y' dimensionless horizontal coordinate, y/1
z,X3 vertical coordinate
aB angle between breaking waves and shoreline
an coefficient in moving wind stress analytical development
a* dimensionless parameter, w/w*
a** dimensionless parameter, c/w*
B constant, T /og
X
xii
"'V Dg
6 impulse function
Ap central pressure index
Lt time increment
AX,Ay grid element dimensions
AT shear stress increment
c dimensionless parameter, mh/l
o, dimensionless parameter, mho /
c dimensionless parameter, mhl/B
n surge height, referenced to MSL
n' dimensionless surge height
e angle between point of observation and line of storm movement
shear stress parameter
\ width of idealized storm
\* dimensionless storm width, .\/
p molecular viscosity
( mn/E; dummy variable
o density of seawater
oa density of air
0,on cnn/2z
o',an' Cnn/2Z
T shear stress
T surface shear stress
T components of shear stress
xTy,
1x and y components of surface shear stress
xq,y
xiii
Tb x and y conoonents of bottom shear stress
x,y
T' Reynolds stress
( latitude
w the earth's angular velocity
V gradient operator
v2 Laplacian operator
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial
Fulfillment of the Requirements for the Degree of Doctor of Philosophy
NUMERICAL CALCULATION OF THE RESPONSE OF COASTAL WATERS
TO STORM SYSTE;IS WITH APPLICATION TO HURRICANE CA'ILLE
OF AUGUST 1722, 1969
By
Bryan Row.ell Pearce
August, 1972
Chairman: Robert G. Dean
Cochairman: Bent A. Christensen
Hajor Department: Civil and Coastal Engineering
A numerical model is employed to describe the water motion in a
coastal region associated with the passage of a hurricane or severe
storm. The model is twodimensional, employs the vertically integrated
or "tidal" equations of motion, and is used to describe the specific
case of Hurricane Camille of August 1'69. Two models are employed, a
large grid with 16 nautical mile grid elements and a small grid model
with 6 nautical mile grid elements. The results of the two models were
not significantly different. The maximum storm surge or maximum
water level above mean sea level was calculated and found to agree to
within twenty per cent with maximum surge heights determined from
high water marks. Calculated surge height was found to be insensitive
to bottom friction coefficients varying from .005 to .02. inclusion
of th? nonlinear convective terms affected the calculated surge
height maximum for Hurricane Camille by less than two per cent and
affected the surge in any grid element adjacent to the coast by less
than five per cent. Twodimensional storm surge plots at different
times are presented.
Hurricane generated currents were calculated and compared to data
taken off the Florida coast. It is concluded that more actual current
data are necessary before hurricane generated currents can be cal
culated with confidence. The current calculations were found to be
sensitive to bottom friction and subject to "windup".
In addition to numerical calculations, basic analytical cases
are covered, including the response of a shelf of uniform depth to
a triangular wind stress distribution moving with constant velocity.
xvi
CHAPTER 1
INTRODUCTION
The problem of predicting the response of coastal waters to
extreme wind systems is becoming increasingly important as the coastal
regions of this country are developed for residential, recreational
and industrial uses. Disasters such as the ones occurring in recent
years at Galveston as a result of Hurricane Carla and along the
Mississippi coast as a result of Hurricane Camille demonstrate that
action should be taken to minimize the tragic consequences of such
events. Prediction of the vulnerability of particular areas to flood
ing in similar storms would be of considerable interest to local
residents and governments as well as insurance underwriters. Moreover,
real time forecasting of storm surge elevation, with confidence, for
an approaching storm would be of great value in planning and imple
menting evacuation and other measures to minimize damage. This is
clearly illustrated by Hurricane Agnes (1972). Agnes was a hurricane
of large areal extent but relatively low intensity. Little attention
was paid to surge forecasting, and as a result the residents of some
areas on the west coast of Florida were surprised by a storm tide 3 to
8 feet above normal. Damage in the TampaSt. Petersburg area was
estimated at about $20 million. Undoubtedly some of this could have
been prevented if adequate real time surge forecasts had been avail
able.
2
The design of permanent facilities at or near the shoreline
must include consideration of the maximum and minimum tides that may
occur as well as the wave action. In the operation of nuclear power
plants, for example, one consideration is the availability of
sufficient cooling water to ensure an orderly shut down prior to a low
tide occurence below the level necessary for plant operation. Thus,
a proper study of the situation would establish a maximum elevation at
which the inlet facilities could be placed.
Various approaches are used in solving this type of problem.
The purely analytical approaches are limited to idealized wind
stresses, pressure distributions, and topography. Lamb [1] considers
the case of "free" long waves in an infinite canal of uniform depth.
In addition to gravity, small disturbing forces are considered to be
acting on the fluid, which vary only by a small fraction of their value
within distances comparable to the depth. This leads to an in
homogeneous wave equation similar to that of Section 3e of the text.
Lamb then finds a solution for a pressure distribution that travels
with unchanged form at constant velocity. The presence of boundary or
terminal conditions greatly complicates the solution. The case of a
finite shelf is considered in the text.
Reid [2] presents the approximate response of a sloping shelf
due to a triangular shear stress distribution traveling directly on
shore. The results, obtained at the shoreline, are based upon the
linear onedimensional wave equation. The solutions are obtained
using a graphical technique and the method of characteristics. He
considers many different cases of fetch lengths and storm speed and
summarizes his results in graphical and tabular form.
Further analytical work in the area of storm surges has been
and will be done, yet the complexities of the actual forcing functions
and coastal topographies lead one to realize that another approach is
necessary, if realistic solutions to actual problems are to be
obtained. Hansen [3,4] is one of the pioneers of storm surge cal
culations. He calculated surges arising from storms occurring in the
North Sea, using the linearized equations of motion. That is, he
ignored surge height with respect to depth and ignored the convective
terms. He laid the foundation, however, for the present work in
numerical surge calculations. Many contemporary investigators still
use, in some form, the staggered integration scheme he devised.
Platzman [5] analyzes the case of a pressure distribution or
squall line moving over Lake Michigan causing a surge of more than
six feet on the Chicago lake front. The volume transport or ver
tically integrated equations of motion were used in central difference
form in conjunction with a twodimensional bathymetry and a moving
pressure disturbance. Platzman's calculated surge is about one
half of the observed surge and he concludes that this is probably due
to the omission of wind stress.
Miyazaki [6] treats the Gulf of Mexico as a closed basin,
assuming the surge height to be zero at the mouth between Yucatan,
Cuba and Florida. This procedure eliminates the difficult problem
of specifying lateral boundary conditions in open water. The grid
he uses is variable and finest at points of interest in the coastal
regions. He includes bottom stress in a form given by Reid [7] which
is a function of both flux, q, and of surface shear stress, Tn
Marinos and Woodward [8] calculate storm surges using the
assumption of bathystrophic flow (currents parallel to the bottom
contours). They arrive at a "quasionedimensional" scheme which
considers only onedimensional depth profiles but includes the
Coriolis effect. The bathystrophic procedure was originally developed
by Freeman, Baer and Jung [9].
Reid and Bodine [10] present a detailed numerical model for
the prediction of surges in Galveston Bay, Texas. The bay is treated
as a nearly closed system which, again, eliminates the problem of
open water, lateral boundary conditions. Measured tidal data at
the entrances to the bay are used as an input to the model and flow
through these areas is allowed. The model includes bottom friction
and allows for overtopping and flooding of adjacent lowlying areas.
Harris and Jelesnianski [11] and Jelesnianski [12,13,14,15]
have presented a series of papers concerning storm surge calculations.
The numerical technique follows that of Platzman. In the paper
"Numerical Computations of Storm Surges with Bottom Stress," Jel
esnianski [14] includes bottom stress in his calculation, using an
eddy viscosity and a slip coefficient. Jelesnianski notes that for
a storm not moving too slowly, the bottom friction term is not im
portant and may be omitted, but for slow storms or storms moving
parallel to the coast bottom friction must be included. In the paper
"Bottom Stress TimeHistory in Linearized Equations of Motion for
Storm Surges," Jelesnianski [15] abandons the slip coefficient and
maintains the eddy viscosity with the assumption that it is constant
in both the vertical and horizontal. Good agreement is obtained
between surge height data and the calculations. He notes, however, that
"currents are only of casual interest compared to the slope or
height variations." In general this is true, yet, it is becoming
apparent that the currents are very important in terms of beach
erosion and forces on various ocean structures. The development of a
more generally applicable dissipation mechanism would decrease the
sensitivity of the current calculations thus decreasing the necessity
for "calibrating" a model which may mask errors and deficiencies in
the model by appropriately selecting bottom friction coefficients.
The volume transport equations are formulated in terms of
flux components and do not describe the velocity distribution near the
bottom, which is related to the bottom friction. Jelesnianski [15]
deals with the problem of bottom friction by providing bottom shear
stress as a convolution integral of the surface shear stress function.
In this technique he linearizes the equation of motion by making the
assumption that the surge height is very small compared to the depth,
n < < h. In doing so he precludes flooding in lowlying areas and
calculations for very shallow waters. In this report it is antici
pated that the modeling of flooding in lowlying areas and surge
calculations in very shallow water will both be desirable, thus a
quadratic bottom friction function which allows the use of the non
linear equations is used. This is equivalent to making the assumption
that the vertical velocity profile is maintained in a near steady
state condition. Future development of the model will rely on measure
ments to provide insight into the processes involved.
Surprisingly little work has been done in measuring hurricane
phenomena. The storm surge data available are primarily from high
water marks. Tidal records are predominantly from enclosed bays and
rivers making comparison to the numerical model difficult. The only
known current data for Hurricane Camille were taken by Murray [16,17]
on the Florida Panhandle. These currents, measured in the surf zone,
were due to waveinduced mass transport as well as wind shear stress and
water slope. Clearly more data of this type are needed.
The difficulty of reaching offshore areas in a hurricane makes
manned observation of these phenomena nearly impossible, so alternate
measures must be found. The appropriate instrumentation must be
either in place continuously waiting for a storm or it must be placed
in the path of the storm prior to the development of adverse sea
conditions in that area. In the first case, the instrumentation sites
must be periodically serviced with the knowledge that the wait for a
hurricane may be a long one. In both cases, there is the problem of
instrument recovery. An inexpensive underwater gage and integrating
current meter have been developed and are described in the report.
CHAPTER 2
STOPRI SURGE CALCULATIONS THEORY
a. Volume Transport Equations
The equations of motion and continuity as used in the develop
ment of the storm surge model are developed qualitatively in this
section with the details of the integration over depth presented in
Appendix A. The coordinate system and nomenclature are presented in
Figure 21. The origin of the coordinate system is arbitrarily chosen
at the shelfbreak and is approximately centered in the region of
interest. The x axis is directed toward shore, y is oriented
approximately parallel to the coast and z is positive upwards. The
NavierStokes or momentum equations for fluid motion are the governing
dynamic equations:
+ u + u wau 2w(sin)v 1 a + 72u (21a)
v+ uv + vv + + 2w(sin)u 2v (21b)
at ax ay az P ay p
aw aw dw )w
+ U% + v~ + =  + V2w g (21c)
at ax ay az p ,z P
where u, v and w are the components of the water velocity in the x, y
and z directions respectively, t is time, p is the water density, p is
pressure, u is molecular viscosity, u is the angular velocity of the
Elemental Water Column
Showing Shear Stresses
and Surface Pressure, P7
Mean Sea Level
(MSL)
FIGURE 2I
PERSPECTIVE OF NEARSHORE
SHOWING COORDINATE SYSTEM
NOMENCLATURE
A z
Y A
Origin
Bottom
ZONE
AND
Earth, } is the latitude, v2 is the Laplacian operator v2 =
dX
a2 a2
+ + and the fluid has been assumed incompressible. The
continuity equation is:
au + av aw (2d)
ax ay az1d)
Equations (2la) through (21c) portray explicitly only the laminar
processes. In this development, the turbulent effects are dominant,
especially the turbulent or Reynolds stresses. Although it is
impossible to describe the turbulent motion in detail, it is useful
to separate the effects of the mean and fluctuating velocity components.
By describing the variables in Equations (21) as the sum of a time
mean value and a fluctuation and then taking a time mean, the turbulent
equations of motion result. Thus setting u = u + u', v = v + v'
w = w + w', and p = p + p' leads to the xequation of motion
au + au 2u au 1 ap + ,!2 a'r+ au'v' + au'w') (22)
+ u + v + w ...u(.+ + (22)
at ax ay az p ax p ax ay az
where
T
UI I u'iu dt
and uj, u2 and u. indicate u', v' and w' respectively. The terms
in parenthesis are the turbulent stresses and T is a period large
compared to the period of a turbulent fluctuation but small compared
to any change in u or v. It is assumed here that terms of the form
7 v2 u are negligible when compared to the turbulent stresses. Thus
the turbulent equations to be used are:
au au + vu au 2w(sin)v B 2 1 xx + xz)
at ax + y az P ax ax 2y a3
av + v V v +p _l a=v aT'
+ + v + w+ 2w(sin,)u 1P a + y v + z) (23b)
at ax 3y 3z p ay P ax ay aZ
a 'I a ' a Iz
aw aw W a w_ 1 1 xz z z
+ u + v + az z g +4 x + + ) (23c)
at ax ay TZ p aZ p ax ay 3Z
where it is now understood that u, v, w, p, T' indicate time mean
xixj
values. Tix the Reynolds stress, is represented by puu.1 u where
;ixj x1J
puiu =PU'.. This is explained in detail in reference [18]. The
I3j 3 1
continuity equation in turbulent form is unchanged
au + av aw (23d)
ax ay az
Frequently a Boussinesq approximation is used for the expression
in parenthesis in Equation (22) which maintains the original form of
Equation (21). The equations remain intractable either in the form
(21)or(23). Since the velocity gradients in the vertical are antici
pated to be, in most cases, an order of magnitude larger than the
velocity gradients in the horizontal, the lateral stresses are
neglected. Thus it is assumed that T' = T' = 0. The normal stress
aT' xy yx
terms of the form xixi are assumed small compared to the surface
axi 1 apn
pressure term  and the hydrostatic pressure. To
further simplify Equations (23) the assumption is made that the
vertical accelerations are small, or as a consequence that the pressure
distribution in the vertical is hydrostatic. Thus,
p(x,y,z,t) = Pr(x,y) + y(n(x,y,t)z)
(24)
where n(x,y,t) is the free surface displacement from mean sea level
(MSL), p (x,y) is the pressure at the free surface, and y = pg is the
specific weight of water.
The Equations (23) are now integrated over the vertical from
h to n where h is the z coordinate of the bottom. This procedure
together with Equation (24) and the free surface and bottom boundary
conditions leads to the vertically integrated or tidal equations of
motion. The details carried out in Appendix A are straightforward
and with approximations to be defined later lead to the following,
using Leibnitz's rule,
qx qx qx + !L sqx D n gD+1 ) (2Sa)
at D ax D ay y p a aX p nx (25a)
x x
qy + x + + 2w(sing)q gD + 1 T ) (25b)
at D ax D aByx p y pny b
aqx anq
x a + 0 (25c)
ax ay at
where the subscripts, n, and b indicate the subscripted quantities are
to be evaluated at the surface and bottom respectively. D, qx and q
are defined as follows:
D(x,y) = h(x,y) + n(x,y,t)
1n(x,y,t)
qx(x,y,t) = u(x,y,t) dz
h(x,y)
rln(x,y,t)
q (x,yt) = v(x,y,t) dz
h(x,y)
It should be noted also that the convective terms of the form
q aq
x as written are only approximate; the actual terms resulting
from the integration are of the form
n
ay J uv dz,
a h
but are written in the differential form in analogy to Equations (21).
The convective terms will be neglected in the further development of
the model until considered again in Appendix B. In summary, the
equations of motion in the final form to be used are:
aq xD apri Dr, 1
aq 2(sin.)qy  gD a + (T T ) (26a)
at y pax ax p x
3qy + 2~(sin.t.)q D a gD 1 + 1 ( Tb ) (26b)
at x 0 ay ;y p n y
The form of the continuity equation is:
aq aq
a ay n 0 (26c)
ax ay at
b. Assumptions and Boundary Conditions
The equations of motion were developed in the preceding
section. To specify a particular problem, various conditions
must be imposed on the boundaries of the model. On the water
surface a forcing function consisting of a normal pressure and a
shear stress is applied. At the bottom, a shear stress or friction
dependent on the water velocity is imposed. In addition conditions
on the surge height and water transport are necessary at the areal
boundaries of the model.
The shear stress at the free surface is calculated by the
method proposed by Van Dorn [19]. The shear stress is based on the
wind velocity and is given as follows:
T p klUlUx
f x
(27)
T = kpklUU
ny
where U = U + U and U and U are the x and y components of the
V y x y
wind velocity, respectively, at a reference elevation of approximately
10 meters. The constant k has been determined by Van Dorn as the
following function of velocity:
k U < Ucr
k = (28)
S+ 2(l Ucr/U)2 ; U Ucr
where
kI = 1.1 x 106
k2 = 2.5 x 106
Ucr = 14 knots = 23.6 feet/second
There are uncertainties in applying Van Dorn's results to hurricane
winds. Van Dorn's data are generally for wind speeds low in comparison
to the wind speeds of a large hurricane. Although there is no proof
that the results are directly applicable, they are used here because
the associated measurements showed little scatter and the results
are reasonably consistent with other data sources. Wu [20] has pub
lished more recent information concerning the wind stress coefficient,
k. Figure 22 summarizes available data. Van Dorn's k, calculated
from Equation (23), is also plotted in Figure 22 to provide a com
parison. Since the hurricane model to be discussed requires a shear
stress field as input, a better shear stress model, if it becomes
available, can be easily incorporated into the surge model without
changes in concept or approach. The wind velocity input will be
discussed in Section 2c.
The treatment of the bottom friction, rb Tb is a much
x y
more difficult problem. In making storm surge calculations the
simplest assumption is to consider b = Tb = o. This is employed
x y
in many developments of storm surge calculations, especially the
earlier ones. Under appropriate conditions this is a justifiable
assumption. For the case of a rapidly moving storm traveling per
pendicular or nearly perpendicular to the shore it is a good assump
.003
0
e o.oe 0
/
/
.002 0
0 oLegend
/ Vo n Dorn's k
S.001 Oceanic Measurements
S0
) 0 Note:
.c Oceanic Measurements
From Reference 20]
S Ii
5
Wind
10
Velocity,
15 20 25
U, (Meters/Second)
FIGURE 22
COMPARISON OF VAN
WITH OCEANIC WIND
MEASUREMENTS
DORN'S k
STRESS
tion since there is little time for "windup" and the friction
effects are small due to low water velocities [13].
A better assumption in some cases is that of a linear or
linearizedd" bottom stress of the form fLqx, f qy where fL is a
linearized friction factor. This is the type of friction used, for
example, by Dean and Pearce [21]. A linearizedd" friction term is
generally necessary if friction is required and any type of analytical
solution is desired. The bottom shear stress employed in this model
is the DarcyWeisbach expression and is the most accurate of the
three,
pflulux pflulu
b 8 Tb 8
x y
or in terms of volume transport
Pflqlqx Pfjqlqy
Tb = 8D ; Tb y 8DP
x y
where u = JUx2 + Uy2 q = /q2 + q 2 and f is the quadratic
bottom friction factor. In incorporating this steady state rela
tionship into the model it is assumed that the law is valid for the
propagation of long waves of the period of interest. This is
equivalent to assuming that the velocity variation over depth is
nearly the same for slowly varying flows as for steady state.
Although this technique has limitations, especially in deep water,
it is the best method available that can be used in conjunction with
the nonlinear equations to allow for flooding and to permit cal
culations in very shallow water. Another important consideration is
that in deep water the friction is relatively less important than in
shallow water since the water velocities are in general much smaller.
As an illustrative example, the response of a water column to an
impulsive wind is depicted in Figure 23.
Returning to the problem specification, the boundary conditions
are illustrated in Figure 35. At the seaward boundary, the wind
setup is assumed zero as a result of considering an infinite depth
at that point. The barometric tide still exists at the seaward
boundary and is given by n(x,y) = 1.15(p p (x,y)) where pc designates
the barometric pressure at a distance large compared to the storm
dimension, where r, is in feet, and the pressure p is in inches of
mercury.
At any landsea boundary the flux normal to that boundary is
zero, q = 0, where the subscript n indicates the component normal
to that boundary.
Most difficult from a physical standpoint are the lateral
boundary conditions. They are, however, necessary to carry out the
calculations. Since no real physical limitations exist, conditions
must be chosen that least disturb the model and therefore allow a
realistic solution to be obtained. A further condition is then
imposed that the boundaries are at a distance large compared to a
characteristic dimension of the storm. This means, in effect, that
the boundary conditions on the lateral boundaries become a second
order effect if they are far away. The boundary conditions used at
the lateral boundaries of this model are: the flow in the x direction
is equal to zero, qx = 0, and the surge height n is set equal to the
barometric tide at the lateral boundaries.
18
Direction of Wind Velocity 
7.
_!1  _ __J
 .
it i"7 2"f<''
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
I =0
I5
10
20
30
40
50
Notes:
Wind
Depth
0
Woter Velocity,
FIGURE 23
sec.
sec.
sec.
sec.
sec.
sec.
sec.
eociy 0; t<0
velocity 9.8 ft/sec, t 0 
= 0.48 Feet
0.1
u, Feet/Second
0.2
TRANSIENT WATER VELOCITY
DISTRIBUTION IN A CLOSED CHANNEL
DUE TO IMPULSIVE WIND (FROM
REFERENCE [22])
I .'I r
(
I
I
I
1
I
.1
A
:I
I
i
I
I
I
ii: '
Y
01
0.
I
The initial conditions are simply a quiescent system. The surge
heights can be initially zero or in barometric equilibrium: this makes
little difference since the storm is initially far from the region
of interest. Initial barometric equilibrium is used in this study.
More importantly the storm is assumed to start at a distance offshore
large enough that the effects due to model startup are negligible.
The effects of various starting points can be seen in Figures 45 and
410, where the surge height is relatively insensitive but the currents
exhibit windup.
c. The Hurricane Model
The actual storm surge model uses T T and p as input.
x y
Since the measured wind shear stress is in general not available,
Van Dorn's technique is used to convert U and U to T and T as
x y
discussed in Section 2b. Actual wind fields may be applied or an
idealized model used to represent the storm. In this section an
analytical model is used to represent wind velocity and pressure fields.
Since the actual surge model is the primary objective of this study,
the hurricane model will be presented and used without further comment,
except for a short note on the numerical technique in Chapter 3.
The hurricane representation used is not unique and a different
one could be substituted into the numerical model. The hurricane
model used is that presented by Wilson [23]. The important parameters
of an idealized hurricane are: the central pressure index CPI or
equivalently the central pressure anomaly Ap = pm p0 (where pm is
the normal pressure or the pressure at a large distance from the
hurricane center and po is the central pressure of the hurricane), the
distance from the hurricane center to the point of maximum winds R,
and the forward velocity V of the hurricane. These parameters Ap,
R, and V are reported by the National Weather Service when character
izing a storm, thereby providing a method to approximate a hurricane.
The isobars are assumed circular and the pressure distribution at
the surface is taken as pR/r
the surface is taken as p = p + Apep The streamlines are
n o
assumed circular for a stationary storm. The equations presented
include the effect of a constant translational velocity, V,
superimposed on the storm as shown in Figure 24. Using Wilson's
notation:
UG = Uc( J2 + 1 
where
V U
Y 1/2( + U)
c g
and
Uc = Ap R eR/r
P r
Pa
Ap R R/r
U =
g 2w(sinT)
V' = V(sine)
The subscript G refers to gradient wind, c to cyclostrophic
wind, and g to geostrophic wind. pa is the density of air, V the
forward velocity of the storm, e is defined in Figure 24 and 2w(sing)
is the Coriolis parameter as used in Equations (21). The gradient
wind UG is the desired parameter. These winds are calculated at a
point high enough above the surface to allow the assumption of zero
friction. To reduce these to the 10 meter elevation a reduction factor
is needed. The factor depends on latitude, divergence of wind
direction from the isobars and other factors. In this study a
constant factor of 0.83 was used. See references [24,25].
UG = Gradient Wind Velocity
SV = Forward Velocity
Storm Center o Stor
of Storm
FIGURE 24 DEFINITION SKETCH FOR STORM TRAVELING AT CONSTANT VELOCITY
The winds of Hurricane Camille are evaluated using this method.
The hurricane position as a function of time is given in Figure 25.
The central pressure p and the peripheral pressure pm, are given as
functions of time in Figure 26a. Figure 26b depicts radius to
maximum winds R, and forward speed of the storm V, as functions of
time. Figure 26c gives maximum sustained wind speeds in knots as
a function of time. These data are from the National Weather Service
[26,27]. Figures 27a through 27c present comparisons, at different
times, of the surface wind fields for Hurricane Camille as developed
by the National Weather Service [27] with the wind fields calculated
by the above model.
Notes: I
Times are GMT
Dotes ore in A
FIGURE 25 STORM TRACK FOR HURRICANE
CAMILLE, AUGUST 1618, 1969
Legend:
940
o
S920
90
0900
o
c
(
. . 0  0 
0600
1200 1800
Time (GMT)
FIGURE 260
PRESSURE CHARACTERISTICS OF
HURRICANE CAMILLE AS A
FUNCTION OF TIME
(DATA FROM REFERENCE [26])
I 20 
2 odius to Maximum Winds
E _, 0odiu 0 Mx u m Wiuds
Er
o C I0. 
"^ i Forward Velocity 15
*0 / 5 /
Y
o / \ 
10 10
0000 0600 1200 1800 0000 0600
Aug 17 Time (GMT) Aug 18
FIGURE 26b RADIUS TO MAXIMUM WINDS AND
FORWARD SPEED OF HURRICANE
CAMILLE AS FUNCTIONS OF TIME
(DATA FROM REFERENCE [26])
 Centrol Pressure
00 Peripheral Pressure
Nole:
Storr Landfall at 0430 GMT, Aug 18
_ ,
oou
0000
Aug 17
)
0000
Aug 18
1010
0600
L.
o
030
Q)
,
L
020a
I
0)
1l
o)
S125
U,
'o
c
Aug 171
E 120 2
D
E
0
115
0000
Aug 17
FIGURE 26c
0600 1200 1800 OC
Time (GMT) A
MAXIMUM WIND SPEED IN
CAMILLE AS A FUNCTION
(DATA FROM REFERENCE
)00 0600
ug 18
HURRICANE
OF TIME
[26])
cn
'.0
4
0
cn
C>
C)
cJ1
4
U
L
0
0
0
C.J
0
rL.
O.J
 O
C)
I i4
J
.I U

 ,a
0
1 u
Ue
>
> rn
27
cn
I>
"/ "
o . co
z CD
01
m 0O
C) CO
if
C,
OL
II /2
t/ C
.O / '
, t,,, l ,,
0 (D
\ 0"O"\ o o0 /
C O
co o L)
n,,
00
9 I / \\ (
r o\ co
s.
In \0 \ 1 \*
0 0 0
0 I I.
0 0 0 CL..
\ 'o^ \fe\   ^
mL L
0 0
0 0
o II
LI L1
I
e '
EJ I
00
C!
d. Analytical TreatmentSteady State, Without Bottom Friction
i. Problem Definition
Having defined the problem and established equations of motion
and boundary conditions, it remains to find a solution for a particular
storm and affected coastal region. Before presenting the numerical
approach necessary to attain a solution of complicated prototype con
ditions and geometry, several simplified analytical approaches will be
considered. These will serve three basic purposes. First, they
allow highly simplified calculations for basic engineering purposes.
Second, they allow intuitive investigations into the physical processes
involved. Third, they produce solutions with which to check the
numerical approaches.
Figure 28 is a definition sketch for this section. The present
treatment includes only motions in the xdirection; Equations (25a)
and (25c) therefore reduce to Equations (29a) and (29b)
respectively.
aq q aq D ap
ax +x x D gD a n + 1 ( ) (29a)
5at D x p xx ;x p n b
x x
x +an 0 (29b)
ax at
Bottom
FIGURE 28 NOTATION FOR ONEDIMENSIONAL PROBLEM
The solution is now independent of y. Two cases will be considered,
the uniform depth case h(x) = h = constant and the case of h(x) = h
mx or uniform slope. In both cases it will be assumed that the wind
stress on the surface is uniform. For the steady state case, the
case of no net transport in the xdirection, q = 0 and therefore the
shear stress ,b = 0, regardless of f, the bottom friction coefficient.
v.
(Note: it is assumed here that there is zero velocity over the entire
q aq
depth.) For this problem the D (x) term is zero and the atmospheric
D dX
pressure gradient  is assumed small in relation to the other terms.
ax
ii. Uniform wind stress Acting Over a Shelf of Uniform Depth
In the first example h(x) is uniform (=h) as shown in Figure
29.
z
9^  1L :,7 P. X
Origin
h x= x
CD
FIGURE 29 ONEDIMENSIONAL SHELF WITH UNIFORM DEPTH
The further assumption is made that the shear stress is uniform and
constant at the water surface z = n. The problem can be stated as
follows:
T
a2n x
gD  + 0
ax p
aq
qx at 0
Since D = h + n
pg(h+n) an T (210)
ax n
Rewriting Equation (210) leads to
T
aX 7 x (211)
ax pg(h+[r)
At x 0, the boundary condition for r is determined by noting that
seaward of that point h m and therefore
an 0 as h ,
ax
Thus n = constant as h , <, and n = 0 is chosen as the boundary con
dition at the shelf break, x = 0. Solving for n yields
T
r(h + n/2) = x
pg
or
n = h + h2 + 2x
pg
Evaluating r using Van Dorn's constants, completes the solution.
nX
An even simpler procedure is frequently the basis for storm surge
calculations. Assuming that the storm surge, n, is small compared to
the undisturbed depth (i.e. h + n , h) leads to the solution
T
nx
pg
Designating xc as the fetch or shelf width and h as depth leads to
the storm surge at the coastline x = xc
T
nx Fetch
x=Xc pg Depth
Figure 210 provides a simplified procedure for making these
calculations and is based on the latter solution. As an example,
assume a depth of 50 feet for a 50 mile shelf width (or fetch) and
a wind speed of 50 knots, then the storm surge calculated from
Figure 210 is 3.23 feet.
s
r. 0
0
LT
(\J
O
O
(1aj) L '10SP!H a6jnS
34
iii. Uniform Wind Stress Acting Over a Shelf of Uniform Slope
As a slightly more realistic situation, the same assumptions
are made as in the previous case except now a sloping bottom
(or h(x) = ho mx) is considered as shown in Figure 211.
T"q
77x
X/ I
xKxK
FIGURE 211
ONEDIMENSIONAL SHELF WITH UNIFORMLY SLOPING BOTTOM
From Equation (211)
(h + rj) ar x constant =
ax og
and noting that h(x) = h mx
0
aD (h+n) an + (h+) (h+n) m (h+n)
ax ax ax ax
DdD
thus dx = Timd
at x = 0 to obtain
Integrate and apply the boundary condition n = 0
1 4
= ln (1  )
(212)
where
mn
mh
mh
mhi
=o
To simplify the calculations obtain n at the coastline x = xc
Equation (212) then reduces to
1 
S= In ( ) (213)
Solutions to Equation (213) at x = xc are provided in Figure 212.
As an example consider a shelf with xc = 50 nautical miles, ho = 80
feet, and hI = 20 feet, this shelf has the same width and the same
average depth as in case ii. For the same 50 knot wind, n at the coast
is 3.53 feet. Heuristically, one would expect a larger n for the
case of the sloping shelf.
Consider the location on the sloping shelf which has the same
depth as the uniform depth shelf. For a point x towards shore
there is a change in depth ah and conversely for Ax away from shore.
Thus
+ constant constant 'h h 2 2
Ar, [  + ( 0 (Ah)3]
constant constant [ +h Ah2h)3]
n h.h. (F) + 0 (Ah) ]
hF fah FI
w
I
S0
N ,0 o
0\ 0
\ \n
0
C\%J 0
C)
"J
0
LLi
)\ (/)
z
S I0
5 X \ \
JD
0 \
4 CL
L
37
and the total over 2Ax is
n = An + An = constant [2 + 2(Ah + 0 (Ah)4]
It can be seen, then, that a shelf of uniform depth will, in all
cases, produce a smaller surge than a shelf of constant slope with
the same average depth.
e. Analytical Treatment Moving Wind Stress and Pressure Field
i. Problem Formulation
The purpose of this section is to contribute to the development of
a simple procedure for the analysis of the dynamic effects involved
in storm surge calculations. This analysis uses the onedimensional
tidal equations and the development will be for an arbitrary forcing
function such as a hurricane. Also the bottom profile must be in an
analytical form simple enough to allow a solution to the equations of
motion. Conceptually, the approach is as follows: find the storm
surge n(x,t) resulting from an impulsive forcing function representing
the wind shear stress, pressure gradient or both. If the impulse
occurred at time t' and position x', the resulting surge at (x,t)
could be written n(x,t) = G(x,t; x',t'); this is sometimes referred
to as the Green's function. Having found the general form of the
solution G(x,t; x',t') for a simple offshore bathymetry one can then
construct the solution for an arbitrary forcing function. This may
be obtained by formal integration or the results may be obtained
numerically.
The development in this section will be for the special case of
a shelf of finite width and uniform depth as shown in Figure 213
and a triangular shear stress distribution as shown in Figure 215
For the equations of motion and continuity, consider Equations
(25) where qx is abbreviated as q, nx as in, and bx as Tb,
\ x
39
9 + gD an D aP b) (214a)
at D ax ax pax p (214a)
S+ = 0 (214b)
The convective acceleration terms will be neglected and the assumption
n < < h is made, and since D(x,t) = h + n(x,t) = h, Equations (214a)
and (214b) reduce to
aq(x,t) + gh an(x,t) h apn 1
at + gh ,ax p (x,t) + (n(x,t) Tb(x,t))
aq(x,t) + an(x,t) = 0
ax at
Ignoring the bottom friction term Tb and letting (1T (x,t) 
h ap (x,t) n
 ax ) = F'(x,t) we have:
p ax
a+ gh F= F'(x,t) (215a)
aqt an
ax a 0 (215b)
The variables are described in Figure 213, where q is the volume
transport per unit width. Differentiating Equation (215a) with
respect to x and Equation (215b) with respect to t and subtracting,
produces the inhomogeneous wave equation
2n gh 2r aF'(,t) F(x,t) (216)
Tt= TX 7 x
T (x, t)
7
Origin ]
h
q(x,t)
SII== Il l_ i _ _
CO
FIGURE 213 IDEALIZED SHELF BATHYMETRY
ii. Solution
First a solution to the homogeneous problem F(x,t) = 0 is
obtained. Using separation of variables, let n = X(x)T(t)
For the spacial part of the problem, we obtain
gh = a2X (217)
where an2 is a constant as yet undetermined. The solution to
Equation (217) is
X = cl sin  x + c cos x (218)
Referring to Figure 213 note that the depth is infinite at x = 0,
and the surge must be zero at that point. At x = z, there can be no
flow through the wall and this leads to the condition = 0 at
x = 9..
The condition n = 0 at x = 0 determines
C2 0 and I lan (cos n ) T(t) = 0
ax x=Z fgh 'gh
leading to = (2n + 1) The general case for the solution
becomes a sum of the solutions in Equation (218). Thus the solution
is taken to be in the form
n E a (t) sin nx (219)
n=1,3,5"*
Referring to the homogeneous form of Equation (216), obtain with
Equation (219) the following:
32a (t) nx C2 n2 a innj =
S {at sin + C2 ) a (t) sin 2 = 0
n=1,3,5 
From orthogonality, obtain
a "(t) + n2 a (t) = 0
where on = 1 and c = g
The two homogeneous solutions are of the form a (t) = sin ant,
Cnj
cos ont. Writing F(x,t) = Z B (t) sin " and B (t) =
n n=1,3,5 n
SF(x,t) sin 2x dx for odd n, and using the method of variation of
parameters, arrive at the following for the inhomogeneous equation;
with the lower limit determined by the fact that a (0) = 0.
sin oat o t c)s o t C(t
an (t) B=(t') cosant'dt' n n B(t )sinot'dt'
0 0
t B (t )
a (t) = B sinan(tt')dt'
r a(t) sin n = B (t')sina (tt')dt'sin 2
n= ,3,5 n=1,3,5 n 0
and finally
n(x,t) = Z F(x',t')sino (tt')sin n dt'dx'sin (220)
n=1,3,5 n. O 0
Letting F(x',t') = ,6(x'x ) 6(t't ) one then has an impulse function
located at x = x and t = to where n represents the Green's function
or response to this impulse function, and : is a constant representing
the strength of the instantaneous surface shear stress function
n(x,t) = G(x,t; Xo ,t) =
C O. 2 z tx n x
n=l
n= 1,3,5 0
G(x,t;x ,to) = ? sine, (tto) ssin sin x (221)
S n=1,3,5 n
It is now possible to construct a solution by the superposition of an
appropriate series of impulses in space and time.
For some hurricane motions, a more convenient form of the Green's
function is obtained by using an impulse function which travels with
constant velocity and remains unchanged during its traverse over the
43
shelf. Assuming the forcing function to always start at x = 0, to
travel at constant velocity C, and that the time that it crosses
x = 0 is arbitrary at to, the impulse function is:
F(x,t) = K<(x C(tt )) 0 < x < z
Without loss of generality let t = t + to which yields
F(x,t) = K6(x Ct) 0 < x < R
Thus for the case of a suddenly increased shear stress as shown in
T
Figure 214, K where T is the amount of the increase in
p 0
shear stress.
P
C0
P
C
x
I ar r 8 (xCt) = F(x, t)
P ax P
0C
Sx
FIGURE 214 IMPULSIVE FORCING FUNCTION MOVING AT CONSTANT SPEED
For the case of an arbitrary shear stress distribution as shown in
Figure 215a, locally LT = x0
ax
FIGURE 215a
ARBITRARY FORCING FUNCTION
For a triangular shear stress distribution as shown in Figure 215b
I a.
P ad
+O
r I
FIGURE 215b TRIANGULAR FORCING FUNCTION
I~
>C
;
.r"'/p
KO 
1 r
r: becomes  Ux, and the following is obtained.
p dX
F(x',t') = 1  %(x'Ct') dx = AT 6(x'Ct')
p dX p
The resulting Green's Function can then be written in integral form as
follows.
J 2 ar j Ct't' nx' n7Tx
G(x,t;O,O)= (x'Ct)sin(tt)sindt'dx'sinx
n=1,3,5 n p 0n
Note that this is really the solution for an incremental element having
a shear stress jump of AT. Since 6(x'Ct') 1 0 only when x' = Ct'
G(x,t;0) =
t
2AT nnCt' nux
E 2A sinn(tt')sin 2 dt' sin 2
n=1,3,5 Pn 0
Letting an = = nV n = o
and integrating yields
2AT t n rx
G(x,t;O) = po(Z [o' sinot a sin'ot] sin 2x
n=1,3,5
Note that this is the solution for a solitary continuous impulse
function traveling at speed C; at time t > */C the forcing function
crosses the shelf, and the solution becomes periodic in time.
Letting c* = C/c, c = /gh and o' = o c* one arrives at
G(x,t;) = +2T [c*sinatsin(ac*t)] sin nx
n=1,3,5 pz(c*zl
with the constraints t < r/C and c* / 1. When t > z/C
G(x,t;0) =
2 fAT k/C
n=1,3,5 p
o
sina(tt') sin n "Ct' d sin n7x
22 21
where a' i/C = nT
2
n7T
cos = 0
nl
sin = (1) 2
2
G(x,t; 0) thus becomes
G(x,t;0) =
S 2Airsin(nnx/2)[+cs t
n=I*3,5 )P(c*) +c*sn
n=1,3,5
n1
 coso(ti/C) (1) 2 ]
a'
with the constraints t > z/C and c*2 / I. If AT is written as r dx'
dX
and it is noted that x = Ct then AT = C 1dt. Letting t = t t
yields G(x,t; t ) where to represents the time at which *(xC(tt ))
crosses the point x = 0. The desired result is then the sum or integral
of the above solutions
n(x,t) = I
G(x,t; t ) dt
The shear stress profile of a hurricane can be approximated by a
triangular function [2]. The derivative of this is two rectangular
pulses as indicated in Figures 215b. The response can be found
using Equation (222) and remembering that G(x,t; t ) is represented
by two functions depending on whether t is less than or greater than
E/C. The solution is straightforward but unfortunately tedious.
n = 1,3,5
(222)
The solution is simplified by putting n into a dimensionless form.
Letting
T = 4k/c
= aT = 2nTr
X* = X/Z
t* = t/T
allows n' to be represented as Equation (223).
l'(x*,t*) = n(xt) (2n)3
(2K 'T3/z)
(223)
The appropriate variables and constants for the triangular stress
distribution are illustrated in Figure 216, where tI has been
arbitrarily set equal to zero for convenience,
C aT
P ax
Ko
K0
tI
. X
FIGURE 216
DEFIINITIO;, SKETCH FOR SHEAR STRESS DISTRIBUTION
and are defined as follows:
C 3a
K ax
p ax
48
ti = 0 tl* = 0
t2 2C
St2 8c*
t3 C 3* 4c*
X+
K = r ;
I =
o
0 < t < 
T < t
Co 2C
2C 'o < C
<' = ; < t
C 0
The solution must be presented as six different forms, each for a
different time span. They are as follows:
(1) In the first case the storm has completely passed by the
position, x = , thus
,1+e **+l
t > C t* >
and the solution for n is
n'(x*,t*) =
nl to=t3
sin(n7x/29) [sin((tt )./C) () 2 +c*coso (tt)]
n=1,3,5 ) n
By defining the contents of the brackets as H (to), the response becomes
n'(x*,t*) = r sin(n x/2) [H 3) H ) 2H )
n=l1,3,5 nF (c l) [HI(t3) + Hl(tl) 2HI(t2)1
where tl, t2, t3 represent values of to and are defined in Figure 216.
(2) In case 2 the storm is completely on the shelf 0 < x < , thus
A/C < t < R/C T I
and n can be written as follows
n'(x*,t*) =
t=t3
sin(ncx/2) [+c*cos( n(tto)) L cos(c*.n(tto))]
n=1,3,5
t=tl
Defining H2(to) as the contents of the brackets, as above, leads to
n'(x*,t*) = sin(nnx/2) [H
n=,3,5 nJ(c*zl) [H2(t3)+H2(tl) 2H2(t2)]
n=1,3,5 nc 1
The time periods from 0 < t < A/C and ,/C < t < have not yet been
covered since they introduce the additional complication of having
only part of the forcing function in the region of interest. The
solutions for these remaining time periods are:
(3)
z + A/2 < + A 1 + \*/2 1 + I *
C t C 4c* < < 4c*
n' (x*,t*) =
sin(nnx/22)
z n(nc ) [H2(t3)H2(t/C)+H1(t/C)2H (t2)+H (tl)
n=1,3,5
(4)
2 + A/2 1 1 + .0/2
r/C < t < ? + 2 1< <1 + /2
t C c
sin(nux/2z) [H (t)
n in x/2)l [H2(t3)2H2(t2) + H2(tt/C)Hl(tz/C)+H(tl)]
n=1,3,5,...
A < t <  < t*
2C c 8c < 4c*
0 < t < X/2C 0 < t* < 8*
8c*
ri'(x*,t*) n=
n=1,3,5,...
sin(nnx/2.) H
n (c* 1) [2t ) H2(t)
n'(x*,t*) has now been described for all t*. It can be seen that n*
can be described as a function of three dimensionless quantities
c*, X*, and t*.
Figure 217 shows profiles of n' on the shelf for various values
of t* and particular values of ,\ and c*.
Figure 218, perhaps of more interest, shows n'ax, the
maximum storm surge that occurs as a result of the storm passage as a
function of both c* and \*.
To determine the accuracy of both theory and calculations, an
independent approach was utilized to obtain equivalent dimensionless
surge heights. The onedimensional tidal equation of motion and
equation of continuity in finite difference form, as developed in
sin(nrx/2g.)
n'(x*,t*) = n(c ) H2(t)+H2(tl)2H2(t2)]
n= I,3,5,... *
0
0 0 z
0
O ()
co
rr
Q)
0
r W
S\ / '
0
0
LLL
,o,)
0
z D
00
*o Xo &
uIn
z
O < L
d 0 .0
o 0 6
11 0
u0 0
0
soun AJDJliqj 't 'Iqb!lH binS i ~
.5
4
.3
2
o
.2 .4
FIGURE 218
DIMENSIONLESS
MAXIMA AS A
DiMENSIONLESS
VELOCITY
2 1.4 1.6 18 20
SURGE HEIGHT
FUNCTION OF
STORM WIDTH AND
Section 3a, were used to obtain a numerical solution to this problem.
The equations were linearized (D h) to agree with the assumptions
in this section, and the bottom friction neglected. The "no flow"
boundary condition qx = 0 was prescribed at the coast with a boundary
condition r, = 0 assumed at the shelf break. Note that the qx = 0
condition is comparable to = 0. Figure 219a depicts a shelf of
constant depth divided into 20 onedimensional grid elements, as used
for the calculations. The "storm" or triangular shear stress dis
tribution as illustrated in Figure 219b was moved over the shelf
waters at a constant speed, C. A slight approximation was used here.
Once the edge of a shear stress step encounters a new grid element,
the stress was assumed constant over that element for a time At = Ax/C,
or the time for the leading edge to traverse the grid element. The
numerical results were put into dimensionless form as in Equation (223)
and compared to the analytical results yielding Figure 220.
54
Notes:
(I) r Lineorly Exlropoloted to Coast
(2) Boundary Conditions
(o) q =0 at Coast
(b) 77 0 in Grid (I)
   l.,!
h I,
I _
,Shear Stress Distribution
Moves at Constant
i/ \ Velocity C.
C z
77 =. 0
I i I I i I I
I I
. I I I
20 19 18 17 16 15 14 13 12 II 10 9 8 7
20 Grid Elements
x
FIGURE 219a DEFINITION SKETCH
CALCULATIONS
To
'/2 To
'/6 r
I I
I I I I
I I I i
S I I
6 5 4 3 2 1=I
 [VIL
Ax~
T h* CO
FOR NUMERICAL
FOR NUMERICAL
Idealized Distribution
. I Discrete Approximation
/ I I I
"
........ l
./ "I i I
/III' x
FIGURE 219b FINITE DIFFERENCE MODEL OF
TRIANGULAR SHEAR STRESS DISTRIBUTION
_ __
I
.I f I I I I I I I I I
.1 .2 3 .4 5 .6 .7 8 9 1.0
r (Analytic)
FIGURE 220 COMPARISON OF NUMERICAL AND
ANALYTICAL SOLUTIONS
CHAPTER 3
THE NUMERICAL MODEL
a. Onedimensional Numerical Model
In most prototype situations it is unrealistic to make the simpli
fying assumptions necessary for a successful analytical treatment of the
storm surge problem. To incorporate such features as an irregular coast
line, irregular bathymetry, islands, and arbitrary wind stress patterns,
a numerical approach becomes necessary. The actual numerical approach or
algorithm is frequently referred to as a numerical model. Models may be
roughly categorized as one, two, and threedimensional. One and two
dimensional models are presented in this chapter.
The onedimensional model considers motion in only one direction,
e.g. along a perpendicular to the coastline. It is presented here for two
reasons: its simplicity allows physical insight into the problem and it
provides a check of the twodimensional model, which can approach a
onedimensional model as a limiting condition.
The development of the onedimensional numerical model begins by
referring to Equations 26.
at 20(sing)qy  gD an + ) (26a)
at P *x yx p ny b)
x. x
y + 2w(sing)qx = D r gD n+ ( ) (26b)
t P 3y 6y 0 n y b y
aqx q n
l+ + (26c)
ax ay at
Considering only motion in the xdirection leads to Equations (3la) and
(31b)
S_ gD a b ) (3la)
at pX x p xI b
aqx +n
X +
xa +t = 0 (3lb)
ax at
pfqxlql
Using the quadratic friction law, Tb 8D Equations
x
(3la) and (3lb) can be written as
aq x D n n TnX fqx qx
 gD 3n + ~ (32a)
at p aX aX p 8D
aq
 = 0 (32b)
Dx at
These are the equations to be used for the onedimensional calculations.
The finite difference integration scheme used for the numerical
calculations is similar to that of Reid and Bodine [10] and Verma and
Dean [28]. Figure 31 provides a schematic diagram of a typical coastal
region showing the coordinate system and defining the variables in use.
The vertical coordinate, z, is oriented positively upward and referenced
to mean sea level, MISL; the horizontal coordinate x is oriented toward
shore with x = 0 occurring at the shelf break. In the computations n,
h, T, and p are considered uniform over the grid width LN,and centered
at the midpoint of the segment. The quantity, qx, denotes the flux in
volume per unit width through the plane representing the seaward bound
ary of a grid element.
C
I C
O
z
0
UI
O
z
J
LL.
Y
0
I 
0
z
U
wo
0B
..........  
For clarity the variables are placed into an indicial form where
qx(i+l,n) indicates the volume transport per unit width through the
plane section at the location x = i Lx and at time nt, i.e. flow into
the grid designated by index i+l. Similarly, rj(i,n+1/2) indicates the
value of n for the i grid element. It is assumed to be centered at
x = (i+1/2)Ax and to be evaluated at time t = (n+1/2)At.
Following the scheme in Appendix A, the equations are placed in
finite difference form as Equations (33a) and (33b)
,n l) ,n) + t [D(i n+ /2)+D(i,n+1/2)
q (i,n+1) [qx 1,n) + at { 2,x
x
{[n(il,n+1/2) r(i,n+l/2)] g + [p ((il,n+1/2) p (i,n+1/2)]}
+ at [T (il,n+1/2) + r (i,n+1/2)]] (33a)
2p r1
Sx x
where
Fx= 1 t f
F = 1+ 2 Iqx(i,n)I[D(il,n+1/2)+D(i,n+l/2)]2
The continuity Equation (32b) provides, in finite difference form,
the equation for the free surface,
n(i,n+3/2) = n(i,n+1/2) + ' [qx(i,n+l) qx(i+l,n+l)] (33b)
To illustrate the calculation procedure assume that all values of
qx are available up to and including t = nAt. Beginning at the
seaward boundary, calculate qx's at time (n+l)at using Equations (33a).
Next, calculate n's at (n+3/2)at using Equation (33b), then advance the
pressure and shear stress to (n+3/2)At and calculate qx's at time
(n42)At, and so forth. In effect, averages over small segments are being
made in both space and time to approximate the true solution. This re
sults in a type of leapfrog method, since qx's at time (n+l)At are used
to extend the solution from n(i,n+1/2) to n(i,n+3/2) and similarly
n and D at (n+3/2) are used to obtain qx(i,n+2). The process is con
tinued until a desired time is reached, such as a surge maximum at the
coast.
The seaward boundary condition, the barometric tide, is applied
to grid (1,n), where the wind setup is zero and h(l,n) o. This con
dition is written as n(l,n) = 1.15 [p. p (l,n)] where n is in feet and
p. is the undisturbed pressure in inches of mercury far from the storm.
The landward boundary condition of no flow requires qx(ic,n) = 0. The
initial conditions used will require all qx's = 0 and the grid to be in
barometric equilibrium. The only remaining data required for the calcu
lations are the pressure and shear stress distributions as functions of
time. Figure 32 shows the movement of such a system at constant
velocity and with constant form. The required wind and pressure data are
extrapolated linearly from the storm model described in Section 2c.
The model also allows for the "flooding" or "draining" of a grid
element. Consider the element (ic,n) in Figure 31, if n(icl,n)
<h(ic) then the system is as shown, if n(icl,n) > h(ic) then the grid
element (ic) must be flooded or activated. In practice when activating
a grid element an instability can occur which results in switching be
tween the activated and inactivated states. Also remember that at the
moment of grid activation D = 0 causing trouble with the calculations at
that point. To avoid these problems, the following steps are taken.
L." .
'Q~ ) (If
) 0. ) L
IJ)
Once activated the grid element is constrained to remain activated fo.r a
specified number of time steps. Five time steps has been found satis
factory in the present model. In addition to the above, D is set equal
to .01 ft. in the justflooded grid element prior to the calculations for
the next time step to prevent zero divide problems. This introduces a
small error into the system. Since flooding occurs infrequently with
respect to the number of time steps this error is considered acceptable.
To calculate a water velocity u from qx an average depth over ad
jacent grids is taken and u(i,n) is approximated as
I
u(i,n) = 2qx(i,n) [D(il,n+1/2) + D(i,n+1/2)] (34)
Since u(i,n) is not used in the calculations, the approximations in
Equation (34) are not considered significant. A more accurate value,
centered at time n, which must be calculated one time step after the
qx's is given by Equation (35).
u(i,n) = 4qx(i,n) [D(il,n1/2) + D(il,n+1/2) + D(i,n1/2)
1
+ D(i,n+1/2)] (35)
It is assumed that the hurricane or pressure and shear stress dis
tributions approach normal to the coast. Numerical stability for the
nt 1/2
onedimensional model requires that < (gD ) where D is the
Ax max max
maximum total depth expected anywhere in the system. Thus given a bottom
profile and, shear stress and pressure, as a function of time, the surge
can be calculated using this reasonably simple technique. This model is
particularly useful in quantifying effects of the numerical model that are
not inherent in the equations of motion but result as an effect of dis
63
cretization. The model is also easily checked and thus provides a basis
for comparison with other models which are not easily checked.
b. Development in Two Dimensions
The development in two dimensions is analogous to that
dimension and therefore much of the detail will be omitted.
that the twodimensional model considers motion in the x and
refer to Equation (26), and employing the quadratic bottom
_q_ an D P +qq
x 2mw(sin,)q + gD n D ;P x f
t y X P X p CD7
aq dp ri fjqq
y + 2w(sinP)qx + gD r D '1 y
at x dy p y P 8D
dq aq
x + V. ; = 0
;X Dy at
in one
Recalling
y directions,
friction obtain
(36a)
(36b)
(36c)
These equations are now in forms which can be applied and where no
other approximations need be made. Expressing Equation (36a) in finite
difference form as described in Appendix B leads to the following for the
xequation.
qx(i,j,n+l) = [qx(ijn) + 2w(sino)qy (i,j,n)&t
+ At { (iI,j,n+1/2) + T (i,j,n+1/2)}
x x
+ {D(i1,j,n+1/2) + D(i,j,n+1/2)) { [pn(i ,j,n+1/2)
P (i,j,n+1/2)] + g [n(i1,j ,n+1/2) n(i,j,n+1/2)]}] (371)
where
F = 1 + { [qx(ij,n]2 + [qy(i,j,n)]2}1/2
[D(il,j,n+1/2) + D(i,j,n+1/2]2 (37b)
Equations (37a), (37b), (38a), and (38b) are adapted using the
form given by Reid and Bodine [10], with changes in notation.
Following the same procedure for the yequation yields
q(i,j,n+) [qy(i,j,n) 2w(sin) q (i,j,n)At
y y
2pA { (i,j,n+1/2) + (i,j1 ,n+1/2)}
y y
At 1
+ x {D(i,j,n+l/2) + D(i,jl,n+l/2)} { [pn(i,jl,n+1/2)
pn(i,j,n+1/2)] + g[n(i,jl,n+1/2) n(i,j,n+1/2)]}] (38a)
where
fAt 1f2
F = 1 + ~ {[qx(i,j,n)]2 + qy(i,j,n)]2} 1/2
[D(i,j,n+1/2) + D(i,jl,n+1/2)]2 (38b)
Here, Ax and Ay are the grid spacings as shown in Figure 33. Figure
33 also defines the variables in Equation (37a) and the sense of
operation of a grid element.
The values for the variables, n, p h, D, T are assumed to be
centered in each grid as shown in Figure 33. The values are, after
that, assumed uniform over the entire grid element for the calculation
at that time step.
Since qx(i ,j,n) and q (i,j,n) occur at different points it becomes
66
c cI
Lw
C.O
 x .i
or
3 c f c L
U5 )

0 I
c c2 LL

ow w
S oo
C "': I O
000
a EE
E5 .e 2
. Q
0 "L
0 E0 O
W V3 C
L. 0
, )<> J
. ....t Ll
a o o
Sr a. .
"r')
S E Eu
" ." 3"
.C Lo
necessary to obtain spatially averaged values q (i,j,n) and q (i,j,n)
for use in the Coriolis term and the friction terms Fx, F Thus,
q,,(i,j,n) and qx(i,ji,n) are considered to be at the same point and
J
qx(i,j,n) and q (ij,n) are considered to be at the same point. This
is illustrated more clearly in Figures 34a and 34b. q Ci,j,n)
is defined as
qy(ij,n) = [qy(i,j,n) + q (il,j,n) + q (i,j+l,n) + q (il,j+l,n)]
and qx is expressed by
1
q (i,j,n) = [qx(i,j,n) + qx(i+l,j,n) + qx(i+l,jl,n) + q (i,jl,n)]
The equation of continuity, Equation (36c), becomes in finite
difference form,
n(i,j,n+3/2) = n(i,j,n+1/2) + A [qx(il,j,n+l)qx(i,j,n+l)]
+ [q (i,jl,n+l) qy(i,j,n+l)] (39a)
The total depth D is defined by
D(i,j,c) = h(i,j) + n(ij,E) (39b)
where is an arbitrary time increment.
To carry out an actual computation, a grid network is set up
as illustrated in Figure 35.
The coordinate convention is the same as Figure 21 with x
positive toward shore and ypositive to the left while facing the
shoreline. The boundary conditions are illustrated in Figure 35 and
are described as follows: 1) the fluid flux across any sealand
COMPONENTS OF y(i, j, n)
t
qx (+ j)
(i, j n)
qy(i, j, n)
qC(i, j, n)
( iI, j)
qx(i+1, jI)
q (i, j 1, n)
t
qx(i, j1, n)
(i1, jl)
FIGURE 34b
COMPONENTS OF (i, j, n)
FIGURE 34a
69
r7 Centered in Seaward and Lateral Boundary Rov's is
Equated to the Sum of Static Borometric Tide
and Astronomical Tide
I
qs Normal to Land Boundaries
are Equated to Zero
FIGURE 35
SCHEMA OF GRID SYSTEM AND
BOUNDARY CONDITIONS EMPLOYED IN
TWO DIMENSIONAL CALCULATIONS
L.
S"0
o
Ju
L5c
aI UJ
boundary is equal to zero; 2) the boundary condition at the seaward
boundary is given as n equal to the barometric tide plus the astro
nomical tide at that point; 3) at the lateral boundaries qx = 0;
4) also at the lateral boundaries n is set equal to the barometric
tide plus the astronomical tide. Furthermore it will be assumed that
the lateral boundaries are far enough away from the region of interest
that the lateral boundary conditions have little effect. It should be
reemphasized that the lateral conditions are those conditions arrived
at empirically, which allow calculations without instabilities or other
problems with the model. Initially the system is considered quiescent
and in barometric equilibrium. This condition is of little effect
providing the storm begins sufficiently far offshore. In this model
the bathymetry and shoreline are arbitrary and the presence of islands
is allowed.
Currents in the coastal waters are of interest for several
reasons, among them: 1) to aid in making and evaluating future current
and surge calculations, 2) at the coastline hurricanegenerated
currents are of interest as a cause of beach erosion, 3) hurricane
generated currents may be instrumental in producing extreme forces on
ocean structures, and 4) currents are instrumental in dispersion and
pollution studies. Thus, it is of further interest to describe the
scheme used to calculate the water velocity at a specified grid
element. As in the onedimensional case the total depth at time n1/2
is used to calculate u(i,j,n) and v(i,j,n). In order to maintain a
consistent approach, an average in both the x and y directions is
taken. This leads to the following value of u(i,j,n) and for v(i,j,n)
as described in Figure 36.
Note:
Shaded Areas
Weighted Twice
Figures Indicate Grid Elements
Involved in Calculating Weighted
Total Depth for Water Velocity
Calculations at Grid (i,j)
FIGURE 36
FILTERS USED TO CALCULATE WATER
VELOCITIES
( (i,) (j1)
(il, j~ l) jl 1 )
(i+l, j) (i+1, j1)
(i,j) 4 (i, j1)
q y(i, j,n)
(i j) (iI, jI)
72
u(i,j,n) = 8qx(i,j,n) [2D(il,j,n)+2D(i,j,n)+D(i,j+l,n)+D(il,jl,n)
+D(i ,j1 ,n)+D(il ,j+l ,n)]'
v(i,j,n) = 8q (i,j,n) [2D(i,j1l,n)+2D(i,j,n)+D(i1l,j,n)+D(i1l,j1l,n)
+D(i+1 ,j,n)+D(i+1 ,j1 ,n)]'
Adjacent to the landsea boundary it is necessary to modify this
procedure because of the no flow boundary condition.
c. Effects of Discretization and Checks of the Numerical Model
In any model of this type it is desirable to obtain some
estimate of the effects of the grid size, time increments, and lateral
boundary conditions on the model output. A check of the programming
is also necessary. One way to accomplish this is with a series of two
or more programs which although of different concept or design, converge
to identical solutions for some limiting case. This does not neces
sarily prove all are correct but provides confidence in the results.
A onedimensional time dependent program should give results
similar to those of a twodimensional time dependent program if the
radius scale of the hurricane is large enough with respect to some
length scale characteristic of the system. A comparison of this type
has been carried out for an extremely large storm, with the results
and storm characteristics as shown in Figure 37. One would, in
general, expect the onedimensional program to over predict water
levels far from the coast since a crosssection of the storm is used
and this effect extends laterally to infinity, but in the twodimen
sional case the pressures and stresses die off with distance from the
center line. Near the coast the Coriolis effect and the twodimensional
effects dominate. The onedimensional plot shown in Figure 37a repre
sents a cross section of maximum surge heights. The bottom was a
linear slope approximated by a finite difference form.
10
Legend:
 One Dmensional Model
Two Dimensional Model
No!es:
Surge Heights Plotted ore Maximum
Values Occuring During Passage of Storm
V=10 Knots f = 02
R=40 Nauticol Miles ,

AP= 432 Inches of Mercury
Bottom Slope =.0005 .
 
.. ... . .. I
100 90 80 70 60 50 40 30 20
Distance From Coast (Nautical Miles)
FIGURE 37a COMPARISON OF MAXIMUM SURGE HEIGHTS
CALCULATED BY THE ONE AND TWO
DIMENSIONAL MODELS
S80
S60
C()
* 40
5
)
/
',
Legend:.
xComponent of Wind Velocity
 Surface Pressure
Note:
30 >
U
(U
0
29 L
28
&
20
Profile Token Through Point of Maximum Winds 0
R,V, and LP c"(e Given in Figure 370
0 L  I 1 '27 U)
100 80 60 40 20 0 20 40 60 80 100
Distance, Along xCoordinale, From Storm Center (Noutical Miles)
FIGURE 37b PROFILES OF STORM CHARACTERISTICS
USED IN COMPARISON
I
''
\
.\
,\
.\
.\
r
S\
'\
,\
'\
Ii
I0
I
120
100
+120
+80
U')
" +40
rI
, 40
5
80
120
Y Component of
Wind Velocity
00 80 60 40 20/ 0 +20 +40 +60 +80 +100
Distance Along XCoordinole
from Storm Center (Nautical
Miles)
YCOMPONENT OF WIND VELOCITY
FIGURE 37c
As a further comparison, a time independent or steady state
model was developed as a further check of the computer modeling. It
aqx
is one dimensional and time independent in that is equal to zero.
L is
qx must then be equal to zero and the bottom friction term is
eliminated. The numerical technique consists of a RungeKutta
starter and Adam's Method predictor [29]. The stepwise bottom profile,
with slope .0005 as before, was included and flooding allowed, also
as before. The storm characteristics for the onedimensional time
independent case are shown in Figure 37b, c, with the storm center
placed 10 nautical miles offshore. Note that although the storm is
stationary it has the characteristics of a storm traveling with
v = 10 knots. The numerical results and the steplike bottom profile
are illustrated by the dotted line in Figure 38.
In order to compare these results with those of the one
dimensional timedependent routine the timedependent routine was
run with identical storm characteristics and then stopped when the
storm center was 10 nautical miles offshore. In this case the
bottom friction, f, was set equal to zero and as a result the solution
oscillated. The shaded boxes in Figure 38 represent the range of
oscillations. The stepwise approximation to the uniform slope,
.0005, was again used in this case.
How well does the 10 nautical mile steplike bottom profile
model a true linear slope? This question can be answered as follows.
The time independent model was run with all variables as before, with
the single exception of a uniformly sloping bottom. These results are
plotted as the solid line in Figure 38. In deep water the results are
t'
L o 0
') 7
01)
U)ciJ
co 0 0
GJ ( 
~0
t o o
Sc _c J
2 U )0
c
I c I a
w1 1
SI I
\ \ O
UJ
O
I 0
C E W \ cr
w0 6
c : 0
o O
0 00
S i LUjo
I .
~
_I I Lii
c 0
L o I
(iOj) 3BJns uJ1S
0l~j 0EJn 0 0I
nearly identical. In the shallow region the differences become
noticeable. At the risk of over generalizing, this stepwise grid
network is permissible for values of D = h + n greater than 20 feet;
for values less than 20 feet the errors can become significant.
In the twodimensional problem another "errordimension"
exists, that of the effect of placement of lateral boundaries. This
is, in addition to, the effects of discretization. The lateral bound
aries of the model are an artifice and if they are not a sufficient
distance from the zone of high wind stress then their proximity will
influence the calculated storm tides and velocities. In order to
evaluate this effect some basis for comparison must be available. A
twodimensional analytical treatment of storm tides on continental
shelves is available from Dean and Pearce [21].
Consider the situation illustrated in Figure 39 in which a
"strip" wind field is uniform within a total width, w, and zero out
side of this width. The Coriolis effect is omitted and the bottom
shear stress is expressed as proportional to the transport components
(i.e. linearized). Only shear stresses in the horizontal plane are
considered, i.e. adjacent water columns are considered to be laterally
uncoupled by stress terms; it may be anticipated that this will result
in a discontinuity in qx, the transport component in the xdirection.
The windstress in the region jyj < w/2 drives the water onshore and
results in a storm tide field which is a maximum at (x' E x/Z = 1.0,
y' = y/k = 0). The qx field in this region is a maximum at the
shelf break and decreases to zero at the coastline. It is clear, on an
intuitive basis, that the narrower the wind field relative to the
shelf width, the smaller the storm tide will be. For a very
of Active
Stress
Elevation
FIGURE 39
IDEALIZED SHELF
WITH UNIFORM D
RESPONSE PROBLEM
EPTH
Shelf 
Break
broad wind field, the storm tide at the shelf should approach the value
corresponding to a onedimensional condition. For wind fields which
are laterally limited, there will be a gradient in the storm tide
field such that the storm tide approaches zero as y' + The
surface gradient and wind stress fields will cause an onshore transport
in the wind stress region, an offshore transport outside of this
region and a y component of transport (" longshore transport") which
is directed outward from the region of onshore wind stress. Response
in terms of dimensionless n, qx and qy are presented in Figures
310a, b, c.
If the (artificial) lateral boundaries in the numerical model are
not a sufficient distance from the zone of high wind stress, the
proximity of the model boundaries will influence the calculated storm
tides and velocities. In order to evaluate this proximity effect,
numerical model computations were carried out for the shelf geometry
shown in Figure 311. The friction and total water depth terms in
the numerical model were linearized in order that the differences in
quantities calculated by the numerical model and the analytical
solution should be indicative only of artificial effects introduced
by the numerical model. The results of these calculations are de
scribedin the following paragraphs.
The first set of calculations was carried out for a ratio of
shelf width to onehalf the wind field breadth (c /(w/2)) of 1.83
and a grid size, Lx, which is variable. The lateral extent of the
shelf, incorporated into the model is denoted w* and the ratio of
wind field width to shelf extent, a*(i.e. a*=w/w*). The parameter
y
t
32
21
C.
o05
^ o\' ~~0~
IL
0 0.2 0.4 0.6 0.8 1.0
FIGURE 310a
ISOLINES OF DIMENSIONLESS STORM
TIDE, ANALYTICAL SOLUTION
Noes
(I) See Figure 39 for
Definition Sketch
(2) Storm Tide Symmetric
About y'= 0
(3) Results Shown Apply
for a= /(w/2) = 2.0
0
( L
tt,
Y
t
3 Notes:
(I) See Figure 39 for /
Definition Sketch
(2) q'x is Symmetric
About y'= 0 A
(3) Results Sho,vn Apply
For a = /(w/2) = 2.0
2C
I 7
/
Lq;
din d S l r e s s /
0 
ILJ / / A ,
0 0.2 0.4 0.6 0.8 1.0
FIGURE 3lOb ISOLINES OF DIMENSIONLESS TRANSPORT
IN xDIRECTION, ANALYTICAL SOLUTION
I
t
2
I
I
L5/
I
Notes:
(I) See Figure 39 for
Definition Sketch
(2) q' is Antisymmetricol
About y'= 0
(3) Results Shown Apply
for a = e/(w/2) = 2.0
0
mn
0.2 0.4 0.6 0.8 1.0
FIGURE 310c ISOLINES
OF DIMlENSIONLESS
IN yDIRECTION, ANALYTICAL
TRANSPORT
SOLUTION
&'J
<