UFL/COELTR1127
NUMERICAL SIMULATION OF MEAN CROSSSHORE
CURRENTS: A STREAM FUNCTION APPROACH
by
Albert E. Browder
Dissertation
2000
NUMERICAL SIMULATION OF
MEAN CROSSSHORE CURRENTS:
A STREAM FUNCTION APPROACH
by
ALBERT E. BROWDER
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2000
ACKNOWLEDGMENTS
I wish to express my sincere appreciation and utmost respect to the chairman of my
graduate committee, Robert G. Dean, for his guidance, trust, and friendship. Working for
someone who grants you the latitude to pursue your own interests and always teaches you
something at each stage of the process has been the most rewarding aspect of my
dissertation work. For providing me the experience of teaching and always being there to
answer my "real quick" questions, I wish to thank committee member Robert J. Thieke.
To the other members of my graduate committee, who have provided me great insight and
guidance in my education and the completion of this document, I express my gratitude.
I am fortunate to have had many dedicated classmates and friends in this endeavor,
almost too many to list. Thanks go to Mark and Eric, qualifier coaches extraordinaire, for
their support and friendship. To Becky, I express my thanks for keeping me on track and
being my friend. Jamie, Justin, Kevin, Bill, Matt, Sean, Dave, Jon, Heather, Beth, and
Chris deserve thanks for all the encouragement, and distractions necessary to complete this
task. For teaching me the basics, I thank Jack Kinder. Thanks are in order for all at Olsen
Associates for their understanding and encouragement. If I have left anyone out, it is by
no means intentional, just a sign of how fortunate I am to have so many good friends.
I am ever grateful to my parents, Susan and Larry, for their never ending support.
Finally, none of this would have been possible without the love and undying support of my
wife, Keeley. This is as much yours as it is mine.
TABLE OF CONTENTS
page
ACKNOWLEDGMENTS ............................................ ii
LIST OF TABLES ................................................. vi
LIST OF FIGURES ................................................... vii
KEY TO SYMBOLS ................................................. xii
ABSTRACT ........................................................ xv
CHAPTERS
1 INTRODUCTION ............................................ 1
1.1 Objectives and Rationale ...................................... 3
1.2 Present Approach to Return Flow Modeling ....................... 4
2 LITERATURE REVIEW ......................................... 8
2.1 Motivating Literature ......................................... 9
2.2 Return Flow Models ......................................... 13
2.3 Turbulent Eddy Viscosity .................................... 18
2.4 Return Flow M easurements ................................... 25
3 DERIVATION OF GOVERNING EQUATION AND NUMERICAL
MODEL DEVELOPMENT ............................... 29
3.1 Derivation of Governing Equation .............................. 30
3.2 Development of OneDimensional Model ........................ 37
3.2.1 Boundary Conditions for OneDimensional Model .......... 40
3.2.2 Solution Scheme for OneDimensional Model ............. 41
3.3 Development of TwoDimensional Return Flow Model ............. 43
3.3.1 Transformation of the Partial Differential Equation ......... 43
3.3.2 Discretization of the TwoDimensional PDE .............. 45
3.3.3 Boundary Conditions for TwoDimensional Model ......... 48
3.3.4 Solution Scheme for TwoDimensional Model ............. 49
4 BOUNDARY CONDITIONS AND TRANSITION
REGION TREATMENT .............................. ..... 52
4.1 W ave Height Transformation .................................. 53
4.2 Determination of Return Flow Rate ........................... 58
4.3 Specification of the Surface Shear Stress Due to Waves ............. 62
4.4 Specification of Bottom and Sidewall Tangential Boundary Conditions 64
4.5 Matching Conditions Across the Transition Region ................. 68
5 VALIDATION OF NUMERICAL MODEL ........................ 75
5.1 OneDimensional Analytic Solutions .......................... 76
5.1.1 SteadyState Simulations .............................. 76
5.1.2 Variable Eddy Viscosity Simulations Steady State ......... 78
5.1.3 TimeDependent Simulation ........................... 81
5.2 TwoDimensional Model Validation ............................ 87
5.2.1 Comparison to OneDimensional Results ................. 87
5.2.2 Comparison to Analytic Series Solution ................. 88
6 COMPARISON OF NUMERICAL MODEL TO
EXPERIMENTAL DATA ................................ 93
6.1 Introduction of Model Output and RMS Errors .................... 94
6.2 Simulation of Cox and Kobayashi (1997) Experiments ............. 104
6.3 Simulation of Ting and Kirby (1994) Experiments ................ 114
6.4 Simulation of Nadaoka and Kondoh (1982) Experiments ........... 124
6.5 Simulation of Okayasu and Katayama (1992) Experiments .......... 130
6.6 Simulation of Roelvink and Reniers (1995) Experiments ........... 139
6.7 Simulation of Smith et al. (1992) Experiments ................... 147
6.8 Comments ............................................... 154
7 IMPLICATIONS FOR CROSSSHORE SEDIMENT
TRANSPORT AND FUTURE STUDY ...................... 156
7.1 Variation in Predicted Bottom Shear Stress and Shear Stress Profiles .. 158
7.2 Transition Region Treatment ................................. 163
7.3 Onshore Bottom Stresses and Velocities Seaward of the Break Point .. 169
7.4 Selection of Wave Theory and Phase Speed Definition ............. 173
7.5 Selection of Vertical Distribution of Eddy Viscosity ............... 176
8 SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS
FOR FUTURE WORK ................................. 182
8.1 Summary ................................................ 182
8.2 Conclusions .............................................. 184
8.3 Recommendations for Future Work ............................ 187
APPENDICES
A NUMERICAL MODEL DEVELOPMENT ......................... 189
A.1 Development of OneDimensional Model ....................... 190
A.1.1 Boundary Conditions for OneDimensional Model ........ 192
A.1.2 Solution Scheme for OneDimensional Model ............ 195
A.2 Development of TwoDimensional Return Flow Model ............ 196
A.2.1 Transformation of the Partial Differential Equation ........ 196
A.2.2 Discretization of the TwoDimensional PDE ............. 201
A.2.3 Boundary Conditions for TwoDimensional Model ........ 207
A.2.4 Solution Scheme for TwoDimensional Model ........... 213
B FORTRAN CODE LISTING RF PSI.F90 ........................ 215
REFERENCES ................................................ 247
BIOGRAPHICAL SKETCH ........................................ 251
LIST OF TABLES
Table page
6.1 Parameters to be Investigated/Adjusted in Numerical Model
for Comparison to Experimental Data. ....................... 94
6.2 Errors in Example Simulation Shown In Figures 6.1 and 6.2 for
Cox and Kobayashi (1997) Laboratory Experiment .............. 99
6.3 Conditions for Return Flow Measurements Simulated by Numerical Model 100
7.1 Comparison of normalized RMS errors calculated from applying the eddy
viscosity distributions in Figure 7.12 for the experiment of Cox and
Kobayashi (1997). ...................................... 181
LIST OF FIGURES
Figure page
1.1 Rendering of waveinduced return flow in the surf zone. ................. 2
1.2 Definition sketch of return flow problem ............................. 5
2.1 Illustration of the effect of the quadratic mixing length expression
employed by Reid (1957). ................................. 22
2.2 Comparison of mean horizontal velocity profiles as reported by Cox and
Kobayashi (1997) and Ting and Kirby (1994). .................. 28
3.1 Definition sketch for model development. ........................... 31
3.2 Definition sketch for onedimensional model.......................... 39
3.3 Example output from onedimensional return flow model ............... 42
3.4 Definition sketch of coordinate transformation between the physical, xz,
domain and the computational, )7, domain .................... 44
3.5 Solution cells for the fourthorder finite difference equation in the xz domain
and the transformed f domain. ............................. 46
3.6 Flowchart for twodimensional numerical model ...................... 51
4.1 Stream function wave theory shoaling table for normally incident waves ... 54
4.2 Demonstration of wave transformation model for two laboratory studies. ... 56
4.3 Demonstration of the treatment of the transition region in the simulation of
the data set ofTing & Kirby (1994). ......................... 71
5.1 Definition sketch of onedimensional problem ......................... 77
5.2 Comparison of 1D numerical model to analytic solution for steadystate
conditions. .......................................... 78
5.3 Comparison of 1D numerical model results to those of Reid (1957) for the
case ofQ= 0 ................. .......................... 80
5.4 Comparison of numerical model results to superimposed analytic solutions
to the diffusion equation. .................................. 84
5.5 Time history of the forcing term, A, in Eq. (5.10). ...................... 85
5.6 Investigation of time step dependence on the onedimensional semi 
implicit numerical model. ................................. 86
5.7 Comparison of 2D numerical model to 1D analytical results of Eq. (5.4)
computed at each vertical section. ............................ 89
5.8 Comparison of numerical model results to Fourier Series solution, Eq. (5.17). 92
6.1 Example contour and velocity vector plot of laboratory experiments by Cox
and Kobayashi (1997). .................................... 95
6.2 Example of cross shore variation in wave height, radiation stress,
volumetric transport, and vertical profiles of horizontal velocity
for the experiments of Cox and Kobayashi (1997)................ 96
6.3 Comparison of cross shore variation in wave height, radiation stress,
volumetric transport, and vertical profiles of horizontal velocity
for the experiments of Ting and Kirby (1994) (e= 0.001 m2/s). ..... 102
6.4 Comparison of present radiation stress transition region patch to the wave
surface roller model of Dally and Brown (1995) for the data set
of Ting and Kirby (1994) (e= 0.001 m2/s) .................... 103
6.5 Comparison of simulations of the data set of Cox and Kobayashi (1997). .. 105
6.6 Comparison of predicted radiation stress and volumetric transport for the
experiments of Cox and Kobayashi (1997)..................... 107
6.7 Comparison of normalized RMS errors, 8, for the experiment of Cox and
Kobayashi (1997) employing spatially variable eddy viscosity fields. 110
6.8 Comparison of vertical profiles of horizontal velocity at station #4 of the
Cox and Kobayashi (1997) data set. ......................... 111
6.9 Comparison of vertical profiles of horizontal velocity at station #4 of the
Cox and Kobayashi (1997) data set for a zero bottom shear stress
boundary condition. ..................................... 112
6.10 Predicted flow field for the data set of Cox and Kobayashi (1997). ....... 113
6.11 Velocity profiles corresponding to the flow field in Figure 6.10 predicted
for the data set of Cox and Kobayashi (1997). ................. 114
6.12 Comparison of simulations for the data set of Ting and Kirby (1994). ...... 115
6.13 Comparison of predicted radiation stress and volumetric transport for the
experiments of Ting and Kirby (1994) ....................... 117
6.14 Comparison of normalized rootmeansquare errors, ,, for the experiments
of Ting and Kirby (1994) employing spatially variable eddy
viscosity fields. ........................................ 119
6.15 Comparison of vertical profiles of horizontal velocity at Station #2 (offshore
of the break point) of the Ting and Kirby (1994) data set for a noslip
bottom boundary condition ................................ 120
6.16 Comparison of velocity profiles at Station #2 (offshore of break point) of the
Ting and Kirby (1994) data set. ............................. 121
6.17 Predicted flow field for the data set of Ting and Kirby (1994) ............ 122
6.18 Velocity profiles corresponding to the flow field in Figure 6.17 predicted for
the data set of Ting and Kirby (1994)......................... 122
6.19 Comparison of simulations of the data of Nadaoka and Kondoh
(1982, Case 1). ........................................... 125
6.20 Comparison of predicted radiation stress and volumetric transport for the
experiments of Nadaoka and Kondoh (1982, Case 1). ........... 126
6.21 Comparison of normalized rootmeansquare errors, f, for the experiments
of Nadaoka and Kondoh (1982) employing spatially variable eddy
viscosity fields. ........................................ 128
6.22 Predicted flow fields for the experiments of Nadaoka and Kondoh (1982). 129
6.23 Velocity profiles corresponding to the flow field in the lower plot of Figure
6.22 predicted for the data of Nadaoka and Kondoh (1982) ....... 130
6.24 Comparison of predicted RMS wave height, radiation stress, volumetric
transport for the random wave experiments of Okayasu and
Katayama (1992)......................................... 132
6.25 Comparison of simulations of the data set of Okayasu and
Katayama (1992)......................................... 133
6.26 Comparison of normalized rootmeansquare errors, 8, obtained for the
experiments of Okayasu and Katayama (1992) by applying spatially
variable eddy viscosity fields. .............................. 135
6.27 Predicted flow field for the experiments of Okayasu and
Katayama (1992)......................................... 136
6.28 Comparison of simulations applying random waves averaged over time and
applying a monochromatic wave field ofRMS height ........... 137
6.29 Comparison of simulations of the data set of Roelvink and Reniers (1995)
from the Delta flume experiments in Delft, The Netherlands ...... 140
6.30 Demonstration of the effect of the choice of surface eddy viscosity. ...... 141
6.31 Comparison of predicted radiation stress and volumetric transport for the
experiments of Roelvink and Reniers (1995) .................. 143
6.32 Comparison of normalized rootmeansquare errors, A, obtained for the
experiments of Roelvink and Reniers (1995) by applying spatially
variable eddy viscosity fields. .............................. 145
6.33 Flow field predicted for the experiments of Roelvink and Reniers (1995). 146
6.34 Comparison of simulations applying timeaveraged random waves versus
monochromatic waves of RMS height......................... 146
6.35 Comparison of simulations of the data set of Smith et al. (1992).......... 148
6.36 Comparison of predicted radiation stress and volumetric transport for the
experiments of Smith et al. (1992). ......................... 149
6.37 Comparison of normalized rootmeansquare errors, f, obtained for the
experiments of Smith et al. (1992) by applying spatially variable
eddy viscosity fields. .................................... 152
6.38 Comparison of simulations applying timeaveraged random waves versus
monochromatic waves of RMS height ....................... 153
7.1 Crossshore variation in mean bottom shear stress for the experiment
of Cox and Kobayashi (1997) ............................... 159
7.2 Crossshore variation in mean bottom shear stress for the experiment
of Cox and Kobayashi (1997) ............................... 161
7.3 Vertical profiles of velocity, shear stress, and eddy viscosity at two measuring
stations of the experiment of Cox and Kobayashi (1997) ......... 162
7.4 Predicted flow field for the experiment ofNadaoka and Kondoh
(1982, Case 1). ........................................... 165
7.5 Comparison of predicted and measured profiles of vertical velocity for the
experiment of Nadaoka and Kondoh (1982, Case 1). ............. 166
7.6 Predicted flow field for the data set of Cox and Kobayashi (1997). ....... 167
7.7 Comparison of predicted and measured vertical mean velocities from the
experiments of Cox et al. (1995)/Cox and Kobayashi (1997). ..... 168
7.8 Schematic of the effect of increasing the magnitude of the offshore directed
shear stress applied in the shoaling region. .................... 171
7.9 Predicted flow field in the offshore region for the experiments of Roelvink and
Reniers (1995)............................................ 172
7.10 Contours of nondimensional wave induced momentum flux (radiation stress)
determined from the 40 cases tabulated for Stream Function Wave
Theory by Dean (1974). .................................. 174
7.11 Comparison of predicted velocity profiles applying offshore and onshore
directed surface stresses. .................................. 177
7.12 Comparison of nondimensionalized vertical distributions of eddy viscosity. 179
7.13 Comparison of predicted velocity profiles applying the eddy viscosity
distributions plotted in Figure 7.12. .......................... 181
KEY TO SYMBOLS
A# Steady state term coefficients describing the transformation of coordinates between
the xz domain and the ,7 domain (# 1 through 14, units vary).
A, Crosssectional area of wave surface roller (m2).
B# Time dependent term coefficients describing the transformation of coordinates
between the xz domain and the rq domain (# 1 through 5, units vary).
C Wave celerity (m/s).
Cj Solution coefficients in Fourier Series analytic solution (s').
CL# Steady state term coefficients arising from discretization process on the solution
column of cell (# 1 through 5, no units).
CR# Steady state term coefficients arising from discretization process off the solution
column of cell (# 1 through 16, no units).
DL# Time dependent term coefficients arising from discretization process on the
solution column of cell (# 1 through 3, no units).
DR# Time dependent term coefficients arising from discretization process off the
solution column of cell (# 1 through 6, no units).
D Total water depth (in Reid, 1957, m).
E Total wave energy per unit surface area (N/m2).
ERMs Rootmeansquare error of each predicted velocity profile (m/s).
g Acceleration due to gravity (m/s2)
H Wave height (m).
H, Significant wave height (m).
Ho Deepwater significant wave height (m).
h Still water depth (m).
h+f Total water depth (still water depth plus setup/setdown, m)
i, j Horizontal and vertical cell coordinates, respectively, in numerical solution.
k Wave number (27t/wavelength).
ko Von Karman's Constant (= 0.4).
K Dimensionless wave decay coefficient.
L Wavelength (m).
L Mixing length (turbulence length scale, m).
L, Deepwater wavelength (m).
m Ratio of bottom to surface shear stress.
n Time step counter in numerical solution.
Q Wave induced volumetric transport (m3/s/m).
S, Radiation stress, onshore flux of waveinduced onshore momentum (N.m/m2).
t Time (s).
T Wave period (s).
U, W Horizontal, vertical components of mean velocity (m/s).
Ub Tangential bottom velocity (m/s).
U. Streaming velocity (tangential velocity at upper edge of bottom boundary
layer, m/s).
u, v Horizontal, vertical components of wave induced, orbital, velocity (m/s).
u w' Horizontal, vertical components of turbulence induced velocity (m/s).
x Horizontal coordinate in physical domain (m).
z Vertical coordinate in physical domain (m).
a Wave number in Fourier Series analytic solution (m).
xiii
p Nondimensional rootmeansquare error in each velocity profile.
F Ratio of stable wave height to local water depth.
At Time step in numerical solution (s).
Ax Horizontal spatial step in physical domain (m).
Az Vertical spatial step in physical domain (m).
Atj Vertical spatial step in numerical solution (no units).
Ae" Horizontal spatial step in numerical solution (m).
e Turbulent eddy viscosity (m2/s).
it Setup/setdown of mean water surface from still water surface elevation (m).
r7 Vertical coordinate in transformed computational domain (= z/(h+~i), no units).
0 Degree of implicitness in numerical solution.
K Ratio of wave height to local water depth at breaking.
v Molecular kinematic viscosity (m2/s).
e Horizontal coordinate in transformed computational domain (m).
p Density of water (freshwater or saltwater, kg/m3).
a Wave frequency, (27t/T, s'l).
1Tb Viscous shear stress in the adirection on the bplane.
rZb Turbulent shear stress in the adirection on the bplane
Zb Bottom shear stress (N/m2).
r, Shear stress at instantaneous free surface (N/m2).
r, Shear stress at mean water surface (N/m2).
fr Stream function (s').
V Del operator (total derivative operator).
I
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
NUMERICAL SIMULATION OF
MEAN CROSSSHORE CURRENTS:
A STREAM FUNCTION APPROACH
By
Albert E. Browder
August 2000
Chairman: Robert G. Dean
Major Department: Civil & Coastal Engineering
A numerical model applying a stream function formulation is developed and tested
to simulate the twodimensional, waveinduced return flow in the nearshore region,
commonly referred to as the undertow. The model is developed using the stream function,
r, to represent the waveinduced return flow, and tests are conducted of the effects of
various wave theories and boundary conditions on the flow. Solution of /produces the
mean velocity in the horizontal and vertical directions, providing a description of the
complete flow field in the twodimensional vertical plane perpendicular to the shoreline.
The governing equation, a fourthorder partial differential equation in f, is not
specific to waveinduced flows; therefore, the details of the wave forcing are introduced
entirely through the boundary conditions. In this way the model provides a convenient
platform to evaluate different combinations of applied boundary conditions. The model
includes the ability to study spatial variations in the turbulent eddy viscosity. The present
model also accommodates monochromatic or random wave fields and arbitrary bottom
profiles, allowing the modeling of flows over bartrough beach profiles or any other
irregular profile geometry.
Comparison to measurements of waveinduced return flows demonstrates that the
model is able to predict the vertical structure and magnitude of the return flow to within
25% to 30% (typical) for most applications. If the return flow rate is properly predicted (or
fit), the model can predict the vertical structure of the horizontal velocity to within 10%
(typical). Simulation of six laboratory and field experiments demonstrates the capability
of the model and highlights the challenges of return flow modeling. The simulations
indicate that the most important parameter in return flow modeling is the proper prediction
of the waveinduced volumetric transport. Further improvement of the predictive capability
of the model inside the surf zone stems from the application of a verticallyvarying eddy
viscosity field based on the turbulence generated by breaking waves. Comparison of the
vertical structure of the return flow seaward and landward of the break point illustrates
differences in the appropriate boundary conditions and eddy viscosity distributions that
should be applied in return flow modeling.
CHAPTER 1
INTRODUCTION
Characterizing fluid flow in the surf zone presents the substantial challenge of
identifying many distinctly different processes and attempting to quantify the effects and
interactions of each. While the flow field in the nearshore region is at least visually
dominated by the passage and breaking of waves on a time scale of seconds to minutes, the
mean flows created by these incoming waves also play an important role in the nearshore.
Additionally, currents introduced by other phenomena, such as tides or winds, may
contribute to the flow field in any nearshore region. On a shorter temporal scale, the effects
of turbulence must also be considered. The combination of these factors presents a
daunting task to the coastal engineer charged with describing or modeling these flows in
order to design coastal structures or predict changes to the nearshore morphology due to
waves.
This study focuses on one aspect of the flow field in the surf zone, that being the
mean flow in the crossshore, twodimensional vertical plane (2DV), as rendered in Figure
1.1. The desire to understand and accurately predict the mean flow patterns in the
nearshore region stems primarily from the interest in predicting changes in the beach
profile induced by waves. The mean waveinduced crossshore and longshore currents act
to transport sediment not only by placing stress directly on the seabed, but also by
transporting the sediment mobilized and suspended by the passage of other currents,
individual waves, or groups of waves (or all three and more in some cases). In the case of
2
the mean crossshore flow, mobilization and/or transport of sediment by the mean wave
induced current obviously plays a role in shaping the crossshore profile and may play a
strong role in the formation of longshore bars in the profile.
The mean crossshore current is frequently referred to as the "undertow" but is more
accurately labeled the return flow and shall be discussed as such in this study. The term
undertow has been used frequently in the media in conjunction with the occurrence of rip
currents, hence the actual definition of the undertow is somewhat vague while the return
flow and rip currents are two distinct entities. The return flow always exists in the surf
zone, whereas rip currents appear under various conditions but are not always present.
' ... , ,
7 7
...  5.LB. 2000
Figure 1.1 Rendering of waveinduced return flow in the surf zone. The return flow is
the twodimensional flow in the vertical plane perpendicular to the shoreline. The
relationship of the return flow to longshore currents and rip currents is illustrated.
(~A
relationship of the return flow to longshore currents and rip currents is illustrated.
3
1.1 Objectives and Rationale
The objective of this study is to develop a numerical model of the mean flows in
the nearshore region via a stream function formulation and to validate the approach. The
driving interest in the development of a stream function return flow model is the desire to
investigate various individual aspects of the mean flow. These aspects include the effects
of wave transformation, the application of various combinations of boundary conditions,
and the effect of spatially variable turbulent eddy viscosity, e, in the nearshore region.
In addition, it is of interest to investigate the differences between onedimensional
and twodimensional formulations of the return flow profile, the role of various
formulations for mass transport and shear stress at the water surface, and the possible role
of the equilibrium of the moments of the wave induced forces across the beach profile.
DyhrNielsen and Sorensen (1970) first hypothesized that the balance of moments at a
given vertical section may indeed play a role in the formation of longshore bars. While not
much attention has been paid to this possibility, the present model provides a convenient
means to demonstrate the idea and its possible alteration of the return flow profile,
primarily through the simple application of additional boundary conditions in the model.
Dally (1980) presented an analytic formulation using the stream function to model
the vertical variation of the horizontal momentum of fluid dominated by uniform eddy
viscosity. The governing equation Dally used is the biharmonic equation, which is also
used to describe creeping flows dominated by molecular viscosity. Dally presented a series
solution for the horizontal momentum in a channel subjected to a surface stress. The
passage of the breaking (or shoaling) wave provided the surface shear stress applied at the
mean water surface; flow conservation was achieved by specifying the values of the stream
function at the boundaries. The fourth condition specified the bottom velocity (= 0).
4
Dally reported that the resulting flow field included only weak flows away from the
point of application of the surface shear stress. The model presented by Dally was
developed for a flat bottom and was not coupled horizontally to shear stresses in adjacent
cells. The use of the stream function in a fourthorder PDE, however, prompted the interest
in creating a twodimensional model that would accommodate an arbitrary beach profile
and would attempt to horizontally couple the flow.
1.2 Present Approach to Return Flow Modeling
Figure 1.2 illustrates the basic problem formulation. LonguetHiggins and Stewart
(1962) derive the balance of waveinduced radiation stress and the subsequent uniform
pressure gradient caused by the water surface slope (setup or setdown). Many researchers
have discussed that these two forces, while balanced in a depthaveraged sense, are not
balanced at every elevation throughout the water column. The gradient in water surface
slope acts uniformly throughout the water column, while the radiation stress is greater at
the surface than the seabed. This imbalance drives an offshore flow along the lower
portions of the water column. This flow is termed the return flow.
Using the scalar stream function, r, in a 2D return flow model is appealing for a
number of reasons. Solving for produces the velocity components in both directions,
whereas many models solve only for the horizontal velocity and do not address the vertical
velocity. This allows for the solution of the entire 2DV velocity field across the profile,
providing a complete description of the predicted mean flow patterns in the surf zone. Use
of the stream function produces one scalar variable, to be solved in one equation (the
crossdifferentiated, time averaged momentum equation). Computationally, the solution
for two variables (each velocity component) adds an additional complexity to the problem.
Figure 1.2 Definition sketch of return flow problem. The return flow is created from the
vertical imbalance of the wave induced radiation stress in the onshore direction and the
uniformly distributed pressure gradient in the offshore direction that arises from the
slope in the water surface.
The crossdifferentiation process also removes the pressure terms from the
equations, eliminating the need to describe the pressure field in the model. Modeling the
problem in two dimensions allows for an investigation of the coupling of stresses from one
vertical section to the next, as opposed to considering the vertical sections separately. This
becomes a particular point of interest in the immediate vicinity of breaking. Noncoupled
computations tend to produce a strong discontinuity in the description of the bottom shear
stress near the breakpoint. When applied to profile evolution models, the discontinuity can
lead to large errors in the predicted behavior of the profile.
The application of the boundary conditions becomes a particularly interesting and
potentially beneficial issue. The governing equation is a 2D, fourthorder PDE in t, hence
four boundary conditions are needed in each direction. In the vertical, this provides the
ability to prescribe the boundary conditions at the mean water surface and at the seabed
while also explicitly dictating that mass (or flow rate) be conserved. This fact provides
additional information not available from many previous return flow model formulations.
6
In the horizontal, the benefit of the inclusion of four boundary conditions is not as
significant, although it provides a convenient means for the transfer of fluid into and out
of the model domain, as might be the case for an external, nonwave driven current.
One drawback to having the additional boundary conditions is then the need to
properly describe these boundary conditions, such as the specification of the velocity or the
shear stress at the seabed. The governing equation contains no information that specifically
relates to a wavedriven flow, thus all the information regarding the wave forcing must be
supplied through the boundary conditions. While the model provides a convenient
approach to evaluate these options, application of specific values for these boundary
conditions may not be supported by the existing knowledge of the processes near the
seabed or the surface.
In addition to the challenges of establishing boundary conditions, another problem
that modelers of waveinduced return flow face is the characterization of the flow
immediately following breaking. In this area, commonly referred to as the transition
region or outer surf zone, the wave undergoes a significant transformation of shape but does
not necessarily experience a significant decrease in its momentum. Researchers have
studied this problem and at present a common technique is to describe the wave in the surf
zone as a turbulent bore, in which the wave form contains a volume of water in front of it
referred to as the surface roller (Svendsen, 1984a).
While the surface roller technique developed by Svendsen appears to apply well in
the inner surf zone, the details of how this roller forms in conjunction with the conservation
of momentum across the entire surf zone are not well defined. It will be demonstrated that
momentum is not conserved between the offshore region and the inner surf zone
description of Svendsen (1984a), as well as other traditional wave theories (e.g., linear
7
wave theory). Dally and Brown (1995) introduce a technique to describe the redistribution
of momentum and the formation of a wave surface roller in the transition region that does
satisfy momentum conservation. In this work, the method of Dally and Brown (1995) is
applied and compared to a simple attempt to "patch" the offshore and inner surf zone
descriptions by applying a simple conservation of momentum argument between the
offshore wave theory and the inner surf zone theories (such as the Svendsen roller model).
It will be shown that the conservation of integral wave properties is imperative in creating
a smoothly coupled model (more so than the use of a 2D versus 1D formulation).
Inclusion of an arbitrary beach profile geometry and a spatially variable turbulent
eddy viscosity field precludes an analytic solution of the 2D problem. Therefore a two
dimensional numerical model is developed to investigate the aforementioned issues. The
2D numerical model is compared to simple flow situations to check the results against
existing 2D analytic solutions and other return flow models.
The inclusion of a timedependent component of the mean flow also provides some
advantages. The time dependency allows for an investigation of the effects of a random
wave field on the return flow since the effects of previous waves are maintained (the model
tracks the 'history' of the flow). The ability to incorporate random waves is beneficial
since it tends to produce bottom velocities and shear stresses that are more uniform and
smooth across the profile, which lends itself to the successful coupling to a sediment
transport module. Additionally, the ability to model a changing wave and flow field could
ultimately lead to the simulation of storm induced changes in a beach profile, where the
wave climate builds then diminishes over several hours.
CHAPTER 2
LITERATURE REVIEW
A significant amount of literature exists that discusses the mean flow field in the
surf zone. Johnson (1919) qualitatively describes the nature of many currents present in
the nearshore region, and in describing 'hydraulic currents' due to waves he notes that there
is a net shoreward movement of water that occurs in shorewarddirected oscillatory waves
that must result in the accumulation of water along the coast. Johnson further states:
It is clear that the water piled up against a shore in the manner just described must
escape, thereby producing more or less continuous 'hydraulic currents.' If the escape
is seaward, along the bottom, we have the current known as the undertow (p. 104);
Johnson clearly describes the physical mechanisms of interest in the present work and their
general cause. A more rigorous explanation of the generating mechanisms of the return
flow is provided in the introductory chapter based on the balance of the waveinduced
radiation stress and the resultant setup or setdown of the water surface, as described by
LonguetHiggins and Stewart (1962) (Figure 1.2).
LonguetHiggins and Stewart (1962, 1964) present the mean horizontal momentum
balance, depthintegrated and timeaveraged, for the passage of waves over a sloping
bottom. As presented graphically in Figure 1.2, the balance equation is:
s pg(h+i)) Y, (2.1)
ax ax
where S, is the wave induced radiation stress (the flux of onshore directed wave
momentum in the onshore, x, direction), iris the mean water surface elevation, its gradient
in x represents the water surface slope, and fb is the mean bottom shear stress. The mean
bottom shear stress in Eq. (2.1) is frequently labeled as a minor contribution to the
momentum balance, and LonguetHiggins and Stewart discarded the term in their
discussion. In a depth integrated sense, the radiation stress and the water surface slope
balance, but as a function of depth (z), the terms are unbalanced, which creates the return
flow in the surf zone.
The objective of this work is to describe the development and application of a
stream function formulation to model the mean flow field in the surf zone. With that goal
in mind, it is useful to review and describe previous efforts that have prompted the
development and use of this particular formulation, as well as the data that exist that allow
its demonstration.
Return flow modeling depends strongly on the selection of boundary conditions.
In each case presented the choice of boundary conditions is discussed. The application of
the boundary conditions to the numerical model and the individual choice of values for
each condition are discussed in Chapter 4. Similarly, a variation in the eddy viscosity field
in the nearshore region can affect the return flow. These possibilities are demonstrated
throughout the present work, and their contributing researchers are noted accordingly.
2.1 Motivating Literature
DyhrNielsen and Sorensen (DN&S, 1970) provide a review of nearshore processes
that sets the stage for much of the recent work done in this area. DN&S qualitatively
describe the balance of wave thrust and pressure gradient and the resultant onshore directed
10
current along the seabed outside the break point. They also describe the balance inside the
surf zone (Figure 1.2) where energy losses due to wave breaking lead to offshore directed
currents at the seabed. DN&S hypothesize that this combination of currents contributes to
the transport of sediment toward the break point and leads to the development of shore
parallel sand bars.
DN&S present the balance of forces in terms of moments which must be balanced
in the water column. The authors conclude that the passage of broken or unbroken waves
creates an unbalanced moment within the water column. According to DN&S, this
unbalanced moment must be taken up at the bed, resulting in onshore and offshore directed
bed stresses seaward and landward of the surf zone, respectively. The possibility of
onshore directed stresses caused by wave shoaling seaward of the break point provides one
motivation for the development of the present model, which provides a means of clearly
demonstrating the effects of the balance of moments and the possibility of onshore directed
transport in the shoaling region.
The work of Dally (1980) provides the motivation for the formulation applied in the
present study. Dally developed a numerical model to study beach profile evolution (Dally,
1980, and Dally & Dean, 1984). In the development of the profile evolution model, the
author derived an expression for the mean return flow profile based on the primitive
conservation of fluid momentum equations. The expression was then coupled with a
description of the sediment suspension concentration profile to produce a sediment
transport model in the surf zone. The return flow model developed by Dally will be
discussed elsewhere in this chapter.
In the same work, Dally (1980) presents a brief attempt to model the horizontal
momentum of fluid in a flat channel under an applied surface stress. In this effort, Dally
11
timeaveraged and crossdifferentiated the equations of motion to produce the biharmonic
equation (with some algebraic manipulation):
1 + 2 4* +  0 (2.2)
X 4 X 2aZ 2 aZ4
This equation can also be used to describe creeping flows in viscous fluids, and while the
analogy with creeping flows does not hold in the present physical situation, the existing
mathematical solutions from creeping flows are useful in this work to verify the numerical
solutions. The present model builds on the work of Dally (1980) by developing a similar
fourthorder partial differential equation (PDE) in terms of the stream function that
considers arbitrary geometries and allows a spatially variable eddy viscosity.
Further motivation for the present work arises from the differences between various
existing return flow models. The secondorder model used by Dally is derived from the
time averaged equations of motion and includes the wave troughtocrest region, a process
that produces a second order equation in the mean velocity, U. This is accomplished by
employing linear wave theory, in which the troughtocrest region is collapsed to the Mean
Water Level (MWL), consistent with small amplitude wave theory. Other modelers, e.g.,
Svendsen (1984b) and Stive and Wind (1986), truncate the model domain at the trough
level so as not to invoke linear wave theory. The inclusion of the troughtocrest region in
the derivation incorporates details of the wave momentum flux as an internal forcing
mechanism in the development (as opposed to applying this flux as a boundary condition).
The other approach, that of truncating the domain at the wave trough level, leaves the
momentum contributions in the troughtocrest region as an option for a surface boundary
condition, and thus somewhat separates the flow description from the wave motion.
12
In the latter approach, a detail of the flow is lost. The need to specify two boundary
conditions, without the details of the wave troughtocrest momentum explicitly included
in the governing equation, results in the loss of one of the three pieces of information
generally used to describe the problem (flow rate, upper boundary condition, and lower
boundary condition). In these cases, previous investigators have generally specified the
conservation of flow through the section and either the bottom condition (e.g., Svendsen,
1984a) or the surface shear stress (e.g., Stive and Wind, 1986). The other, unspecified
condition is found as a product of the derivation.
The differences in these two approaches leads to confusion with regard to the
boundary conditions and the application of a boundary condition versus a constraint on the
flow. In the present development, the governing equation is derived from the equations of
motion and is somewhat generic (little wave information is evident in the final equation).
Waverelated flow characteristics are introduced strictly through the boundary conditions.
In this manner the effects of various aspects of the wave passage are elucidated by the use
of the fourthorder PDE. The four boundary conditions needed (in the vertical) will be
applied at the surface and the bottom to dictate a surface boundary condition, a bottom
boundary condition, and the conservation of the flow through any vertical section.
In the well known work of LonguetHiggins (1953), the author develops a fourth
order stream function expression for profiles of the mass transport velocity under water
waves. In the 'conduction solution' developed, LonguetHiggins introduces the application
of the same four boundary conditions that will be applied in the vertical coordinate in the
present model. The author discusses the use of the first and secondderivative boundary
conditions (velocity and shear stress, respectively) where the chosen values correspond to
matching conditions at boundary layers at both the free surface and the bottom. It is noted
13
that the velocity profiles developed by LonguetHiggins (1953) apply to uniform eddy
viscosity flows in relatively deep water (compared to the wave height), and the specific
boundary conditions applied are not expected to hold in the present surf zone application
(e.g., LonguetHiggins applied a bottom streaming velocity condition, which, as will be
shown, is not well supported by experimental measurements in the surf zone).
It is the efforts of LonguetHiggins, DyhrNielsen and Sorensen, and Dally that
provide the basis and motivation for the present numerical modeling investigation of the
mean flow in the surf zone. It is the goal of the present work to combine these efforts into
a model that more completely describes the mean flow processes in the nearshore region.
2.2 Return Flow Models
Many wellknown return flow models are developed directly from the two
dimensional momentum equations and produce a second order differential equation in U
by introducing a turbulent eddy viscosity formulation. Bearing that in mind, it is useful to
investigate the full 2D Reynolds equation of turbulent flow in the xdirection, where the
horizontal velocity components u and w have been decomposed into individual mean ( U,
W) fluctuating ( u, w), and turbulent ( u w ') components:
+ [ (U + u7 + 2 i] + [ (U + u7 + u)(W + w~ + w ) ]
at ax az
1 Bji a a
1 + [ 2va(U++u/)] (2.3)
p ax ax ax
a aa
+ [ v( (U++u) + (W+i+w + )]
az az ax
Generally, in the literature the longperiod fluctuation of the mean flow is neglected
along with many of the wave induced fluctuating components. Additionally, it is
14
customary to neglect any turbulent normal stresses. The turbulent shear stresses, known
as the Reynolds stresses, can be replaced by a turbulent shear stress representation
employing the eddy viscosity, e. The introduction of the turbulent shear stress
representation then produces a second order equation in U, resulting in the need for two
boundary conditions for solution:
( w) = 2 + _[E (2.4)
x x x dz 8z
The reader will note that the depthintegration of Eq. (2.4) produces the original
balance shown in Eq. (2.1), where the term on the lefthand side of Eq. (2.4) leads to the
horizontal gradient in wave radiation stress, S,, and the integration of the second term on
the right hand side leads to terms for the surface and bottom shear stress, which are then
frequently neglected. Eq. (2.4), or a similar form, is the foundation for many return flow
models including Dally (1980), B6rek9i (1982), Svendsen (1984b), and Stive and Wind
(1986), among others. These investigators applied Eq. (2.4) with a variety of boundary
conditions and wave descriptions. A detailed derivation of the Reynolds equations and the
substitution of the eddy viscosity representation is included in Chapter 3 as part of the
development of the governing equation used in the present model.
Dally (1980) and Dally and Dean (1984) develop an expression for the return flow
based on Eq. (2.4). Dally includes the wave troughtocrest region in the solution domain,
retaining the wave momentum flux in that region as an internal forcing mechanism in the
governing equation. In linear wave theory, this region is considered to be collapsed to the
MWL. This assumption results in the apparent application of a surface shear stress without
the application of a formal boundary condition. The shear stress at the instantaneous free
15
surface and the velocity at the bottom are then applied as the two boundary conditions. The
equation, which contains the term describing the gradient in water surface slope, is
constrained such that the net flow through the vertical section, Q, equals the volumetric
transport of the waves (i.e., the Stokes drift). This constraint eliminates the water surface
slope term from the equation. The resultant expression for the return flow with uniform
eddy viscosity was found to be:
U(z) = gh 8H 2 l 
8E ax 8 h 2 h 8 (2.5)
+ 3(b )2 1 3Q ) 1
2 h 2 2h h
where the coordinates are defined in Figure 1.2. Dally included the specification of a non
zero bottom velocity to provide a means, if desired, of matching the velocity in the interior
of the flow to a boundary layer edge velocity, after LonguetHiggins (1953).
The first term on the righthand side of Eq. (2.5) arises from the wave momentum
flux in the troughtocrest region. The use of linear wave theory raises a unique point about
the Dally (1980) model. Applying linear wave theory extends the limits of the problem to
the free surface. By doing so, Dally is able to include the effect of the radiation stress in
the trough to crest region as a shear stress on the upper limit of the model domain, without
applying a boundary condition to do so.
The work of Bdrek9i (1982) follows the same approach of collapsing the trough
tocrest region down to the MWL. Subsequent researchers have questioned the use of
linear wave theory and the application of a zero surface shear stress condition at the free
surface in the surf zone setting. The application of a zero surface shear stress at the
instantaneous free surface, when in fact the free surface is assumed to be coincident with
16
the MWL via linear theory, does create a degree of confusion, as evidenced by subsequent
discussion of this model (e.g., Svendsen, 1984a; Dally and Dean, 1986).
Dean (1995) employed a similar approach to the Daily (1980) model, but included
the possibility of a wind stress. Dean applied the horizontal momentum balance
(essentially Eq. (2.4) with the left hand side set equal to zero) with an applied surface shear
stress (the wind, r,) as one boundary conditions and a set bottom boundary velocity (=0).
The modified expression from Dean follows:
U(z) h [ 2T E 3(2 l() 1i + 3Q[ (2 1 (2.6)
pe ax 8 h 2 h 8 2h h
where the flow rate through the section is dictated by substituting for the water surface
slope in Eq. (2.4). In this way, Dean constrains the solution to include the bottom and
surface boundary conditions and conserve the flow through the vertical section. The first
term in Eq. (2.6) includes the wind stress, denoted as r,, as well as the wave momentum
flux in the troughtocrest region.
Hansen and Svendsen (1984), and Svendsen (1984b) employ a similar approach to
return flow modeling in which the net flow and bottom velocity are specified as boundary
conditions. In these instances the flow domain is truncated at the wave trough level, and
these investigators leave the troughlevel boundary condition unspecified. The discussion
in these models then focuses on the appropriate choice of the bottom velocity. Svendsen
et al. (1987) discuss the need to include the details of the flow near the bottom, where a
matching technique is applied to a boundary layer representation. Hansen and Svendsen
(1984) note that in several laboratory measurements, the measured bottom velocities seem
to agree very well with the magnitude of the waveinduced streaming velocity, but not in
17
direction. The direction of the measured bottom velocity in many return flow experiments
is offshore, not onshore as a streaming velocity argument would suggest (see Chapter 6).
Stive & Wind (1986) chose to adopt a different boundary condition in the
application of Eq. (2.4). The authors truncate the domain at the wave trough level and
argue that since the imbalance that drives the return flow is focused in the troughtocrest
region, it is more appropriate to specify the troughlevel shear stress, along with the
specification of the flow rate, rather than the bottom velocity. The details of the velocity
at the seabed, which is a product of the solution and is generally nonzero, are not
discussed.
Svendsen et al. (1987) revisited the original development of Svendsen (1984b) to
focus on the matching of the bottom boundary layer velocity and shear stress to the near
bottom velocity and shear stress in the body of the fluid. This matching produces a much
smaller eddy viscosity in the boundary layer, where a zero velocity at the bottom is
required. In a related effort, Putrevu and Svendsen (1993) extended the analysis of
Svendsen et al. (1987) to the shoaling region and showed improvement in the agreement
with experimental data by including a streaming velocity effect in the wave induced bottom
boundary layer.
Cox and Kobayashi (1997, 1998) used detailed measurements of the return flow
both in the boundary layer and in the body of the fluid to develop a logarithmic expression
for the velocity in the boundary layer. This solution was matched to a parabolic description
of the velocity in the body of the fluid. The authors also demonstrate that a quadratic
friction equation with a fitted friction factor was able to reasonably predict the resultant
timedependent bottom shear stress. The authors compare their developments to laboratory
data as well as field data (see Section 2.4, Return Flow Measurements, below).
18
2.3 Turbulent Eddy Viscosity
The turbulent eddy viscosity, e, is the diffusive mechanism in the present model.
In order to produce accurate simulations of the return flow, it is obviously necessary to
assign a value to the turbulent eddy viscosity. It will be shown in this work that the
application of a spatially variable turbulent eddy viscosity, particularly in the vertical
direction, produces significant improvement in the agreement between model simulations
and measured data. In this section a review of available literature regarding the selection
of eddy viscosity values is presented for both spatially uniform and variable cases. To
begin, the determination of an appropriate value to use in the uniform eddy viscosity case
is pursued.
The eddy viscosity is a property of the flow situation, which makes its
determination difficult. In reviewing the substantial amount of available literature, an
attempt is made to bracket the possible range of eddy viscosity values to apply in the
numerical model. Furthermore, an effort is made to find a physical basis for applying
spatially variable eddy viscosity values. Given this uncertainty, the possible variations in
eddy viscosity, both in a spatially uniform and spatially varying scheme, are investigated.
Hansen and Svendsen (1984) present a 1D analytic expression for the return flow
in which the eddy viscosity varies exponentially in the vertical direction to represent the
decay of turbulence from the surface down through the water column. The authors produce
expressions eof between O.005ch and 0.Olch, where c is the wave celerity and h is the local
water depth. A typical laboratory breaking wave condition of 0.2 m water depth translates
to eddy viscosity values in the range of e= 0.001 to 0.003 m2/s, assuming shallow water
conditions. Although the authors discuss only laboratory waves, application of the above
expressions to field conditions (1 m water depth), produces values in the range e= 0.016
19
to 0.031 m2/s. It is noted that by assuming shallow water conditions ( c = (gh)" ) in these
expressions, the horizontal variation in eddy viscosity would be expected to depend on h32.
Stive and Wind (1986) discuss the values of e used in both the Hansen and
Svendsen (1984) experiments and their own laboratory measurements. The authors
compute the value of the depthuniform eddy viscosity at each crossshore measuring
profile from the collected data. For the experiments of Hansen and Svendsen (1984), Stive
and Wind report values of between 0.0017 and 0.0008 m2/s in surf zone water depths
ranging between 14.5 cm and 8.6 cm, respectively. For their own experiments, Stive and
Wind report depthuniform eddy viscosity values of between 0.0026 to 0.0005 m2/s in surf
zone water depths ranging between 18.5 to 7.7 cm respectively. Stive and Wind discuss
an analytic approach to determining an appropriate value of Eusing an analogy to the eddy
viscosity in wake flows. The resulting expression, 6 = 102ch, is similar to that proposed
by Hansen and Svendsen (1984). Stive and Wind do not report any information regarding
the vertical variation of the eddy viscosity.
Svendsen et al. (1987) discuss the matching of the velocity profile in the bottom
boundary layer, in which eddy viscosity values are hypothesized to be significantly lower,
to the velocity profile in the body of the fluid. Using the original data reported by Hansen
and Svendsen (1984), Svendsen et al. propose values of e in the body of the fluid of
between approximately 0.001 and 0.003 m2/s for surf zone water depths between 14.5 cm
and 8.6 cm, respectively, in the laboratory. Svendsen et al. discuss the likelihood that the
eddy viscosity in the bottom boundary layer may be as much as 103 times smaller than the
appropriate value in the body of the fluid.
Haines and Sallenger (1994) discuss the application of horizontally variable eddy
viscosity fields in field data from Duck, North Carolina. The authors relate the horizontal
20
variation in the eddy viscosity to the dissipation of the turbulent bore formed by wave
breaking. The eddy viscosity is at its maximum value just landward of the first break point
(on the seaward face of the bar). The value decreases rapidly following breaking, reaching
a minimum equilibrium value before rising slightly at the shore break region. The
magnitudes of the surface eddy viscosity in the cross shore were determined from fitting
the measured field data at each horizontal measuring location. Inspection of the data
presented by Haines and Sallenger also indicate a strong dependence of the eddy viscosity
on h3/2. The range of fitted values of e found by the authors for this particular field
experiment range between 0.0055 and 0.075 m2/s. The authors employ a verticallyuniform
eddy viscosity at each horizontal location.
Reid (1957) describes the mean horizontal velocity profile in bounded and
unbounded channels subjected to surface wind stresses. By linking linearly increasing
mixing length descriptions of the turbulent shear stress from the surface and bottom
boundaries, a parabolic description of the mixing length is found. Using this relationship
in combination with a linear shear stress distribution in the vertical, Reid developed non
dimensional equations describing the velocity profile under a wide range of surface and
bottom stress conditions. The two pertinent equations are:
du du
T(z) = pL2 d (2.7)
dz dz
ko
L(z) = (z+zo)(D+z1 z) (2.8)
D
where ris the shear stress whose value at the bottom is described by a quadratic friction
expression, L is the mixing length for turbulent flow, ko is von Karman's dimensionless
constant (= 0.40), D is the water depth, and zo and z, are characteristic bottom and surface
21
roughness scales, respectively. Reid related the bottom shear stress to that of the surface
shear stress by the coefficient, m (generally, m I < 1 in the example presented by Reid).
Reid (1957) provides one of the only analytic means of validating the numerical
model for flows with nonuniform eddy viscosity. The mixing length, L, and the turbulent
eddy viscosity (used in the numerical model) can be related as follows:
E = L 2lu (2.9)
8z
The variation in the mixing length over depth and hence the eddy viscosity is
explained in Figure 2.1. In the figure, the analogy between the present situation and open
channel flow is studied. In the upper three cases, there is zero shear stress and no
generation of turbulence at the surface. In the first plot, the assumption of a depth uniform
shear stress (an incorrect but occasionally made assumption) and a linearly increasing
mixing length from the bottom leads to a linearly increasing eddy viscosity value from the
bottom upward via Eq. (2.9). It is noted that the use of a linearly increasing mixing length
argument is approximated from the Prandtl development of mixing length above a flat
surface in a semiinfinite fluid.
In the second plot of Figure 2.1, a linear distribution for the shear stress is applied
in which the shear stress varies linearly from the bottom value to zero at the surface. In this
instance the eddy viscosity produced is parabolic, with a maximum at middepth. In the
third plot, a parabolic distribution of the mixing length is applied, assuming that the mixing
length is limited at both boundaries (i.e., no longer a semiinfinite fluid approximation).
In this case the resultant eddy viscosity curve resembles a higher order behavior, with a
maximum that is shifted toward the bottom where the turbulence is being generated. The
fourth plot in Figure 2.1 is the situation modeled by Reid (1957) where turbulence is
0 (a)
0
.... c = f(z2
(b)
= f(z)
L = (z+zo)(h+zlz)
.E = f(Zn) (C)
0 0 0 V
= f(z) L = +zo)(h+zlz) E = f(zn) (d)
Figure 2.1 Illustration of the effect of the quadratic mixing length expression employed by
Reid (1957). The analogy to gravitydriven open channel flow is drawn to illustrate the
different distributions of eddy viscosity that are produced under different assumptions for
the distribution of shear stress and mixing length. The bottom figure depicts the situation
applicable to the present model. The assumption of uniform shear stress shown in plot (a)
is incorrect but is shown for comparison.
23
generated at both boundaries, and the turbulence generated at the surface by the wind stress
is much greater than the turbulence generated by the stress at the bottom.
The fourth example in Figure 2.1 is analogous to the present flow situation, in
which the surface stress generated by breaking waves is believed to be much greater than
the oppositelydirected bottom stress generated by the return flow. This distribution of the
shear stress, sketched in the fourth and final plot in Figure 2.1, produces a similar higher
order variation in the eddy viscosity as seen in the third plot, with the exception that the
maximum value of the eddy viscosity occurs closer to the surface, where the majority of
the turbulence is introduced from the breaking waves. This shift of the curve can also be
deduced from Eq. (2.8) in that the roughness scale at the surface due to the breaking waves
is expected to be much greater than the roughness scale along the bottom. B6rek9i (1982)
discussed the profile shape produced by the Reid mixing length model and investigated the
role of the roughness scales, zo and z,, in Eq. (2.8). Borekgi applied the wave height as an
upper limit on the surface roughness scale and plotted the resultant profiles.
The vertical variation of e sketched in Figure 2.1d is a maximum a short distance
from the surface (approximately 0.2 or 0.3h below the Mean Water Level). The value of
the eddy viscosity quickly decreases to a much lower value at the bottom, perhaps as much
as 102 or 103 times lower than the surface value (Svendsen et al., 1987). The significantly
lower value of eddy viscosity near the seabed acts to allow the fluid to essentially "slip"
over the bottom, acting much like a thin wave/current bottom boundary layer solution. In
the present model, no mathematics is incorporated to simulate a bottom boundary layer.
Thus only the vertical variation in the eddy viscosity can produce this phenomenon in the
model. The variation in also serves to produce much more uniform profiles in the body
of the fluid when compared to the spatially uniform case.
24
The vertical variation in eddy viscosity can in some ways be inferred from
measurements of turbulent kinetic energy. Stive (1984) presents vertical profiles of
turbulent energy under breaking waves in the laboratory. These profiles exhibit a shape
similar to the vertical variation in e found from Reid (1957). The data of Stive (1984)
indicated a peak value of turbulent kinetic energy occurring just below the Mean Water
Surface. Ting and Kirby (1994) present similar profiles of turbulent kinetic energy from
laboratory measurements.
In the horizontal, substitution of Eq. (2.7) into Eq. (2.9) produces a dependence of
e on L3/2 (the shear stress, r, varies linearly in the vertical). Since the length, L, is
dependent on the local water depth, h, the eddy viscosity is expected to depend on h3/2 as
well. This is consistent with the findings of Stive and Wind (1986) and Svendsen et al.
(1987). This is also the general finding when the eddy viscosity is considered in terms of
Froude scaling. Given the Froude scaling time dependence of L'2, the eddy viscosity,
which has dimensions of [L2/T], would be expected to vary as L31. Use of the depth as a
length scale is then an obvious choice.
Inspection of data regarding the variation of ein the vertical and horizontal leads
to several conclusions. The eddy viscosity appears to be strongly related to the water
depth, which can be considered a scale for the turbulent mixing length. Reid (1957)
presented a parabolic mixing length argument in which the mixing length reached its
maximum at middepth. Translation of this mixing length produces a vertical variation in
that depends on the shear stress at the boundaries and varies to a higher order, with the
peak value occurring well above middepth. For simplicity in fitting, and due to the
uncertainty in the predicted value near the region of zero shear stress, a third order
polynomial is fit to the Reid profile to investigate the vertical variation in 6.
25
The horizontal variation in e appears to be related to h3/, among other factors.
Various researchers have proposed values for the eddy viscosity at a particular depth that
take the form Bch, where B is a constant, typically between 0.005 and 0.01. The overall
magnitude of the eddy viscosity may be estimated from this simple expression. Using
typical values, assuming a laboratory depth of 0.1 m produces a typical eddy viscosity
value of 0.001 m2/s. Similarly, applying a field depth of 1.0 m in the expression yields a
value of = 0.03 m2/s. Haines and Sallenger (1994) proposed a horizontal variation in
eddy viscosity (depthuniform) based on turbulent bore dissipation. Their hypothesis was
borne out through field measurements.
2.4 Return Flow Measurements
A significant amount of data exist that describe the return flow under waves. This
section provides a introduction to the existing data and highlights the experimental data that
will be used in the present discussion to test the numerical model. Russell and Osorio
(1958) presented some of the first laboratory measurements of mean velocity profiles,
which were aimed primarily at investigating wave induced drift profiles over horizontal
bottoms and submerged bars. Drift profiles were measured in a Lagrangian manner by
timing neutrally buoyant particles over fixed distances.
More recent laboratory measurements of return flow profiles begin with the work
of Stive (1980), who collected velocity measurements under monochromatic, spilling
waves breaking on a smooth concrete slope of 1:40. Velocity profiles were measured by
Laser Doppler Anemometry (LDA) at a spacing of approximately 1 m across the profile
and included profiles seaward of, at, and landward of the break point. Stive focused on the
profiles well landward of the break point in the organized inner surf zone region.
26
Nadaoka and Kondoh (1982) present laboratory LDA velocity measurements of
return flows under spilling and plunging breakers measured at seven profiling stations. The
waves were monochromatic, breaking on a smooth wooden slope of 1:20. Okayasu and
Katayama (1992) present both monochromatic and irregular wave cases for waves breaking
on plane and barred planar laboratory beaches. Both beach profiles were constructed of
smooth planking of 1:20 slope. The barredbeach case consisted of a bar with 1:20 slopes
on both sides. No monochromatic waves were tested on the barred beach.
Ting and Kirby (1994) present velocity profile measurements for monochromatic,
spilling waves breaking on a smooth laboratory slope of 1:35. Cox and Kobayashi (1997)
present similar data collected in the same wave tank used by Ting and Kirby for
monochromatic, spilling waves breaking on an artificially roughened slope of 1:35.
One set of data are simulated with the numerical model to study the return flow
profile over a moveable bed profile in the laboratory. Roelvink and Reniers (1995) present
return flow data from several experiments conducted as part of the Delta Flume
Experiments conducted at the Delft Hydraulics facility in the Netherlands. At least 11
velocity profiles were collected for each experiment. Specific cases from these experiments
will be investigated and compared to the numerical model results. Field data from Duck,
North Carolina, will also be applied. Smith et al. (1992) present velocity measurements
from the DELILAH series of field experiments.
Four plane bed laboratory experiments will be simulated in detail in the present
work for purposes of model comparison. The periodic wave, plane slope experiment of
Nadaoka and Kondoh (1982), the random wave, barred slope experiment of Okayasu and
Katayama (1992), and the periodic wave, plane slope experiments of Ting and Kirby (1994)
and Cox and Kobayashi (1997) will be simulated in the numerical model. The Delta Flume
27
movable bed experiment will be investigated for cases of random waves over (relatively)
fixed but irregularly shaped profiles (Roelvink and Reniers, 1995). Lastly, the DELILAH
Duck field data set of Smith et al. (1992) will be simulated to study field conditions. These
six data sets were chosen primarily for their density of mean velocity measurements across
the entire profile.
Some experimental data in the literature were collected under reasonably similar
conditions. In these cases the differences and similarities in the results is investigated.
Figure 2.2 presents portions of the data collected by Ting and Kirby (1994) and Cox and
Kobayashi (1997). Both data were collected in the same tank under similar wave
conditions with the exception of a smooth bottom (Ting and Kirby) versus roughened
bottom. Inspection of the two experiments reveals that the two cases produced very similar
profiles overall, both in shape and magnitude, with the exception of the second profile of
Ting and Kirby. Some slight differences in the shape of the individual profiles is observed.
While the data of Ting & Kirby do not extend into the bottom boundary layer, as
do the data of Cox and Kobayashi, it is clear that if the older dataset extended through the
boundary layer the agreement (at least visually) would improve. The agreement between
these experiments is discussed only in terms of overall or visual similarity due primarily
to the discrepancies in measuring locations. The numerical model will simulate these two
experiments (and others). In these simulations a more quantitative measure will be applied
to determine what model parameters and boundary conditions are 'better' or 'worse.'
The agreement between the simulations and the experiments, or between the
simulations themselves, is the indicator of what input conditions produce the best answers.
The laboratory data certainly have various degrees of error that affect the result. This is
particularly true close to the surface, where the data collection techniques (such as LDA)
are affected by aeration and dropout of the signal above the trough level. Although these
data do exhibit some scatter and some individual problems, the quality of the return flow
data available in the literature is generally considered quite good and is used to validate the
present model.
0)
.1'
f.u
a)
5 6 7 8 9
Distance from beach toe (m)
a)
.4,
SE
) j
.0
4ZU
5 6 7 8 9
Distance from beach toe (m)
Figure 2.2 Comparison of mean horizontal velocity profiles as reported by Cox and
Kobayashi (1997) and Ting and Kirby (1994). These two datasets were collected under
very similar conditions in the same wave tank. The primary difference in the experiments
is the bottom surface conditions (roughened versus smooth bottom). The wave conditions
in the Cox and Kobayashi experiments includes a 2.2 s wave period (periodic waves) with
a 17.1 cm breaking wave height over a rough bottom. The Ting and Kirby data set used
a 2.0 s wave period and had a breaking wave height of 16.3 cm over a smooth bottom.
Both conditions are classified as spilling breakers. Not all profile lines are shown.
CHAPTER 3
DERIVATION OF GOVERNING EQUATION
AND NUMERICAL MODEL DEVELOPMENT
This chapter presents the derivation of the governing equation and its subsequent
discretization for the present numerical model of the return flow in the surf zone.
Developments for both onedimensional (lD) and twodimensional (2D) versions of the
model are presented (both options are contained in the same program). For the 2D model,
the governing equation is modified to transform the arbitrary domain based on a real beach
profile shape to a rectangular, orthogonal domain for solution. Details of the
transformation are provided. The performance of the model depends strongly on the
selection of boundary conditions. The range of boundary conditions to be applied is
discussed at length in Chapter 4. The derivation of the governing equation and its
discretization are quite lengthy, particularly for the 2D case. The full discretization,
including a listing of all coefficients used, is contained in Appendix A. A complete listing
of the computer code written in FORTRAN90 is found in Appendix B.
The modeled flow field extends vertically from the mean water surface to the
seabed and horizontally from the shoreline out well beyond the surf zone. The governing
equation is obtained through crossdifferentiating the twodimensional NavierStokes
equations for turbulent flow (e.g. the Reynolds equations). The resulting thirdorder partial
differential equation (PDE) in terms of velocity is then recast in terms of the stream
function, ultimately producing a fourthorder PDE that retains the generalization of time
30
dependence and a nonuniform eddy viscosity field. This turbulent eddy viscosity, e, is
considered a variable quantity in both the crossshore (x) and vertical (z) directions.
This development, simplified by the assumption of uniform eddy viscosity and
steady state, produces the wellknown biharmonic equation. The biharmonic equation can
also describe creeping flows of constant molecular viscosity. In the present application the
addition of time dependence and variable eddy viscosity significantly alters the governing
equation, making analytic solutions obtainable for only the simplest of geometries and
viscosity fields and thus prompting a numerical approach to the problem.
3.1 Derivation of Governing Equation
To begin, consider the twodimensional equations of motion given in Equations
(3.1) and (3.2) and defined as shown in Figure 3.1. Dally (1980), among others, showed
that by time averaging Eqs. (3.1) and (3.2) over a long period, neglecting the acceleration
terms on the left hand side of the equations, and crossdifferentiating the equations, the
biharmonic equation ( V4 f = 0) for the mean flow can be obtained by assuming uniform
viscosity (Dally simply replaced the molecular viscosity with the eddy viscosity).
au au 2 auw 1 ap 1 dx x
+ + = ( ( + Z) (3.1)
at ax 9z p dx p dx 9z
aw + w + gZ + + )( (3.2)
at dx dz p dz p dx dz
direction of
wave propagation
VSWLx
still water depth, h mean water depth, h+i
seabed
Figure 3.1 Definition sketch for model development (setdown of the MWL shown).
The shear stress terms in Eqs. (3.1) and (3.2) are defined as follows:
6u 6w
r = T,. = PV(a + a
az ax
6u
rx = 2pv( )
ax
6w
,z = 2pv(a)
dz
(3.3)
(3.4)
(3.5)
Dally (1980) derived the biharmonic equation by neglecting acceleration terms,
directly substituting the eddy viscosity for the molecular viscosity in Eqs. (3.3) through
(3.5) and crossdifferentiating the time average of Eqs. (3.1) and (3.2). This approach
produces the same PDE, but is somewhat misleading because the contribution to the mean
flow field from the turbulent eddy viscosity arises from the nonlinear advective acceleration
terms in Eqs. (3.1) and (3.2), the terms that Dally considered negligible. A more rigorous
approach, applying the classic Reynolds equations for turbulent flow, traces the
introduction of the turbulent eddy viscosity into the problem and provides a clearer picture
32
of the individual components that contribute to the variation of the eddy viscosity across
the surf zone. The reader is referred to Schlichting (1968) for more information regarding
the Reynolds equations.
The derivation of the Reynolds equations for turbulent flow from the NavierStokes
equations (Eqs. (3.1) and (3.2)) requires the decomposition of the velocity terms u and w
into steady and fluctuating components:
u = U + u + u' w = W + w + w' (3.6)
where U and W represent the mean flow components in the horizontal and vertical
directions, respectively, u and iv represent the oscillating wave orbital motions, and u 'and
w 'represent the fluctuating turbulent components of the flow. Substitution of these
decompositions into Eqs. (3.1) and (3.2) and subsequent time averaging (denoted by the
overbar) over a representative wave period, such as the peak or modal period, yields:
8(U+17+u') a
O+ [U2+ 2U + 2Uu' + 2u' + a2 +
at ax
+ [ UW + U + + Uw + iW + u'W + aw + w + u + u'w ]
z __(3.7)
1 dji a a8
[ 2v(U+,+u') ]
p x ox ax
+ [ v( (U+a+u ) + (W+w+w') )]
dz dz ax
9(W+W +w ') W a 1 
at + [ UW +UVw +Uw +7W +u W +aw + + u7W / +u '/ ]
at 9x
+ [ W + 2W" + 2Ww/ + 2ww + 42+w2 g (3.8)
9z p dz
a a a a a /
+ x[ v( (U+U+u') + (W+W+w ) )] + [ 2v(W+w+w ) ]
9x 9z 8x 9z 9z
33
Many of the product terms in Eqs. (3.7) and (3.8) time average to zero. Following
the classic eddy viscosity model first introduced by Boussinesq, some of the nonlinear
terms are related to the mean flow as "apparent" shear stresses (denoted by the primes),
which are similar in form to the viscous shear stresses in Eqs. (3.3) through (3.5):
/ / a  + aw
x = TZ, = pu'w pa paw' puw = pe( + ) (3.9)
az ax
u = p 2 p2 2pu = Zp( ) (3.10)
ox
/ /2 aw
tZZ =pw p2 2pw = 2p( ) (3.11)
az
These definitions of the apparent shear stresses differ slightly from the classic
definitions due to the inclusion of the wave orbital motion terms. Careful inspection of
Eqs. (3.9) through (3.11) indicates which terms have been treated as zero after time
averaging and which terms have been retained in Eqs (3.12) and (3.13). From Eq. (3.7) the
last three components of the xgradient term have been retained, along with the last four
components of the zgradient term (Eq. (3.12)). From Eq. (3.8) the last four components
of the xgradient term and the last three components of the zgradient term have been
retained (Eq. (3.13). The remaining components have been taken as zero after time
averaging. This assumption neglects the acceleration of the mean flow (terms such as the
time average of U2) and seems reasonable in the inner surf zone, and even more reasonable
outside the surf zone. In the outer surf zone, or transition region, this assumption is more
suspect, but for present purposes these terms are neglected.
34
The terms retained in Eqs. (3.9) through (3.11) represent various contributions to
the apparent stress from turbulence and wave motion. The first terms retained, those
containing only u 'and/or w represent the classic Reynolds normal and shear stresses in
turbulent flow. The terms involving only u and/or iv represent wave shear and normal
stresses (the reader will note from the discussion in Chapter 2 that the gradient in time
averaged values of u and ii represent the radiation stress (LonguetHiggins and Stewart,
1962)). In the present development, no details of wave induced flows are included and
these terms are incorporated into the turbulent eddy viscosity representation.
The nonlinear terms may or may not be significant in various parts of the surf zone.
Inside the surf zone, it is assumed that the turbulent Reynolds shear stresses are the primary
source of turbulent diffusion, generated by the breaking waves. Outside the surf zone, the
importance of the Reynolds stresses is unclear and the wave shear stresses may play a
greater role. These terms were investigated by Rivero and Arcilla (1995), who
demonstrated that the timeaverage of quantities such as ui~ is not zero for waves shoaling
on sloping beds. The magnitude of these wave stresses is not well known, however, their
inclusion in the model provides a plausible mechanism for applying turbulent diffusion
seaward of the break point.
The cross product terms, such as the time average of u 'ii, represent the possible
stresses arising from the transport of turbulence by the wave orbital motion. Again, the
importance of these terms relative to the Reynolds stresses is unclear, but they are included
in the apparent shear stress terms for completeness. It is important to realize that
quantifying the eddy viscosity is a difficult task, given that eis both a property of the fluid
and the flow. However, the need to quantify each of the nonlinear terms in Eqs. (3.7) and
(3.8) is relieved by grouping them into the eddy viscosity form of Eqs. (3.9) through (3.11).
35
Substituting the apparent stresses with the viscous shear stresses into Eqs. (3.7) and
(3.8) and dropping the terms that equal zero after time averaging over a representative wave
period produces the Reynolds equations for the mean flow:
U 1 9pi 1 9 9U aU
+ [ pv + 2pe ]
at p ax p ax ax ax
(3.12)
1 au 9U 6w
+ [ pv + PE( + ) ]
p z az az ax
1w 1 1 [ 168 W + U 6W
S=g + 1[ p + pE(p + )]
at p az p ax ax P z ax
(3.13)
+ [ pv + 2pE ]i
p z 6z 6z
In general, the viscous terms in Eqs. (3.12) and (3.13) are much smaller than the
turbulent eddy viscosity terms (v << e). These terms will be neglected. Cross
differentiation of the two equations produces the following third order PDE, in which the
pressure terms have canceled and g is assumed to be invariant in x:
6U 6U 6W
w 2(2E) 6 2(E ) (E )
a 9U aw a x 9z ax
( )=
at 9z ax 9xaz 8z2 9z2
(3.14)
82(e ) Eaw) (2e)
9z ax 9z
+ z
9x2 ax2 axaz
A stream function notation is adopted for the mean flow such that:
U W a (3.15)
9z ax
which modifies Eq. (3.14) to read:
a( 2 ) a 2 a2(e1
a2(E ) a (E
a 8)2 qj a2 aZ 2 aX 2
a( + 4j) = xz 2 + x2
dt z 2 x 2 axaz az2 9z2
(3.16)
a2a a2 2 a2) a2( )
(E) a2(2E 2E
8z 2 x2 axaz
ax2 8x2 axaz
It is recognized that Equation (3.16) is very similar to the equation describing the transport
of vorticity due to turbulent diffusion. The term on the left hand side of Eq. (3.16) is the
time derivative of the Laplacian of /r, which is equal to c, where c represents the vorticity
in the fluid ( 7V2' = w).
The application of the 2D model will be presented in terms of uniform and non
uniform eddy viscosity. With uniform eddy viscosity, Eq (3.16) is simplified subscriptt
notation for all derivatives is adopted at this point):
(IV, + IlJ), = E( ixl + 24ijx,, + 1.z. ) (3.17)
and the additional assumption of steady flow (over a long period) produces the biharmonic
equation:
4j. + 2iJ .il + 4 = 0 (3.18)
As mentioned, Equation (3.18) is commonly known as the biharmonic equation, and can
be used to describe the creeping flow of viscous fluids. While the analogy to creeping flow
37
is not particularly relevant, the steady state version of the equation developed herein is
identical to the biharmonic equation if the turbulent eddy viscosity is uniform, thus the few
existing analytical solutions for creeping flow situations can be used to test the present
model in simple geometries.
Dally (1980) presents an analytic series solution to Eq. (3.18) and attempts to model
the horizontal fluid momentum resulting from a surface shear stress applied over a finite
length. The results of the analysis by Dally, while limited, provide the impetus for the
development of the present approach to return flow modeling. One distinct advantage of
this approach is that crossdifferentiation has removed the necessity to describe the pressure
variation in the flow. Inclusion of a spatially variable eddy viscosity field and the
numerical solution provides the opportunity to model the return flow more realistically.
The partial derivatives in Eq. (3.16) are expanded, maintaining that e= e(x,z). Eq.
(3.19) represents the partial differential equation to be modeled in the xz domain:
..zt + El = EIJ + 2E*Ixxzz + E*.
+ 2E, E + 2E + 2Ex, + 2e. vz (3.19)
+ (e e z)xx + 4Exi ix + (EZ Ex)iJ
3.2 Development of OneDimensional Model
For a simple investigation of the vertical variation of the horizontal velocity in the
return flow profile, Eq. (3.19) is simplified and solved numerically. The 1D partial
differential equation in z is:
iJl = el + 2 EI* + eZ1 E.
(3.20)
38
Eq. (3.19) is discretized using central difference approximations for the spatial gradients
in the body of the model domain (to second order in Az). The time derivative is modeled
as a forward difference approximation, and is represented in a variable, semiimplicit
fashion by the inclusion of the timeweighting term, 0. The model domain is illustrated
in Figure 3.2. The fourthorder derivative in z yields five terms after discretization; upon
collection of all f terms, the following finite difference equation results, involving five
different points in the domain at two successive times, where n is the time axis (n is the
time at the previous time step, n+1 is at the new time step). In the 1D model the vertical
origin is taken at the seabed.
[, 21); + 1, + 11 =
SAt [(e AzEl,_2 + (4E, + 2AzE +(Az)2 Z).Ii 1
(az)2
+ (6, 2(Az)2 E,) + (4E 2A. + (z)2 +1 (/ AEz +2 +1 (3.21)
+ (10) [( AzE4_;2 + (4e1 + 2Aze + (Az) 2E. ,_1
(Az)2
+ (6E6 2(A,)2 E)IJ,+ (4E 2AzE6 + (Az)2 E)l,+1 + (E AZE6),,I,]n
+ [*y 1 2,1 + *,, ,]"
The value of Ocan be modified to run the model in a fully implicit fashion (0=1),
a fully explicit fashion (0= 0), or any intermediate value. Setting 8= 0.5 is equivalent to
using the well known scheme introduced by Crank and Nicholson (1947). The variable
implicitness is included to alleviate the time step restrictions of the process modeled, which
is similar to modeling the diffusion ofvorticity. In a fully explicit scheme, the time step
is limited to the following value:
At (AZ)2 (3.22)
26
z
direction of
wave propagation
applied surface stress, T
grid point N (w = Q) j+2
j+1
Mean water
j1 depth, h+Ti
applied bottom velocity ub) j2
or applied shear stress ('b) grid point 1 (9 = 0)
seabed
Figure 3.2 Definition sketch for onedimensional model. The figure illustrates the four
applied boundary conditions and fivepoint solution cell.
which is quite restrictive, resulting in long run times in a 2D model. The use of a semi
implicit or fully implicit scheme with 8> 0.5 allows for a much larger time step to be
employed, one closer to the order of the physical forcing mechanism in the model, with
unconditional stability and convergence (Smith, 1985). It does not, however, guarantee
convergence to the correct answer for any time step value. In the present case, the physical
time step is controlled by the eddy viscosity, not the wave period. The appropriate choice
of time step is discussed in the model validation portion of this work (Chapter 5).
It is noted that higher values of 0 introduce artificial numerical diffusion
(smoothing) into the modeling process. For this reason the value of 0 is kept at
approximately 0.5 to gain the benefit of the relaxed time step requirement but not introduce
excessive numerical diffusion. Since it is not the goal of this work to exercise the details
of numerical modeling, the reader is referred to any number of numerical modeling texts,
such as Hoffman (1989) and Smith (1985). In the models, Ois set to 0.51, and the model
is exercised at a range of time steps to assure consistency in the final answers.
40
3.2.1 Boundary Conditions for OneDimensional Model
Application of Eq. (3.21) requires the specification of four boundary conditions
(Figure 3.2). As with any numerical model, the selection and implementation of the proper
boundary conditions represents a substantial fraction of the modeling task. Recent
literature regarding return flow modeling is dominated by the discussion of the proper
choice of boundary conditions, and the present model is no exception. The detailed
selection of boundary condition values is found in Chapter 4, and the mechanics of
applying the boundary conditions in the model are found in Appendix A. The boundary
conditions applied herein are those described by LonguetHiggins (1953) for establishing
shear stress, velocity, and stream function value conditions at the boundaries.
In the 1D model, four boundary conditions and one initial condition are required
to produce the desired vertical profile of the stream function. Two of the boundary
conditions are simply the specification of the value of the stream function at the seabed and
at the Mean Water Level (MWL). The difference in these two conditions dictates the flow
rate through the model domain, thus assuring that flow is conserved based on the computed
volumetric transport, Q, of the wave being modeled. Generally, the value of fis set to zero
at the seabed and the surface value is specified as / = Q, where Q is the forward
volumetric transport of the wave and Q is the flowrate through the model domain (this
results in positive values, generally producing offshore return flows based on Eq. 3.15).
The third condition is the application of a surface shear stress caused by the passage
of the waves. Outside the surf zone this stress may be the result of the gradient in wave
height due to shoaling, while inside the surf zone the shear stress is created by the wave
height gradient caused by breaking. In terms of fr, the shear stress condition is equivalent
to a secondderivative condition at the surface. The final boundary condition is applied at
41
the seabed and the option is given of specifying either the tangential bottom velocity or
tangential shear stress condition. The bottom velocity is specified by the first derivative
of /with respect to z.
The initial condition for the model is generally that of a quiescent condition, or
"cold start," where the values of 0fin the body of the domain are set to zero. The surface
boundary conditions then provide the driving force to the problem. As an alternative, an
initial estimate of the vertical profile can be provided, either from an assumed profile or
from the results of another simulation. This would provide a "hot start" to the model and
lead to a steadystate solution in fewer time steps. A steady state solution in the 1D model
is quickly determined and provision of a "hotstart" condition is not needed. In the time
dependent 2D cases, a good initial estimate is beneficial.
3.2.2 Solution Scheme for OneDimensional Model
The stream function values across the model domain at the n+1 time step are solved
simultaneously after creating a 5band diagonal matrix ofN2 equations using Eq. (3.21).
Because the values of fkat the surface and bottom are given by the boundary conditions,
the matrix and vectors have the dimension N2. The boundary conditions are input directly
to the matrices and vectors. The program employs two subroutines from the Numerical
Recipes suite of routines (Press et al., 1986). The first routine (BANDEC) performs a
lowerupper (LU) decomposition on the matrix. These matrices are then sent to the second
subroutine (BANBKS) which performs a back substitution operation to solve for the vector
of values. In the 1D model, only one performance of each routine is required at each
time step. Other routines are available to perform the same operations, such as the IMSL
routine LSARB.
42
In timedependent simulations, the values of the stream function at each point are
averaged over many waves to produce the final mean stream function profile, such as to
simulate the return flow produced by a random wave field. Intermediate values of the
solutions are stored in order to observe the time scales at which the profile evolves due to
changes in the input wave parameters. The computer program proceeds to process the
values of stream function to produce as output a list of the horizontal velocity and the shear
stress, both as the mean value and at any desired intermediate time steps. An example of
the model output for uniform eddy viscosity and monochromatic waves is shown in Figure
3.3 for input values matching those used in the return flow discussion of Dean (1995). The
computer program is provided in Appendix B.
shear stress (N/m2)
40 30 20 10 0 10

0.8 
0.6 I /
0.4 
0.2 
0.4 0.3 0.2 0.1 0 0.1
velocity (m/s, onshore positive)
Figure 3.3 Example output from onedimensional return flow model. The solid lines
present the velocity and shear stress distributions assuming a zero bottom velocity
boundary condition. The broken lines present the velocity and shear stress distributions
assuming a bottom shear stress boundary condition defined to be equal in magnitude and
opposite in sign of the applied surface stress (m = 1.0). Input values for these examples
were chosen to match those used by Dean, 1995, and the zero velocity boundary condition
case matches Dean's example (H= 0.78 m, h = 1.0 m, e= 0.04 m2/s, z; = 7.9 m2/s).
43
3.3 Development of TwoDimensional Return Flow Model
Building on the 1D model, Eq. (3.19) is revisited and applied in the 2DV
nearshore region. The same boundary conditions used in the 1D model apply to the 2D
model, with the addition of four more boundary conditions in the horizontal direction, i.e.
sidewall boundary conditions. Since Eq. (3.19) is written for the physical, xz coordinate
system, a coordinate stretching scheme is implemented to allow for the modeling of
arbitrary beach profiles. This coordinate stretching simply nondimensionalizes the vertical
coordinate, but lengthens the derivation substantially. For clarity, only the pertinent
equations are given in this chapter. Appendix A contains the full development.
3.3.1 Transformation of the Partial Differential Equation
Figure 3.4 illustrates the procedure used to transform the physical geometry of the
nearshore region into an orthogonal computational domain. Following Hoffman (1989) the
vertical coordinate is nondimensionalized by the local water depth (Eq. (3.23)). This
simple gridstretching technique allows for the modeling of any profile (Figure 3.4). As
contrasted to the 1D model, the 2D model will take as the vertical origin the Still Water
Level to maintain a horizontal datum since setup and setdown effects will be included, thus
grid point #1 corresponds to the surface and grid point #Nj applies to the seabed.
rl = z/h (3.23)
Transformation of Eq. (3.19) from the xz domain to the e]q domain simply
involves repeated applications of the chain rule for derivatives, as demonstrated below:
+ ll (3.24)
az az a 8z ar
0
. 4
E
12
6
50 100 150 200
xaxis, onshore distance (m, typical)
I
250
.0
4'
1 onshore gridpoint # N
5axis, onshore gridpoint #
Figure 3.4 Definition sketch of coordinate transformation between the physical, xz,
domain and the computational, er7, domain.
The approach throughout this derivation will be to separate the timedependent
terms from those that comprise the steady state solution. This produces several more terms
and coefficients but makes the operation of the model much clearer. Each term in frin Eq.
(3.19) is expanded; including the eddy viscosity terms. Carrying the chain rule through the
fourthorder PDE produces additional lower order derivatives of rdue to the coordinate
45
transformation required. The transformed PDE from Eq. (3.19) is:
[ Bly + B2*E + B3,,, + B4*E + B5i ]
Al_ + A2Ji + A39* + A4~n + A5I111 (3.25)
+ A6*EEE + A7%Ij + A8+* + A9gPT
+ AIOlIV + A7lllJ + A12*Tl + A13*E + A14*
where the left hand side represents the time dependent terms and the right hand side
represents the steady state solution. Each of the 19 coefficients in Eq. (3.25) contains the
transformation details of each derivative. As an example, one of the more lengthy
coefficients is given below; all the coefficients are defined in Appendix A. The values of
Ex, 9 E,, e _, and c that appear in the coefficients of Eq. (3.25) also must be expanded to
account for the coordinate transformation via Eqs. (3.23) and (3.24). The value of
coefficient A8 from Eq. (3.24) is given to demonstrate the number of terms involved:
A8 = [ e(12,xxx + 6,r 2 + 2SxT2 + 4xxrTI + 8xTrix
+ 8x% r + 8 1xr + 4lx2rl2 + 2 +r. 2 + 6r 1 2 + 12lzz) (3.26)
+ Ex(6(x2 + 2Exr12 + 4r,1r) + E(4 + 2,2 + 6+r2) ]
3.3.2 Discretization of the TwoDimensional PDE
Eq. (3.25) is now discretized for the numerical solution of fon the interior points
of the grid (recall that the values of the stream function along the edges are specified as
boundary conditions, as in the 1D model). Figure 3.5 compares the discretized solution
cell for the xz domain versus the solution cell for the transformed e? domain. The
transformation process adds an additional eight points to the finitedifference equation.
Eq. (3.25) is discretized using central difference spatial approximations of fr with
a forward difference approximation in time. As in the 1D model, the equation includes
z
j2
qF
j1
j+
j+1
j+2 I 
i2 i1 i i+1 i+2
x
V
j2
j1
j
j+1
j+2
i2 i1 i i+1 i+2
Figure 3.5 Solution cells for the fourthorder finite difference equation in the xz
domain and the transformed { j7 domain. The five circled points on column i represent
the solution points in the set of simultaneous equations. Rowj=1 corresponds to the
surface row.
a variable degree of implicitness in the solution, represented by the value of 0. Again,
efforts are made to separate those terms arising from the time derivative from those arising
from the spatial derivatives. In a somewhat schematic way, the equation is discretized as
follows:
[ time dep. terms ]n1 [ time dep. terms ]" =
OAt[ s.s. solution terms ]"1 + (1O)At[ s.s. solution terms ]n
(3.27)
In Eq. (3.27) the time dependent terms on the left hand side are the discretized
forms of the left hand side of Eq. (3.25), while the steady state (s.s.) solution terms are the
discretized forms of the right hand side of Eq. (3.25). The full, discretized equation is
47
given by Eq. (3.28). As written, Eq. (3.28) would produce 21 solution points at the n+1
time step, as opposed to the previous five points in the 1D model. This presents a problem
in establishing a spatially consistent and still manageable matrix of coefficients to solve for
simultaneously. This problem is relieved by moving any points that do not reside on the
ith column to the right hand side and solving the system of equations iteratively at each time
step (refer to Figure 3.5). Appendix A contains the definitions and rationale behind the
assignment of coefficients. A typical steadystate coefficient, corresponding to the [f,i]
point in the solution cell in Eq. (3.28), is presented in Eq. (3.29):
[ DR11P1, .1 + DR2*1,._ + DR3ij,1 1
+DL11 J, + DL2 ,i + DL31*1+, i
+ DR4I_1, i+ + DR511, ,+ + DR6a+1,1 i 1+1
[ DRlj1, ,1 + DR21, ,_, + DR3,+1,,1
+DLll.1, i + DL2n, + DL31*+1,
+DR4*_1, + DRS5i. ,t + DR61iJ+I, ,1 1
= OAt[ CRII_1, 12 + CR21j1, 2 + CR3lJ+1, 12
+ CR4J,2_,_1 + CRI5_, ,_1 + CR6ir,i, + CR71j+,, 1 + CR81Ir/+2, (3.28)
+CL1 _2, + CL211, + CL31P, + CLI,11 + CLS51,+2,
+ CR91_2,i+1 + CR10,_1, ,+1 + CR111 ,,+1 + CR12,+1'+1 + CR13*,+2, +1
+ CR141j1, 1+2 + CR15) 1+2 + CR161j+1, t12 ]n+1
+ (1 O)At[ CRl_J., 12 + CR2, 12 + CR31i+1, 12
+ CR4*12 ,_ + CRSlj l I + CR6j 1 + CR71j+l 1 + CR815+2, 1
+CL11I12, + CL2j_1, i + CL3, t + CL,1+, + CLS1+2, i
+ CR912, i+1 + CRIO_1, t+1 + CR111 ,, +1 + CR121j, +1, + CR13 1+2, i+1
+ CR141111l, +2 + CR151J 1,+2 + CR16*1+1,+2 1]
(A7)4 + A12 (A)4 2A12(ar 2 (3.29)
CL3 = [ 6 (A) + 4A3 + 6A5 2AO1 2A12(Ar)2 ] (3.29)
(Ag)4 (Ag2 (A)2
48
3.3.3 Boundary Conditions for TwoDimensional Model
For the twodimensional model, the boundary conditions are essentially the same
as in the onedimensional case. The value of the stream function along the boundaries is
specified, occupying four of the eight boundary conditions needed. Along the bottom and
side boundaries, Vfis generally set to zero, although inflow or outflow could be specified
in the model by changing the value of r from one grid point to the next (refer to Eq.
(3.15)). The value of f along the surface is specified as before by the modeling of the
wave transformation across the surf zone. At each grid point the value of the volumetric
transport created by the wave passage is used to establish f~at the surface.
The remaining two boundary conditions at the seabed and the surface are also
applied in a fashion similar to the 1D model case. The two remaining conditions along the
sidewall are specified as velocity conditions, with the velocity generally set to zero.
However, in the cases of both velocity and shear stress specifications, the transformation
of coordinates used to achieve a rectangular model domain complicates their application.
As an example, to specify the tangential bottom velocity along a sloping bottom,
the continuity equation must be invoked, along with the chain rule to translate the sloping
bottom to the rectangular solution domain:
U Wh = t, *Ah = Ub h + 1 (3.30)
(3.30)
which when transformed, discretized and solved for the outlying point just below the
seabed produces:
= ()h1 M )1 + + ,(3.3
., ( + ( + A, r5 "
+ A (3.31)
+ +T .gr ,.'+ +
hTl, + 1}, J
49
As the reader might imagine, the construction of a proper boundary condition that
involves a second derivative of fis substantially more involved. Further difficulties arise
at the comers of the model domain, where two separate boundary conditions share common
points in the domain that must be solved for simultaneously. The reader is referred to
Appendix A for the full development of the boundary conditions in the transformed
coordinate system.
3.3.4 Solution Scheme for TwoDimensional Model
The solution scheme for the 2D model is similar to that used in the 1D case. A
simultaneous set of Ni2 equations is solved for the interior values of f at each vertical
section (Figure 3.5). The same subroutines used in the 1D model are employed to solve
the system of equations (BANDEC and BANBKS, Press et al., 1986). This process is
swept horizontally from the offshore boundary to the shoreline for each time step.
As discussed previously, the method only solves for the five points on the ith column
at the n+1 time step. In the complete formulation, there would be 21 points to solve for
(Eq. (3.28)). This presents difficulties in creating a spatially consistent coefficient matrix.
As a concession, the offcolumn points, those with R in their coefficient labels, are moved
to the right hand side of the equation and their respective (values are treated as known.
The points moved from the solution side to the known side are then updated in a series of
iterations for each time step. This entails tracking three different values of fkat each point
in the domain, the old value, the new intermediate value (assumed to be known from the
previous iteration at each time step), and the newly solvedfor value at each iteration.
The passage of waves across the profile represents the driving force in the model
and is entirely described by the boundary conditions (no knowledge of the effect of waves
50
was incorporated in the derivation of the governing equation other than time averaging).
The model begins by determining the setup and setdown of the MWL and the wave height
change across the profile for a chosen typical wave height, providing the initial surface
boundary condition (the height used is generally either the rootmeansquare (RMS) wave
height or the significant wave height, adjusted with an appropriate factor). The model runs
through a series of time steps comprising the duration of each wave period (generally the
modal period for irregular waves), computing the flow field at each step.
The model then changes the surface boundary conditions based on the next
randomly chosen wave height, and runs through the same number of time steps. The
change in the value of fat the surface is assumed to occur over a single time step. This
is clearly erroneous, however, no data were found to describe the lag in flow rate in a surf
zone of changing wave conditions. Obviously the narrower the surf zone width the faster
the return flow rate will adjust to a change in condition. This is thought to be especially
true in a laboratory setting, where the adjustment in setup and return flow is rapid.
The setup and setdown remain fixed throughout the run, thus there is no grid
adjustment or wetting and drying of cells in the model. The model continually updates the
average value of the stream function at each point in the model and periodically can record
the instantaneous values in order to study the time evolution of the flow field. At each time
step, the model compares the updated flow field to the previous flow field after each
iteration and considers the maximum change in the value of fas a cutoff criteria. Figure
3.6 illustrates the operational steps of the numerical model. Validation of the numerical
model operation is discussed in Chapter 5.
PROGRAM FLOWCHART
Read initial wave parameters
and discretize beach profile
Transform domain to orthogonal
geometry, compute spatial derivatives
Perform initial wave transformation,
compute flow rate, mean water
surface elevation, surface shear stress
at each horizontal location
Compute turbulent eddy viscosity field
Compute transformed coefficients
in governing equation and create
coefficient matrix
Ses Update flow
Monochromatic waves? field, ,
no
choose wave
N= N + 1 height, compute
new Q,rs (AsVmax< n
tolerance?
Update flow
yes
esono
t = T?
yes
nN = Nwaves
yes
Compute final products:
y, velocity, shear stress fields
and profiles
AEB
Figure 3.6 Flowchart for twodimensional numerical model.
CHAPTER 4
BOUNDARY CONDITIONS AND
TRANSITION REGION TREATMENT
This chapter describes the selection of boundary conditions for use in both the one
dimensional and twodimensional models. As with any return flow model, the successful
prediction of the mean flow field requires careful selection of these conditions. In the
present model, this issue is increasingly important, since the use of the fourth order PDE
in V/requires the application of four boundary conditions in each direction, as opposed to
the customary two conditions applied in the simpler equation for U. Chapter 2 introduced
the literature regarding the development of return flow models and their boundary
conditions, and Appendix A details the mechanics of applying the conditions in the
numerical model.
In this chapter, the options and resultant values applied for each boundary condition
are discussed. In many cases, there is more than one option available for each condition.
These options will be presented herein and the merits of each will be demonstrated in
following chapters. Further, the chosen values are expected to differ significantly between
the region seaward of the break point and the surf zone. The applicable mechanisms and
equations employed to determine the conditions will also differ. The difficulty of matching
conditions across the entire profile arises in most models of the return flow. The details of
how a wave transforms in shape, transport, and momentum across the transition region
between the break point and the beginning of the inner surf zone, is an area of active
53
interest. Herein, two methods are applied to smoothly connect the region seaward of the
break point to the inner surf zone while maintaining conservation of momentum.
4.1 Wave Height Transformation
In order to determine the two surface boundary conditions it is necessary to describe
the wave height across the entire beach profile. This description includes the shoaling
region from the offshore boundary to the first break point, the breaker zone, and any
intermediate wave reformation and secondary breaker zones within the surf zone.
To begin, the wave height at the offshore boundary is determined. This height can
be given directly if known, or the corresponding offshore wave height can be given. In
either case, the wave height at breaking, if known, is matched by adjusting the input wave
condition. The wave period is provided by the user. In the offshore region, the wave is
marched across the model grid toward the shoreline until breaking is achieved, as
determined by a usersupplied breaking criteria (such as x= H/h = 0.8). The waves are
shoaled using Stream Function Wave Theory (SFWT, Dean, 1974). This selection was
made in order to produce steeper gradients in the wave height near breaking, which
provided better overall agreement with the laboratory measurements discussed in this study.
Furthermore, the nonlinear stream function theory was chosen over linear shoaling because
it better approximates the wave profile in the shoaling region, producing narrower waves
with higher peaks and flatter troughs (and ultimately lower and more realistic values of
volumetric transport and radiation stress, Sx).
Stream function theory is a numerical theory, and does not produce simple analytic
expressions for shoaling calculations. For this reason, the graphs and lookup tables
produced by Dean (1974) were digitized into an interpolation table in the numerical model
54
to provide the necessary shoaling information. Figure 4.1 provides a reproduction of the
stream function shoaling table of Dean (1974) for the case of normally incident waves, as
would be the case for this 2D investigation.
Inside the surf zone, the wave transformation model of Dally et al. (1985) is
employed to describe the wave breaking and possible reformation process in the surf zone.
This model compares the difference between the local wave energy flux and the stable
wave energy flux that would exist at that local water depth. Similar to the wellknown
constant criteria to describe wave breaking, Dally et al. employ an empirical constant, F,
to provide a stable wave criteria, where is the ratio of the stable wave height to the local
* i0
h/Lo
Figure 4.1 Stream function wave theory shoaling table for normally incident waves.
Figure reproduced from Dean (1974).
Legend
Wove Steepness H/Lo
 For Given Deep Water
" Wove Steepess Ho/Lo
:*** Limit of Linaor Wove
Theory Applicabilty For_
1% in HLo
. I
 4 f
55
water depth. The resultant equation to describe the change in wave height is:
a [H2(h)1/2] K [H2(h)1/2 F2(h)5s2] (4.1)
ax h
where h represents the mean water depth and K is a dimensionless decay coefficient. In the
present model Eq. (4.1) is solved numerically to accommodate arbitrary beach profile
shapes, beginning at the first break point as determined from the offshore shoaling routine.
The model marches the finite difference version of Eq. (4.1) across the surf zone, breaking
the waves, reforming them, and rebreaking them.
Dally et al. provide bestfit values for the stable wave criteria, T, and the
dimensionless decay coefficient, K, of approximately 0.4 and 0.15, respectively, for beach
slopes milder than 1/20. For steeper slopes, the authors recommend that K be set to 0.2.
The authors note that this ambiguity in the choice of K presents some unwanted empiricism
into the model overall, but more definitive data to determine its value is not available.
By employing the wave height transformation model of Dally et al. (1985), more
realistic wave height profiles can be reproduced, including examples with barred beach
profiles, where wave reformation and secondary breaking occurs. A comparison of wave
height profiles for two laboratory cases, one on a planar beach and one on a barred beach
(Okayasu and Katayama, 1992), is shown in Figure 4.2.
One primary advantage of using the Dally et al. model is that it better describes the
sudden decrease in wave height that occurs immediately upon breaking, whereas a simpler
model, based only on a depth limiting criteria, would predict a more gentle decrease in
wave height. The sharp gradient in wave height following breaking leads directly to
increased values of surface shear stress, which will be shown to noticeably affect the shape
measured
predicted, K = 0.15
predicted, K = 0.20
0
0.1
0.2
0.3
0.4
2 0 2 4 6 8 10 12 14
onshore distance from beach toe (m)
VSWL
20
2 1 0 1 2 3 4 5 6 7 8
onshore distance from beach toe (m)
Figure 4.2 Demonstration of wave transformation model for two laboratory studies.
The upper plot (A) compares the measured wave height transformation from the
experiments of Cox & Kobayashi (1997) for a planar beach of 1/35 slope and a period
of T = 2.2 s. The lower plot (B) illustrates the predicted transformation for the random
wave spectrum (the RMS value is plotted) used by Okayasu & Katayama (1992) for a
barred profile of 1/20 slope. The period in the lower plot was 1.14 s. Seaward of the
break point, Stream Function Wave Theory is used for wave height transformation,
whereas landward of the first break point, the model of Dally et al. (1985) is employed.
E
S"o
OL
:J w
.1J
0.)=
S,
E
1>
o
4,
I
57
of the return flow profile. In the barredbeach example of Figure 4.2, the model is able to
characterize the reformation and rebreaking of the wave following the trough.
In the lower plot of Figure 4.2, the random wave spectrum used by Okayasu and
Katayama (1992) was modeled by running individual waves over the profile at a period
equal to the spectral peak (1.14 s). The timedependent feature of the model is applied to
produce the average flow field over a large number of randomly chosen waves. The model
assumes that the wave height distribution is described by a Rayleigh distribution. The
Rayleigh distribution is integrated to produce a cumulative distribution function that varies
from zero to one. A random number generator is then applied to select the series of waves
to model. Typically, 300 waves are run to produce an average flow field. The first few
waves are neglected to avoid averaging "spinup" conditions into the flow field.
The use of the stream function to model the return flow implicitly includes the
effects of pressure. While the pressure terms were eliminated by crossdifferentiation, their
effect is still represented in the physical problem. The effect of pressure is manifest in a
change in the water surface slope, thus resulting in setup or setdown. To accommodate this
in the model, the setup and setdown are calculated at each column and the depth is adjusted
accordingly. The values of setup and setdown are calculated from linear theory. Setdown,
calculated from the offshore boundary to the first break point, is determined as follows:
H2k
1loffshore = 8sin ) (4.2)
8sinh(2kh)
where k represents the wave number (2 z/L). The setup is found from a depthlimited wave
height criteria starting from the knowledge of the setdown at the breakpoint:
oonshore = break + 0.186(hbreakh)
(4.3)
58
4.2 Determination of Return Flow Rate (Selection of pat the Surface)
The first two boundary conditions to be determined govern the flow rate through
each vertical section. This rate, Q, represents the depthintegrated return flow and is equal
in magnitude to the waveinduced volumetric transport, or Stokes Drift, directed onshore
(and assumed to exist above the model domain). In this way, the combination of the flow
through the vertical section of the model domain and the flow assumed to pass above each
vertical section sum to zero, thus conserving volume in the surf zone. Determination of the
volumetric flow rate at each point across the model domain begins with the prediction of
the wave height transformation, described in the previous section.
The flow rate establishes the difference between the value of the stream function
at the bed and at the surface (recall the definition of the stream function from Eq. (3.15)).
For convenience, the value of the stream function at the bed is arbitrarily set to zero, thus
making urface = Q. Outside the surf zone, two choices are presented. The first option is
to apply linear theory to the value of the volumetric flow rate. Linear wave theory predicts
the transport to be:
E
9 (4.4)
pc
where E is the wave energy (= epcgH2) and C is the wave celerity (=V(gh) in shallow
water (Dean and Dalrymple, 1991). An alternative selection would be to choose the
volumetric flow rate that results from a higher order wave theory, such as SFWT. The
development of the solution for Q based on SFWT is quite lengthy and no simple analytic
expression is available, therefore the value of Q chosen is determined from the lookup
tables presented by Dean (1974).
59
To compare, for a wave condition of H= 0.78 m, T= 8 s, and h = 1.0 m, linear
wave theory (Eq. 4.4) predicts a return flow value of 0.238 m3/s per meter of crest width
(this example is the same example previously presented from Dean, 1995). Stream
function theory, on the other hand, predicts a much smaller value of Q in this case, 0.073
m3/s/m. In this particular example the difference in predicted values is quite extreme (the
stream function value is only 30% of the linear value). In general the stream function
values will be lower, owing to the nature of the wave profile predicted, which is much
narrower and more peaked than the sinusoidal form predicted by linear theory.
The decision to choose either the linear theory or SFWT is based on the nature of
the wave height transformation from offshore up to breaking. In general, linear shoaling
produces smaller wave heights at breaking relative to SFWT, and the gradient in wave
height over a differential distance is likewise less. Overall comparison of the two theories
to the laboratory measurements of return flow investigated herein (Chapter 6), both in wave
height profile and return flow magnitude, suggests that SFWT provides better agreement
across the range of experiments. Thus SFWT is coded into the numerical model as the
theory used to shoal waves from the offshore boundary to the first break point.
Inside the surf zone, several options are available. The wave transformation model
provides the wave height at each cross shore position. From that wave height, linear theory
and SFWT can be applied directly for Q. The appropriateness of using these theories by
themselves (i.e., with no consideration given for transition region or roller effects) is
somewhat suspect. For comparison and for use with a roller model, these theories will be
applied to determine their agreement with measurements.
An alternative theory is to adopt a wave height profile resembling a sawtooth, then
develop the expression for Q in the same way the linear value is obtained. In the surf zone,
60
the shallow water assumption is applied, in which the horizontal velocity is assumed to be
uniform over depth and given by:
a = li (4.5)
where 7 in this case represents the instantaneous elevation of the water surface (the wave
height profile). Integration of Eq. (4.5) from the seabed to the water surface and
subsequent time averaging produces the volumetric transport based on the time average of
the square of J, just as linear theory predicts. However, for a sawtooth wave, the sinusoidal
water surface profile is replaced by a linear profile given by:
r = H (L/2 < x < L/2) (4.6)
L
which when squared and averaged between L/2 and L/2 produces a value of Q that is 2/3
the value of Q predicted by linear theory.
Some investigators report that inside the surfzone the value of the volumetric
transport from linear theory is lower than the actual transport due to an additional
contribution from the volume contained in the surface roller fronting the broken wave crest
(Svendsen, 1984a, and Stive & Wind, 1986). Stive and Wind, 1986, report the increased
value of volumetric transport as:
Q = (g/h)1/2H(h +) (4.7)
10
61
While the representation of Stive and Wind is intended to better typify a borelike wave
profile in the surf zone, for the example conditions discussed above, the predicted value of
Q is 0.250 m3/s/m, which is not substantially larger than the linear value (roughly 5%).
At this stage there are now four options for predicting the value of Q inside the surf
zone. If the same example conditions are applied in the surf zone, the range of Q now
varies from the SFWT value (0.073 m3/s/m) to the value predicted by Stive & Wind (1986),
0.250 m3/s/m. The sawtooth wave theory predicts a rate of 0.159 m3/s per meter of crest
width. Lacking further information to refine the estimate, the range of options will be
examined, both to produce the best fit to the measured data, and to compare the validity of
each. It will be shown in a later section that the prediction of the transport value plays a
very important role in the prediction of the return flow and is one of the primary objectives
in simulating the flows in the transition region.
Dally and Brown (1995) present a wave surface roller model that combines SFWT
and a surface roller to predict the momentum flux and volumetric transport throughout the
surf zone. SFWT is applied to describe the properties of the organized wave motion that
persists after breaking. The surface roller component begins to grow after initial breaking
and increases the predicted volumetric transport as well as the momentum flux (or radiation
stress). The use of SFWT inside the surf zone allows that portion of the description of the
parameters to vary smoothly from offshore to onshore. The roller contributions grow
smoothly from the break point and then diminish in the inner surf zone.
Definition of Phase Speed. In all the theories applied herein, the wave theory used
does not include the effect of the return flow on the wave transformation. This is known
as the first definition of phase speed. The SFWT tables presented by Dean (1974) assume
that the wave is propagating without lateral boundaries, so that no external or return flow
62
affects the wave field. Inclusion of the effect of the return flow on the wave
transformation, known as the second definition of phase speed, would entail the
reformulation of the SFWT to include the effect, a procedure described by Dalrymple and
Dean (1975) for uniform currents. This was not included in the present work, where waves
propagate over a varying current. The effect of the return flow on the wave transformation
is expected to be of secondary importance relative to the application of the various
boundary conditions discussed herein. As a check, the wave celerity, C, (first definition
of phase speed) associated with an 8 s linear wave in 1.0 m shallow water is 3.1 m/s. The
expected depth averaged return flow magnitude for such a wave is roughly 0.24 m/s, less
than 8% of the wave celerity. The ratio of U/C alters the apparent period of the wave and
thus adjusts the predicted volumetric flow rate and radiation stress.
4.3 Specification of the Surface Shear Stress Due to Waves
The third boundary condition to be specified is the mean surface shear stress caused
by the passage of the waves. This stress, while actually a force distributed between the
trough to crest region of the waves, is modeled as a concentrated stress acting at MWL, as
described by Dally, 1980. The force arises from the gradient in wave momentum flux
between the trough and the crest. Referring to Figure 1.2, the momentum flux in this
region accounts for onethird of the total momentum flux (via linear theory), and is
generally of greater magnitude than the force due to the uniform pressure gradient. The
surface stress is applied in the numerical model as a secondderivative boundary condition
in z.
a ai)
p( ) (4.8)
az az
63
The value of the mean surface stress, z, is found from the gradient in the wave
induced radiation stress in the troughtocrest region, denoted Sx,. In shallow water, linear
theory produces S, 1 = E/2 (Dean & Dalrymple, 1991), which leads to the following
equation for the applied surface stress outside the surf zone:
d(E12) 1 9H2 E ah
S= pg (4.9)
ax 16 ax 4(h+~i) ax
For a negative bottom slope (decreasing water depth in the positive x, onshore, direction),
Eq. (4.9) produces an offshore directed stress at the surface. Inside the surf zone, where
waves are breaking and decaying in height, r, is found by assuming that the wave height
is related to the local water depth by a constant (i.e. H = Ah). With S, equal to E/2, is
assumed to be:
1 g2h 2 E ah (4
T, pgK (4.10)
16 ax (h+lj) ax
which produces an onshore shear stress at the surface when the bottom slope is negative.
As with the flow rate, researchers have suggested that the additional volume of
water carried forward by the surface roller contributes an additional stress to the water
surface. Svendsen (1984a) suggests an alternate representation for r,:
1 Ar h aH H2
S( + _)pg (4.11)
16 H2L ax
The second term in parentheses in Eq. (4.11) represents the effect of the surface roller,
whose cross sectional area is Ar (c 0.9H2) (Svendsen, 1984a). The application of Eq. (4.11)
64
compared to Eq. (4.12) generally applies a greater stress to the water column. The surface
roller model of Dally and Brown (1995) also produces the radiation stress at each point
across the domain. From that description, the gradient can be computed to determine the
surface shear stress. Lacking better information on the division of the radiation stress
between the troughtocrest region and the body of the fluid in nonlinear waves, it is
assumed that onethird of the computed radiation stress from this and every other model
applied represents the surface shear stress.
It will be shown that the presence of a surface shear stress in most surf zone
conditions drives the peak of the offshore return flow to a lower point in the water column
and reduces the offshore directed return flow near the surface. In some instances it may
cause the opportunity for flow reversal in the upper portion of the water column (onshore).
The measurement of surface shear stress due to breaking waves is extremely difficult;
information on this subject is inferred from velocity measurements but is affected by the
aerated nature of the surf zone, particularly in the upper portions of the water column.
4.4 Specification of Bottom & Sidewall Tangential Boundary Conditions
The final boundary condition to be specified is the tangential condition along the
bottom and sidewalls. Here the options are to assign a velocity or a shear stress condition.
Generally, the model domain is established such that the offshore boundary is well offshore
of the break point. At that position, the mean velocities in both directions are assumed to
be very small, thus the tangential condition along the sidewalls is chosen to be a zero
velocity (noslip) condition. It is noted that there is little data to support applying any other
condition. At the onshore boundary, the domain has been extended nearly into the swash
zone; assignment of a vertical noslip condition at this point is considered reasonable.
65
Along the bottom, the modeled PDE does not expressly consider the presence of a
bottom boundary layer, induced by the oscillatory motion of the waves and/or by the return
flow. Given that fact, it may be desirable to specify the bottom velocity as either zero or
that velocity corresponding to the velocity at the edge of the boundary layer. For generality
in the model, the 1D model bottom velocity boundary condition is specified as follows:
Ub = ()b (4.12)
8z
As mentioned, the basis for selection of a bottom velocity boundary condition other
than zero is fairly limited. It would seem appealing to apply the velocity at the edge of the
boundary layer, also known as the streaming velocity, at this stage (LonguetHiggins,
1953). The streaming velocity is given by:
3 H20k
U= H (4.13)
16 sinh2kh
where a= 2)rT, with Tbeing the wave period. However, it will be shown in following
chapters that very little data exist to support an onshore directed velocity at the bottom.
This is particularly true in the surf zone, where strong offshore directed bottom velocities
are frequently measured. Hansen and Svendsen, 1984, noted that although the direction
was incorrect, the magnitude of the bottom velocity predicted by Eq. 4.13 provided good
agreement with experimental data. Hansen and Svendsen provided no supporting
hypothesis for reversing the direction of this velocity.
It is equally easy to apply a zero velocity condition in Eq. (4.12), and much more
defensible since it is required that the velocity ultimately go to zero at the bottom through
66
some boundary layer. The difficulty in using a zero velocity at the bottom is that the
mathematics of the governing PDE are not intended to describe a boundary layer, and in
the case of uniform eddy viscosity the resultant velocity profile is parabolic and pinned at
zero at the bed. This condition does not always correspond well with measurements, as
will be shown in a later chapter. It is possible that applying an adjusted eddy viscosity
value at the bottom, in which the eddy viscosity is significantly reduced, may alleviate the
strict condition at the bottom, making the resultant velocity profile appear to "slip" over the
bottom. A variable eddy viscosity, however, introduces an entirely different set of
challenges; these will be discussed in a subsequent chapter.
The other alternative considered is the specification of the bottom shear stress at
each vertical section. Eq. (4.14) provides this condition, in which the mean bottom shear
stress has been defined as a fraction, m, of the mean surface shear stress following the work
ofReid (1957). The generality of a varying eddy viscosity in the numerical model has been
respected:
T b P(Eb ) (4.14)
"c = m', = Pz z) (4.14)
az dz
For the case of surface stress along a bounded channel, such as in the present beach
profile situation, the quantity m is generally negative, indicating the bottom stress is
directed opposite to the surface stress, as might be expected inside the surf zone for the
mean flow in the absence of wind. Reid (1957) discusses the value ofm for steady, wind
driven flows in bounded channels. Applying the example conditions from Dean (1995) (Q
= 0.238 m3/s/m, r, = 7.9 N/m2, and h = 1 m), the ratio of bottom to surface shear stress is
predicted to be approximately 0.09 to 0.16, depending on the choice of roughness height
67
(reference Reid, 1957, Figure 11 and Figure 12). It is interesting to note that the shear
stress value given by Dean for the surface shear stress due to momentum transfer from
breaking waves ( 7.9 N/m2) corresponds roughly to the stress generated by winds just
reaching hurricane force (>121 km/h, per Reid (1957)). This illustrates the strong effect
breaking waves have on the surf zone flow.
Reid (1957) developed expressions for the velocity profiles in channels due to
applied wind stresses by applying the quadratic friction law:
Tb = P ave lUve (4.15)
wherefis the DarcyWeisbach friction factor and Uave is the depthaveraged value of the
return flow. Applying Eq. (3.46) to Dean's example yields a friction factor of roughly 0.1,
which approaches the upper limit for friction factors applied in openchannel and pipe
flows.
Diegaard et al. (1991) measured the mean and timevarying shear stress under
regular waves breaking on a 1:30 beach slope. Few details are given regarding the wave
height transformation; the wave height information provided is listed as regular waves of
10 to 13 cm height (presumably on the horizontal portion of the wave tank or at the toe) and
1.5 to 3.0 s period, similar in scale to the experiments discussed above. The authors present
data indicating mean bottom shear stresses of 1 to 2 N/m2 for spilling breakers, and up to
3 N/m2 for plunging breakers. The breaker heights are not presented, however, using crude
assumptions and basic linear theory, Eq. (4.10) would predict a surface shear stress of
roughly 3.8 N/m2, producing a ratio of stresses of between 0.26 and 0.53 for the spilling
breakers. SFWT would produce a higher wave height and a higher shear stress, but a lower
68
volumetric transport and hence net return flow rate, leading to lower values of the ratio of
bottom to surface shear stress.
Overall, the surface shear stress due to breaking waves is expected to be
substantially higher than the resultant mean bottom shear stress, perhaps by as much as ten
times (per Reid, 1957). It is noted, however, that the return flow may exist in the absence
of a significant surface shear stress, such as in the region just seaward of breaking. In those
areas, the ratio of the mean bottom to surface shear stress may be greater than one, or may
have no physical meaning (as in the case of zero surface shear stress). The analogy is
drawn to open channel flow on a gently sloping surface, where the bottom shear stress is
nonzero, but the surface shear stress is zero. In application to the model, the existing
literature suggest that in the surf zone, where the return flow is driven by breaking waves,
the mean bottom shear stress is not expected to exceed 0.3 to 0.5 times the surface shear
stress. Those values are believed to be the extreme, high estimate, and an estimate of 0.1
or less may in fact be more appropriate. Again, however, guidance as to the selection of
the mean bottom shear stress or the ratio of stresses is limited.
4.5 Matching Conditions Across the Transition Region
The present model is driven from the surface boundary conditions where shoaling
and breaking waves exert stresses on the water surface and transport mass onshore which
must be returned offshore through the model domain. Techniques for applying these
boundary conditions are numerous, as demonstrated in the previous sections. Most of these
techniques, however, do not consider the continuity in certain wave properties, such as
radiation stress, across the entire profile. For example, the model of Svendsen, 1984a,
begins at the end of the transition region, some distance landward of where breaking is
69
initiated. The expressions developed for radiation stress and mass transport result in large
discontinuities at the break point when linked to the corresponding offshore descriptions
of these properties.
Discontinuities in the mass transport are not as great a concern as are discontinuities
in momentum. The surface roller represents an additional volume of water that can be
simply recirculated within the surf zone. This volume of water could originate from
outside the surf zone at the beginning of a wave event or from alongshore. Regardless, the
fact that there is a recirculating volume of water that does not leave the surf zone is not
problematic so long as conservation of mass (volume) is achieved at each vertical section.
The increase in momentum represented by the surface roller does present continuity
problems, however. Just prior to breaking, there is a certain amount of momentum in the
wave available to be transferred into the surf zone. This momentum must be distributed
across the surf zone. Generally this momentum is transferred to the water column via the
gradient in radiation stress seen in the decreasing wave height due to breaking. This is
principally manifest in the surface shear stress and set up of the MWL in the surf zone.
Depending on the representation used to describe this surface shear stress (which is
determined in practice by the reduction in wave height across the surf zone), integration of
these shear stresses from the shoreline seaward to the break point may result in a total
required momentum input to the surf zone that is greater than the momentum input
determined from the seaward side (i.e. from SFWT).
Figure 4.3 illustrates this problem and one simple solution adopted to span the
transition region and maintain continuity of momentum. The figure plots the value of
several wave quantities computed across the surf zone. The data set of Ting and Kirby is
chosen to demonstrate the problem. In the upper plot, the variation in wave height across
the surf zone is shown. SFWT is linked at the break point to the wave breaking model of
Dally et al. (1985). This variation is, by default, fairly smooth and does not exhibit any
substantial discontinuities (except for the cusp at the break point itself).
The second plot in Figure 4.3 illustrates the variation in radiation stress across the
surf zone, dependent on the choice of description used. At the break point, the decision
must be made as to which theory to apply to describe the necessary properties. The plot
shows the discontinuity that arises from any choice of theory sawtoothh, linear, or surface
roller). The surface roller approximation produces the greatest discrepancy. A similar
discontinuity at the break point is seen in the third plot of the volumetric transport.
There is no mechanism to generate the additional radiation stress, thus the switch
in theories presents a problem with conservation of momentum. This problem, the
transition region problem, arises in any return flow model, generally because the majority
of data collected in the surf zone exhibit integral wave properties that are greater than
predicted by higher order wave theories, such as SFWT. Svendsen (1984a) noted this
discrepancy in the volumetric transport measured at vertical sections in laboratory surf
zones and was prompted to develop the surface roller theory. Svendsen, however, did not
discuss the issue of conservation of momentum from the break point landward.
One ad hoc approach to this dilemma is to determine the distance required for the
radiation stress to return to the level of the radiation stress just offshore of the break point.
This distance may be considered to be the length of the transition region. In application to
the periodic wave experiments discussed in this work, this distance is on the order of 6 to
11 water depths at breaking. This appears to be in reasonable agreement with the
measurements and discussion of transition lengths (e.g. Nadaoka and Kondoh, 1982).
Svendsen (1984a) describes how breaking waves experience a transition region in which
t. E
o U
2 E
'Y
0
.>_
UE
Surf Zone Conditions
SLinear Theory
 Roller (Svendsen)
Sawtooth Profile
Patch w/ Roller
7
5 6 7 8 9
Distance from beach toe (m)
Figure 4.3 Demonstration of the treatment of the transition region in the simulation of
the data set of Ting & Kirby (1994). The expressions for radiation stress inside and
outside the surf zone are investigated to determine the length necessary for the surf zone
value to match the incident value seaward of breaking. This distance is used as a
folding length to patch the offshore descriptions of radiation stress and volumetric
transport to the inshore descriptions.
72
the wave is reshaped from a nearly breaking, steep, asymmetric shape just offshore of the
break point to a fully developed turbulent bore. Laboratory data cited by Svendsen indicate
that in this region there is little momentum loss, as evidenced by the lack of setup of the
water surface. Given the lack of momentum loss, the technique of matching the radiation
stresses and then patching the offshore description to the onshore description gains some
degree of merit.
The offshore and onshore descriptions are patched together in the following fashion.
The distance across the transition region is determined by inspecting the radiation stress
determined from the surf zone theory for the point where the radiation stress matches the
radiation stress just seaward of the break point (b.p.). This distance, A/is then used to patch
the two descriptions of radiation stress such that the stress varies consistently from the
break point shoreward. The patching process is demonstrated in Eq. (4.16):
x
(Sxx = (Sxxner ( (Sxx)b.p.+,nner (Sxx)b.p.stream function )e K (4.16)
where the values denoted by inner could be determined from any theory, such as linear
wave theory or a surface roller approximation.
Similarly, the offshore and onshore values of the volumetric transport, Q, are
patched to produce a consistent variation. The variation in Q is shown in the third plot of
Figure 4.3. The transition distance from the radiation stress plot produces a slight increase
in the transport just landward of breaking, but still reduces the predicted value of Q
substantially compared to the unpatched case. Inspection of the measured values of Q
indicates that the maximum value of transport occurs some distance landward of the
published break point (the velocity profiles were extrapolated to the bottom and to the
73
MWL). Published literature suggest that in this region the wave is reforming itself into the
surface roller form described by Svendsen, while only transferring a small amount of
momentum to the water column. The patch described herein reproduces this phenomenon.
Using the patched description of the radiation stress, the surface shear stress can be
computed and thus varies smoothly across the profile as well.
Further evidence of the transition region is seen by inspection of the velocity profile
immediately landward of breaking in Figure 4.3. This profile suggests that the strong
return flow profile characteristic of inner surf zone areas is not present. In fact the profile
at the break point exhibits a shape typical of offshore velocity profiles. It is apparent from
the plot of volumetric transport that the use of the higher value of transport (e.g. from linear
or surface roller expressions) would lead to a significant overestimate in the transport
value. The merged curve reproduces the measured transport values quite well, with the
exception of the fifth profile (x 8.5 m). No explanation for the transport measured at that
station is provided. A coefficient could be introduced to modify the characteristic distance,
/, applied in Eq. (4.16), however, no guidance for such a coefficient exists, except for a
fitting to existing data, which would introduce an unwanted degree of empiricism.
Dally and Brown (1995) present a surface roller model in which the dissipation of
wave energy is assumed to be derived from the shear at the interface between the organized
wave form (predicted by SFWT) and the turbulent surface roller fronting the wave. The
dissipation is represented by contributions from both the organized wave motion and the
roller, which first grows in size then dissipates across the surf zone. The initial growth of
the roller limits the decrease in total radiation stress, and provides a plausible means of
increasing the volumetric transport landward of the break point. Dally and Brown provide
expressions for the momentum flux and the volumetric flux of the roller. The volume flux
74
of the roller is then added to the volume flux of the organized wave form and thus produces
somewhat higher predicted values of transport across the transition region.
In this work, all simulations incorporate some attempt to patch the offshore region
to the inner surf zone by providing a treatment of the transition region. This is
accomplished either through the characteristic patchlength approach or through application
of the Dally and Brown (1995) roller model. In the first instance, the characteristic length
will be determined from linear and sawtooth theories as well as the roller model of
Svendsen (1984a).
CHAPTER 5
MODEL TESTING & VALIDATION
Before proceeding to the modeling of actual surf zone flows, the model will be
applied in simple flow situations for which analytical solutions exist in order to validate its
performance. While the inclusion of an arbitrary domain geometry and a spatially variable
eddy viscosity field precludes an analytic solution to the full problem, certain solutions are
available for simpler cases, such as simple onedimensional flows and twodimensional
flows of uniform eddy viscosity and horizontal bottom geometry. Using these simple
examples, the models may be tested and validated to develop confidence in the results
produced from the more complicated simulations.
This chapter presents test cases designed to demonstrate the proper performance of
the numerical model developed in Chapter 3. The steadystate 1D simulations are
compared to solutions of horizontal flows subject to various boundary conditions. The
timedependent 1D flow field is compared to solutions of the diffusion equation, a second
order, time dependent PDE in terms of the horizontal velocity, U. The present, governing
fourthorder equation in fris essentially equivalent to the equation describing the diffusion
of vorticity, and solutions in should produce equivalent solutions in U when converted.
This is quite advantageous, since solutions to the time dependent fourthorder equation,
even in one dimension, are difficult to obtain for the desired boundary conditions.
For the 2D case, obtaining an analytic solution for any timedependent flow in any
geometry has proven difficult. In light of this problem, the steadystate problem is
76
investigated. Solutions for the steady state problem in two dimensions follow those
solutions of the biharmonic equation. A Fourier series solution is developed to compare
to the results of the numerical model. Additionally, for extremely simple flow situations,
vertical sections taken well away from sidewall boundaries and areas of spatial gradients
in the horizontal are expected to behave very similarly to the onedimensional solutions.
5.1 OneDimensional Analytic Solutions
5.1.1 SteadyState Simulations
To begin, an analytic solution to the steady, onedimensional, fourthorder equation
in $'is found. The PDE is:
S= 0 (5.1)
Solutions of the form l= e" are sought. Application of this form to Eq. (5.1) yields four
repeated roots (=0), resulting in the following general solution:
I* = Az3 + Bz2 + Cz + D (5.2)
The boundary conditions to be investigated are illustrated in Figure 5.1 and are listed as:
(0) =(h) =0 z(0)=0 z(h)= (5.3)
PE
Application of the four boundary conditions in Eq. (5.3) results in the following
equation for the steady state, analytic solution:
/ 3 2/ 2 T1h 2 )
_(1D ,, =  z + _ 2 ()2(1 ) (5.4)
4h 4 4 h h
1
Figure 5.1 Definition sketch of onedimensional problem. The four boundary
conditions are shown along with typical profiles of fand the resulting velocity
distribution. Note that the schematic velocity profile results in zero net flow through the
vertical section, consistent with the applied boundary conditions of 0(0) = /(h) = 0.
where z'= r/peand the origin is taken at the seabed. This solution, when converted to the
corresponding horizontal velocity, U, produces the profile of Dally (1980) and Dean (1995)
described in Chapter 2 (Equation (2.5)). Figure 5.2 plots the results of the simulation under
the conditions in Eq. (5.3) for the example values of Dean (1995). The plot shows the near
perfect agreement of the simulation with the analytic solution in this simple example.
In Figure 5.2, the resulting velocity and shear stress profiles demonstrate the
parabolic and linear shapes, respectively, that one would expect from taking successive
derivatives of Eq. (5.4). This example was chosen to produce no net flow through the
vertical section. In the surfzone, this situation could be replaced by some nonzero flow
value, such as the waveinduced volumetric transport, simply by adjusting the surface
boundary condition in Eq. (5.3) to the desired value, e.g. Of(h) = Q.
 Analytic solution
 Numerical Solution
0 0.002 0.004 0.006 0.008
stream function (m2/s)
shear stress, r (N/m2)
4 0 4 8
0.02 0 0.02 0.04
velocity, U (m/s)
Figure 5.2 Comparison of 1D numerical model to analytic solution for steadystate
conditions. The chosen example simulates flow in a bounded channel driven only by a
surface shear stress. The example conditions are similar to those given in Dean (1995)
(r= 7.9 N/m2, e= 0.04 m2/s, Q = 0.0).
5.1.2 Variable Eddy Viscosity Simulations Steady State
Reid (1957) presented velocity and shear stress profiles in a bounded channel (i.e.,
Q = 0) subject to a surface stress using a quadratic expression for mixing length and hence
turbulent diffusion. As discussed in Chapter 2, Reid combined a parabolic description of
the turbulent mixing length, L, to the quadratic friction law to determine the resultant linear
shear stress profile. The description of the mixing length employs coefficients describing
the relative roughness of the bed and the water surface (Eq. (2.8)). As discussed in
previous chapters, limited experimental results are provided to guide the selection of these
characteristic lengths and the ratio of the bottom and surface shear stresses.
0.06
To demonstrate the ability of the numerical model to properly accommodate
variable eddy viscosity fields, the case of zero net flow was simulated. The results of the
simulation are illustrated in Figure 5.3. The example of Reid is nondimensionalized by
the surface shear velocity (I(r/p)), therefore, the surface shear stress applied in the
previous example by Dean (1995) is employed to present the results in dimensional units.
For the case of m = r/z, < 0, the model employs two different expressions, joined at the
elevation where the return flow value is a maximum (the point where the shear stress is
zero, oTJ/& = 0). Figure 5.3 also plots the profile predicted by Reid (1957). The profile
produced is much more uniform in the middle of the domain, suggestive of the turbulent
flow present. The resulting parabolic mixing length produced in this example is shown in
the first lower plot (b).
Using the relationship between mixing length and eddy viscosity presented by
Eq.(2.9), the eddy viscosity distribution suggested by the example in Figure 5.3 is plotted
in Figure 5.3(c). The difficulty of the use of two equations to fit the profile, along with the
indeterminate value of e that results where oJ/cG = 0, is apparent in the plot near the
bottom. An exact fit of predicted from the method of Reid creates sharp gradients in the
derivatives of the eddy viscosity distribution that do not allow solution in the numerical
model. To reproduce the Reid velocity profile, the distribution of ewas fit to a smooth
thirdorder polynomial as shown in Figure 5.3(c). The velocity profile produced by the
numerical model is in very good agreement with the analytical approach of Reid (1957).
1
0.8
S0.6
S0.4
0)
0.2
0
0.4 0 0.4 0.8 1.2
velocity (m/s)
I i
0.04 0.08
mixing length (L/h)
1
0.8
, 0.6
0
m 0.4
0.2
0
0 0.002 0.004 0.006 0.008
kinematic eddy viscosity (m2/s)
Figure 5.3 Comparison of 1D numerical model results to those of Reid (1957) for the
case of Q = 0. The upper plot depicts the analytic and numerical solutions to the
horizontal velocity profile (a). Reid employs a parabolic mixing length theory (b). In
order to reproduce the predicted eddy viscosity seen in plot (c), it was necessary to fit
the predicted curve as a 3rd order polynomial. The difficulty in fitting the predicted
eddy viscosity from Reid is apparent in plot (c), where the two expressions used by Reid
are fit at the point in the profile where cU/c& = 0, thus producing an indeterminate
condition.
1.6 2
1
0.8
0.6
0.4
0.2
_(b)
. .
0.12
81
Note that the eddy viscosity distribution predicted from Reid (1957) peaks at a
depth approximately equal to 0.6 times the local water depth. In surf zone flows, it is
expected that this peak may be shifted further toward the surface. In laboratory
experiments with spilling breakers, Ting and Kirby (1994) present the measured vertical
distribution of the timemean horizontal turbulent velocity. These distributions exhibit a
similar shape to the eddy viscosity distribution shown in Figure 5.3, however the peak of
the mean turbulent velocity measured in the laboratory under breaking waves appears to
be as high as 0.7 to 0.8 times the local water depth. Stive (1984) presents similar
distributions of turbulent intensity under breaking waves in which the peak value is shifted
toward the surface.
5.1.3 TimeDependent Simulation
The ability of the model to reproduce timedependent behavior is now investigated.
In this instance the onedimensional equation in becomese:
= E zzz (5.5)
where for simplicity the eddy viscosity has been assumed uniform. The boundary
conditions of Eq. (5.3) are retained along with an initial condition of quiescence (f/(z, 0) =
0). In the time dependent case, analytical solutions to the fourth order PDE of Eq. (5.5)
were unobtainable, however, it is recognized that in terms of velocity, U, the model should
reproduce results from simulations of the secondorder diffusion equation:
S 2 (5.6)
at 8z'
82
Solutions to Eq. (5.6), typically Fourier sine series solutions given the finite
physical domain, are readily obtainable. Carslaw and Jaeger (1986), provide two solutions
that can be superimposed to produce the timedependent version of the steadystate
example presented in Figure 5.2. To begin, the solution for the horizontal velocity profile
that develops from a fixed bottom velocity (=0) and a specified surface shear stress, r, is
found:
00 (2n+l)n 2t
Sh8h (1)" (2n+l)xr 4) ]t
U, =  z + L ( sin( z) e 4h (5.7)
PE T2 pE n=0 (2n+1)2 2h
with boundary and initial conditions:
U(z,O) = 0 U(0,t) = 0 U1z(h,t) (5.8)
PE
The second solution needed is similar to the solution that results from modeling the
diffusion of heat within a slender rod with internal heat generation. The PDE in this case
is:
U A(z,t) (5.9)
at 8z2
In the present physical problem, the additional term, A, represents the combination of the
pressure and momentum flux gradients, which would result in the slope of the water
surface. One solution to Eq. (5.9), for the case of a uniform forcing term, Ao, is:
Aoh 2 2 32 (1)" (2n(2l) E n ,
Sz 2 [ cos( n1)z) e (5.10)
2 h 2 n n=0 (2n+1)3 2h
subject to the conditions:
U2(z,0) = 0 U2(0,t) = 0 U2(2h,t) = 0 (5.11)
Note that in the formulation for U2 by Carslaw and Jaeger, symmetry has been employed
to solve the equation in the domain 0 z 2h, thus producing a situation in which the
shear stress at the surface (z = h) is zero. The superposition of the two solutions, U, and U2,
between 0
The first solution, Eq. (5.7), applies to the case of zero internal forcing (A = 0), such
as the case of wind blowing across an unbounded channel. Application of Eq. (5.10) then
constrains the solution to achieve the desired flow rate. In Eq. (5.7), the flow rate at infinite
time is dictated by the shear stress, which controls the slope of the velocity profile. At
infinite time the velocity profile would be linear with a slope of r/pe and the
corresponding Q would be rh2/(2pe) per unit channel width (see Figure 5.4).
The 1D simulation is compared with solutions from Eqs. (5.7) and (5.10) in Figure
5.4, where the constant term Ao in Eq. (5.10) is adjusted at each time, t, to produce the
desired flow rate through the section. The pressure gradient that produces the desired Q
must change with time. In keeping with the previous example, the case to be investigated
will be the "spinup to" and spindown from cases of the example in Figure 5.2.
Figure 5.4 plots the results of the simulation to bring the initially quiescent domain
to the steadystate condition seen in Figure 5.2. The figure plots the two separate solutions
from Eqs. (5.7) and (5.10) at t = 1 s. For reference the steady state result of each solution
and their sum is shown. The numerical result at the corresponding time agrees extremely
well with the analytic solution in this case. The numerical model is run with a time step
of 0.1 s and a spatial resolution of 1/24th the water depth, h (h = 1 m in this example).
