Numerical simulation of mean cross-shore currents

UFL/COEL-TR1127

NUMERICAL SIMULATION OF MEAN CROSS-SHORE
CURRENTS: A STREAM FUNCTION APPROACH

by

Albert E. Browder

Dissertation

2000

NUMERICAL SIMULATION OF
MEAN CROSS-SHORE 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 One-Dimensional Model ........................ 37
3.2.1 Boundary Conditions for One-Dimensional Model .......... 40
3.2.2 Solution Scheme for One-Dimensional Model ............. 41
3.3 Development of Two-Dimensional Return Flow Model ............. 43
3.3.1 Transformation of the Partial Differential Equation ......... 43
3.3.2 Discretization of the Two-Dimensional PDE .............. 45
3.3.3 Boundary Conditions for Two-Dimensional Model ......... 48
3.3.4 Solution Scheme for Two-Dimensional 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 One-Dimensional Analytic Solutions .......................... 76
5.1.1 Steady-State Simulations .............................. 76
5.1.2 Variable Eddy Viscosity Simulations Steady State ......... 78
5.1.3 Time-Dependent Simulation ........................... 81
5.2 Two-Dimensional Model Validation ............................ 87
5.2.1 Comparison to One-Dimensional 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 CROSS-SHORE 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 One-Dimensional Model ....................... 190
A.1.1 Boundary Conditions for One-Dimensional Model ........ 192
A.1.2 Solution Scheme for One-Dimensional Model ............ 195
A.2 Development of Two-Dimensional Return Flow Model ............ 196
A.2.1 Transformation of the Partial Differential Equation ........ 196
A.2.2 Discretization of the Two-Dimensional PDE ............. 201
A.2.3 Boundary Conditions for Two-Dimensional Model ........ 207
A.2.4 Solution Scheme for Two-Dimensional 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 wave-induced 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 one-dimensional model.......................... 39

3.3 Example output from one-dimensional return flow model ............... 42

3.4 Definition sketch of coordinate transformation between the physical, x-z,
domain and the computational, -)7, domain .................... 44

3.5 Solution cells for the fourth-order finite difference equation in the x-z domain
and the transformed f- domain. ............................. 46

3.6 Flowchart for two-dimensional 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 one-dimensional problem ......................... 77

5.2 Comparison of 1-D numerical model to analytic solution for steady-state
conditions. .......................................... 78

5.3 Comparison of 1-D 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 one-dimensional semi -
implicit numerical model. ................................. 86

5.7 Comparison of 2-D numerical model to 1-D 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 root-mean-square 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 no-slip
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 root-mean-square 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 root-mean-square 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 root-mean-square 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 time-averaged 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 root-mean-square 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 time-averaged random waves versus
monochromatic waves of RMS height ....................... 153

7.1 Cross-shore variation in mean bottom shear stress for the experiment
of Cox and Kobayashi (1997) ............................... 159

7.2 Cross-shore 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 non-dimensionalized 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 x-z domain and the ,7 domain (# 1 through 14, units vary).

A, Cross-sectional area of wave surface roller (m2).

B# Time dependent term coefficients describing the transformation of coordinates
between the x-z 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 Root-mean-square 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 wave-induced 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 Non-dimensional root-mean-square 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 a-direction on the b-plane.

rZb Turbulent shear stress in the a-direction on the b-plane

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 CROSS-SHORE 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 two-dimensional, wave-induced return flow in the nearshore region,

commonly referred to as the undertow. The model is developed using the stream function,

r, to represent the wave-induced 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 two-dimensional vertical plane perpendicular to the shoreline.

The governing equation, a fourth-order partial differential equation in f, is not

specific to wave-induced 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 bar-trough beach profiles or any other

irregular profile geometry.

Comparison to measurements of wave-induced 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 wave-induced volumetric transport. Further improvement of the predictive capability

of the model inside the surf zone stems from the application of a vertically-varying 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 cross-shore, two-dimensional vertical plane (2-DV), 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 wave-induced cross-shore 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 cross-shore flow, mobilization and/or transport of sediment by the mean wave-

induced current obviously plays a role in shaping the cross-shore profile and may play a

strong role in the formation of longshore bars in the profile.

The mean cross-shore 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 wave-induced return flow in the surf zone. The return flow is
the two-dimensional 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 one-dimensional

and two-dimensional 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.

Dyhr-Nielsen 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 fourth-order PDE, however, prompted the interest

in creating a two-dimensional 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. Longuet-Higgins and Stewart

(1962) derive the balance of wave-induced radiation stress and the subsequent uniform

pressure gradient caused by the water surface slope (set-up or set-down). Many researchers

have discussed that these two forces, while balanced in a depth-averaged 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 2-D 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 2-DV 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

cross-differentiated, 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 cross-differentiation 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. Non-coupled

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 2-D, fourth-order 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, non-wave 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 wave-driven 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 wave-induced 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 2-D versus 1-D formulation).

Inclusion of an arbitrary beach profile geometry and a spatially variable turbulent

eddy viscosity field precludes an analytic solution of the 2-D problem. Therefore a two-

dimensional numerical model is developed to investigate the aforementioned issues. The

2-D numerical model is compared to simple flow situations to check the results against

existing 2-D analytic solutions and other return flow models.

The inclusion of a time-dependent 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 shoreward-directed 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 wave-induced

radiation stress and the resultant set-up or set-down of the water surface, as described by

Longuet-Higgins and Stewart (1962) (Figure 1.2).

Longuet-Higgins and Stewart (1962, 1964) present the mean horizontal momentum

balance, depth-integrated and time-averaged, 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 Longuet-Higgins 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

Dyhr-Nielsen 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

time-averaged and cross-differentiated 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

fourth-order 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 second-order model used by Dally is derived from the

time averaged equations of motion and includes the wave trough-to-crest 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 trough-to-crest 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 trough-to-crest 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 trough-to-crest 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 trough-to-crest 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).

Wave-related 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 fourth-order 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 Longuet-Higgins (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, Longuet-Higgins 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 second-derivative 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 Longuet-Higgins (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., Longuet-Higgins 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 Longuet-Higgins, Dyhr-Nielsen 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 well-known 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 2-D Reynolds equation of turbulent flow in the x-direction, 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 + -[ 2v-a(U++u/)] (2.3)
p ax ax ax
a aa
+ -[ v( --(U++u) + -(W+i+w + )]
az az ax

Generally, in the literature the long-period 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 depth-integration of Eq. (2.4) produces the original

balance shown in Eq. (2.1), where the term on the left-hand 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 trough-to-crest 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 Longuet-Higgins (1953).

The first term on the right-hand side of Eq. (2.5) arises from the wave momentum

flux in the trough-to-crest 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-

to-crest 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 trough-to-crest 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 trough-level 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 wave-induced 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 trough-to-crest

region, it is more appropriate to specify the trough-level 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 non-zero, 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

time-dependent 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 1-D 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 depth-uniform eddy viscosity at each cross-shore 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 depth-uniform 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 = 10-2ch, 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 vertically-uniform

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 non-uniform 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 semi-infinite 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 mid-depth. 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 semi-infinite 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+zl-z)
.E = f(Zn) (C)

0 0 0 V

= f(z) L = +zo)(h+zl-z) E = f(zn) (d)

Figure 2.1 Illustration of the effect of the quadratic mixing length expression employed by
Reid (1957). The analogy to gravity-driven 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 oppositely-directed 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 mid-depth. 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 mid-depth. 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 (depth-uniform) 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 barred-beach 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 one-dimensional (l-D) and two-dimensional (2-D) versions of the

model are presented (both options are contained in the same program). For the 2-D 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 2-D 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 cross-differentiating the two-dimensional Navier-Stokes

equations for turbulent flow (e.g. the Reynolds equations). The resulting third-order partial

differential equation (PDE) in terms of velocity is then recast in terms of the stream

function, ultimately producing a fourth-order PDE that retains the generalization of time

30

dependence and a non-uniform eddy viscosity field. This turbulent eddy viscosity, e, is

considered a variable quantity in both the cross-shore (x) and vertical (z) directions.

This development, simplified by the assumption of uniform eddy viscosity and

steady state, produces the well-known 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 two-dimensional 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 cross-differentiating 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 cross-differentiating 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 Navier-Stokes

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 x-gradient term have been retained, along with the last four

components of the z-gradient term (Eq. (3.12)). From Eq. (3.8) the last four components

of the x-gradient term and the last three components of the z-gradient 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 (Longuet-Higgins 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 non-linear 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 time-average 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 2-D 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 cross-differentiation 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 x-z 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 One-Dimensional 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 1-D 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, semi-implicit

fashion by the inclusion of the time-weighting term, 0. The model domain is illustrated

in Figure 3.2. The fourth-order 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 1-D 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)
+ (1-0) [( 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
j-1 depth, h+Ti

applied bottom velocity ub) j-2
or applied shear stress ('b) grid point 1 (9 = 0)
seabed

Figure 3.2 Definition sketch for one-dimensional model. The figure illustrates the four
applied boundary conditions and five-point solution cell.

which is quite restrictive, resulting in long run times in a 2-D 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 One-Dimensional 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 Longuet-Higgins (1953) for establishing

shear stress, velocity, and stream function value conditions at the boundaries.

In the 1-D 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 second-derivative 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 steady-state solution in fewer time steps. A steady state solution in the 1-D model

is quickly determined and provision of a "hot-start" condition is not needed. In the time-

dependent 2-D cases, a good initial estimate is beneficial.

3.2.2 Solution Scheme for One-Dimensional Model

The stream function values across the model domain at the n+1 time step are solved

simultaneously after creating a 5-band diagonal matrix ofN-2 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 N-2. 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

lower-upper (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 1-D 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 time-dependent 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 one-dimensional 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 Two-Dimensional Return Flow Model

Building on the 1-D model, Eq. (3.19) is revisited and applied in the 2-DV

nearshore region. The same boundary conditions used in the 1-D model apply to the 2-D

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, x-z coordinate

system, a coordinate stretching scheme is implemented to allow for the modeling of

arbitrary beach profiles. This coordinate stretching simply non-dimensionalizes 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 non-dimensionalized by the local water depth (Eq. (3.23)). This

simple grid-stretching technique allows for the modeling of any profile (Figure 3.4). As

contrasted to the 1-D model, the 2-D 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 x-z 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

1-2

-6

50 100 150 200
x-axis, onshore distance (m, typical)

I

250

.0
4'

1 onshore grid-point # N
5-axis, onshore grid-point #

Figure 3.4 Definition sketch of coordinate transformation between the physical, x-z,
domain and the computational, e-r7, domain.

The approach throughout this derivation will be to separate the time-dependent

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

fourth-order 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 Two-Dimensional 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 1-D model). Figure 3.5 compares the discretized solution

cell for the x-z domain versus the solution cell for the transformed e--? domain. The

transformation process adds an additional eight points to the finite-difference equation.

Eq. (3.25) is discretized using central difference spatial approximations of fr with

a forward difference approximation in time. As in the 1-D model, the equation includes

z

j-2
qF

j-1

j+

j+1

j+2- I- -

i-2 i-1 i i+1 i+2

x
V

j-2

j-1

j

j+1

j+2

i-2 i-1 i i+1 i+2

Figure 3.5 Solution cells for the fourth-order finite difference equation in the x-z
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 + (1-O)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 1-D 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 steady-state 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, 1-2 + CR21j1, -2 + CR3lJ+1, 1-2
+ 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.-, 1-2 + CR2, 1-2 + CR31i+1, 1-2
+ CR4*1-2 ,-_ + CRSlj l I- + CR6j -1 + CR71j+l -1 + CR815+2, -1
+CL11I12, + CL2j_1, i + CL3, t + CL,1+, + CLS1+2, i
+ CR91-2, 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 Two-Dimensional Model

For the two-dimensional model, the boundary conditions are essentially the same

as in the one-dimensional 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 1-D 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 Two-Dimensional Model

The solution scheme for the 2-D model is similar to that used in the 1-D case. A

simultaneous set of Ni-2 equations is solved for the interior values of f at each vertical

section (Figure 3.5). The same subroutines used in the 1-D 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 off-column 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 solved-for 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 set-up and set-down 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 root-mean-square (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 set-up and return flow is rapid.

The set-up and set-down 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 two-dimensional 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 two-dimensional 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 user-supplied 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 non-linear 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 2-D 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 well-known

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 best-fit 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
.1-J

0.)=
S-,

E

1>
o-

4-,

I

57

of the return flow profile. In the barred-beach 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 time-dependent 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 "spin-up" 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 cross-differentiation, 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 depth-limited wave

height criteria starting from the knowledge of the setdown at the breakpoint:

oonshore = break + 0.186(hbreak-h)

(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 depth-integrated return flow and is equal

in magnitude to the wave-induced 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 surf-zone 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 bore-like 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 one-third 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 second-derivative 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 trough-to-crest 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 trough-to-crest region and the body of the fluid in non-linear waves, it is

assumed that one-third 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 (no-slip) 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 no-slip 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 1-D 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 (Longuet-Higgins,

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', = -P-z --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 Darcy-Weisbach friction factor and Uave is the depth-averaged 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 open-channel and pipe

flows.

Diegaard et al. (1991) measured the mean and time-varying 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

non-zero, 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 patch-length 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 one-dimensional flows and two-dimensional

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 steady-state 1-D simulations are

compared to solutions of horizontal flows subject to various boundary conditions. The

time-dependent 1-D 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

fourth-order 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 fourth-order equation,

even in one dimension, are difficult to obtain for the desired boundary conditions.

For the 2-D case, obtaining an analytic solution for any time-dependent flow in any

geometry has proven difficult. In light of this problem, the steady-state 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 one-dimensional solutions.

5.1 One-Dimensional Analytic Solutions

5.1.1 Steady-State Simulations

To begin, an analytic solution to the steady, one-dimensional, fourth-order 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 -)
_(1-D ,, = -- z + _- 2 ()2(1 -) (5.4)
4h 4 4 h h

1

Figure 5.1 Definition sketch of one-dimensional 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 surf-zone, this situation could be replaced by some non-zero flow

value, such as the wave-induced 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 1-D numerical model to analytic solution for steady-state
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 non-dimensionalized 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

third-order 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 1-D 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 time-mean 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 Time-Dependent Simulation

The ability of the model to reproduce time-dependent behavior is now investigated.

In this instance the one-dimensional 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 second-order 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 time-dependent version of the steady-state

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 1-D 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 "spin-up to" and spin-down 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 steady-state 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).

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
57	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file