UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 NUMERICAL MODELING OF WAVES AND WAVEINDUCED CURRENTS IN THE NEARSHORE ZONE By YANG ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 PAGE 2 2 2008 Yang Zhang PAGE 3 3 To my mom, dad and my wife PAGE 4 4 ACKNOWLEDGMENTS I am sincerely thankful to my advisor a nd supervisory committee chairman, Dr. Andrew Kennedy, for his guidance, support a nd assistance. I also thank Dr Alex Sheremet, Dr. Tom Hsu, Dr. Renwei Mei, and Dr. Clint Slatton, for their participation on my supervisory committee, and for their assistance. I extend my thanks to th e entire Coastal and Oceanographic Engineering Program faculty, staff and graduate students for th eir help during my study at the University of Florida. I thank my wife for always being there for me with love, encouragement, and support. Finally, I thank my parents, who have been behi nd me every step of the way, providing their unconditional love and support. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS..............................................................................................................4 LIST OF TABLES................................................................................................................. .........8 LIST OF FIGURES........................................................................................................................9 ABSTRACT..................................................................................................................................20 CHAP TER 1 INTRODUCTION.................................................................................................................22 1.1 General Introduction....................................................................................................... .22 1.2 Background......................................................................................................................23 1.2.1 Simulation of Nearshore Waves............................................................................ 23 1.2.2 Simulation of Nearshore Waveinduced Currents................................................. 25 1.3 Study Motivations and Objectives................................................................................... 27 2 REDUCED SPECTRAL WAVE MODEL........................................................................... 32 2.1 Introduction............................................................................................................... .......32 2.2 Literature Review of Wave Models.................................................................................33 2.2.1 Phase Resolving Wave Models............................................................................. 33 2.2.2 Spectral Wave Model............................................................................................ 35 2.2.2.1 Wave spectrum............................................................................................... 35 2.2.2.2 Evolution of wave action balance equation.................................................... 37 2.3 Development of the Reduced Spectral Wave Model....................................................... 42 2.3.1 Introduction to Wave Spectral Moment Equations............................................... 42 2.3.2 Governing Equations of the Reduced Wave Spectral Model................................ 46 2.3.2.1 Closure assumptions and evaluation of integrals........................................... 48 2.3.2.2 Improvement to the original governing equations.......................................... 53 2.3.2.3 Spectral wave breaking model........................................................................ 54 2.3.3.4 Effects of frequency/directional sp reading on wave radiation stresses..........57 2.3.3 Numerical Implementation.................................................................................... 61 2.3.3.1 Temporal differencing....................................................................................61 2.3.3.2 Spatial differencing.........................................................................................63 2. 4 Preliminary Tests......................................................................................................... ...64 2.4.1 Onedimensional Evolution of an In itial Sharp Bump in Deep Water.................. 64 2.4.2 Wave Shoaling over a Onedimensional Sloping Beach.......................................68 2.4.3 Shoaling and Refraction of Spectral Waves over a Planar Beach......................... 70 2.4.4 Transformation of Irregular Waves over a Beach with Periodic Rip Channels.... 73 2.4.5 Physical Effects of Current on Waves................................................................... 74 2.5 Discussion and Summary................................................................................................79 3 STEADY TWODIMENSIONAL NEAR SHORE CIRCULATION M ODEL.................. 100 PAGE 6 6 3.1 Introduction............................................................................................................... .....100 3.2 Formulations............................................................................................................... ...102 3.3 Numerical Procedures....................................................................................................106 3.3.1 Discretization Approach......................................................................................106 3.3.2 Grid and Variable Arrangement.......................................................................... 107 3.3.3 Solution Algorithm: SIMPLE Algorithm............................................................ 109 3.3.4 Discretization of Governing Equations............................................................... 112 3.4 Boundary Conditions.....................................................................................................116 3.4.1 NoSlip Wall Boundary Condition...................................................................... 116 3.4.2 Periodic Boundary Condition.............................................................................. 117 3.4.3 Symmetric Boundary Condition.......................................................................... 118 3.4.4 Inlet Boundary Condition.................................................................................... 118 3.4.5 Outlet Boundary Condition................................................................................. 119 3.4.6 Prescribed Pressure Boundary Condition............................................................120 3.5 Numerical Tests............................................................................................................ .122 3.5.1 Pure Diffusion in a Rectangular Basin................................................................ 122 3.5.2 LidDriven Cavity Flow...................................................................................... 125 3.6 Waveinduced Alongshore Currents on a Planar Beach................................................ 126 4 APPLICATIONS OF THE WAVE/CIRCULATION MODEL ..........................................141 4.1 Waveinduced Alongshore Currents over Barred Beach............................................... 141 4.1.1 Introduction......................................................................................................... 141 4.1.2 Comparison to DUCK94 Fi eld Experiment Data................................................ 143 4.1.3 Effects of Frequency/ Directional Spreading....................................................... 148 4.2 Shear Stability of Alongshore Currents......................................................................... 149 4.2.1 Introduction......................................................................................................... 149 4.2.2 Linear Stability Analysis..................................................................................... 153 4.2.3 Numerical Simulation of Nonlinear Shear Instability......................................... 157 4.2.3.1 Unsteady circulation model.......................................................................... 157 4.2.3.2 Simulation results.........................................................................................159 4.3 Simulation of Rip Currents............................................................................................ 169 4.3.1 Steady Rip Currents............................................................................................. 171 4.3.2 Stability of Rip Currents......................................................................................180 4.3.2.1 Experimental plan.........................................................................................182 4.3.2.2 Results.......................................................................................................... 184 4.3.3 Summary..............................................................................................................191 5 CONCLUSIONS AND RECOMME NDATIONS.............................................................. 225 5.1 Conclusions....................................................................................................................225 5.2 Discussion and Recommendations for Future Work..................................................... 228 APPENDIX A DERIVATION OF MOMENTS OF TH E GAUSSIANSHAPED SPECTRA .................. 231 PAGE 7 7 B EVALUATION OF THE INTEGRALS IN GOVERNING EQUATIONS ....................... 234 C DISCRETIZATION OF TERMS IN THE SHALLOW WA TER EQUATIONS.............. 239 D TWO NUMERICAL ISSUES FOR TH E CIRCULATION MODEL: UNDERRELAXATI ON TECHNIQUE AND MOMENTUM INTERPOLATION METHOD ...... 245 LIST OF REFERENCES............................................................................................................250 BIOGRAPHICAL SKETCH......................................................................................................259 PAGE 8 8 LIST OF TABLES Table page 31 Effect of the underrelaxa tion factor on num ber of itera tion required for convergence and effect of MMIM on the converged computational results....................................... 130 41 Experimental conditions and parameter se ttings for shear instability sim ulations. ak is the bottom friction coefficient, M the lateral momentum mixing coefficient, S the directional bandwidth of incident waves, maxV the maximum longshore velocity of the base mean longshore currents............................................................................... 193 42 Experimental conditions and parameter set ting s for rip current instability simulations. WFinteraction stands for whether wave forcing is computed with the wavecurrent interaction being accounted for; 0 is the angle of wave incidence, ak the bottom friction coefficient, M the lateral momentum mixing coefficient................................. 193 PAGE 9 9 LIST OF FIGURES Figure page 11 Waveinduced longshore currents..................................................................................... 31 12 Rip currents............................................................................................................... ........31 21 Comparison of quantities (existing in the propag ation velocities) approximated using Taylor expansion with exact valu es at various water depths............................................ 82 22 Ratios of the exact radiation stress co mponents to approxim ated ones for Gaussian spectra with varying frequency bandwidth.......................................................................83 23 Ratios of the exact radiation stress com ponents to approxim ated ones for PM spectra with varying peak frequency............................................................................................. 83 24 Ratios of the exact radiation stress co mponents to approxim ated ones for Gaussian directional spectra with vary ing directional bandwidth....................................................83 25 Evolution of an initial sharp bump with incident monochromatic waves of 3 s over constant water depth 30 hm 1st order upwind solutions (left), QUICK solutions (right)................................................................................................................................84 26 Evolution of an initial sharp bump with incident Gaussain wave spectrum with frequency bandwidth.........................................................................................................84 27 Evolution of an initial sharp bump with incident Gaussain wave spectrum with frequency bandwidth 0.2HzS ....................................................................................85 28 The 1D bathymetry with a slope of 1:50 resting on a flat bed with 10 m water depth... 85 29 Shoaling of monochromatic waves on th e 1D bathym etry (Figure 28). Numerical results (solid line), analytical solutions (circles). Nonbreaking waves assumed for comparison purpose.......................................................................................................... 86 210 (a) The JONSWAP wave frequency spectrum and its approximations using the Gaussianshaped spectra: Gaussian spectrum 1 has the same peak frequency 0.25Hz as the JONSWAP spectrum and the frequency bandwidth /0.11mS ; Gaussian spectrum 2 has the same mean frequency and moment 2,0 E as the JONSWAP spectrum. All the three spectra give the same RMS wave height of 1m. (b) Comparison of analytical results for s hoaling of the JONSWAP spectrum, the two Gaussianshaped spectra and modeled resu lts for the two Gaussian spectra. Nonbreaking waves assumed for comparison purpose............................................................ 86 211 Shoaling of irregular waves with various frequency bandwidths on the 1D bathym etry (Figure28). (a) 0.1HzS ; (b) 0.2HzS ; (c) 0.3HzS .....................87 PAGE 10 10 212 Effects of frequency bandwidth on th e wave shoaling over the 1D bathymetry (Figure28). ................................................................................................................... ....87 213 Comparison of wave height and direction between num erical re sults and analytical solutions for monochromatic wave propagation over alongshore uniform 2D bathymetry (water depth profile is the same as Figure 28)............................................. 88 214 Transformations of spectral waves over the s loping beach (Figure 28). Incident Gaussian waves are defined: 2/4m 0S /6m and various 00010,20,30 S ............................................................................................................... 88 215 Comparison of computational results with analy tical solutions for offshore wave spectra: /6m 0.2Hz S o10S ....................................................................... 89 216 Bottom contours of the beach bathym etry with periodic rip channels. ............................ 89 217 Contour plots of computed wave height and wave direction f or normally incident monochromatic waves. 0.5 Hm and 4 Ts at the offshore boundary......................... 90 218 Contour plots of computed wave height, m ean direction and directional bandwidth for irregular waves: 0.5rmsHm 4mTs o0m 0S and o10S at the offshore boundary............................................................................................................. 90 219 Contour plots of computed wave height, m ean direction and directional bandwidth for waves: 0.5rmsHm 4mTs o0m 0S and o20S at the offshore boundary........................................................................................................................... 91 220 Comparison of rmsH along transect s1 and s2 (red lines) for various directional bandwidths: o0S o10S and o20S ................................................................... 91 221 Contour plots of computed m and S for waves: 0.5rmsHm 4mTs o0m 0.3HzS and o0S at the offshore boundary............................................................ 92 222 Comparison of rmsH along transect s1 and s2 (red lines) for various directional bandwidths: 0S 0.1HzS 0.2HzS and 0.3HzS ..................................... 92 223 Contour plots of computed wave hei ght and wave direction f or monochromatic waves: 0.5 Hm 4 Ts and o30 at the offshore boundary................................... 93 224 Contour plots of computed wave height, m ean direction and directional bandwidth for waves: 0.5rmsHm 4mTs o30m 0S and o10S at the offshore boundary........................................................................................................................... 93 PAGE 11 11 225 Contour plots of computed wave height, m ean direction and directional bandwidth for waves: 0.5rmsHm 4mTs o30m 0S and o20S at the offshore boundary........................................................................................................................... 94 226 Comparison of rmsH along transect s1 (black lines ), s2 (red lines) for various directional bandwidths: o0S o10S and o20S ................................................ 94 227 Comparison of m along transect s1 (black lines ), s2 (red lines) for various directional bandwidths: o0S o10S and o20S ................................................ 95 228 Contour plots of computed m and S for waves: 0.5rmsHm 4mTs o30m 0.3HzS and o0S at the offshore boundary............................................................ 95 229 Monochromatic waves obliquely incident on a longshore current over constant deep wat er 100 hm The offshore wave conditions are 00.61Hm 4 Ts and 0 022.5 (a) longshore velocity distribution; (b) wave number along xdirection, red dash line indicates yk; (c) wave direction; (d) normalized wave energy distributed along xdirection............................................................................................. 96 230 Unidirectional waves with various fre quency bandwidths obliquely incident on a longshore c urrent over constant depth. The offshore wave conditions are the same as Figure 229 but with various frequency bandwidths: 0S 0.1mS 0.2mS and 0.3mS ................................................................................................................. 96 231 Waves with various directional bandwi dths obliquely incident on the alongshore current over constant depth. The offs hore wave conditions are the same as Figure 229 but with various di rectional bandwidths: 00.15.73Srad 011.46S and 017.19S ........................................................................................................................ 97 232 Monochromatic waves obliquely incident on the alongshore curre nt over a planar beach with the slope of 1:20. For co mparison pu rpose, nonbreaking waves are assumed............................................................................................................................. 97 233 Monochromatic waves with paralleling currents over cons tant water depth of 100m The offshore wave conditions are 00.61Hm 4 Ts (a) crossshore velocity distribution; (b) wave nu mber along xdirection; (c ) normalized wave energy............... 98 234 Waves with various frequency bandwidths with paralleling currents over deep water. The offshor e wave conditions are the same as Figure 233 but with various frequency bandwidths: 0S 0.1mS 0.2mS and 0.3mS ....................................... 98 PAGE 12 12 235 Waves with various directi onal bandwidths with parallel ing currents over deep water. The offshor e wave conditions are the sa me as Figure 233 but with various directional bandwidths: 00.15.73Srad 011.46S and 017.19S .................... 99 236 Monochromatic waves with paralleling currents over a plan e beach with the slop e of 1:20. The offshore wave conditions are 00.61Hm 4 Ts .......................................... 99 31 Typical nonorthogonal 2D c ontrol volum e and notations............................................ 131 32 Types of variable arrangement: collo cated (left) and staggered (right).......................... 131 33 Typical 2D control volume and auxiliary nodes used................................................... 131 34 Control volume with neighboring CVs........................................................................... 131 35 Example of the periodic boundaries............................................................................... 132 36 Sketch of Poiseuille flow (left) and Couette flow (right). .............................................. 132 37 Comparison of the velocity profile between m odeling results and analytical solutions for Poiseuille flow; Comparison of norma lized velocity distribution with various values for Couette flow between modeling results and analytical solutions.................. 132 38 Sketch of the rectangular basin with symmetry boundaries. .......................................... 133 39 Forcing distribution and an alytica l solutions for the pure diffusion problem. Stream function contours, vorticity cont ours and the velocity field........................................... 133 310 Comparison of u and v at a sample section betw een modeling results and the analytical solutions.......................................................................................................... 134 311 Computational error as functions of grid num ber (uniform Cartesian grids used); the model is expected to be 2nd order accuracy since 2nd order discretization schemes were used........................................................................................................................ 134 312 The nonorthogonal grid used and velocity com parison along a monitored grid line between computed and exact values............................................................................... 134 313 Computational error as functions of nu m ber of grid (nonorthogonal grid used).......... 135 314 Liddriven square cavity flow geom etry and boundary conditions................................ 135 315 Velocity field, vorticity contour, and com parison of velocities along the two geometric centerlines for Re100 (80 uniform grids used)................................... 136 316 Velocity field, vorticity contour, and com parison of velocities along the two geometric centerlines for Re1000 (120 uniform grids used)............................. 136 PAGE 13 13 317 The 150 nonuniform Cartesian grids used for the Re5000 .............................. 137 318 Velocity field, vorticity contour, and com parison of velocities along the two geometric centerlines for Re5000 (150 nonuniform grids used)..................... 137 319 Number of iterations re quired to reach convergence using various velocity underrelaxation factors (Re100 the 30 uniform grids used)........................................ 138 320 Planar beach with a slope of 1:20................................................................................... 138 321 Computed results of (a) energy dissipa tion due to waver breaking, (b) viscosity coefficient (c) wave horizontal orbital veloc ity near water bottom, (d) bottom friction coefficient ...................................................................................................... 138 322 Effect of the bottom friction factor ak on the alongshore curren t distributed along the crossshore distance........................................................................................................ 139 323 Effect of the lateral mixing coefficient M on the alongshore current............................ 139 324 Alongshore momentum balance fo r com putational results with 0.25M and 0.03akm ...................................................................................................................... 139 325 Comparison of wave properties with and without including eff ects of wavecurrent interaction. (a) norm alized wave height ve rsus crossshore dist ance; (b) wave angle variation in the crossshore direction.............................................................................. 140 326 Comparison of alongshore current profiles between neglecting wavecurrent interaction case and in cluding interaction case............................................................... 140 41 Sea bottom elevation relative to mean sea level versus crossshore distance on October 14, 1994 at DUCK, NC and selected crossshore locations of pressure sensors and /or current meters.......................................................................................... 194 42 Predicted alongshore current versus cro ssshore distance without roller model and including the rollers with vari ous wavefront slope values: 0.02 0.05 0.1 and 0.2 ........................................................................................................ 194 43 Computational results of (a) Radiation stress x yS, (b) Wave forcing in crossshore direction.......................................................................................................................... 194 44 Modeldata comparison: (a) wave height rmsH; (b) alongshore velocity; (c) bathymetry used.............................................................................................................. 195 45 Alongshore momentum balance..................................................................................... 195 PAGE 14 14 46 Effects of frequency spreading on the modeled (a) rmsH, (b) x yS and (c) alongshore velocity. Bathymetry, offshore rmsH, mean period and direction are the same as in Figure 44........................................................................................................................ 196 47 Effects of directional spreading on the m odeled (a) rmsH, (b) x yS and (c) alongshore velocity. Bathymetry, offshore rmsH, mean period and direction are the same as in Figure 44........................................................................................................................ 196 48 The artificial longshore cu rrent profile of the study case Bowen and Holm an (1989), a constant water depth of 1 m is used............................................................................. 197 49 Results of linear stability calculations fo r the alongshore velocity profile (Figure 48) in term s of growth rate and propagation ve locity versus alongshore wave number....... 197 410 Spatial form of shear wave associated with the most unstable linear mode (1 00.026km ). (a) stream function pattern; (b) fl uctuation velocity field; (c) total velocity pattern, the shear wave has been scal ed so that its peak magnitude equals the peak of the mean alongshore current.............................................................................. 198 411 Bathymetry for a barred beach using the form ulation given by Slinn et al. (2000), which is an approximate fit to topogra phy measured at Duck, North Carolina, on October 11, 1990............................................................................................................. 198 412 Linear stability calculations for wa veinduced alongshore curre nts over the barred beach for varying bed ap parent roughness, 0.01akm 0.02akm and 0.03akm (a) Mean velocity profiles, (b) bottom fric tion coefficient, (c) growth rate versus alongshore wavenumber, (d) propagation velo city versus alongshore wavenumber..... 199 413 Predicted steady alongshore velocity () Vx for various experimental conditions..........199 414 Time series of velocitiesu, v at ,(420,100) x ymm with 0.03akm and 0.5 M for different values of wave directional spreading S from top to bottom 0000,10,20 S(Run name E01, E02, E03)................................................................... 200 415 Snapshots of the contour plots of vorticity at tim es separated by 0.25 h for 0.03akm and 0.5M. Top panels for 00 S bottom panels for 020 S............ 201 416 Alongshoreaveraged alongshore velocity (,) vxt (m/s) and the alongshoreaveraged perturbation kinetic energy density 221 (,) 2 uvxt (m2s2) as a function of the acrossshore coordinate x and time t for 0.03akm ,0.5 M and 00 S ............... 201 PAGE 15 15 417 Comparison of timeand alongshoreaveraged (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) between simulations using different time step sizes to chec k whether the computations results are independent to the time step size.................................................................................... 202 418 Comparison of timeand alongshoreaveraged (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) between simulations using different grid sizes to check whethe r the computations results are independent to spatial resolution......................................................................................................... 202 419 Timeand alongshoreaverag ed (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) for 0.03akm 0.5 M and 0000,10,20 S .............................................................................................................. 203 420 Timeand alongshoreaveraged acrossshore momentum flux 2hu for 0.03akm and 0.5 M and 0000,10,20 S .............................................................. 203 421 Time series of velocitiesu, v at ,(420,100) x ymm with 0.01akm and 0.5 M for different values of wave directional spreading 0000,10,20 S ............... 204 422 Snapshots of the contour plot of vorticity at tim es separated by 0.25 h for 0.01akm and 0.5 M (a) for 00 S (b) for 020 S .............................................................. 204 423 Alongshoreaveraged alongshore velocity (,) vxt (m/s) and the alongshoreaveraged perturbation kinetic energy density 221 (,) 2 uvxt (m2s2) as a function of the acrossshore coordinate x and time t for 0.01akm ,0.5 M and 00 S ................ 205 424 Timeand alongshoreaverag ed (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) for 0.01akm 0.5M and 0000,10,20 S .............................................................................................................. 205 425 Timeand alongshoreaveraged acrossshore momentum flux 2hu for 0.01akm and 0.5 M and 0000,10,20 S .............................................................. 206 PAGE 16 16 426 Snapshots of the contour plots of vorticity at tim es separated by 0.25 h for 0.3M. Top panels for 00 S and 0.03akm ; middle panels for 010 S and 0.0235akm ; bottom panels for 020 S and 0.012akm ........................................ 206 427 Timeand alongshoreaverag ed (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) for 0.5 M and varying 0000,10,20 S .............................................................................................................. 207 428 Time series of velocities u, v at ,(420,100) x ymm with 0.01akm and 00 S for different values of mixing coefficient from top to bottom 0.5,0.3,0.2 M.............................................................................................................. 207 429 Snapshots of the contour plot of vorticity at 10ht for 0.01akm 00 S and varying lateral mixing coefficient M ............................................................................. 208 430 Timeand alongshoreaveraged perturbation kinetic density 221 () 2 uvx (m2s2) for 0.01akm, 00 S and varying lateral mixing coefficient M .............................. 208 431 Timeand alongshoreaveraged acrossshore momentum flux 2hu for 0.01akm, 00 S and varying lateral mixing coefficient M .................................... 209 432 Comparison of timeand alongs horeaveraged alongshore velocity v (dashed lines) with the steady current V(solid lines) for 0.01akm 00 S and varying lateral mixing coefficient M .......................................................................................... 209 433 Time series of velocities u, v at ,(420,100) x ymm with 0.5M and 00 S for different bottom roughness from top to bottom 0.03,0.02,0.01,0.005akm .......... 210 434 Contour plots of vorticity at tim es separated by 0.25 h for 0.005akm and 0.5 M 210 435 Timeand alongshoreaverag ed (a) alongshore velocity v (m/s) and (b) perturbation kinetic density 221 () 2 uvx (m2s2) for 0.5M 00 Sand varying bed roughness ak ............................................................................................... 211 436 Timeand alongshoreaveraged acrossshore momentum flux 2hu for 0.5M, 00 S and varying bed roughness ak ........................................................................... 211 437 Rip channel bathymetry with 0.95eM ........................................................................ 212 PAGE 17 17 438 Contour plot of the computed wave height and wave direction. (a) withou t the wavecurrent interaction; (b) with the wavecurrent interaction. Parameters are 00 and 01 Hm in deep water, 0.8 M and 0.08akm ......................................................... 212 439 Computed velocity fields: (a) without th e wave current interaction; (b) with the wavecurrent interaction. Parameters are 00 and 01 Hm in deep water, 0.8 M and 0.08akm................................................................................................ 213 440 (a) Acrossshore prof iles of uvelocity at 100 y m the center line of rip channel. (b) Longshore velocity profiles of uvelocity at 100 x m 10 m seaward of the longshore bar. Parameters are 00 and 01 Hm in deep water, 0.8M and 0.08akm...................................................................................................................... 213 441 Computed mean water level (wave set up) relative to the still water level. .................... 214 442 The uvelocity profiles of rip curren t with different values of bottom roughness ak .... 214 443 Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent of rip current, with the bottom roughness ak ......................................................................................... 214 444 Effects of the incident deep water wave height on (a) W ave setup along the center line of rip channel; (b) Acrossshore profiles of the uvelocity at 100 y m................. 215 445 Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent of rip current, with incident wave heights.............................................................................................. 215 446 Variation of bar/rip channel in the lo ngshore coordinate with changing elliptic modulus e M Values are shown for 0,0.50,0.75,0.95,0.99eM with rip channels monotonically decreasing in width e M increases. The solid curve corresponds to 0.95eM, which is the value for most of the test cases................................................ 216 447 The uvelocity profiles of rip curr ent for different rip chann el widths........................... 216 448 Longshore profiles of uvelocity at 100 x m different rip channel widths.................. 216 449 Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent and the width of rip current, with relative rip channel width................................................................ 217 PAGE 18 18 450 The uvelocity profiles of rip current for different values of the wave frequency bandwidth S .................................................................................................................. 217 451 The uvelocity profiles of rip current for different values of the wave directional bandwidth S .................................................................................................................. 217 452 Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv with the wave di rectional bandwidth S .................................................................................................................................... 218 453 Convergence tests for the test case E01. The timeaveraged vorticity fields for different grid sizes and f or different time step sizes....................................................... 218 454 Comparison of acrossshore velocity profiles along the rip channel center f or different grid sizes and time st ep sizes for test case T01................................................ 219 455 Convergence tests for the test case E08. The timeaveraged vorticity fields for different grid sizes and f or different time step size......................................................... 219 456 Comparison of acrossshore velocity profiles along the rip channel center f or different grid sizes and time st ep sizes for test case T08................................................ 220 457 Time series of velocities u and v at the location of ,100,302.5 x ymm (a) for T01, (b) for T02, (c) for T03, (d) for T04....................................................................... 220 458 Snapshots of vorticity contours at various tim es. (a) simulation results of T02; (b) simulation results of T04................................................................................................ 221 459 Contour plots of timeaveraged acrossshore mass flux hu, timeaveraged acrossshore momentum flux huu, and timeaveraged 221 2TKEuuvv for test cases. (a) for T01, (b) for T02, (c) for T03, (d) for T04........................................... 221 460 Time series of velocities u and v for simulation results with offshore incident angles 0000001,2,3,4,6,8 around the rip channel center, at the monitored location of ,100,302.5 x ymm ............................................................................................. 222 461 Snapshots of contour plot of vorticity at various tim es for the experiment cases of 02, 04 and 08............................................................................................... 223 PAGE 19 19 462 Timeaveraged vortici ty, acrossshore mass flux hu, and timeaveraged turbulent kinetic energy 221 2TKEuuvv for simulation cases with various wave incident angles................................................................................................................ 224 PAGE 20 20 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL MODELING OF WAVES AND WAVEINDUCED CURRENTS IN THE NEARSHORE ZONE By Yang Zhang December 2008 Chair: Andrew B. Kennedy Major: Coastal and Oceanographic Engineering Ocean waves and waveinduced currents are tw o of the most important components of nearshore hydrodynamics. This study focuses on transformations of waves and waveinduced currents in the nearshore area through numerical modeling efforts. The first portion of this dissertation documents the development of a redu ced wave spectral model, which describes the evolution equations of wave moments instead of individual wave spectral components as in the full spectral wave models (e.g., SWAN and WAVEWATCH III). This wave model is computationally much cheaper than full spectral models and more accurate than monochromaticbased wave models for wave spectra with finite wave frequency and dire ctional bandwidths, but relatively compact shapes with clear peak dir ections and frequencies. The model can provide information on some of the most important wa ve properties including RM S wave height, mean frequency, mean direction, and frequency and directional spreading. Based on the wave model, formulations for the radiation stresses including the effects of the frequency and directional spreading are derived. The radiation stress co mponents can be overestimated or underestimated by over 50% with the monochromatic approximati on when waves have finite frequency and directional bandwidths. The simulation results showed that the model simulates shoaling and PAGE 21 21 refraction of irregular waves with Gaussianshaped spectra very well for simple bathymetries, and reasonably well for complex bathymetries such as planar beach incised with periodic rip channels. Following that, a 2D finite volume st eady circulation model based on the nonlinear shallow water equations is de veloped using the pressurecor rection solution method and validated by comparing to analytical solutions and to the available benchmark solutions, and the circulation model is then coupled with the reduced wave spectral model. Using the coupled models, waveinduced steady longshore currents are simulated over planar beach bathymetries. For the approximately longshore uniform barred beach at DUCK, both the wave field and the longshore currents are accurately simulated by th e models compared against the DUCK94 data. To simulate the time domain solution of nears hore currents, an unsteady version extended from the steady circulation model is developed. Very low frequency motions, shear waves or shear instabilities of longshore currents, are modeled and compared with results of a linear stability analysis. The steadystate waveinduced rip curr ents and circulations are simulated and the relative importance of various wave and curren t parameters is evaluated. Numerical results suggest that the wavecurrent in teraction produces a negative forcing on the wave forcing, hence decreasing the rip current intensity. Effects of the wave frequency and directional spreading on mean properties of waveinduced longshore cu rrents and rip currents are studied, and the simulation results suggest that the wave freque ncy and especially dir ectional bandwidths have significant impacts on both the mean and the instab ility properties of the longshore currents. For rip currents, however, the impacts are much le ss significant because refraction and breaking are the dominant processes in this case. Finally, num erical experiments on rip currents show that wavecurrent interaction has si gnificant impacts on the unsteady flow and the instabilities are sensitive to a small variation range of the angle of wave incidence. PAGE 22 22 CHAPTER 1 INTRODUCTION 1.1 General Introduction In the twentieth century, coasta l regions have become one of the most important areas on the earth for their economic vita lity and crowded population (Book man et al., 1999). About half the U.S. population now resides within 50 miles of the coastline and the coastal population continues to grow (Dean and Da lrymple, 2002). More recently, re creation and tourism at beaches have become more important economically. During the past several decade s, coastal engineering has received increasing attention with the objec tives of obtaining a bett er understanding of natural coastal processes and developing more effective coastal protection (NRC 1999). Ocean waves and waveinduced currents are tw o of the most important components of nearshore hydrodynamics. Waves can be seen in all sizes and forms almost in every body of water open to the atmosphere, depending on magn itude of the forces acting on the water (Dean and Dalrymple, 1998). Ocean waves usually consist of wind waves and swell waves. Wind waves are generated locally by wind blowing over the water surface. Swell waves are seas that have moved away from the area in which they we re formed. Swell can be waves from a distant storm or combination of waves from different storm systems. As waves propagate from deep water to shallow water, wave transformations, such as shoaling, refraction, diffraction and dissipation, could take place. As water continues to become shallower, waves eventually become unstable and break at some depth, dissipating ener gy in the form of turbulence in the surf zone (Lippmann et al., 1996). The break ing of waves transfers momentum into the water column, thus generating currents (LonguetHigg ins and Stewart, 1964). Depending on the incident wave angle and the specific bathymetry, the waveinduced currents can be either alongshore currents (oblique incident waves on beaches without major alongshore variations) or nearshore PAGE 23 23 circulation cells (shorenormal incidence over complex bathymetries). Since turbulent fluctuations in the surf zone lif t many sediment particles from th e seabed and mix them with the fluid, nearshore currents will carry sediments and drive sediment transport, and eventually shape the evolution of the beach. Historically, a tremendous effort has been invested in predicting nearshore waves and waveinduced currents. A brief lit erature review in simulation of nearshore waves and waveinduced nearshore currents is given in next sec tion. The motivations and objectives of this study are described in Section 1.3. 1.2 Background 1.2.1 Simulation of Nearshore Waves Ocean waves have significant impacts on human activities, which has inspired scientists, researchers and engineers to study and to unders tand them. Water wave theory has a very long history that has been over one and half centuries, and many water wave theories have been developed. Two categories can be roughly classified for all these wave theories: phase resolving wave models, and spectral wave models. The m ildslope wave models and Boussinesq wave models both belong to the phase resolving wave model. Spectral wave models are based on wave action conservation and treat th e wave field as a stochastic phenomenon, belonging to the phasedaveraged model category. To reduce the complexity, monochromatic wave s (waves that only have one frequency and one direction) have been traditionally assumed in developing wave theories in early years. A number of monochromatic wave th eories have been developed, and these include the relatively lately developed mildslope wave models. Berkhoff (1972) first introduced the mildslope equation, which assumes a slowly varying bathym etry and can be used to study the effects of PAGE 24 24 combined refractiondiffraction. Since its debut, th e mildslope equation has been widely used in simulating wave transformations over complex bathymetries and received a great amount of efforts to improve its performance and to extend th e range of its validity. These efforts include developing the parabolic mildslope equations wh ich are much easily to be solved than the original elliptic equation (Radder, 1979; Boo ij, 1981; Kirby, 1986), accounting for the effect of ambient currents (Kirby and Dalrymple, 1984), incl uding the effects of the nonlinearity of waves (Kirby and Dalrymple, 1986; Hedge 1987; Li, et al., 2003). The monochromatic wave models based on mildslope equation have been developed and some of them are widely used now such as REF/DIF1 (Kirby and Dalrymple, 1994). Like the mildslope equations, the Boussinesqtype equations describe evolution of surface elevation and Boussinesq wave models belong to phase resolving model too. Since the in troduction of standard Boussinesq equations for variable water depth by Peregrine (1967), tremendous effort has been devoted to the development of extended equations with improved the dispersion and non linearity properties (e.g. Madsen, et al., 1991, Nwogu, 1993, Wei et al., 1995, Gobbi et al., 2000 ). Numerical models based on these equations have been developed and have been shown to give good predictions in comparison with field data when applied within their range of va lidity. One popular Boussinesq wave model is FUNWAVE, which was developed by Kirby et al. (1998). Monochromatic wave theories provide insights about the na tures of waves and build the theoretical foundation for studying th e real irregular ocean waves. In practice, ocean engineers have often utilized the monochromatic wave models to investigate the transformation of ocean waves by using a representative monochromatic wave to approximate the irregular waves. However, many researchers have shown that such an approximation may result in significant errors due to vast dissimilarity in the refract iondiffraction patterns from monochromatic and PAGE 25 25 irregular wave fields even if the wave ener gy spectra are narrowly distributed (Goda, 1985; Vincent and Briggs, 1989; Panchang et al., 1990). Ocean waves have irregular wave heights and periods and the sea surface is continuously varying, which make it infeasible to use a determin istic approach to describe sea surface. On the other hand, statistical propertie s of the surface like mean wave height, mean wave periods and directions appear to vary slow ly in time and space compared to typical wave periods and wave lengths. The surface displacement in the ocean can be seen as the sum of many harmonic waves that are statistically independent. Since the Fourier series approach was applied to random waves, the statistical theory of ocean waves has developed rapidly in past decades. In spectral wave models, waves generally are described with di rectionfrequency (or other forms) variance spectrum (,,,)Ffxt which characterizes the statistical properties of sea surface, or wave action (,,,)/NfxtF Evolution of wave spectrum is described by the spectral action balance equation. Several such models are cu rrently popular and have been applied to open oceans, shelf seas and nearshore regions, including WAM (WAMDI, 1988), WAVEWATCH III (Tolman, 1991; Tolman and Chalikov, 1996), SWAN (Holthuijsen et al., 1993; Ris et al., 1994; Booij et al., 1999; Ris et al., 1999). 1.2.2 Simulation of Nearsho re Waveinduced Currents Waves become unstable and begin to break as they propagate into shallow nearshore areas, dissipating energy in the form of turbulence in the surf zone Water piles up onshore due to breaking, leading an increase of the mean water level (setup). This gradient in mean water level is balanced by the gradient of the excess flow of momentum due to the presences of waves, i.e. radiation stress. The gradients in surface wave level and momentum fluxes give rise to the generation of waveinduced nearshore circulatio n. Two typical wavedriven flow patterns are PAGE 26 26 alongshore currents (LonguetHiggin s, 1970, among many others) a nd rip currents (Shepard et al., 1941, among many others), see Figure 11 and 12. Alongshore cu rrents are often generated by obliquely incident waves breaking on a relatively alongshore uniform beach. Celllike circulations can occur when waves are nearly at normal incidence. Often described as narrow, strong seawarddirected currents that extend from the inner surf zone out through the breaking line, rip currents are integral parts of these circulation cells fed by the longshore currents converging nearly at the shoreline. The introduction of the concep t of radiation stress by Longue tHiggins and Stewart (1964) has made it possible to simulate waveinduced curre nts numerically. In past decades, a variety of current models have been developed and used to study the waveinduced currents in coastal areas. Early models simultaneously developed by LonguetHiggins (1970), Bowen (1969) and Thornton (1970) were reduced to one dimens ion and used to predict the waveinduced alongshore currents. Thornton and Guza (1986) de veloped a model that takes into account irregular waves, and successfully reproduced th e measured alongshore curre nt profile on a planar beach. For alongshore uniform barred beaches, models (Church and Thornton, 1993; Feddersen et al., 1998) predicted longshore cu rrent profiles with maxima not over the bar crests, which differ from the observations. Ruessink et al. (2001) demonstrated that including the rollers in the wave forcing greatly improved the predictions. Waves (nearly) normally incident on alongshore nonuniform bathymetries typically leads to twodimensional horizontal fl ow pattern, i.e., cellular circulat ion or rip current. A variety of twodimensional models have been developed in past decades, early studies including Bowen (1969), Noda (1972), Liu and Mei (1976), and Ebersole and Dalrymple (1980), among others. These models are highly approximated for simplicity, for example, neglecting the nonlinear PAGE 27 27 convective terms and/or neglecting wavecurrent in teraction. More recently, more sophisticated modeling effects have been repo rted. Haas et al. (1998) conduc ted numerical simulation of a laboratory rip current system a nd found that the offshore exte nsion of rip currents can be significantly reduced by consideri ng wavecurrent interaction. Chen et al. (1999) utilized Boussinesq equations to simulate the same la boratory set up. Yu and Slinn (2003) conducted a numerical study of rip currents on a barred beach with gentle sinusoidal alongshore variations. Studies of 3D variations of n earshore currents has also been conducted by many researchers in past years, including Van Dongeren et al. (1994), Putrevu and Svendsen (1999), Haas et al. (2003), and Newberger and A llen (2007), among others. 1.3 Study Motivations and Objectives Currently, transformations of ocean waves in coastal areas are simulated either using the full spectral wave models or using the monochr omatic wave models. Traditional spectral wave models have full directional a nd spectral description of the waves from offshore to onshore and do not need a priori assumption of the spectral shape. However, in spectral wave models the wave spectrum needs to be solved in spectral sp ace, besides in the temporal and spatial space, which makes them very timeconsuming. The complexity of these models, meanwhile, makes them difficult to use. Recently, methods using monochromatic wave m odels to compute dir ectional spectral sea state, based on discretiza tion of the offshore wave spectrum in to individual directionfrequency components, have been developed. Theoretical background of this kind of methods can refer to Izumiya and Horikawa (1987). Numerical models based on this idea have been developed including Panchang et al (1990) and Grassaa (1990), among others. Chawla et al (1998) using this approach extended the monochromatic wave m odel REF/DIF to a spectral version. In these PAGE 28 28 models, the individual wave components are co mputed separately using monochromatic model, and the stochastic characteristics of the spectrum are obtained by assembling the wave components by linear superposition. This type of methods can give good results in certain circumstances, however, it is unable to predict wavewave interaction due to the assumption that individual wave components propagate independently. Besides, intuitively, one can tell this approach is not computationally efficient becau se it needs to go through the whole computation process for every single wave component. Many beaches are dominated by narrowbanded wa ves or swell waves. In this case, using the full spectral wave model is expensive and un necessary if a reduced spectral wave model, which could be order of magnitudes faster and also gives quite good re sults, is available. So it is worthy of efforts to develop a relatively simple spectral wave model as an alternative of the traditional spectral models under certain circumstan ces. One of main purposes of this study is to develop a reduced wave spectral model describing the evolution of wave moments instead of wave action, which is the fundamental variable in the full spectral wave models. Extended from the concept of moments of a fr equency spectrum (Holthuijsen, 2007, among others), the general wave moment (,)nEtx could be defined by: 0(,) ,(,;,)nn E twFtdd xx (11) where ,nw are designated weighting f unctions, and the subscript n is the index of weighting function and wave moment. With th e definition of wave moment, the evolution equations of wave moments can be obtained by integrating the standard wave action balance equation multiplied by the weighting functions ,nw over angular frequencies and directions. Compared with conventional sp ectral wave models, in the present model, the dependent PAGE 29 29 variables are wave moments, and any wave moment (,)nEtx is only a function of time and geographical coordinates. As discussed above, a momentsbased wave mode l is an approximation to the full spectral wave models, but it could be computationally mu ch cheaper. To solve the wave action balance equation (the basic equation used in full spectral wave models), at every computational location and time the wave spectrum is usually discretized into tens of frequency bins and direction directional bins. For instance, 30 bins in frequency and 36 bins in direction are typical discretizations, which will result in 1080 components. Combined with dimensions in space and time, it is very computationally expensive to solve the wave spectrum completely. Since moments of a wave spectrum are independent from the specific energy distribution over frequencies and directions, equati ons of wave moments do not need to be solved in the spectral domain. Instead, several equations of wave mome nts will be solved only in space and time (for unsteady problems). Therefore, w ithin its applicative range, a mo mentsbased wave model that gives sufficiently accurate result s with much less computational e ffort is very useful. On the other hand, the momentsbased wave model solves wave moments, which can contain information on other important wave parameters than the RMS wave height, mean frequency, and mean wave angle provided by using the monochromaticbased models, for example, wave frequency bandwidth and directional bandwidth Consequently, the momentsbased model is more accurate that the monochromaticbased models. In addition, the monochromatic approximation tends to do a very poor job in calcu lating the wave radiation stresses (see Section 2.3.3.4, chapter 2), especially when the wave spectrum have finite frequency bandwidth and directional spreading. When incorporated with a circulation model as the wave driver, the momentsbased model can give a much better estimate of wave radiation stresses. PAGE 30 30 As for modeling nearshore waveinduced curr ents, most of current models are timedependent and developed by finite difference methods, and some of them are very computer intensive and require very small time steps to reach steadystate solutions. For a stationary wave field, a steady waveinduced current field usually re sults. In this case, a steady circulation model that is robust and computationally much less tim e consuming will be very useful in practical sense. Up to date, to our know ledge steady current models exclusively developed for nearshore waveinduced circulation are not widely ava ilable. Ruessink (2001) developed a 1D steady model to study the waveinduced longshore curren ts. For complex bathymetries, a 2D model is definitely necessary. In this study, a 2D steady circulation model based on the nonlinear shallow water equation is developed usi ng the finite volume method. The model used the pressurecorrection method solving the govern ing equations, which is suitable for steady flow problems, and is robust and efficient. Another main goal of this study is utilize th e newly developed wave model and circulation model as research tools to i nvestigate various aspects of n earshore waves and waveinduced current. Objectives of this research can be summarized as follows To develop and validate the mome ntsbased spectral wave model To develop and validate a robust and efficient steady 2D nearshore circulation model Using the wave model to investigate the e ffects of wave freque ncy and directional bandwidth on wave transformations and on the wave radiation stresses Incorporate the two models to reproduce the waveinduced alongshore current data, to investigate the effects of wave frequency and directional bandwidths on waveinduced longshore currents and rip currents To develop a unsteady circulat ion model to investigate the very low frequency motions of nearshore currents, including shear instabi lities of alongshore curr ents and rip current instability To investigate effects of wavecurrent in teraction on waves as well as on rip currents PAGE 31 31 Figure 11. Waveinduced longshore currents. Figure 12. Rip currents. PAGE 32 32 CHAPTER 2 REDUCED SPECTRAL WAVE MODEL 2.1 Introduction Tremendous effort has been devoted to near shore wave modeling. Two types of wave models can be roughly classified: the phaseresolving models and phaseaveraged models (spectral models). The former models include the two well known and in tensely investigated categories: mildslope wave m odel and Boussinesq wave model. The spectral wave models are developed based on wave action conservation. In the early years of nearshore wave m odeling, the coastal engineers typically approximated the offshore irregular sea state by a representative monochroma tic (e.g., significant) wave. Great achievements in developing the m onochromaticbase wave models have been earned in past years, and many coastal and ocean engineering projects ha ve benefited greatly by these achievements. However, many researcher s such as Goda (1985) (using an analytical approach), Vincent and Briggs (1989) (by experimental study), and Panchang et al. (1990) (by numerical study) have shown that such an approximation may result in significant errors even if the wave energy spectra ar e narrowly distributed. In the past several decades, extensive research efforts have been invested in developing spectral wave models and several models are now available and widely used including WAM (WAMDI, 1988), WAVEWATCH III (Tolman, 1991; Tolman and Chalikov, 1996), and SWAN (Holthuijsen et al., 1993; Ris et al., 1994; Booij et al., 1999; Ri s et al., 1999). These models entail full directional and spectral description of the ocean waves and do not need a priori assumption of the spectral shape. These models are doing very good jobs in predicting the sea state either in open oceans or in coastal areas The complexity of these models, meanwhile, can make them difficult to use and generally time consuming. PAGE 33 33 Many beaches are dominated by narrowbanded wave spectra or swell waves. In this case, a full spectral wave model may be unnecessary if a simplified spectral wave model, which could be order of magnitudes faster and also gives qu ite good results, is available. In this study, evolution equations of several wave moments of the wave en ergy density spectrum for wave spectra with relative compact shap e with a clear peak frequency a nd direction are derived. Based on these equations, a reduced spectral model is developed, which gives better results than the monochromaticbased models by accounting for th e directional and frequency bandwidth of a wave spectrum. For traditional spectral wave models, the wave spectrum (or wave action) needs to be solved in spectral space, besides in the temporal and spatia l space. On the other hand, the model in this study does not need solve the wave spectrum in spectral space, instead only several wave moments need to be solved. Consequently, the model is computationally more efficient in comparison with spectral wave models. 2.2 Literature Review of Wave Models 2.2.1 Phase Resolving Wave Models Propagation of water waves over irregular bathymetries involves many processes, including shoaling, refrac tion, diffraction and energy dissipati on. Airy wave theory includes the phenomena of shoaling, refraction and energy di ssipation, but not wave diffraction. The combined refractiondiffraction equation, whic h includes diffraction a nd refraction phenomena explicitly, is introduced first by Berkhoff (1972). It is also known as mildslope equation since in deriving the equation it was assumed that the bo ttom slopes are mild. Radder (1979) developed a parabolic mildslope equation, which is easy to deal with the down wa ve boundary and to be solved numerically too. However, the parabolic model is restricted to waves propagating within 45 of the assumed direction. More recently, Booij (1981) and Kirby (1986) have improved the PAGE 34 34 parabolic model and extended th e range of validity of the mo del. The combined refractiondiffraction wave model REF/DIF is developed based on the parabolic mildslope equation (Kirby 1986). Chawla et al. (1998) suggested that rand om sea state can be discretized into wave components and wave components are calculate d simultaneously using the parabolic wave model, and the statistical properties can then be evaluated afterwards. Based on the nondispersive linear long wave theory, Boussinesqtype wave equations have been developed including the lowest order effects of nonlinearity and frequency dispersion. Boussinesq wave equations have provided a sound basis for studying wave transformation in coastal regions since their debut. Peregrine (1967) derived the Boussinesq equations for variable water depth for the first time using the depthaver aged velocities as the independent variables, which are called the standard B oussinesq equations. The standard Boussinesq equations only work well for shallow water areas and weak nonl inearity cases due to the assumption of weak dispersion and weak nonlinearity. Numerous ef forts have been invested to improve the dispersion and nonlinearity properties of the standard Boussinesq e quations. Madsen et al. (1991), Madsen and Srensen (1992) and Nwogu (1993), am ong others, have derived extended forms of Boussinesq equations. Madsen et al improved the dispersion property by introducing additional small value terms, which are governed by the co nstraint of obtaining th e best possible linear dispersion relation. Nwogu used a different method of derivation: the velocity at certain water depth was used as the dependent variable instead of the depth averaged velocity, and the choice of the representative depth was again determin ed by the goal of obtaini ng the best dispersion relation. Both of their equations have been prove d to be valid for water depths from relatively deep to shallow water. To enable the Boussine sq equations to simulate strong nonlinearity, Wei et al. (1995) derived the fully nonlinear Boussi nesq equations by adapting the approach of PAGE 35 35 Nwogu (1993). The equations retain the full representation of fluid kinematics in nonlinear surface boundary condition terms by not weak nonlin earity. Gobbi et al. (2000) further derived a new set of equations by using a ne w dependent variable which is defined as a weighted average of the velocity at two distinct water depths. 2.2.2 Spectral Wave Model 2.2.2.1 Wave Spectrum The sea surface elevation at one location can be described as superposition of many wave trains: cos()iii itat (21) where is the sea surface elevation; ia,i and i are the amplitude, angular frequency and phase of the ith wave component respectively. Ocean wa ves are chaotic and in practice it is impossible to specify the sea state in complete detail. One usually uses a statistical description, in which the probability of finding a particular sea state is considered. A good approximate description is provided by the covariance function of the sea surface elevation. tt (22) The Fourier transform of the c ovariance function is defined as the wave spectrum, or the variance density spectrum 2(),if F ftted (23) The physical meaning of wave spectrum is the density function specifying the distribution of the variance of the seasurface elevation ove r the frequencies of the harmonic components. Also, the integral of wave spectrum over all wave components is propor tional to total wave PAGE 36 36 energy per unit area. From Equation (23), the va riance of sea surface elevation can be given by (Holthuijsen, 2007): 2 0Ffdf (24) and the potential wave energy per un it surface area can be written as 21 2Eg (25) In cases of no ambient currents, the energy (variance) of a wave package is a conserved quantity. With the presence of currents, the en ergy or variance of a spectral component is no longer conserved since there will be some mean momentum tran sfer of waves done by currents (LonguetHiggins and Stewart, 1961, 1962). Fortunately, the wave action density /AE ( is the intrinsic angular frequency) is still conserved (Whitham, 1965). Therefore, the wave action density spectrum /NF is usually chosen as the fundamental quantity to be solved in spectral wave models. When the random seaelevation is treated as a stationary, Gaussian pr ocess, then all its statistical characteristics can be expressed in terms of moments of the variance density spectrum F f (Holthuijsen, 2007): 0 n nmfFfdf, ...,2,1,0,1,2,...n (26) Obviously, the zerothorder moment is equal to the variance: 2 0 0mFdf (27) Compared with Equation (25), it is clear that the zeroth moment of a wave spectrum is simply proportional to wave energy. PAGE 37 372.2.2.2 Evolution of Wave Action Balance Equation Wave energy density spectrum and wave acti on density spectrum have been introduced above. Willebrand (1975) derived the evolution equa tion of wave action for variable water depth and ambient currents. There is an alternative way to derive the wave actio n balance equation. In this section, we first briefl y introduce the derivation by following Willebrand, followed by using the alternative method. It was noted by Willebarnd (1975) that the conservation of wave action holds for every single wave component separately 0nxknnNN t (28) where, 2,2/nnnNta x x is the horizontal gradient operator. n is the absolute frequency, which has the following form nnnkku (29) Wave number vector and wave frequency are linked by conservation of waves and by the dispersion relation too: 0nxnt k (210) 2tanhnnngkkh (211) Willebrand (1975) argued that the action de nsity could be introduced in the form: ,, ,d n nAtdNtkkxkx (212) where the superscript dkon the summation sign indicates that the sum is taken only over values of n for which ,ntkx lies in the range dk around the fixed wavenumber k. Instead of the PAGE 38 38 index n, the initial values ,0nnP kx of the wavenumber vectors are used to identify different wave components: ,,nnNNPt x (213) ,,nnPt kkx (214) Without losing generality, assume nP is independent of x and then Equation (2 12) is equivalent to: ,,,,,,,,nnnnAtdAtDPtdPNPtdP kxkkxx x (215) where, dP is a volume element in P space corresponding to dk in kspace, and ,,n n nDPt P k x is the Jacobian of Equation (214). From Equation (215), ,,,,nnnNAtDPt kxx (216) The derivative of Equation (210) with respect to nP gives: ,, ,,0n nxn nDPtDPt tk xx (217) Substituting Equation (216) into (28), comb ined with Equation (2 17), it follows that ,, ,,0n nxn nAtAt t kx kx k (218) It should be noted that the derivatives with respect to x and t are taken for fixed value of nP Introducing derivatives for fixed ,nt kx, Equation (218) can be further written as 0n nnn n xx n nnnAA A tt k kk k kkk (219) Again from Equation (210) it can be shown that PAGE 39 39 nnn xnxn nt kk k k (220) Substituting Equation (220) into (219) and dropping indices, we finally obtain 0 A AA t kxkx (221) Using the following trivial identity xkkx (222) Equation (221) can be rewritte n in the conservative form: ,, ,, ,, 0gFtFtFt t xk xkx kx kx cu (223) This form conserves the total wave action and also holds in other spec tral coordinates, for example in term of k An alternative way, which may be more stra ight and easier to understand, to derive the wave balance equation is given as follows. Fo r simplicity, the 1D equation including ambient currents is presented and extension to 2D is st raightforward. Starting wi th the 1D wave energy balance equation, 0gNcuN tx (224) in which, /NE is the wave action. 21 8rmsEgH is the total wave energy. For a wave spectrum, N can be regarded as 0NSkk (225) PAGE 40 40 where / SkFk is the wave action density spectrum It should be pointed out that ,,,SkSxtkxt and special care is needed when computing /St and /Sx. Substituting Equation (225) into (224), we arrive at 0gSkkcuSkk tx (226) Thus for any individual component 0gSkkcuSkk tx (227) Expanding the first term in Equation (227) Sk kk SkkSkkSkk tt tk t (228) The second term can be rewritten as gg gcuSkkSkkcucuSkk xxx (229) In which, Sk kk SkkkSk Sk x xkxx (230) Plugging Equations (228)(230) into (227), it gives g g gggSkkcuSkk tx Sk kk SkkSkkSkkcu ttktx Sk Sk kk cukcukcuSk x kx x (231) In which, gg gkkk SkcuSkSk kcukcu txt xx (232) PAGE 41 41 According to conservation of waves, the underlined terms in Equation (2 32) is equal to 0. So, substituting Equation (232) back into (231) and divided by k we arrive at 0gg gSk kSk SkcuSkSkcucu tx xxkkt (233) Since the target equation we are aiming at is 0guh SkcuSkkSk txkxhx (234) Now the problem becomes to prove th at the last term in Equation (233) is equal to the last term in Equation (234). It can be shown that ggSkSk uhk kSkSkcucuk u kxhxx kxkx (235) Using the conservation of waves law again, Equation (235) can be rewritten as ggSkSk uh kk kSkSkcucu kxhxx kxkt (236) Substituting Equation (236) into (234), we obtain the 1D wave action balance equation. Extension to 2D equation is straightforward but lengthy, and the detailed derivation will not be reproduced here. In practice, several modes of wave spectrum have been used. The spectrum,,, Fftx is used in WAM model (WAMDI, 1988) The wavenumberdirection spectrum ,,, F ktx is used in WAVEWATCHIII (1999). The spectrum ,,, Ftx distributing wave energy over radial frequencies and directions is used in SWAN model. The different modes of wave spectrum can be switched using strai ghtforward Jacobean transformations: 2 ,,,gk FfFkFk fc and PAGE 42 42 1 ,,,gk FFkFk c (237) The advantage of the wavenumberdirection spectrum Fk is its invariance characteristics with respect to physics of wave growth and decay for variable water depths. The more traditional directionfrequency spectrum F or Ff is physically easy to understand and may be more convenient to deal with in th e numerical computation processes. 2.3 Development of the Reduced Spectral Wave Model In previous sections, the defi nition of wave spectrum and the physical meaning of the wave spectrum have been discussed. Two methods to derive the wave action balance equation were given. Wave spectra in different spectral c oordinates were briefly reviewed and relations between different wave spectrum forms were discussed. In this section, we will discuss the incentive to develop a momentsbased reduced spectral wave model and the underlying idea for the model. The detailed derivation of governing e quations of the model will be shown, followed by several numerical tests. 2.3.1 Introduction to Wave Spectral Moment Equations The wave action balance equation is the basic equation for all traditional spectral wave models. In this study, it will be the starting point for developing the reduced spectral model. As mentioned above, various modes of wave spectrum have been used as dependent variables in spectral models. The spectrum,,, F tx distributing wave energy over radial frequencies and directions will be used in this study. As discussed above, wave action density (,;,)/ NtF x is the always conserved quantit y rather than the energy densityF. The conservation of wave action equation can be written as (Hasselmann et al., 1973) PAGE 43 43 totS FFFF cc t gcU (238) gchc hs U Uk (239) 1h c khmm U k (240) where is twodimensional geogra phical gradient operator; c and c are propagation velocities in spectral spaces space and space, respectively. h is water depth, s is the coordinate in the direction of and m is the coordinate perpendicular to s. k is wave number vector and k the modulus of k. (,) uv U is the current velocity vector. In Equation (238), the second term on left hand side describes the pr opagation of wave ener gy in twodimensional geographical space. The third term stands for the effects of shifti ng of the intrinsic frequency due variations in water depth and currents. The forth term represents depthinduced and currentinduced wave refraction The right hand side term totS is the source/sink term which presents all physical processes of wave generation, dissi pation, nonlinear wavewave interaction. To solve the wave action balance equation numerically, the con tinuous wave spectrum (,) F is always approximated at a finite number of discrete frequencies and directions. In addition, the evolution equation ha s to be discretized into com putational spatial grids to be solved numerically. Thus, a tradition spectral mode l is actually a four dimensional model besides time step marching for unsteady cases. Obviously, it is computationally heavily laden, especially when accuracy is required which means high resolutions both in spatial and spectral spaces are necessary. The tremendous computational costs may be greatly reduced if full discretization of the wave spectrum into spectral space can be avoi ded. This actually is able to be achieved by looking at the evolution of mo ments of wave spectrum instead of wave spectrum itself. PAGE 44 44 With the concept of moments of a frequency sp ectrum discussed in Section 2.2.2.1, for any frequencydirection spectrum(,;,)Ft x, wave moment (,)lEtx corresponding to an arbitrary weighting function ,lw could be expressed as 0(,) ,(,;,)ll E twFtdd xx (241) where, the subscript l is the index of weighting function and wave moment. With the definition of wave moment, the evolution equations of wa ve moments can be obtained by integrating the standard wave action balance equation multiplied by the weighting functions ,lw over angular frequencies and directi ons. The dependent variables in the derived equations are wave moments, and any wave moment (,)lEtx is only a function of time a nd geographical coordinates. Therefore, the governing equations do not need to be solved in spectral do main and consequently the computational burden can be greatly relieved. There may be many possible forms one can choose for the weighting function. A good weighting function should be eas y to treat mathematically, and more importantly the resulting wave moments can provide information about the most important wave parameters such as the RMS wave height, mean frequency, mean wave angle, and wave frequency and directional spreading etc. In this study, the weighting function nime is used and the corresponding moment ,nm E is given by 0(,) (,;,)nim nmEteFtdd xx (242) where, ,0,1,2... nm are the indices of moments. It should be noticed that when1n and 0m Equation (242) reduces to the definition of 0,0 E which is the variance of surface elevation 2 as discussed in Section 2.2.2. With the de finition of wave mome nts given by Equation (2 PAGE 45 45 42), the mean angular frequency m is given by 1,00,0/EE, which is less dependent on highfrequency noise of the wave spectrum compared to the mean zerocrossing period (Holthuijsen, 2007). According to the definition proposed by LonguetHiggins(1975), the frequency bandwidth parameter is given by 1/2 2 0,02,01,0/1SEEE Mean wave direction and directional bandwidth can be computed using moments 0,m E and their specific formulations will be given in late sections. To derive the ,nm E moment equation, we start with the wave action ba lance equation. Multiplying Equation (238) by 1nime and then putting it inside all the derivatives by using the product rule of derivative, but dropping the ri ght hand side source te rm temporarily for convenience, we will get 1 + 1 0nim nim nim nim nim nimeF eFceF t ceFnceFimceF gcU (243) Integrating Equation (243) over frequencies and directions, comb ined with the definition of moment ,nm E gives: 00 1 00 0 + 1 0nim nim nm nim nim nimEe F d dc e F d d t ceFddnceFdd imceFdd gcU (244) Physically, since (0)0 F and ()0 F for sufficiently small n 00nimceFdd (245) PAGE 46 46 00nimceFdd (246) The wave spectrum moment equation for ,nm E is reduced to 1 00 01 0nim nim nm nimE eFddnceFdd t imceFdd gcU (247) This equation, actually, is not truly a closed system since wave spectrum (,;,) Ft x is still in the equation. If we can find a way to approximate the three integrals in Equation (247) in terms of moment nmE and/or moments of other orders, a system of evolution of moments of spectrum can then be established. Details in establishmen t of such a system will be described in next section, including the closure assumptions that w ill be made and methods used to evaluate these integrals. 2.3.2 Governing Equations of the Reduced Wave Spectral Model Theoretically, the m oments of wave spectrum can be taken to any order as necessary. However, the values of higherorder moments are rather sensitive to noise in the highfrequency range of the spectrum (Holthuijs en, 2007) and problems with convergence may be encountered for too many moments involved. Generally, mo re information can be obtained from more moments. At the same time, however, more mome nts means a more complicated system. So in practice, the number of moments used should be determined by the specific wave spectrum and by the accuracy requirement of problem interest ed. In this study, a system involving five moments, 0,01,02,00,1,,,REEEE and 0,1 IE is developed. This fiveparameter system is chosen on purpose because it is simple but still contains information on: wave RM S height, mean period, mean direction, frequency bandwidth and directional ba ndwidth. The system itself can be used to PAGE 47 47 replace the traditional spectral wave models to some degree, especially when welldefined wave spectra involved, such as narrowbanded Gaussiant ype spectra. What is even more important is that this system can serve as a foundation for a more comprehensive model that is able to simulate accurately the realistic wave spectra. Referring to Equation (247), it is straightforward to give th e governing equations of the five parameter system as follows 00 000 c EFddFdd t gcU (248) 10 0020 EFddcFdd t gcU (249) 2 20 0030 EF d dc F d d t gcU (250) 01 0000iiic E eFddeFddiecFdd t gcU (251) It should be noted that Equation (251) is a complex equation and consists of two equations: the real equation and an imaginary eq uation. It is a five parameter system with five equations and seems to be a closed system. However, we still need to eliminate the (,;,) Ft x from the equations and to make them only in terms of wave moments. As we all know, wave group velocity gc is a function of wave frequency and direction: cos,sin,gchgc. Ambient current is a known quantity and only a function of geographical location and time:(,) t UUx. Propagation velocities in spectral spaces are functions of coordinates of all dimensions: ,,,cctx, ,,,cctx. Now the problem becomes how those integrals in the govern ing equations can be expressed in terms of wave moments. Essentially, once this problem is addressed, the system is closed. PAGE 48 482.3.2.1 Closure Assumptions and Evaluation of Integrals To close the system introduced above, the inte grals have to be able to be expressed in terms of wave moments. In this section, we will list the assumptions required and give detailed description of evaluatio n of the integrals appearing in th e wave moment equations. The first assumption needs to be made is as follows Assumption 1: frequencydirection spectrum (,) F is separable and can be written as 0,0(,) ()FEMD where 2 0,0/8rmsEH is the total wave energy. M and () D are the frequency spectrum and directional spectrum, respectively. They are normalized so that 01Md and 1Dd The real ocean wave spectra are rarely separable, but this provides a springboard to evaluate the unsolved integrals. Now, we look at the in tegral over frequency first. According to the definition of moments of the wave spectrum density it is clear that the integra nds of integrals in Equation (247) should be able to be expresse d as polynomials of the frequency such that the integrals are possible to be evaluated in terms of mo ments. However, neither of integrands gc, c, c is a polynomial of The simplest way to address this pr oblem may be using the Taylor Series approximation, specifically, taking the Taylor ex pansion of each velocity about the mean frequency m For example, the group velocity g c can be approximated by 2 2 21 2m m mgg gg m mcc cc (252) PAGE 49 49 Substituting this approximated g c into an integral, for instance, 0 n gcMd the integral now can be evaluated in terms of ,0n E 1,0nE,, etc. As known, using a higher order of Taylor series is more accurate. However, it also result s in more moments involved and may make the system very difficult to solve. After a negotiation between the two conflicti ng factors, the Taylor expansion up to 2nd order will be used in this study. To avoid significant errors introduced by using the Taylor expansion approximation, the frequency spectra are required to be narrowbanded. As discussed above, besides the signif icant wave height, mean frequency and mean direction, the directional bandw idth and frequency bandwidth may be the other most important wave parameters highly interested in ocean wa ve simulations. The frequency bandwidth can be defined as 1 2 2 1,00,02,01/ SEEE, a larger value of S corresponding to a broader frequency spectrum. Combined with the Taylor expansion approximation, the integral over frequency can be translated into terms of wave moments for a precisely defined spectru m such as a Gaussiandistributed spectrum, JONSWAP spectrum a nd the PM spectrum, among others. It is straightforward to calculate the mean frequenc y and frequency bandwidth if a Gaussianshaped spectrum is assumed owing to its formulation: 211 ()exp 2 2mM S S (253) in which, m is the mean intrinsic frequency; S represents the standard deviation (bandwidth) of the spectrum in frequencies. Advantages of using Gaussianshape wave spectra include: (1) they are precisely defined by spectra l parameters of interest. (2) It is relatively easy to evaluate integrals. In this model, the Gaussianshaped sp ectrum is assumed although other spectra can be used too. Therefore, we introduce the second assumption as follows PAGE 50 50 Assumption 2: narrowbanded Gaussianshaped frequency spectra are assumed and /0.25mS to ensure the narrowbanded limitation. Ocean wave spectra generally contain a combin ation of local windgenerated seas with swell waves, and their shapes are usually irregular and might be very co mplex. It is also obvious that no perfect Gaussianshape wave spectra can be se en in real sea. While, for some coasts where local wind is not the main driver for the genera tion of ocean waves, narrow banded wave spectra or swell waves dominate. Approximation of this kind of spectra using Ga ussianshape spectra may be valid. Moreover, it is instructive to build a system of moments even a crude assumption of the wave spectra. In additi on, the realistic spectra can be well approximately by combining suitably several (or some number) Gaussianshape spectra. A comprehensive model based on this idea can be developed. In the moment equation nmE higher order moments concer ning the directional spectrum ,1 nmE and ,2 nmE appear (see Appendix B). To make a closed system, additional higher order moments are required to be expressed in terms of ,nm E and/or lower moments. Investigation of possible relations be tween these moments using different spect ra suggests that there is no general rule that can be drawn, and the relation depe nds largely on the specifi c type of spectra. For simplicity, the current model is limited to the specif ic type of directional spectra as shown in the following assumption Assumption 3: Gaussianshaped directional spec tra are assumed. That is, 211 ()exp 2 2mD S S (254) PAGE 51 51 where, m is the mean wave direction. S is onesided the directiona l spread, with a larger value implying broader spectrum. Combination of equati ons (253) and (254) provides the description of the target spectra in this model. With the target wave spectra, calculation of the general moment nmE of the spectrum density according to the definition gives 221/2 1 0,0 /2 0 ,! 2 1/2, =0,2,4... !()! 2 0, otherwisemn nk nkk m mSim k nmn Sn kn k E Eee knk (255) where 1/2 nk is the gamma function. Detailed derivations refer to Appendix A. Substituting corresponding values of n and m into Equation (255), it is straightforward to obtain the following relations 1,00,0/m E E 22 2,00,0/mSEE 2/2 0,10,0mSi E Eee (256) Now it is very clear that this system, describi ng evolution of five moments of the spectrum density, can also be described by five spectral parameters: total variance0,0E mean angular frequency m frequency bandwidth S mean direction m and directional bandwidth S With the assumptions and approximations disc ussed above, integrals in Equation (247) can be evaluated in terms of wave moments. Before solving the integrals, we first look closely at the spectral velocities (see Equation (239) and (240)). The wave ray coordinatessand m involved in spectral velocities can be given by a coordinate transformation with a clockwise rotation of the geographic coordinate system x and y. The relation between these two coordinate systems is: cossin sincos sx my (257) PAGE 52 52 From the linear dispersion relation 2tanh gkkh we can have 11 n sinh22 k hkhh (258) where h is the water depth and /gncc is the ratio of group velocity and wave phase velocity. Then propagation velocity in space can be rewritten as 2211 cossincos sin 2uuvv chnn hx y x y U (259) Similarly, the propagation velocity in space is 2211 sincoscossincos sin 2ghhuuvv ccc hxyyx yx (260) In velocity c n is a function of and /2gcc in c is also a function of These two terms need the Taylor expansion appr oximation as the group velocity 2 2 21 2m mmmmnn nn (261) For convenience, let /2gcccc Then, 2 2 21 2m mmmmcc cc cccc (262) Several derivatives used in Taylor series approximations are listed ()111 1gg gcc knc n n kkcknk (263) 2 2 2221'1'''1gg g gcc nknnnn kcknc (264) 1 2gc ccc 2 22 2221 2gc cc c (265) 1 1 c n nk (266) PAGE 53 53 2 22 2'1 1 g nknn c nkc (267) 22 sinh2sinh2tanh2 nhkh n kkhkhkh (268) 223 2 2 244(cosh2tanh2sinh2sec2) '' sinh2tanh2 (sinh2tanh2)nhkhkhkhkhhkh n kkhkhkhkh (269) Errors will be introduced by using Taylor expa nsion approximation, so it is necessary to check the significance of errors. Comparison of g c /gcc and /2gcc between approximation using Taylor expansion and exact values are conducted and results are plotted in Figure 21. The wave period of 5 s is used for water depths ra nging from deep, intermediate to shallow water. Figure 21 suggests that all the three quantities approximated using Taylor expansion are reasonably accurate within a not too broad frequency range. Approximations are very good in shallow water and approximated factor n and /2gcc match the exact values well too in deep water. For frequencies fairly off the mean value, group velocity is not very well approximated in deep water and in intermediate water depth al l approximations have maximum error about 10%20%. For relatively narrowbanded wave spectra cas e which is interested in this study, wave energy belonging to where approximation is not good enough is only a very small fraction of total energy. Overall, it is reasonable to use the Taylor expansion approximation. Substituting these Taylor series approximations (to 2nd power) into the integrals, it is lengthy but not difficult to evalua te them. Detailed derivation and the final governing equation of this reduced wave model refer to Appendix B. 2.3.2.2 An Improvement to the Original Governing Equations Recall Equation (256), 1,00,0/mEE and 1/2 2 2,00,0/mSEE With the narrowbanded limitation /0.25mS the frequency bandwidth S is usually much smaller than the PAGE 54 54 mean angular frequency m This makes S very vulnerable and numerical errors in m which are not of significance can greatly contaminate S Actually, this problem was first observed through some numerical experiments. For init ial condition problems, numerical spurious oscillations were seen and the system is not st able. A more wellconditi oned set of equations is required to make the system more stable and more accurate. One simple way to achieve this is to redefine2,0E while other moments remain unchanged: 2 2,0 0mEF d d (270) Correspondingly, the 2,0E equation becomes 2 2,0 00 02 3 0mmm mmEF d dF d d t F cd d ggcU cU (271) where, 2 2,00,0 E ES (272) From the expressions, the slightly modifi ed system is more wellconditioned because 2,0E is only a function of 0,0 E and S and numerical errors in m will not have significant impact on S any more. Numerical experiments also show a better performance and be tter stability of the modification (not shown). 2.3.2.3 Spectral Wave Breaking Model For irregular waves, the spectral and directi onal details of depthinduced wave breaking are not well understood and little is known about how to simulate numerically. Fortunately, the total energy dissipation due to depthinduced brea king has been studied and can be well modeled based on different theories, such as the borebased breaking model (Battjes and Jassen 1978, PAGE 55 55 Rattanapitikon 1998) and the random wave tran sformation model by Thornton and Guza (1983). Laboratory observations indicate that the shape of a wave spectrum is conserved as waves propagate cross a simple beach profile, based on which Eldeberky and Battjes (1996) developed a spectral version of breaking model of Battjes and Jassen (1978). This breaking model is used in SWAN and the formulation is ,,tot br totD SF E (273) where brS is the source term due to wave breaking. totD is the rate of total energy dissipation due to wave breaking and totE is the total wave energy, defined by 0totEFdd (274) Lim et al. (2003) incorporated Rattanapi tionkon (1998)s breaking model and Thornton and Guza (1983)s breaking model into the SW AN and evaluated their performance by comparing with experimental data, and concluded that model of Thornton and Guza (1983) yielded the smallest errors. In the present m odel, the breaking model of Thornton and Guza will be incorporated, and the formulation will be us ed is a slightly modified version by Whitford (1988), based on a best fit to field data. 2.5 2 3 3 13 11 16rms rms tot pHH DfBC hh (275) where, 11tanh80.99rmsH C h (276) in which p f is the peak frequency, the coefficient 0.8 B (measuring the intensity of wave breaking) and 0.42 PAGE 56 56 To develop a wave breaking model fit to the current reduced wave spectral model, similar operations apply to the breaking dissipation term ,brS Multiplying Equation (273) by nime and integrating over frequencie s and directions, we can have 00 nim nim tot tot br nm tot totDD eSdd eFddE EE (277) In which, 0,0tot E E The root mean square wave height in totD can be evaluated by 0,0 E with the following relation 0,08rmsHE (278) According to Equation (277), the break ing dissipation terms associated with 0,0E ,1,0E 2,0E 0,1E equation are as the following, respectively 0 brtotSddD (279) 0 brtotmSddD (280) 2 2 0 mbrtotSddDS (281) 01 0 00i br totE eSddD E (282) Thus, the dissipation due to breaking has been built into this wave model by including the breaking terms, given by equations (279) to (282). The energy dissipation rate totD is infinitesimally small when no breaking happens but increases to finite value when waves start breaking. It may be worthy men tioning that no assumption has been made in developing this breaking model. Therefore, this breaking mode l essentially does not limit to narrowbanded regular spectra. PAGE 57 572.3.3.4 Effects of Frequency/Directional Spreading on Wave Radiation Stresses Wave radiation stress gradients define the wave forcing, which drives the currents in the surf zone and leads to shorel ine setup (LonguetHiggins and Stew art, 1964). Accurate estimates of the radiation stresses are required in nears hore circulation modeling. For simplicity, radiation stresses computed using the wave statistics (rmsH peak frequency and mean direction) have been used to drive the nearshore circ ulation models (e.g., Church a nd Thornton, 1993; Ruessink et al., 2001). For random waves, however, the wave radia tion stress is also a function of frequencydirectional bandwidth and the monochromatic (o r narrowbanded) approximation can result in significant errors in radiation stresses (Battjes, 1972). Fedders en (2004) compared the exact radiation stresses with approximations for an em pirical wave spectrum, and showed that the component x yS could be overestimated up to 60% using the narrowbanded approximation, which is roughly consistent with the combined Duck94 and SandyDuck observations. The present wave model solves the frequenc y and directional bandwidth in the entire computational domain besides the wave height, m ean direction and mean frequency. Therefore, it can provide a much more accurate estimate of the radiation stress than the monochromatic approximation by including the effect of fr equency and directional spreading. For a monochromatic (single frequency and single direction) wave of height H, the radiation stress can be written as (LonguetHiggins and Stewart, 1964) 21 cos1 2xxSEn 21 sin1 2yySEn 1 sin2 2xySEn (283) where, 21/8 EgH and /gncc For random waves, Battjes (1972) proposed a similar formulation but in terms of wave energy density spectrum Ff To be consistent with the notations used before, the formulation can be written as follows PAGE 58 58 2 01 cos1, 2g xxc SgFdd c (284) 2 01 sin1, 2g yyc Sg Fdd c (285) 0sincos,g xyc Sg Fdd c (286) where g c and c are the group and phase ve locity, respectively, and th ey are functions of the redial frequency Assuming F is separable: ,()() FMD (287) where () M is the function of fre quency distribution, and2 0,0 0() /8rmsMdEH. () D is normalized so that ()1 Dd With Equation (287), the ex pressions for radiation stress components become: 2 0,0 01 ()1()cos () 2g xxc SgMdDdgE c (288) 2 0,0 01 ()1()sin () 2g yyc SgMdDdgE c (289) 0()()sincos ()g xyc SgMdDd c (290) Frequency Spreading To examine the effect of frequency spread ing, the unidirectiona l waves are assumed temporarily, i.e. () D is a delta function. Then Equation (288) becomes 2 0,0 01 (1cos) () ()2g xxc Sg MdgE c (291) PAGE 59 59 For any given () M and water depths, the integral in the above equation can be computed numerically. The amount of errors introduced by using the monochromatic approximation, xx monoS, is evaluated by the ratio 2 0,0 0 2 0,01 (1cos)/()() 2 1 /()(1cos) 2g xx xx mono gccMdE S S Ecc (292) Ratios for the other two components can be derived in the same way and will not be reproduced here. For analysis pu rpose, two types of empirical frequency spectra are examined. One is the Gaussiantype distri bution around the mean frequency m with onesided bandwidth S The other is the Pierson and Moskowta (1964) spectrum: 4 25 4exp1.25 2pgf f Ff f (293) where defines the total wave energy, p f is the peak frequency. For Gaussian frequency spectra, ratios of ex act radiation stresses to approximated ones for intermediate water depth are computed, and they bias from 1 and decr ease monotonously with increasing the frequency spreading bandwidth (see Figure 22), i.e. radiation stresses are overestimated. At the fairly broad spectra /0.5mS ratios are all larger than 0.91. Thus, for Gaussianshaped frequency spectra, the frequenc y bandwidth does not have strong influences on the radiation stresses. For the PM spectra, the ra tios are calculated for varying the free parameter p f and all radiation stress components are underest imated in this case (Figure 23). However, the errors are not sign ificant either. Although a realistic frequency spec trum can by no means be a Gaussian distribution or PM spectrum, this st rongly suggests that the frequency spreading does PAGE 60 60 not change the wave radiation stresses significantly and the monochromatic approximation in calculating radiation stresses could be reasonable for most cases. Directional Spreading In examining the effects of directional spre ading on radiation stress, the assumption of 0,0()mME is made. Substituted into Equation (290), it gives: 0,0()sincosxySgEDd (294) If () D is a delta function, it is straight forward to show that the ratio for x yS is as simple as 2/exp(2)xyxy monoSSS (295) where S the standard deviation of spectrum in directions. The above equation is valid for any nonzero mean wave angle. If 00 both xyS and mono xyS are equal to zero. Expressions for x xS and yyS are 22 211 11cos2 22 1 ()(1cos) 2S xx xx monone S S n (296) 22 211 11cos2 22 1 ()(1sin) 2S yy yy monone S S n (297) If assuming a very small wave angle (i.e. cos1 ) and shallow water, the ratios become 22/(1)/3S xxxx monoSSe and 22/2S yyyy monoSSe respectively. Figure 24 shows how the ratios vary with the wave directional bandwidth S For shallow water and small nonzero PAGE 61 61 mean wave angles, x xS and x yS will be overestimated by using the monochromatic approximation, whereas yyS is underestimated. When the directional spreading 034 S 18%, 51% overestimation of x xS and x yS respectively are observed, and yyS is underestimated by over 50%. 2.3.3 Numerical Implementation The reduced wave model derived above consists of five equations that are nonlinear and intricately coupled. It may be extremely difficult to apply an im plicit scheme if not impossible. So in this study an explicit propagation scheme is used. 2.3.3.1 Temporal Differencing The simplest method for unsteady problems is the explicit Euler method, but it only has the first order truncation error in time and requires very small step sizes to keep stability. The RungeKutta type formulas are widely used b ecause of their good stabil ity characteristics. Another advantage of the RungeKutta methods is that they have a very loose restriction on the step size. The 2nd order formula and 4th order formula are two most widely used ones, which has been proved to be very stable. In the following, application of them to governing equations will be described in detail. For convenience, we rewrite the governing equation in a simple way and here we only write down the 0,0E equation: 0,00,01,02,00,1(,,,)0 EfEEEE t (298) The 2nd order RungeKutta formula can be treated as two steps. St ep 1 calculates the estimators using explicit Euler approximation. Step 2 computes the new time level values using estimators. Step 1: PAGE 62 62 1/2 0,0 0,0 0,01,02,00,1 *(,,,) 2nn nt EEfEEEE (299) The same implementation should be made si multaneously on other moment equations. For simplicity, we will not show all of them. Afte r all moments are computed at time level 1/2n we arrive at Step 2: 1 1/2 0,0 0,0 0,01,02,00,1*(,,,)nn nEEtfEEEE (2100) The 4th order RungeKutta formula is more accurate than the 2nd order. Meanwhile it requires more steps to do and is more expensive. It has four steps, and each step calculates an intermediate estimator. The estimator is then used in the following step to compute the new estimator. The intermediate estimators are interd ependent and the order given in the following must be followed. Step 1: 1/2 0,0 0,0 0,01,02,00,1 *(,,,) 2nn nt EEfEEEE (2101) Step 2: 1/2 1/2 0,0 0,0 0,01,02,00,1* **(,,,) 2nn nt EEfEEEE (2102) Step 3: 1 1/2 0,0 0,0 0,01,02,00,1** *(,,,)nn nEEtfEEEE (2103) Step 4: 1 1/2 0,0 0,0 0,01,02,00,1 0,01,02,00,1* 1/2 1 0,01,02,00,1** 0,01,02,00,1*(,,,)2(,,,) 6 2(,,,)(,,,)nn nn nnt EEfEEEEfEEEE fEEEEfEEEE (2104) PAGE 63 632.3.3.2 Spatial Differencing Several spatial differencing schemes have been built into the model, including the upwind differencing (UD) scheme, the cen tral differencing (CD), the 2nd order upwind scheme and 3rd order upwind method. Upwind (upwave) differencing method can be simply illustrated by 1 ii x x (2105) The upwind scheme is simple, robust and efficien t, and it has been widely applied in CFD calculation. The scheme accounts fo r the direction of the flow, so transportivene ss is naturally built into the formulation. Solution of upwind me thod is bounded, hence no wiggles occur in the solution. The upwind scheme is based on the backward differencing formula so the accuracy is only first order on the basis of the Taylor Series truncation e rror. The main drawback of the upwind theme is that it causes distribution of tran sported properties to become smeared in some types of problems, such as initial value problem (a sharp bump at the beginning of propagation in a computational domain; results a nd figures will be shown in next section). The resulting error has a diffusionlike appearance and is usually referre d to as false diffusion, or dispersive error. Using higher order discretizati on can minimize the dispersion e rrors. Higher order schemes involve more neighboring points and reduce the e rrors by bringing in a wider influence. The central differencing is 2nd order accuracy proved to be unconditionally unstable when used in conjunction with explicit Euler method for convec tive problems. And it is also found that the central differencing theme does not possess the transportiveness property under certain circumstances. The 2nd and 3rd order upwind themes are higher or der schemes which also retain the transportiveness property. For uniform grids, the 2nd order upwind scheme is as follows PAGE 64 64 121 34 2iiixx (2106) The 3rd order upwind scheme can be written as 121 337 8iiiixx (2107) This scheme may be viewed as a parabolic co rrection to linear interpolation for the cell face values, hence its also called the wellknown QUICK (Quadratic Upstream Interpolation for Convective Kinetics). Here the 2nd and 3rd schemes may be combined into a single expression: 11 121 134 2ii iiixx (2108) where, 0 yields to the 2nd upwind theme; 0.75 the QUICK; 1 the central differencing theme. 2. 4 Preliminary Tests The present model was tested by analytically solving the energy conservation equation or the wave action balance equation fo r several simple test cases: (1) 1D evolution of an initial sharp bump in deep water, (2) shoaling on over a 1D planar beach, (3) wave transformations over a 2D beach. 2.4.1 Onedimensional Evolution of an Initial Sharp Bump in Deep Water It is logical to test a new model first for the simple cases and then for complex ones. In this study, a onedimensional version of the wave model was first developed. Evolution of an initial surface sharp bump in deep water is simulated to test the validity of physical description of governing equations and performa nce of numerical schemes. For the 1D problem, the governing equations can be reduced as follows with the assumption of no ambient currents for simplicity, PAGE 65 65 2 0,00,02,0 21 0 2m mggEcEcE tx (2109) 2 1,02,0 1,0 1,0 2,0 2 0,01 0 2m m mgggEE EcEcEc txE (2110) 2 2 2,0 2,02,02,0 2 0,03 20 2m m mgg m gE cc EcE E tx Ex (2111) (a) Monochromatic Waves First, we look at the monochromatic wave case. For monochromatic waves, the 1D unsteady wave energy conservation equation for constant deep water can be simply written as 0 4 gT EE tx (2112) where, 21 8 E gH is wave energy, g is the acceleration of gravity; T is wave period. Equation (2112) is a onedimens ional unsteady advection equation. The exact transient solution of the equation over time t is as simple as ,, 0 4 gT ExtExt (2113) This test is aiming to check the models capability in simulating a pure convective problem. Consider an initial sharp Gaussian profile travels over a flat bottom with water depth 30 hm Wave period is 3 seconds, and group velocity2.342/gcms The computational domain was chosen long enough such that no influences of wa ve conditions at either boundary can reach to the bump profile, and then the conve ction characteristics of the pr oblem can be clearly observed. Simulations are conducted using different numerical schemes discussed in Section 2.3.3. Figure 25 shows the comparison of the exact soluti ons, first order upwind solutions and QUICK solutions. As discussed above, the first order up wind scheme introduces d iffusion error, and the PAGE 66 66 initial profile is spread out considerably and th e peak amplitude has decreased significantly but nowhere has the solution exceeded the initial bounds. In contrast, QUICK method has a very good agreement with the analytical solutions. For the present reduced spectral wave model, the monochromatic wave case actua lly can be regarded as an extreme wave spectrum, 0 S Equations (2109)(2111) are esse ntially equivalent to one anot her, and all of them can be reduced to the wave energy conserva tion equation given by Equation (2112). (b) Irregular Waves For monochromatic waves, the analytical solu tion of the advection equation can be easily obtained as discussed above. For irregular waves, the wave energy conservation equation can not be used any more. Instead, the wave action ba lance equation should be referred. For 1D, no current and constant deep wate r and ignoring source terms, the wave action balance equation is reduced to ,, ,,0 4 gT FxtFxt tx (2114) With the target Gaussian spectra 2 0,011 ,,,exp 2 2mFxtExt S S (2115) The initial surface bump is given by0,0(,0) Ex At any time t the analytical solution of Equation (2114) is 2 0,011 ,,(/4)exp 2 2mFxtExgT S S (2116) And then it is straightforward to compute0,0E m and S PAGE 67 67 The reduced model uses spectral parameters or wave moments as dependent variables. This implies that it has been assumed that an input wave spectrum at offshore boundary will remain as a Gaussian distribution at any time a nd at any computational grids. However, wave components with different frequencies tend to separate from each other, known as the wave dispersion effect. Consequently, th e shape of a wave spectrum may change significantly as wave transformations, such as shoaling, refraction, diffraction and dissipation, take place. Two simulations are presented out of many experiment al runs: shorter time simulation (the total time is 600 s) for the spectrum with bandwidth 0.1Hz S longer time simulation (2000 s) for the broader spectrum with bandwidth 0.2Hz S For the first simulation, the predicted results agree very well with the exact solutions, the wave energy at va rious times simulated almost perfectly as well as the mean fr equency and bandwidth well predicte d, as shown in Figure 26. The last panel gives the comparison of the real wave spectrum with the predicted spectrum at a monitored location at 600ts. Difference between the predicted spectrum and the real spectrum is minor, implying that in this case the input wave spectrum roughly remains as a Gaussian distribution. For the second test, the wave energy distribution is well modeled while spurious oscillations begin to appear in the mean freque ncy and especially in the bandwidth (Figure 27). By checking the spectrum at 5200 x m and 2000ts obviously, the real spectrum has double peaks and can by no means be approximated by a Gaussian distribution. This significant disagreement in spectrum is believed to be th e main reason for the poor simulation of frequency bandwidth, in addition to the errors introduced by using Taylor series expansion. PAGE 68 682.4.2 Wave Shoaling over a Onedimensional Sloping Beach Shoaling of monochromatic waves and irregu lar waves was simulated on a planar beach with the slope of 1:50. The bat hymetry consists of a flat bottom connected to a constant slope shown in Figure 28. In this test, no energy di ssipation was assumed for comparison purpose. (a) Monochromatic waves The offshore boundary condition is given by gene rating steady 4 s, monochromatic waves. The 1D steady energy conservation equati on (no dissipation) can be written as 0gcE x and wave energy can be easily computed analytically using the wave energy and group velocity at the offshore boundary. Figure 29 shows the co mparison between computational results with exact solutions. The model predicts wave shoaling very accurately for monochromatic waves. (b) Irregular waves In irregular wave case, shoaling of waves with the mean period of 4 s and varying frequency bandwidth is predicted. Recall the onedimension, no currents wave balance equation: 0gFcF tx where(,,) FFt x is the wave energy spectrum density. For stationary incoming waves at offshore boundary condition (wav e input), it is essent ially a steadystate problem and the time differencing term / Ft can be dropped. The analytical solution of the equation 0gcF x can be used to test the perf ormance of the model. Since both g c and F are known at the offshore boundary, the energy density F at every location can be computed analytically. It is worth pointing out that the way to deriving the analytical solutions is valid for any spectral shape. Field measurements of ocean wave spect rum have shown that the onedimensional frequency spectrum appears to roughly have a universal shape: the JONSWAP spectrum PAGE 69 69 (Hasselmann et al., 1973). The current wave model uses the Gaussianshaped wave spectra, which are different from JONSWAP spectra to some degree. Therefore, it is necessary to examine the differences between using JONSWAP spectra and using Gaussian spectra as an approximation. Figure 210(a) shows the co mparison of the JONSWAP wave frequency spectrum and its approximations using two different Gaussianshaped spectra. The first Gaussian spectrum use the peak frequency of the JONSWAP spectrum as its mean (or peak) frequency and the frequency bandwidth /0.11mS is chosen to fit the JONSWAP. The second Gaussian spectrum has the same mean frequency as the JONSWAP spectrum and its frequency bandwidth is given so that the moment 2,0E is identical to that of the JONSWAP spectrum. All the three spectra give the same RMS wave height of 1 m. Comparison of analytical results for shoaling of the JONSWAP spectrum, these two Gaussianshaped spectra and the modeled results for the Gaussian spectra is shown in Figure 210(b). Nonbreaking waves are assumed for comparison purpose. For the first Gaussianshaped spectrum, the maximum difference between the modeled wave heights and analytical values is onl y about 1%. The maximum difference between analytical solutions for the JONSWAP spectrum and modeled results using the first Gaussian spectrum is about 3.6% at the sh allowest water depth, and the RMS difference is less than 1%. For the second Gaussian, the analytical resu lts are more different between the JONSWAP spectrum and the Gaussian spectrum due to the vast discrepancy in the spectral shapes. In addition, since the second Gaussian spectrum is too broad and the assumption of relatively narrowbanded spectrum does not hold here, the model did not ma tch the analytical results accurately. Figure 211 shows the comparison between mode ling results and analytical solutions. For 0.1Hz S the model predicts the spatial variation of energy and the mean frequency variation PAGE 70 70 accurately but slightly overestimates the frequency bandwidth. As wave spectra become broad, the model is less accurate but overall giving reasonable results. For 0.3Hz S the frequency bandwidth is poorly predicted and even the wave energy is not predicted very well, suggesting that the model may lose its accuracy or even give bad results if wave spectra are broad. This in fact is what we expect because the narrowbanded wave spectrum is one of the base assumptions for developing the reduced model. Effects of the frequency bandwidth on wave shoaling are also tested (see Figure 212). Genera lly, waves with broader spectra shoal more quickly than narrowbanded waves if there is no any dissi pation and no breaking. This test roughly demonstrated the capability of the reduced model as well as some of its limitations. Overall, the model can simulate shoali ng of waves with Gaussian spectra well, while, considerable errors are observed for broad wa ve spectra. These errors are believed to be introduced because of two reasons. One reason is that some quantities are approximated using the Taylor expansion which is not accurate wh en the wave spectrum is relatively broad. The other reason is that the assumption that an i nput Gaussian spectrum will maintain its Gaussian shape is more likely to break down when the wave spectrum is broad. 2.4.3 Shoaling and Refraction of Spectral Waves over a Planar Beach In Section 2.4.2, shoaling of on a 1D, sloping beach bathymetry with either monochromatic or spectral waves was simulated numerically. In this section, a simple twodimensional test is used to check the capab ility of the model in simulating shoaling and refraction of spectral waves. For slopes with straight bottom contours parallel to the shoreline, the yderivative of any quantity is equal to zero. For monochromatic wave s, conservation of wave energy shows that for any two points of interest, 111222coscosggEcEc According to Snells law, the refraction PAGE 71 71 coefficient is straightforward to compute. Comb ined with the shoaling coefficient, wave energy at any point of interest can be computed exactly using wave conditions at incident boundary. So, the analytical solutions can be easily establishe d for monochromatic waves. For spectral waves, however, the analytical solutions are not as appa rent, and need more efforts. The basic idea to derive the analytical solutions is splitting the total wave energy into many components, and each component corresponds to a dir ection and a frequency and satisfies the energy conservation equation. Each component is solved analytically, and then wave spectrum corresponding to each component can be computed by dividing the bin si zes in frequency and direction. Once having the wave spectra everywhere, the wave moment s can be computed following their definitions. One thing should be aware of is that j varies spatially due to wa ve refraction. Specifically, step 1: ,(,)ijij ijEF (wave energy contained in certain spectral bin ij) step 2: ,cos0igjijcE x (using wave energy c onservation to compute ijE at every grid location) step 3: ,(,)/ijijijFE (obtain the wave spectrum) step 4: compute wave moments and spectral parameters. In this test, the 2D bathymetry used is ydirection uniform and the acrossshore depth profile is the same as the 1D bathymetry (Fi gure 28). Waves of 4 s (mean period if irregular waves), 1 m in height (rmsH if irregular waves) at the offshore boundary, a nd the incident angle (mean direction if multidirectional waves) is 30. Figure 213 shows wave refraction and wave shoaling as the offshore monochromatic waves propagate towards the shoreline. The top panel compares predicted wave heights with exact values. The bottom panel gives comparison of wave direction between simulated values with values given by the Snells law. Clearly, the model predicts wave refraction and wave shoaling very well for monochromatic incident waves. PAGE 72 72 Simulations for waves with a single frequency (i.e. 0 S ) but multiple directions have been conducted. Figure 214 gives the comparison of the simulation resu lts with analytical solutions for waves with various directional bandwidths ooo10, 20, 30 S For relatively narrow directional spreading cases o10 S and o20 S very good agreement between predicted rmsH m S and analytical values is observed. For o30 S as waves approach to the shoreline, the predicted S gradually deviate from analytical values and differences in mean direction and wave height can also be seen. It is also found that wave heights become smaller with the increase in directional bandwidth and bandwidth S itself decreases as waves approach to the shoreline. This is what is supposed to happen because the characteristics of wave refraction suggest that waves tend to travel in th e direction normal to bot tom contours, thus all wave components bend to x direction and consequently bo th the mean direction and the directional bandwidth become smaller. Simulations for incident waves with nonzero frequency bandwidth and nonzero directional bandwidth are also conducted. Model resu lts of waves with 0.1Hz S and o10 S are presented and compared to the analytical results as shown in Figure 215. Acrossshore variations of rmsH m and m are accurately modeled, and noticeable underestimation in S and overestimation in S are observed. Overall, this test demonstrated the models capability to predict shoaling and refraction of irregular waves with certain fr equency bandwidth and directional bandwidth. PAGE 73 732.4.4 Transformation of Irregular Waves over a Beach with Periodic Rip Channels Simulating wave transformation over a beach with periodic rip channels is a great challenge for this model, especially for waves wi th a broad directional spectrum. The bathymetry used in this test is defined by 1/3 8,0.025120exp3/20sin 80 y hxyx x (2117) The bathymetry is similar to that used in N oda (1974). Bottom contours of the bathymetry are shown in Figure 216. Periodic boundary conditi on is applied to the lateral boundaries. Waves of 4 s, 0.5rmsHm at offshore boundary are used throughout this section. Numerical experiments are carried out for two objectives: (1) to check the capability of this model in predicting wave transformations over this relatively complex 2D bathymetry; (2) to investigate the effects of frequency and directional ba ndwidths on wave propert ies. Both the normal incidence and oblique incidence (30 from xdirection) are studied. For each case, monochromatic waves, waves with nonzero frequency bandwi dths and waves with nonzero directional bandwidths are simulated. For normal wave incidence, waves bend from rip channels and consequently wave focusing occurs in the region between rip chan nels (see Figure 217). As the water depth decreases, waves break very str ongly outside rip channels. For in cident waves with directional bandwidth o10 S (Figure 218) and o20 S (Figure 219), the wave height distribution and wave direction are similar to the monochrom atic case. The model predicted very large S in the area near to the shoreline outside rip channels (see Figure 218, 219). The possible explanation is that wave caustics begin to happen and the mo del is not able to give correct values of S Figure 220 gives rmsH along the two transects for cases w ith various values of directional PAGE 74 74 bandwidth. Difference in wave heights due to th e presence of directiona l bandwidth is not of significance but is not negligible either for the case o20 S Modeled rmsH for waves with nonzero frequency bandwidth is very similar to th e monochromatic wave ca se and will not be presented here. Figure 221 gives the distribution of the mean angular frequency and the frequency bandwidth for the case 0.3Hz S It can be seen that both the mean angular frequency m and bandwidth S decreases from offshore to onshore with strong variation outside the rip channels where strong wave refraction happens. Wave height rmsH along the two transects are also compared for different fre quency bandwidths (Figure 222). No appreciable changes are observed with the presence of fre quency bandwidth because wave refraction and breaking are the dominant factors and ove rwhelm others in this case. Another set of experiments is carried out for oblique wave incidence with an angle of 30 from the xdirection. Figure 223 shows the simulated results for monochromatic waves. Similar to the normal incidence case, contour plots of rmsH did not change very much with the presence of directional bandwidth o10 S (Figure 224) and o20 S (Figure 225). This is further demonstrated by comparing rmsH along the two transects for va rious directional bandwidths (Figure 226). The mean wave direction along transect s1 and s2 is only slightly changed by the presence of directional bandwidth (Figure 227 ). Again, the frequency bandwidth does not have considerable impacts on wave field. Figure 228 shows the spatial distri bution of the computed mean frequency and frequenc y bandwidth. 2.4.5 Physical Effects of Current on Waves In addition to the topographical variation, ambient currents also lead to wave transformations. Several general conclusions have been drawn about effects of currents on waves PAGE 75 75 including: when the currents and waves are in the same direction, waves will be lengthened and the wave height will be decreased. However, waves are shorten and steepened by the current in opposite direction, and strong e nough currents can completely bl ock the propagation of waves, resulting in wave breaking. In this section, numerical si mulations of wave transformati ons due to artificial alongshore currents were conducted, and for monochromatic waves the modeling results were compared to the analytical solutions. Good agreement between the analytical solutions and the numerical predictions demonstrated the wave models abil ity to accurately predict wave transformations due to ambient currents. The model then was used to investigate effects of frequency/directional bandwidth upon evolution of waves with the pres ence of these artificial ambient currents. Two simple cases, whose analytical solutions are available, were used: waves obliquely incident on an alongshore current, and waves with currents in the same direction. With the purpose of studying effects of currents on waves, the constant wate r depth was used in most case to avoid wave transformations due to topographical variations. For a sloping beach bathymetry, modeling results for monochromatic wave cases only are presented here. As waves propagate over the sloping beach, effects of frequenc y and direction bandwidth are rela tively small compared to the impacts of wave shoaling, refraction and breaking. Waves with Oblique Incidence on Alongshore Currents With the presence of currents, the absolu te angular frequency can be written as 2/ T Uk (2118) and the intrinsic frequency is defined by the linear dispersion relation 2tanh gkkh (2119) PAGE 76 76 For steady waves and current, assuming no energy i nput or energy dissipa tion, the wave action conservation equation is reduced to 0gE UC (2120) If both the currents and the bathymetry are l ongshore uniform, Equation (2120) can be further reduced to cos/0gEC x (2121) Since /0 y and0 k, it is easy to show that 0sinsinykkkconst (2122) where, the subscript represents values at the reference location (often at offshore boundary). From (2121), the wave energy variation due to the effects of current can be found 00 00cos cosg gc E Ec (2123) in which, can be computed according to Equation (2119) and wave direction is calculated using Equation (2122). The artificial alongshore current profile 4()4sin2/ x vxxl was used for this test case. x l is the computational domain in xdirection. Waves of00.61 Hm 4Ts and 0 022.5 propagate through the alongshore currents. Figure 229 shows the comparison of modeled results with the analytical solutions for monochromatic waves on constant water depth of 100 m. Spatial variation of wave number due to the presence of current is sh own in panel (b), and the red dashed line indicates yk which is constant (fro m Equation (2122)). Very interesting conclusions have been drawn about the relation between yk and the minimum k: If min ykk wave rays can PAGE 77 77 penetrate the current after deflection; if min ykk wave caustics will happen and wave rays incident from outside must bend backward afte r touching the caustic edge; for the extreme situation, 0yk which means waves propagate against th e current, waves can be trapped near the current peaks. For more details refer to Meis textbook (1983). As expected, the wave model failed for wave caustics cases and we will limit ourselves to minykk case. Strong wave refraction around the current peaks is observed, as shown in panel (c), and variation in energy along xdirection is plotted in pa nel (d). The wave number, wave direction and energy variation were accurately simulated by the model. To investigate the effects of frequency and di rectional bandwidths, simulations with same current profile and wave conditions but with nonzero S and S were also conducted. Increasing S from 0 to 0.3 Hz, the mean wave number slightly decreases and also reduces a bit around the current peaks (Figure 230). While, the wave energy is not affected significantly, which can be explained by considering that th e effect of the changes in wave number cancels out that of the changes in wave direction. The effects of directional spreading S are demonstrated in Figure 231, and a weaker current with the same shape is us ed for this run to avoid the possible instability problem. For a large directional spreading 017 S the mean wave angle increases significantly around the current peaks and surprisingly the ener gy variation pattern is completely changed. Wave energy increases instead of decreasing for unidirectional waves or waves with sufficiently small directional spreading. Figure 232 shows the simulation results and analytical solutions of the same wave condition but on a planar beach with a slope of 1:20. In this test, no wave breaking was assumed such that anal ytical solutions are available. Waves with Paralleling Currents PAGE 78 78 For the case of waves with paralleling currents, we can assume that waves and currents are in the xdirection, and then the current is given by (),0 uuxv It should be noted that this kind of current requires a vertic al velocity component to satisf y the continuity equation. To illustrate the ability of the wave model to simu late how depthaveraged currents affect wave propagation, however, attention will only fo cus on horizontal components here. For 1D (alongshore uniform) topography, the wave action conservation Equation (2120) can be written as /0gECu x (2124) and the energy variation can be obtained by 00 0 0 g gcu E E cu (2125) In this set of experiments, the artificial velocity 3()0.6245sin(2/) x ux xl was used. Figure 233 shows the comparison of computational results with analytical solutions of wave transformations due to the paralleling currents over the flat bed. Currents in the same direction of waves lengthen the waves and decrease the wave number. On the other hand, opposing currents increases the wave number and s horten the waves (Figure 233(b)). Wave energy increases when waves meet the opposing currents and decrease s when waves come across the downstream currents, see Figure 233(c). Increasing frequency bandwidth S results in the strengthening of wave number variation due to the current, as shown in Figure 234(b) However, the wave energy was not altered much by increasing S from 0.1m to 0.3m As expected, the directional spreading does not ha ve significant influences since the current is parallel to the propagation direction of waves, which can be seen in Figure 235. Figure 236 shows the PAGE 79 79 comparison of modeling results with analytical solutions for monochrom atic waves over the sloping beach bathymetry, and good agreement was obtained. 2.5 Discussion and Summary A new type of wave model, wh ich is different from the trad itional spectral wave models and approximate monochromaticbased models, ha s been developed and several simple tests have been conducted in this chap ter. Unlike the full wave spectral models, this reduced model describes the evolution of wave moments inst ead of individual wave spectral components. Evolution equation for the general moment nmE was introduced by integrating the standard wave action balance equation multiplied by the weighting function, 1 nime over angular frequencies and directions. As a preliminary study in developing this momentsbased wave model, a relatively simple system with moments 0,01,02,01,0,,, E EEE as fundamental variables, which contain some of the most important wave properties including rmsH mean frequency, mean direction, frequency and directional bandwidth is developed. To uniquely define all terms in the governing equations in terms of wave mo ments, several approximations and assumptions were made including 1. The input wave spectrum is a separable frequencydirection spectrum. 2. In order to use the Taylor expansion appr oximation, the input frequency spectrum was required to be narrowbanded. 3. Only the Gaussianshaped spectra are used in this study and thus the higher order moments existed in equation nmE can be computed using the wave spectral parameters including the wave energy 0,0E mean angular frequency m mean direction m directional bandwidth S and frequency bandwidth S PAGE 80 80 4. An input spectrum was assumed to conserve its spectral shape. All these approximations and assumptions introd uce errors and/or limit the application of the model to some degree. It should be pointed ou t that other wave frequency spectra such as the JONSWAP spectra and PM spectra rather than the Gaussianshaped spectra can be simulated using the evolution equations of wave moments although we limit to Gaussianshaped spectra in this model. In addition, theoretically the directio nal spectrum is not necessary to be a Gaussian distribution either. Waves rarely maintain their spectral shapes as they propagate from offshore to the shoreline. Redistribution of wave energy over the spectrum can be caused by many processes including refraction, shoaling, diffraction and nonlinear wavewave interactions. For the directional spectrum, the shap e at certain locations can be ve ry different from that of the incident spectrum at offshore, especially when strong wave refractiondiffraction happens. Reilly and Guza (1991) investigated the spectral transfor mation of a symmetric, singlepeak directional spectrum over a circular shoal using the b ackward ray tracing method. They found that directional spectra at locations be hind the shoal can be very complex and several peak directions were observed. Despite all the approximations and limitations in applicability, for b eaches typically having wave spectra with finite relatively compact shapes with clear peak direction and frequency, the reduced spectral wave model can give better results than monochr omaticbased wave models and is much computationally cheaper than a full spectral wave model. Therefore, this momentsbased model is useful and can be an alternative to th e traditional spectral wave models within its applicative range. Another advantage of this model is that it can directly investigate effects of directional/frequency bandwidth upon wave fields and then the waveinduced currents. In PAGE 81 81 addition, the realistic spectra can be well appr oximately by combining su itably several Gaussianshape spectra. A comprehensive model base d on this idea can be developed. It is also worth pointing out that this model does not include wind input, energy dissipation due to whitecapping and wavewave interaction. Wave reflection and di ffraction can not be simulated by this model either. So, future studies to lift limitations and improve the performance of the model include: 1. Developing a more sophisticated model which ca n deal with more realistic wave spectra, for example, spectra with multipeaks in both frequency and direction 2. Seeking a sound way to approximate higher order moments of the directional spectrum appearing in equation nmE such that the assumption that a directional spectrum conserve its shape can be removed. This is important because the directional spectrum may change its shape vastly under certain circumstances 3. Accounting for wave diffracti on and wavewave interacti on, which may play very important roles in evolution of spectral waves in nearshore areas, especially for complex bathymetries PAGE 82 82 (a) (b) (c) Figure 21. Comparison of quantiti es (existing in the propagation velocities) approximated using Taylor expansion with exact values at va rious water depths. (a) for group velocity g c ; (b) for /gncc ; (c) for /2gcc solid lines represent appr oximated values, dots for exact values. The vertical dot lines mark the position of m PAGE 83 83 0 0.1 0.2 0.3 0.4 0.5 0.9 0.92 0.94 0.96 0.98 1 S/mRatio Figure 22. Ratios of the exact radiation stress components to approximated ones for Gaussian spectra with varying frequency bandwidth. 2/6m 030 and water depth 6 hm are used. /xxxx monoSS(circle), /yyyy monoSS(xmark), /xyxy monoSS(square). 0.05 0.1 0.15 0.2 0.8 0.9 1 1.1 fp(Hz)Ratio Figure 23. Ratios of the exact radiation stress components to a pproximated ones for PM spectra with varying peak frequency. 030 and water depth 6hm are used. /xxxx monoSS(circle), /yyyy monoSS(xmark), /xyxy monoSS(square). 5 10 15 20 25 30 35 0 0.5 1 1.5 2 S(deg)Ratio Figure 24. Ratios of the exact radiation stress components to approximated ones for Gaussian directional spectra with varying directional bandwidth S 01.0 and shallow water are assumed. /xxxx monoSS(circle), /yyyy monoSS(xmark), /xyxy monoSS(square). PAGE 84 84 Figure 25. Evolution of an initial sharp bump with incident monochromatic waves of 3 s over constant water depth 30 hm. 1st order upwind solutions (left), QUICK solutions (right). Analytical solution, Computed results Figure 26. Evolution of an initial sharp bump with incident Gaussain wave spectrum with frequency bandwidth0.1Hz S Wave energy (first panel), mean frequency (second panel), frequency bandwidth (third panel) and the last panel shows comparison between computed and the real wave spectra at2800 x m 600 ts Analytical solution, Computed results. PAGE 85 85 Figure 27. Same as Figure 26 but with 0.2Hz S; The last panel shows comparison between computed and the real wave spectra at5200 x m ,2000 ts Analytical solution, computed results. 0 100 200 300 400 500 600 700 800 900 1000 1100 10 8 6 4 2 0 2 x(m) =>onshoreh(m) 1:50 Figure 28. The 1D bathymetry with a slope of 1:50 resting on a flat bed with 10 m water depth. PAGE 86 86 500 550 600 650 700 750 800 850 900 950 1000 0.9 1 1.1 1.2 Distance(m) =>onshoreE0,0(m2) Numerical results Analytical solution Figure 29. Shoaling of monochromatic waves on the 1D bathymetry (Figure 28). Numerical results (solid line), analytical solutions (circles). Nonbreaking waves assumed for comparison purpose. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.05 0.1 0.15 0.2 0.25 f(Hz)wave spectral density (m2/Hz)(a) JONSWAP spectrum Gaussianshaped spectrum 1 Gaussianshaped spectrum 2 500 550 600 650 700 750 800 850 900 950 1000 0.9 0.95 1 1.05 1.1 1.15 x(m)Hrms(m) (b) analytical resultsJONSWAP analytical resultsGaussian 1 modeled results for Gaussian 1 analytical resultsGaussian 2 modeled results for Gaussian Figure 210. (a) The JONSWAP wave frequency spectrum and its approximations using the Gaussianshaped spectra: Gaussian spectrum 1 has the same peak frequency 0.25Hz as the JONSWAP spectrum and the frequency bandwidth /0.11mS ; Gaussian spectrum 2 has the same mean frequency and moment 2,0E as the JONSWAP spectrum. All the three spectra give the same RMS wave height of 1m. (b) Comparison of analytical results for s hoaling of the JONSWAP spectrum, the two Gaussianshaped spectra and modeled resu lts for the two Gaussian spectra. Nonbreaking waves assumed for comparison purpose. PAGE 87 87 500 600 700 800 900 1000 0.8 0.9 1 1.1 1.2 1.3 Distance(m) =>onshoreE0,0(m2) 500 600 700 800 900 1000 1.562 1.564 1.566 1.568 1.57 1.572 1.574 Distance(m) =>onshorem(Hz) 500 600 700 800 900 1000 0.1 0.1002 0.1004 0.1006 0.1008 Distance(m) =>onshoreS(Hz) (a) Analytical Computed Analytical Computed Analytical Computed 500 600 700 800 900 1000 0.8 0.9 1 1.1 1.2 1.3 Distance(m) =>onshoreE0,0(m2) 500 600 700 800 900 1000 1.54 1.545 1.55 1.555 1.56 1.565 1.57 Distance(m) =>onshorem(Hz) 500 600 700 800 900 1000 0.195 0.2 0.205 0.21 Distance(m) =>onshoreS(Hz) (b) Analytical Computational 500 600 700 800 900 1000 0.9 1 1.1 1.2 1.3 Distance(m) =>onshoreE0,0(m2) 500 600 700 800 900 1000 1.5 1.52 1.54 1.56 1.58 1.6 Distance(m) =>onshorem(Hz) 500 600 700 800 900 1000 0.3 0.305 0.31 0.315 0.32 Distance(m) =>onshoreS(Hz) (c) Analytical Computational Figure 211. Shoaling of irregular waves w ith various frequency bandwidths on the 1D bathymetry (Figure28). (a) 0.1Hz S ; (b) 0.2Hz S ; (c) 0.3Hz S Nonbreaking waves assumed. 500 550 600 650 700 750 800 850 900 950 1000 0.8 0.9 1 1.1 1.2 1.3 Distance(m) =>onshoreE0,0(m2) Figure 212. Effects of frequency bandwidth on the wave shoaling over the 1D bathymetry (Figure28). Lines from bottom to top stand for 0, 0.1, 0.2, 0.3 (Hz) S respectively. Nonbreaking waves assumed. PAGE 88 88 0.9 0.95 1 1.05 1.1 H(m) 500 600 700 800 900 1000 0 10 20 30 40 (degree)x(m) Exact Computed Snell's law Computed Figure 213. Comparison of wave height and direction between num erical results (solid lines) and analytical solutions (circles) for monochromatic wave propagation over alongshore uniform 2D bathymetry (water de pth profile is the same as Figure 28). Nonbreaking waves assumed. 0.85 0.9 0.95 1 1.05 Hrms(m) 10 20 30 m(degree) 500 600 700 800 900 1000 0 10 20 30 S(degree)x(m) Figure 214. Transformations of spectral waves over the sloping beach (Figure 28). Incident Gaussian waves are defined: 2/4m 0 S /6m and various S (Lines from top to bottom represent 00010,20,30 S respectively. Symbols indicate analytical solutions). Nonbreaking waves assumed. PAGE 89 89 0.8 1 1.2 Hrms(m) 1.5 1.55 1.6 m(Hz) 0.19 0.2 0.21 S(Hz) 0 20 40 m(o) 500 550 600 650 700 750 800 850 900 950 1000 0 5 10 S(o)x(m) Figure 215. Comparison of computa tional results (lines) with analytical solutions (circles) for offshore wave spectra: /6m 0.2Hz S o10 S Nonbreaking waves assumed. 0 50 50 511111 51 51 52222.52 52 53333 53.53 5444x(m) =>onshorey(m)s1 s2 0 20 40 60 80 100 120 140 160 0 20 40 60 80 100 120 140 160 Figure 216. Bottom contours of the beach bathym etry with periodic rip channels. Two acrossshore transects are marked by the thick dashed lines. PAGE 90 90 0 20 20 20 30 30 30 30 40 40 40 4 0 50 50 50 50 .50 50 50 50 60 60 60 60 70 70 70 7x(m) =>shoreliney(m)Wave Height(m) 0 50 100 150 0 20 40 60 80 100 120 140 160 2 0 2 01 0 1 01 01 01 01 02 02 0x(m) =>shoreline Wave Direction(degree) 0 50 100 150 0 20 40 60 80 100 120 140 160 Figure 217. Contour plots of computed wave height (left) and wave direction (right) for normally incident monochromatic waves. 0.5 Hm and 4 Ts at the offshore boundary. 0 2 0 20 2 0 30 3 0 30 40 40 40 50 50 50 50 50 5 0 50 60 60 60 6x(m) =>shoreliney(m)Hrms(m) 0 50 100 150 0 20 40 60 80 100 120 140 160 202 01 01 01 0 1 02 02 0x(m) =>shoreline m(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 1 01 01 01 01 01 01 0 2 0 2 02 0 x(m) =>onshore S(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 Figure 218. Contour plots of computed wave height (left), mean direction (middle) and directional bandwidth (right) for irregular waves: 0.5rmsHm 4mTs o0m 0 S and o10 S at the offshore boundary. PAGE 91 91 0.20 20 20 30 30 .30 40 40 40 50 50 50 50 50 50 50 60 60 .60 6x(m) =>shoreliney(m)Hrms(m) 0 50 100 150 0 20 40 60 80 100 120 140 160 2 0 2 01 01 01 01 02 0 2 0x(m) =>shoreline m(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 2 02 02 02 03 03 03 04 040405 05 0 5 0 6 06 06 06 0x(m) =>onshore S(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 Figure 219. Contour plots of computed wave height (left), mean direction (middle) and directional bandwidth (right) for waves: 0.5rmsHm 4mTs o0m 0 S and o20 S at the offshore boundary. 0 20 40 60 80 100 120 140 160 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x(m)Hrms(m) Figure 220. Comparison of rmsH along transect s1 (black line s), s2 (red lines) for various directional bandwidths: o0 S (solid lines), o10 S (shortdashed lines) and o20 S (longdashed lines). PAGE 92 92 Figure 221. Contour plots of computed m (left) and S (right) for waves: 0.5rmsHm 4mTs o0m 0.3Hz S and o0 S at the offshore boundary. 0 20 40 60 80 100 120 140 160 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x(m)Hrms(m) Figure 222. Comparison of rmsH along transect s1 (black line s), s2 (red lines) for various directional bandwidths: 0 S (solid lines), 0.1Hz S (shortdashed lines), 0.2Hz S (longdashed lines) and 0.3Hz S (dashdot lines). PAGE 93 93 0 20 20 20 3 0 30 30 40 40.40 50 50 50 50 50 50 50 50 50 5x ( m ) =>shoreliney(m)Wave Height(m) 0 50 100 150 0 20 40 60 80 100 120 140 160 1 01 01 01 02 02 02 02 02 03 0x ( m ) =>shoreline Wave Direction(degree) 0 50 100 150 0 20 40 60 80 100 120 140 160 Figure 223. Contour plots of computed wave height (left) and wave direction (right) for monochromatic waves: 0.5 Hm 4 Ts and o30 at the offshore boundary. 0 20 20 20 30 30.30 40.4 0.40 50 50 .50 50 50 50 50 5x(m) =>shoreliney(m)Hrms(m) 0 50 100 150 0 20 40 60 80 100 120 140 160 1 0 1 01 01 02 0202 02 0x(m) =>shoreline m(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 88881 01 01 01 01 21 21 2 1 41 41 41 61 61 61 8 1 8 1 8 2 0 2 0 2 0 x(m) =>onshore S(o) 0 50 100 150 0 20 40 60 80 100 120 140 160 Figure 224. Contour plots of computed wave height (left), mean direction (middle) and directional bandwidth (right) for waves: 0.5rmsHm 4mTs o30m 0 S and o10 S at the offshore boundary. PAGE 94 94 0 20 20.20 30.30 30 40.4 0.40 50 50 50 50 50 5x(m) =>shoreliney(m)Hrms(m) 0 100 0 20 40 60 80 100 120 140 160 1 01 01 01 02 02 02 02 0x(m) =>shoreline m(o) 0 100 0 20 40 60 80 100 120 140 160 2 02 02 03 03 03 04 0 4 0 4 0 5 0 5 0 5 0 5 0 5 0x(m) =>onshore S(o) 0 100 0 20 40 60 80 100 120 140 160 Figure 225. Contour plots of computed wave height (left), mean direction (middle) and directional bandwidth (right) for waves: 0.5rmsHm 4mTs o30m 0 S and o20 S at the offshore boundary. 0 20 40 60 80 100 120 140 160 0.1 0.2 0.3 0.4 0.5 0.6 x(m)Hrms(m) Figure 226. Comparison of rmsH along transect s1 (black line s), s2 (red lines) for various directional bandwidths: o0 S (solid lines), o10 S (shortdashed lines) and o20 S (longdashed lines). PAGE 95 95 0 20 40 60 80 100 120 140 160 0 5 10 15 20 25 30 x(m)m(o) Figure 227. Comparison of m along transect s1 (black line s), s2 (red lines) for various directional bandwidths: o0 S (solid lines), o10 S (shortdashed lines) and o20 S (longdashed lines). Figure 228. Contour plots of computed m (left) and S (right) for waves: 0.5rmsHm 4mTs o30m 0.3Hz S and o0 S at the offshore boundary. PAGE 96 96 0 10 20 30 40 50 60 0 2 4 v(m/s)(a) 0 10 20 30 40 50 60 0 0.2 0.4 k(m1) k2 (b) 0 10 20 30 40 50 60 20 40 60 (degree)(c) 0 10 20 30 40 50 60 0.5 1 1.5 E/(E)0x ( m ) (d) Figure 229. Monochromatic waves obliquely incide nt on a longshore current over constant deep water 100hm. The offshore wave conditions are 00.61 Hm 4Ts and 0 022.5 modeling results; o analytical solutions. (a) longshore velocity distribution; (b) wave number along xdirection, red dash line indicates yk ; (c) wave direction; (d) normalized wave ener gy distributed along xdirection. 0 10 20 30 40 50 60 0 2 4 v(m/s)(a) 0 10 20 30 40 50 60 0 0.2 0.4 k(m1)(b) 0 10 20 30 40 50 60 20 40 60 (degree)(c) 0 10 20 30 40 50 60 0.5 1 1.5 E/(E)0x(m) (d) Figure 230. Unidirectional wave s with various frequency bandw idths obliquely incident on a longshore current over constant depth. The offshore wave conditions are the same as Figure 229 but with various frequency bandwidths: (solid curve), 0 S ; (dotted line), 0.1mS ; (dashed line), 0.2mS ; (dashdotted line) 0.3mS PAGE 97 97 0 10 20 30 40 50 60 0 0.5 1 v(m/s)(a) 0 10 20 30 40 50 60 0.2 0.25 0.3 k(m1)(b) 0 10 20 30 40 50 60 20 25 30 35 (degree)(c) 0 10 20 30 40 50 60 0.9 1 1.1 E/(E)0x(m) (d) Figure 231. Waves with various directional bandwidths oblique ly incident on the alongshore current over constant depth. The offshore wave conditions are the same as Figure 229 but with various di rectional bandwidths: 00.15.73 Srad dotted line; 011.46 S dashed line; 017.19 S dashdotted line. 0 10 20 30 40 50 60 0 2 4 v(m/s)(a) 0 10 20 30 40 50 60 0 2 4 k(m1)(b) 0 10 20 30 40 50 60 0 20 40 (degree)(c) 0 10 20 30 40 50 60 0 5 10 E/(E)0x ( m ) (d) Figure 232. Monochromatic waves obliquely in cident on the alongshore current over a planar beach with the slope of 1:20. For comparison purpose, nonbreaking waves are assumed. PAGE 98 98 0 10 20 30 40 50 60 1 0 1 u(m/s)(a) 0 10 20 30 40 50 60 0.2 0.3 0.4 k(m1)(b) 0 10 20 30 40 50 60 0.5 1 1.5 2 (c)E/(E)0x(m) Figure 233. Monochromatic waves with paralleli ng currents over constant water depth of 100m. The offshore wave conditions are 00.61 Hm 4Ts modeling results; o analytical solutions. (a) crossshore veloci ty distribution; (b) wave number along xdirection; (c) normalized wave energy. 0 10 20 30 40 50 60 1 0 1 u(m/s)(a) 0 10 20 30 40 50 60 0.2 0.3 0.4 k(m1)(b) 0 10 20 30 40 50 60 0.5 1 1.5 2 (c)E/(E)0x(m) Figure 234. Waves with various fr equency bandwidths with parall eling currents ove r deep water. The offshore wave conditions are the same as Figure 233 but with various frequency bandwidths: 0 S ; (dotted line), 0.1mS ; (dashed line), 0.2mS ; (dashdotted line) 0.3mS PAGE 99 99 0 10 20 30 40 50 60 1 0 1 u(m/s)(a) 0 10 20 30 40 50 60 0.2 0.3 0.4 k(m1)(b) 0 10 20 30 40 50 60 0.6 0.8 1 1.2 1.4 1.6 (c)E/(E)0x(m) Figure 235. Waves with various directional bandwidths with paralleling currents over deep water. The offshore wave conditions are th e same as Figure 233 but with various directional bandwidths: 00.15.73 Srad dotted line; 011.46 S dashed line; 017.19 S dashdotted line. 0 10 20 30 40 50 60 1 0 1 u(m/s)(a) 0 10 20 30 40 50 60 0 2 4 k(m1)(b) 0 10 20 30 40 50 60 0 5 (c)E/(E)0x(m) Figure 236. Monochromatic waves with paralleling currents over a plane beach with the slope of 1:20. The offshore wave conditions are 00.61 Hm 4Ts modeling results; o analytical solutions. For comparison pur pose, nonbreaking waves are assumed. PAGE 100 100 CHAPTER 3 STEADY TWODIMENSIONAL NEARSHORE CIRCULATION MODEL 3.1 Introduction Waveinduced nearshore currents, including alon gshore currents, circulation cells and rip currents, carry sediments and ar e directly responsible for beac h morphology evolution. Energetic nearshore currents, especially ri p currents, are dangerous for b each swimmers. According to the United States Lifesaving Association, rip currents alone cause more than 100 deaths each year in United States alone, which on average is more deat hs than hurricanes, tr opical storms, lightning and tornadoes combined (Lascody, 1998). In additi on, nearshore circulations make exchange of offshore and onshore water, affecting the crossshore mixing of heat, nutrient and chemical substances in nearshore regions. The importan ce of the nearshore circulation has inspired extensive research efforts by numerous scientis ts and coastal engineer s in past decades. Numerous field observations and laboratory studies have been conducted all over the world and great achievements in improving the understa nding of the mechanism, forcing, and stability of the waveinduced nearshore currents have been accomplished. Mathematical models, which do not have some limitations of the fiel d or laboratory measurements, can study the nearshore circulation systems in more detail and more systema tically. A variety of numerical models studying the nearshore circulation system have been developed since the pioneering work on introducing the concept of waveinduced radi ation stresses for small amplitude waves by LonguetHiggins and Stewart (1964). Among them, early works include Bowen (1969), LonguetHiggins (1970), Thornton (1970), and Noda (1974). More recently, more sophisticated modeling efforts include Van Dongere n et al. (1994), Srensen et al. (1998), Chen et al. (1999), Haas and Svendsen (2000), and Yu and Slinn (2003), among others. PAGE 101 101 Most of nearshore circulation models are ge nerally based on either the Nonlinear Shallow Water Equations (NSWE) (e.g. D ongeren and Svendsen, 1997) or the Boussinesqtype equations (e.g. Chen et. al, 1999), and the former is used more often. Most of these models are timedependent and were developed using the finite difference method, and some of them are very computer intensive and require very small time st eps to reach steadystate solutions. Aside from current instabilities which can happen under certain circumstances steady wave forcing usually leads to steady current motions in cluding longshore currents, rip currents and circulation cells, depending on the wave forcing itself and the specific bathymetry. In this case, a steady circulation model that is r obust and computationally much less time consuming than timedependent models will be very useful in a prac tical sense. Up to date, to our knowledge steady current models exclusively developed for nearshore circulation are rare. In this study, a steady circulation model ba sed on the twodimension depthaveraged NSWE is developed by the finite volume method (FVM). In com putational fluid dynamics (CFD) field, the FVM method has become more and mo re popular owing to its advantages over the finite difference method (FDM) and the finite element method (FEM). For example, FVM naturally satisfies many conservation laws, is co mputationally efficient and can accommodate in any types of grids. Unlike most of the existing models, in this model the pressurecorrection method is used to solve the NSWE, which is ve ry efficient and robust especially for steady problems. The model is developed using the bounda ry fitted structured nonorthogonal grids. Consequently, the model can simulate currents in relatively complex bathymetries. In this chapter, we focus on the development of the 2D steady nears hore circulation model. The numerical procedures, discretization of gove rning equations, treatment of various boundary conditions, and the linear algebraic solvers will be discussed in detail. The model validation and PAGE 102 102 convergence tests are conducted us ing simple computational test s. Applications of this circulation model incorporated with the wave model introduced in the previous chapter in simulating nearshore waveinduced currents are presented in chapter 4. 3.2 Formulations In modeling nearshore twodimensional horizon tal circulations, the 2D nonlinear shallow water equations (NSWE) are often used to desc ribe depth averaged mean flows. The NSWE describe incompressible flow when the ratio of the vertical scale to the horizontal scale is a small value, derived based on the depthintegrated (or depth averaged) NavierStokes equations. Hydrostatic pressure is usually assumed in modeling the NSWE since hydrostatic balance is accurate when the characteristic horizontal length s cale is much larger than vertical length scale. Hydrostatic balance states that gravity balances the pressure gradient in the vertical equation of motion, implying that ve rtical acceleration are negligible: P g z (31) where, the density is constant. This equation implies that the horizontal pressure gradient is independent of z This is the key simplification that unde rlies the 2D shallow water system. For steady problem, dropping the time derivative terms of the traditional shallow water equations, the governing equations of this steady circulation model are written as follows Continuity Equation: 0 huhv xy (32) Momentum Equations: bxx xhP huuhuv F xyx (33) PAGE 103 103 byy yhP huvhvv F xyy (34) where uv are the depthaverag ed velocities in x y directions, h the local water depth, the water density, Pgz is the pressure, the mean water surface displacement relative to the undisturbed free surface, g the gravity, b x and b y denote the bottom friction effect in x and y coordinates respectively. x and y represents the eddy viscos ity (or lateral mixing in oceanic processes) effect. 1 x y xx xS S F x y and 1yxyy ySS F x y are the wave forcing term in x and y direction, respectively. In which x xS x yS ,yyS are wave radiation stress components. It is worth mentioning that the rigidlid approximation has been used in the governing equations given above for simplicity. The rigidlid approximation neglects the effect of surface elevation on the continuity equation. Justificati ons of this approximation in simulating nearshore currents can be found in Allen et al. (1996), zka nHaller and Kirby (1999), among others. Bottom Friction Bottom friction plays a very important role in the waveinduced nearsh ore circulation. In the surf zone, bottom stress could be the major fo rce to balance the wave forcing under certain conditions. It is very difficult to develop a rigorous formula for bottom friction, especially when both waves and currents are present. Jonsson et al. (1974) have proposed a semiempirical formulation for the combined case 1 2bfww (35) PAGE 104 104 where, cwwuu is the vector sum of the current velocity vector and wave orbital velocity vector. The friction coefficient sinwcwffff and is the angle between the total velocity vector w and the wave orbital velocity wu. In this general formulation, information about w f is rarely available and is also difficulty to obtain, which make it extremely difficult to apply directly. Instead, the simple linear da mping formulations are often used in practice by assuming that the incident waves approach the s hore with a small angle a nd the current is weak compared to the wave orbital velocity: 0bfcU uu (36) where, 0U is the amplitude of wave horizontal orbital velocity at the sea bottom. According to the linear shallow water theory, 0U can be evaluated by (Thornton and Guza, 1983) 01 4rmsg UH h (37) Furthermore, the drag coefficient f c is also a major source of uncertainty. Constant f c with the typical value of O(0.01) from offshore to shore line has been used by many investigators (zkanHaller and Kirby, 1999; Slinn et al.2000; Yu and Slinn, 2003). Ru essink et al. (2001) suggested that f c could be parameterized with the Ma nningStrickler equation (Sleath, 1984) 1/30.015a fk c h (38) where, ak is the apparent bed roughness and is assu med to be constant and timeindependent. According to Ruessink et al. (2001), ak is typically in the range of 0.01~0.06m in simulating longshore currents. In the present circulation m odel, bottom friction uses weak current, small incident angle with variable drag coefficient f c given by Equation (38). PAGE 105 105Lateral Momentum Mixing There are three sources of lateral momentum mixing: (a) mixing owing to turbulence induced by wave breaking; (b) mixing due to bo ttominduced turbulence; (c) mixing due to the vertical nonuniformity of horiz ontal velocities. The first two sources are turbulence induced stress and can be modeled by utilizing th e eddy viscosity approach (Sancho, 1997) tS x and thu hu S x x (39) where, t is the turbulence eddy viscosity and a ccounts for both breaking turbulence and bottom turbulence. Putrevu and Svendsen (1999) proposed the following formulation to approximate the momentum mixing due to vertical vari ation in horizontal velocities dhu hu DD x xx (310) They stated that D is roughly proportional to the wa veinduced fluxes in the horizontal directions. Therefore, for waves that are obliquely incident at small angles, the volume fluxes due to waves will most occur in x direction. Under this assumption, we know the dominant dispersive mixing term is x xD term and only this term will remain as a first approximation, which gives 2d xxxxxhu hv DD x xyx d xxxhv D x x (311) Combining the turbulent momentum mixing and the dominant dispersive mixing gives 2xtx xtx xthu hvhu DD x xyxyy (312) PAGE 106 106 2ytx xtthv hvhu D x xyyxy (313) Svendsen and Putrevu (1994) showed that x xD can be much larger thant Therefore, it is convenient and reasonable to remain only the txxD terms, and then the terms are reduced to 2xhuhv x xyx yhv x x (314) The combined eddy viscosity is parameterized as follows 1/3 10/bcUhMh (315) where, 1c is a constant coefficient typically being O(0.01) and 0U is the wave bottom orbital velocity. M is the mixing coefficient and is a cons tant being 0.06~0.5 suggested by zkanHaller and Kirby (1999). bt o tgD represents the energy dissipa tion due to wave breaking, and totD is defined in chapter 1. 3.3 Numerical Procedures 3.3.1 Discretization Approach Three discretization approaches widely used in computational flui d dynamics (CFD) are: Finite Difference Method (FDM), Finite Volume Method (FVM) and Finite Element Method (FEM). FDM is the oldest method for numerical solution of partial differential equations, believed to have been introduced by Euler back to 18th century. FDM is also the easiest one to understand and to use for simple geometries. On e main disadvantage of FDM is that the conservation is not satisfied naturally unless special care is taken. In addition, the restriction to simple geometries is another signi ficant unfavorable feature of FDM. PAGE 107 107 The distinguishing feature of FEM is that th e equations are multiplied by a weight function before they are integrated over the entire domain. The most important advantage of FEM is the ability to deal with arbitrary geometries. The gr ids are easily generated and local refinement can be easily done too. The principle drawback is that the matrices of the algebraic equations are not as well structured as those for regular grids. FVM uses the integral form of the conser vation equations as its starting points. The solution domain is subdivided into a finite number of contiguous control volumes (CVs), and the integrated conservation equations are applied to each CV, as well as to the solution domain as a whole. If we sum equations for all CVs, the global conservation equations will be obtained since the surface integrals of inner CV faces cancel out. Thus, global conservation is built into the method, which is one of its principle advant ages. In addition, the FVM can accommodate any type of grids, so it is suitable for complex geometry. All terms using FVM have physical meaning, which makes it easy to understand a nd implement. FVM is also usually more computationally efficient than FEM. In this study, the FVM is used to disc retise the integral form of the governing equations. 3.3.2 Grid and Variable Arrangement Computational grids can be r oughly divided into two categorie s: structured (regular) grid and unstructured grid. Structured grids consist of families of grid lines with the property that members of a single family do not cross each othe r and cross each member of other families only once. Unlike the structured grids, unstructured gr ids do not have a restriction on the shape of the elements or control volumes nor number of neighbor elements or nodes. Unstructured grids are flexible and can fit an arbitr ary solution domain boundary. However, the matrix of the algebraic PAGE 108 108 equation system is no longer well structured. Conse quently, it is difficulty to find a good solver, and solvers are usually slower than those for structured grids. In order to be able to calculate flows in rela tively complex geometries at the same time to use structured grids, the boundary fitted structured nonortho gonal grids are used in this circulation model. Since the grid lines can follow the boundaries, the boundary conditions are more easily implemented. One disadvantage of such grids is that more terms introduced into the discretized equations thereby incr easing both the difficulty of coding and the cost of solving the equations. Moreover, for a grid with strong nonorthogonality, unphysical solutions may result, and accuracy and efficiency of the algorithm w ill be affected. Nevertheless, this can be circumvented by making the grid as orthogonal as possible in practice. In practice, for a simple geometry the orthogonal grids are favorable since it is easier to generate and the algebraic system associated is wellstructured too. In this model, the principles of discretizat ion and solution methods are developed using the general nonorthogonal grids and they are of c ourse valid for orthogonal grids, which are essentially a special case of nonorthogonal grids. Figure 31 shows a typical nonorthogonal 2D control volume (CV) and the notation used. The CV surface consists of four plane faces, denoted by lowercase letters corres ponding to their direction (,, and wnes) with respect to the central node ( P ). The neighboring CV centers are repres ented by uppercase lett ers corresponding to their directions too. This notation sy stem is used throughout this chapter. There two ways to arrange unknown dependent variables on a grid: staggered arrangement and collocated arrangement, shown as Figure 32. The staggered arrangement makes the stencil for the continuity equation very compact, requiri ng no interpolations in evaluation of mass fluxes. In addition, the pressuregradi ent terms in the momentum equa tions require no interpolation. PAGE 109 109 These two features make the staggered grid appro ach responsive to gridl evel fluctuations, thus uniquely capturing all the resolvable modes and preventing accumulation of spurious energy at the grid level. However, staggeredgrid scheme s become very awkward to generalize to complex geometries. It also has been poi nted out that the use of stagge red grid for complex geometries leads to either high memory re quirements, or inefficient solution methods (Zang et al, 1994). Unlike staggered arrangement, the collocated (nonstaggered) arrangement stores all variables at the same set of grid points and all variables shar e the same control volumes. Generally, using a collocated grid requires lower memory and less computational cost. Moreover, the collocated arrangement has significant advantages in comp licated computational domain, especially when the boundaries have slope discontin uity. The main drawback of collocated arrangement is the occurrence of spurious oscillat ion in pressure field due to the weak coupli ng between the velocities and the pressure. Fortunately, this pr oblem has been greatly alleviated since the improved pressurevelocity c oupling algorithms were develo ped in 1980s. The current circulation model is developed ba sed on collocated variable arrangement. 3.3.3 Solution Algorithm: SIMPLE Algorithm Solving the 2D depthaveraged NavierStokes equations is very complicated by the natures of the equations themselves: the conv ective terms in the momentum equations are nonlinear; dependent variables are intricatel y coupled, and there is no apparent dominant variable for each equation; there is no indepe ndent equation for the pressure, whose gradients contributes to momentum equati ons. Solution of steady equations is even more difficult since many methods such as fractional step methods and some explicit time advance schemes that work well for unsteady problems may not work for steady problems. PAGE 110 110 The pressurecorrection method, which belongs to the segregated algorithm, is an approach widely used for solving the Navierstokes equations for fluid flow problems. The SIMPLE (SemiImplicit Method for PressureLinked Equa tions) algorithm, proposed by Patankar and Spalding (1972), was the first such algorit hm widely used in the literature. The distinguishing characteristic of SIMPLEfamily algorithm is that a pressurecorrection term is introduced to the calculation procedure of each segregated solution step (outer iteration level) to improve the velocity field computed from the linearized momentum equation as the intermediate value. The pressurecorrection e quation is built up by subs tituting the corrected velocities into the discretized continuity equa tion. Once the pressure correction obtained, the velocities solved from linearized momentum equations can be updated and the modified velocities satisfy the mass conser vation for each control volume at each outer iteration level, which is critical for iteration convergence. Now we introduce the SIMPLE algorithm in some detail. At the iteration level m, the discertized momentum equa tions can be written as ,,ii iiuu mmmm P iPliluu lAuAuSSP (316) in which, the subscript P denotes the node at which the partial differential equation is approximated and index l represents neighbor nodes involve d in finite volume approximation, the source term im uS contains all of the terms other th an pressure term, that might depend on the mu and/or other vari ables at the level m. SP is the pressure term, PhP SPdV x Due to the nonlinearity and c oupling of the underlying differential equations, Equation (316) can only be solved using iterative methods. Th e iterations within one time step, in which the coefficient and source matrices are updated, are called the outer iterati on to distinguish them PAGE 111 111 from the inner iterations performed on linear algebraic systems with fixed matrix coefficients. On each outer iteration, the e quations to be solved are: **11 ,,ii iiuu mmmm PiPliluu lAuAuSSP (317) The superscript m stands for solution to be sought. Th e right hand side of the equation is evaluated using the variables at the preceding outer iteration. Rewriting Equation (317) gives, 1 1 ,1i i iiu mm m ul i l m l iP uu PP P SAu hP u AAx (318) where / x represent a discretized spatial derivative. Since the pressure used in Equation (318) was obtained from the previ ous outer iteration, the velocities computed from Equation (318) do not necessarily satisfy the di scretized continuity equation. The discretized continuity equati on may be written as following: 0m i ihu x (319) Now we define the pressure correction P as the difference between the corrected pressure field and the previous pressure field1 mP. Similarly velocity corrections u and vare introduced as small corrections to the intermediate veloci ties computed from th e linearized momentum equations: mmuuu mmvvv and 1 mmPPP (320) Substitution of Equation (320) into the momentum equations (318) yields the relation between the velocity correction and pressure correction as follows ,1i iiu lil l iP uu PPi P Au hP u AAx (321) PAGE 112 112 For the sake of convenience, the first term on the right hand side of the above equations is represented by, iPu : ,,1iiPiP u Pi P hP uu Ax (322) Substituting Equation (322) into the discretized continuity Equation (319), we obtain the following pressure correction equation: *0im ii u iPii i PP PhhP hu hu xAxx x (323) The velocity corrections iu are unknown at this point, so the common practice is to drop them and neglecting this term does not influence the final steadystate soluti on. The above equation is then reduced to a normal Poisson equati on and can be solved using regular methods. Omission of the last term is the main approxi mation of the SIMPLE algorithm and is probably the major reason why the resulting method does not converge very rapidly. Various modified SIMPLE methods have been developed, and they are called SIMPLEfamily methods. The detailed descriptions of these methods can be found in many CFD textbooks such as Patankar (1980), Fletcher (1991) and Ferz iger and Peric (2002) etc. 3.3.4 Discretization of Governing Equations In this section, the finite volume schemes us ed to discretise the governing equations based on the SIMPLE algorithm on a structured nonorthogona l collocated grid is described in some detail. Details of discretizati on of terms in governing equati ons are given in Appendix C. Once the approximations for all terms are co mpleted, rearranging all coefficients and source terms, we can obtain the following alge braic equation systems for momentum equations: ,,, ,,,ii iuu mmm PiPlilu lAuAuQlEWNS (324) PAGE 113 113 where, min(,0)ieE u e Ee PESh Am L (325) min(,0)iwW u w Ww PWSh Am L (326) min(,0)inN u n Nn PNSh Am L (327) min(,0)i s S u s Ss PSSh Am L (328) iiiiiuuuuu pEWNS PAAAAA h (329) 111 111 ,,, ,,, 111 111 ,,, ,,,max(,0)min(,0)+max(,0)min(,0) +max(,0)min(,0)+max(,0)min(,0) +immm mmm ueie eiP eiEwiw wiP wiW mmm mmm nin niP niNsis siP siS e eQmumumumumumu mumumumumumu h S 11 11 mm iiii w w ee ww mm ii ii nsi P ns nn ssuhu huhu S nn huhu huhu SSR S S P nn tt tt (330) Solving the linearized momentum equations using the iterative solver gives the intermediate values of velocities *u and *v which may not satisfy the continuity equation since old mass fluxes and pressure from previ ous outer iteration were used in discretizing the momentum equations. So: *****ewnsPmmmmm (331) where the mass fluxes are calculated using the intermediate velocities *u and *v. P m stands for the residual of the continuity equation. Substituting Equation (331) into Equation (323) leads to the pressurecorrection equation: *, ,,,pp PPllPP lApApmmlEWNS (332) in which, PAGE 114 114 221 1/pu Eeep eAhSA (333) 221 1/pu Wwwp wAhSA (334) 221 1/pv Nnnp nAhSA (335) 221 1/pv Sssp s AhSA (336) The term P m is analogous to P m with uv replacing **, uv see Equation (321) and (331). Since the velocity corrections at this point are not known prior to the solution of the pressurecorrection equation, the P m term is neglected if using SIMPLE algorithm as mentioned above. 3.3.5 Solution of Linear Algebraic Equations In the previous section, we presented a de tailed description of discretization of the governing equation and three sets of algebraic equations were obt ained, two for the momentum equations and one from the pressurecorrection equations. All the three systems have the same form and the algebraic equation for a contro l volume (see Figure 34) can be written as P PWWEESSNNPAAAAAQ (337) Over the whole computational domain, the matrix version of the algebraic system is given by: A Q (338) The coefficient matrix is a fivediagonal matrix, two other grid nodes in e ach direction involved. There are two families of solution techniques fo r linear algebraic equations: direct methods and iterative methods. Simple exam ples of direct methods are Cramers rule matrix inversion, Gauss elimination and LU decomposition. Direct me thods are designed to deal with full matrices and are usually very expensive. Iterative met hods are based on the repeated application of a PAGE 115 115 relatively simple algorithm leading to eventual convergence after certain number of repetitions. For sparse matrices, if each iteration is very ch eap and the number of iterations is small, the iterative methods may be efficient and cost much less than a direct method. Wellknown iterative solvers include Jacobi methods, GaussSeidel methods, Conjugate Gradient methods and others. The iterative method called strong ly implicit procedure (SIP) propos ed by Stone (1968) is used in this circulation model. The SIP method is essentially a version of incomplete lowerupper decomposition method. The ILU methods are based on the idea of usi ng an approximate LU factorization of the coefficient matrix A as the iteration matrix M : M LUAN (339) where L U are the lower and upper ma trix respectively and N is small compared to A Since the convergence is faster if the preconditioner matrix M is a good approximation of coefficient matrix A we wish to select L and U such that M is as good approximation to A as possible. Unfortunately, this method converges slowly. St one (1968) proposed that convergence can be improved by allowing N have nonzero elements on the diagonals corresponding to all nonzero diagonals of LU. And elements on nonzero diagonals of N are chosen so that 0 N Detailed description and the algorithm for fivediagonal matr ices can be found in the CFD textbook Ferziger and Peric (2002). The SIP method is designed for equation syst em arising from disc retization of partial differential equations. Incomplete L U factorization of matrix A need be calculated only once, prior to the first iteration. Us ually, convergence is obtained in a small number of iterations and it always has a smoothing error distribution, wh ich makes it a good solver by itself and a good smoother when the multigrid method is applied. PAGE 116 1163.4 Boundary Conditions Careful treatment of boundary conditions is important for a numerical model to simulate a physical flow correctly and accurately. The fi nite volume method requ ires that fluxes at boundaries either be known or expressed in te rms of known quantities and/or interior nodal values. Since no additional equations are given from boundary conditions, no additional unknowns should be introduced either. The boundary conditions often enc ountered in modeling nearshore circulation have been built in this model. Implementation of those boundary conditions is described in some detail. 3.4.1 NoSlip Wall Boundary Condition The noslip wall is one of most often used boundaries in computational fluid dynamics. Noslip condition means: ii walluu (340) In words, it is saying that at a wall boundary the ve locity of fluid is the sa me as the wall velocity. Meanwhile, the noslip wa ll boundary also implies 0 0sn wallwallvv sn (341) where, n is the local coordinate normal to the wall outwardly; s is the local coordinate parallel to the wall. s v nv are the velocity components in s and n direction, respectively. Since there are no fluxes across the wall, conve ctive fluxes of all qua ntities are zero. The diffusive fluxes should be treated ca refully. The diffusive flux in the u momentum equation at the wall can be written as d uw whu FS n (342) The velocity u can be written as PAGE 117 117 nnssuvxvx n x xn s x xs (343) Substituting Equation (343) into (342), comb ined with Equation (341), we arrive at: 'ss d P w uwssws wwww Pwhvhv FShvxSx nL (344) The quantity 's P hv can be evaluated according to Equati on (C13). The matrix coefficients for u at the wall boundary can be collected accord ing to Equation (344). Treatment of the v momentum equation is the same and will not be reproduced here. At the wall boundary the fluid velocity is alwa ys equal to the wall velocity, so in the process of solving the pr essurecorrection equation the velocity correction u is zero. Consequently, the matrix coefficient for th e pressurecorre ction equation should be set correspondingly. 3.4.2 Periodic Boundary Condition A periodic boundary condition is another often used boundary condition in nearshore circulation models. Periodic boundary conditions are usually used when the bathymetry and flow patter are approximately periodic with certain length scale, such as periodic rip channels. Periodic boundary condition is also used as an approximation sometimes when information at a boundary is not available. Figure 35 shows the sketch of a pair of periodic boundaries. The periodic boundary condition could be regarded that the opposite ends are tied together All quantities at one periodic boundary should be exactly the same as those at the other boundary. Veloci ties at boundaries are interpolated using values of the two neighboring CV centers: 1(,)(,)(1)(,)bn xuvfuvfuv (345) PAGE 118 118 where f is the linear interp olation factor. The mass fluxes a nd convective fluxes at boundaries can be computed using the interpolated veloc ities. The diffusive fluxe s at boundaries are also convenient to be evaluated usi ng the interpolated ve locities. Unlike the wall boundary condition, the pressure corrections at pe riodic boundaries usually are not zero. The pressure correction equations at CVs around the boundaries should be tr eated the same way as those at inner CVs, except the notation of neighboring CVs. 3.4.3 Symmetric Boundary Condition In many flows, there are one or more sy mmetry planes. At a symmetry plane, the following conditions apply: 0nv 0tv n (346) And thus the convective fluxes of all quantities are zero. Also the normal gradients of the velocity components parallel to symmetry plane and of all scalar quantities are zero too. The diffusive flux in the u momentum at a symmetry plane can be written as nn d P w uwnnwn wwww Pwhvhv FShvxSx nL (347) It is similar to the case of wall boundaries, except displacing the normal gradient of shear velocity tv with that of norm velocity nv Treatment of a symmetry boundary condition in the pressurec orrection is quite simple since the mass flux is zero there. 3.4.4 Inlet Boundary Condition At an inlet boundary, all quanti ties are prescribed. In this model, the three fundamental dependent variables ,, uvP need to be known at inlet bounda ries. All the convective fluxes and mass fluxes can be calculated, and the diffusi ve fluxes are approximated using known boundary PAGE 119 119 values of the variables. Onesided finite difference is used to compute the gradients. Since the velocities and the pressure are known at boundaries, both the velocity correction and the pressure correction should be set to zero at boundaries. 3.4.5 Outlet Boundary Condition At outgoing boundaries, the flow properties are usually not known. For this reason, the outgoing boundaries should be select ed as far from the downstream of the region of interest as possible. If the flow at the out going boundaries reaches a fully de veloped state and variation of all quantities ,, uvP in the flow direction is small, it may be reasonable to make the approximation of zero gradients along the boundary grid line. The zero gradient condition on a grid line can be implemented very easily. For example, the outgoing boun dary is located on the east side of domain. The zer o gradient condition gives: ePuu ePvv (348) The condition can be implemented implicitly, and the algebraic equati on for a CV next to boundary can be written as PEPWWSSNNP A AAAAQ (349) The zero gradient approximation is conveni ent to implement. However, it does not guarantee the mass conservation of neither single CVs nor over the computational domain as a whole. If higher accuracy is required, mass cons ervation rule can be used to approximate the velocities at the outgoing boundaries. For a CV at the boundary, mass flux out should be equal to mass flux coming into the CV, which gives: s nw n e eemmm v hS (350) PAGE 120 120 where, n ev is the velocity component normal to the east cell face. The velocity component parallel to east face can be gi ven still using zero gradient approximation. To ensure the global mass conservation, the normal velocity com puted by Equation (350) should be further corrected by: () ()in nn ee out M ass vv M ass (351) where ()in M ass ()out M ass are the total mass flux coming in and going out of the entire computational domain. ()out M ass can be computed approximately using n ev. The velocity at the outgoing boundaries is not corrected by means of pressure corrections. Hence in the discretized pressurecorrect ion equation, the matrix coefficient 0EA 3.4.6 Prescribed Pressure Boundary Condition In some situations, the exact details of th e flow distribution ar e unknown but the boundary values of pressure are prescribed. A common way of dealing with this type of boundary condition is to extrapolate the ve locities at the boundary. Since th e pressure is prescribed, the pressure correction is set to zero at boundary. Unlike the inlet boundary condition, the velocity correction is not zero for this boundary condition. In order to test the implementation of a prescribed pressure boundary condition, two simple flow problems, Poiseuille flow and Coue tte flow, are simulated us ing the current model. The Poiseuille flow is known as the incompressibl e, viscous flow in two paralleling stationary plates. If one plate is moving in its own pl ane, it becomes the Couette flow problem. The geometry and boundary conditions for the two flow problems are shown schematically in Figure 36. The governing equation for these two flow problems is PAGE 121 121 2 21 dudp dydx (352) Prescribed pressure values at the up and down boundaries are given. For Poiseuille flow, the lateral boundary conditions are 00 uy 0 uyh (353) And the analytical solution of the Poiseuille flow is 1 2 dp uyyhy dx (354) For Couette flow, the late ral boundary conditions are 00 uy uyhU (355) where U is the speeding of the moving plat e. The analytical solution is 2yUyU uy hy hh (356) in which, 22 hdp Udx is the nondimensional pressure gradient. When 0 pressure decreases in flow direction and u is positive everywhere; when 0 pressure increases in flow direction and the adverse pressure gradient works against the local viscous forces. Flow will reverse when 1 which indicates the adverse pressure gradient is larger than the viscous forces. In simulating the two flow cases, convection and bottom friction in the circulation model are turned off and the values of pressure at bo th ends are input. Figure 37 shows the comparison between the simulation results and analytical results. It can be seen that good agreement with analytical solutions has been achieved, and c onfidence in dealing the prescribed pressure boundary conditions of this model has been gained through the two simple tests. PAGE 122 1223.5 Numerical Tests 3.5.1 Pure Diffusion in a Rectangular Basin Generally, the possible best wa y to test a newly developed model is to simulate test examples whose analytical solu tions are available, allowing a thorough investigation of models performance, convergence and accuracy. An artificial test, pure diffusion in a rectangular basin, is designed and analytical solutions are available with some simplifications. Neglecting the advection terms and assuming constant water depth, the shallow water equations (32)(33) are reduced to: 0xyuv (357) 2 21 x F Pu u x hh (358) 2 21yF Pv v yhh (359) where x F and yF are the input forcing in xand ydirection, respectively. Using the stream function method, equations (3 57)(359) can be combined to just one equation. In detail, it starts with the relation between the stream function and velocities. u y v x (360) Taking y derivative with Equation (358) gives: 2 2 211 x F Pu u x yhyyhy (361) Taking x derivative with Equation (359) gives: 2 2 211yF Pv v x yhxxhx (362) Combining equations (361) and (362) and el iminating the pressure derivatives, it gives PAGE 123 123 2 21 0y xF F hhyx (363) where, vu x y is the vorticity. Using Equation (36 0), it is straightforward to obtain: 2 (364) Rewriting Equation (363) in st ream function, it becomes: 24 21 0y xF F hhyx (365) For a random forcing, finding th e analytical solutions of E quation (365) may be not easy at all. However, with special forcing and bounda ry conditions analytical solutions are not too difficult to obtain. For example, we assumed the following forcing: 442222 22s i n c o s x xyxyxyxy yAh Fkkkkkkkxky kh 0yF (366) where 2/ x xkL 2/yykL x L and yL are the size of the rectangular basin. We also assumed that all the four si des of the rectangular basin are symmetry boundaries, shown as Figure 38. The analytical solu tions for this designed forcin g distribution gi ven above and symmetry boundaries are not difficu lt to obtain, and they are The velocity field: sincosyxyuAkkxky cossin x xyvAkkxky (367) The stream function: sinsinyxyAkkxky (368) The vorticity: 22sinsin x yxywAkkkxky (369) PAGE 124 124 Figure 39 gives the plots of the input forcing and the analytical solutions on the 49 CV uniform grids. The designed forcing is symmetr ic about the horizontal centerline and antisymmetric about the vertical centerline. This spat ial distribution of forci ng leads to that at a boundary the normal velocity component is equa l to zero. It is the main reason why the symmetry boundary is used for all the four sides of the rectangular basi n. Both the streamline contour and the vorticity contour are point symmetric about the cen ter point. The velocity field composed of four circulation ce lls is point symmetric too. Calculations for the same forcing and bounda ry conditions using th e circulation model were conducted. Figure 310 shows comparison of velocities between numerical solutions and analytical solutions along horizon tal and vertical central lines. The computational velocities match the analytical solutions very well. Since the analytical solutions are available fo r this problem, computational errors can be easily calculated. A systemic and rigorous assessme nt of the computation errors can be made and then credibility of the model can be established. Second order sc hemes are used in discretizing the governing equations. Therefore, the model is expected to be 2nd order accurate. The root mean square values of difference between computat ional and analytical velocity field have been calculated. Computational errors using uniform grids closely follow the theoretical reduction path for second order schemes as grids become finer, as shown in Figure 311. In order to evaluate performance of the mode l on nonorthogonal girds, the same test was conducted on a nonorthogonal grid (Figure 312). Velocity comparison along a sample grid line is shown in Figure 312 (t he right panel). It can be seen that accurate results are achieved on the nonorthogonal grid too. Estimation of computation errors on nonorthogonal grids also shows the second order accuracy as well (Figure 313). PAGE 125 1253.5.2 LidDriven Cavity Flow The incompressible flow in a square cavit y with the top wall moving with a uniform velocity in its own plane is usually called the liddriven ca vity flow problem. This twodimensional flow problem has been used by many authors for testing and evaluating their models. Benchmark solutions for this problem are availa ble in the literature (Ghia et al., 1982). The geometry and boundary conditions are shown sc hematically in Figure 314. The moving lid creates a strong vortex around the center and sequen ce of weaker vortices in the corners. The Reynolds number is defined by: RelidUL (370) where is the viscosity. Calculations are carried out for Re100,1000,5000 For Re100 80 uniform grids are used. The 120 uniform grids are used for the case of Re1000 For the very high Reynolds number case Re5000 a 150 nonuniform Cartesian mesh is used. Figure 315 shows the computational results for Re100 The u velocity profile along the vertical line and v velocity profile along the horizontal line passing through geomet ric center of the cavity show good agreement with the results of Ghia et al. (1982). Computational resu lts and comparison of velocity profiles for Re1000 are shown in Figure 316. Fi gure 318 shows the simulation results of the high Reynolds number flow Re5000 on the nonuniform Cartesian grids (see Figure 317). As seen from Figure 315, 316 and 318, th e center of the primary vortex moves toward the geometric center of the cavity as Reynolds number increases. For Re1000 and Re5000 the secondary vortices at bottom corners become obvious, and a secondary vortex at the top left PAGE 126 126 corner can be apparently seen for the Re5000 case. Velocity profiles along the geometric centerlines for all Re numbers match very well with the results of Ghia (1982). Two important numerical techniques have been adapted in this ci rculation model: UnderRelaxation technique and Momentum Interpolation method. For details, see Appendix D. To test the techniques and to demonstrat e the effects of the modified momentum interpolation method, numerical simulations have been conducted for this cavity flow problem with Re100 The 30 uniform grids are used throughout these calculations. The underrelaxation factor for velocities varying from 0.3 to 0.9 were used in these calculations, and the optimum underrelaxation factor for pressure correction 1.1 P u is used to update the pressure correction at every outer iteration level. Convergence is de clared when the maximum RMS residual of the three governing equati ons is less than 106. Use of the modified momentum interpolation method is expect ed to give unique converged results for different underrelaxation factors. This is well illustra ted by Table 31. The velocities at monitored location are independent of underrel axation factors and unique results are obtained after convergence. Number of iterations required for convergence increases quickly as the velocity underrelaxation factor decreases from 0.9 to 0.3 (see Figure 319). 3.6 Waveinduced Alongshore Currents on a Planar Beach In chapter 2, the development of the presen t reduced spectral wa ve model accounting for currents has been described in great detail. In previous sections of this chapter we introduced the development of the steady 2D nearshore circulati on model, including the formulation, numerics and treatment of various boundary conditions. The model verification was followed using several simple numerical experiments. With the concep t of radiation stress, the wave model and the circulation model can be coupled easily in logical sense. Sp ecial care should be taken in PAGE 127 127 computing radiation stress when irregular waves are involved, as discussed in Chapter 1. The finite volume method (FVM) is used in both the wave model and the circulation model, making it easier to couple the two models. The main ch allenge for coupling the two models lies in the fact the wave model is in time domain while the circulation model is timeinvariant. This, however, does not prevent from coupling the tw o models. For steady wave fields, the wave model, in fact, can be regarded as using a transient method to solve steady problems. The wave field is computed with a pseudo time marching and the final solutions that do not vary with time can be obtained after a sufficient run time. The wave model itself of course can be used to predicted time varying wave fields. The detailed procedure of coupling the wave mo del and the circulation model can be put in the following way. The wave model is originally executed with zero ambient currents, giving an approximate wave field. The radi ation stresses and wave forcing can then be calculated and input to the circulation model as a priori. After the ci rculation model reached to a convergent solution, we rerun the wave model with the newly computed currents with the old wave field as the initial condition. The loop should be repeated until variat ions of both wave properties and current field from last loop are within the tolerance value. It should be noted that this is only a loose twoway coupling. In this section, the coupled wave/circulation model is used to simulate the waveinduced steady alongshore currents over a longshore uniform planar beach. Using this simple bathymetry it is easier to verify the new wave/circulation model, by avoiding the complexity that a more complicated bathymetry may cause. It is also co nvenient to examine effects of bottom friction and lateral mixing effects on currents and influe nces of wavecurrent feedback both on waves and on currents on a relatively simple bathymetry The formulations of the bottom friction and PAGE 128 128 lateral momentum mixing used in this circula tion model were introduced in Section 3.2 and the free parameters ak and M are very important but usually not known as a priori in nearshore circulation modeling. To inve stigate the effect of bottom friction on alongshore currents, simulations with varying bottom fr iction factors are conducted. In addition, a set of simulations with varying lateral mixing coefficient is carried out to study the effects of lateral mixing. The effects of wavecurrent interaction on the alongshore currents are also studied. The artificial planar beach with a constant slope of 1:20 is used here, and the offshore water depth is 3 m (Figure 320). For illustrati on purpose, a small computational domain of 40 m in width, 60 m in length and co arse uniform grids (60) are used. The noslip wall boundary condition is applied to the shor eline boundary, and late ral boundaries are periodic. The offshore wave conditions are0.61Hm, 4Ts and 022.5 Since whether the wave is monochromatic or irregular is i rrelevant to the investigation of effects of the bottom friction, lateral momentum mixing and waveinteraction, regular waves are used here. Effects of frequency/directional spreading on waves and currents will be studi ed separately in chapter 4. The typical acrossshore variati ons of energy dissipation due to breaking, wave bottom orbital velocity, eddy viscosity and bo ttom friction coefficient for 0.2 M and 0.03akm are shown in Figure 321. For a fixed M a decrease in ak increases the magnitude of longshore current v without significantly altering the shape of v velocity profile (Figur e 322). Simulations for varying values of M with the fixed 0.03akm show that lateral mixing smoothes the crossshore distribution of v without shifting the location of maxv (Figure 323). The momentum balance in the longshore direction resulting from a simulation for 0.25M and 0.03akm is PAGE 129 129 shown in Figure 324. The wave forcing is mo stly dissipated by the bottom friction and the mixing process is represented by the eddy viscosity term. As discussed in chapter 2, both the wave hei ght and wave direction can vary significantly as waves propagate through a strong shear current with a considerable angle of incidence. Figure 325 displays the computed wave height and wave direction with and wit hout the influences of wavecurrent interaction. The wa ve height decreases slightly with the presence of longshore currents in the region right before wave break ing starts because of currentinduced wave refraction, which can be seen more clearly in the plot of wave direc tion comparison (see Figure 325b). Wavecurrent interaction, in this case, po stpones the outset of wave breaking. Closer to the shoreline, wave breaking dominates and wave current interaction ha s less influence on wave height. Results also show that the effect of wc interaction on currents is not very significant but is not negligible neither (Figure 326). Generally, taking into a ccount the interaction shifts the velocity profile shoreward and decr eases the current peak slightly. PAGE 130 130 Table 31. Effect of the underrel axation factor on number of it eration required for convergence and effect of MMIM on the converg ed computational results. Monitored velocities Underrelaxation factor uv U(m/s) V(m/s) Iteration required 0.9 0.1140 0.0755 95 0.8 0.1140 0.0755 198 0.7 0.1140 0.0755 332 0.6 0.1140 0.0755 511 0.5 0.1140 0.0755 761 0.4 0.1140 0.0755 1137 0.3 0.1140 0.0755 1764 PAGE 131 131 P E S W Nnw ne sw se we s n nnnn Figure 31. A typical nonorthogonal 2D control volume and notations. u u v v P P,u,v Figure 32. Types of variab le arrangement: collocated (left) and staggered (right). Pw nE eP' E' tsw nw se ne e' s n Figure 33. A typical 2D control volume and auxiliary nodes used. P E S W N Figure 34. A control volume with neighboring CVs. PAGE 132 132periodic boundary periodic boundary1 nx 1 nx Figure 35. An example of the periodic boundaries. u(y) h xyh UP1P2P1P2 Figure 36. Sketch of Poiseuille flow (left) and Couette flow (right). Figure 37. Comparison of the velocity profil e between modeling results (solid line) and analytical solutions (circles) for Poiseu ille flow (left); Comparison of normalized velocity distribution with various values for Couette flow (right), modeling results (solid lines), analytical solutions (circles). PAGE 133 133 symmetry boundary symmetry boundary symmetry boundary symmetry boundaryLy Lx Figure 38. Sketch of the recta ngular basin with symmetry boundaries. Figure 39. Forcing distri bution (top left) and an alytical solutions for the pure diffusion problem. Stream function contours (top right), vorticity contours (bottom left) and the velocity field (bottom right). PAGE 134 134 Figure 310. Comparison of u (left) and v (right) at a sample sec tion between modeling results (solid lines) and the analytic al solutions (circles). Figure 311. Computational error as functions of grid number (uniform Cartesian grids used); the model is expected to be 2nd order accuracy since 2nd order discretization schemes were used. Figure 312. The nonorthogonal grid used (left) and velocity compar ison (right) along a monitored grid line between computed (solid line) and exact values (red points). PAGE 135 135 Figure 313. Computational error as functions of number of grid (nonorthogonal grid used). Ulid=1, V=0 U=0 V=0 U=0 V=0 U=0, V=0 L L Figure 314. Liddriven s quare cavity flow geometry and boundary conditions. PAGE 136 136 Figure 315. Velocity field (left top), vorticity contour (right to p) and comparison of velocities along the two geometri c centerlines for Re100 (80 uniform grids used). Figure 316 Velocity field (left top), vorticity co ntour (right top) and co mparison of velocities along the two geometri c centerlines for Re1000 (120 uniform grids used). PAGE 137 137 Figure 317. The 150 nonuniform Cartesian grids used for the Re5000 Figure 318. Velocity field (left top), vorticity contour (right to p) and comparison of velocities along the two geometri c centerlines for Re5000 (150 nonuniform grids used). PAGE 138 138 Figure 319. Number of iterati ons required to reach convergence using various velocity underrelaxation factors (Re100 the 30 uniform grids used). 0 20 40 60 0 20 40 3 2 1 0 x(m) y(m) h(m) Figure 320. A planar beach with a slope of 1:20. 0 10 20 30 40 50 60 0 50 100 b(N/m.s)(a) 0 10 20 30 40 50 60 0 0.1 0.2 (m2/s)(b) M=0.2 0 10 20 30 40 50 60 0.4 0.6 0.8 U0(m/s)(c) 0 10 20 30 40 50 60 0 0.005 0.01 x(m)(m/s)(d) ka=0.03m Figure 321. (a) energy dissipation due to wa ver breaking, (b) viscosity coefficient (c) wave horizontal orbital velocity near water bottom, (d) bottom friction coefficient M is the mixing coefficient and ak is the apparent bed roughness. PAGE 139 139 0 10 20 30 40 50 60 0.5 0 0.5 1 1.5 2 x(m)v(m/s) ka=0.015 ka=0.03 ka=0.06 Figure 322. Effect of the bottom friction factor ak on the alongshore curr ent distributed along the crossshore distance. 0.2 M is used for the expe rimental runs. 0 10 20 30 40 50 60 0.5 0 0.5 1 1.5 2 2.5 x(m)v(m/s) M=0.0 M=0.1 M=0.2 M=0.4 M=0.6 M=1.0 Figure 323. Effect of the lateral mixing coefficient M on the alongshore current. 0.03akm is used for the experimental runs. 0 10 20 30 40 50 60 6 4 2 0 2 4 6 8 x 103 x(m)Momentum Balance(m/s2) Figure 324. Alongshore momentum balance for 0.25 M and 0.03akm Bottom friction term v (solid line), wave forcing term (dashed line), lateral mixing term (dotted line), and residual (dashdotted line). PAGE 140 140 0 10 20 30 40 50 60 0 0.2 0.4 0.6 0.8 1 1.2 1.4 H/H0x(m) (a) 0 10 20 30 40 50 60 0 5 10 15 20 25 x(m)(o) (b) Figure 325. Comparison of wave properties with (s olid lines) and without (dashed lines) effects of wavecurrent interaction. (a) normalized wave height; (b) wave angle. 0.25M and 0.03akm are used. 0 10 20 30 40 50 60 0.5 0 0.5 1 1.5 2 x(m)v(m/s) w/o wavecurrent interaction with wavecurrent interaction Figure 326. Comparison of alongshore current profiles between neglecting wavecurrent interaction case (dashed line) and incl uding interaction case (solid line). 0.25M and 0.03akm are used. PAGE 141 141 CHAPTER 4 APPLICATIONS OF THE W AVE/CIRCULATION MODEL In previous two chapters, we introduced the reduced spectral wave model and the steady circulation model, and the coupling of the tw o models was briefly discussed, followed by application of the coupled wave/circulation mo del to waveinduced long shore currents over a planar beach. In this chapter, we focus on furthe r applications of the wa ve/circulation model to nearshore currents, including longshore curren ts over barred beaches and rip currents. The DUCK94 field experiment is simulated and compar ed to measured data. Another typical current pattern in the nearshore area, rip current, is also studi ed using the coupled models. In addition, the effects of frequency/directional spreading on the wave field and wa vedriven currents are investigated. Oscillations in velocity field with long peri ods have been observed in field measurements, laboratory experiments and numerical studies for longshore currents (OltmanShay et al., 1989; Allen et al., 1996; zkanHaller an d Kirby, 1999; Slinn et al., 2000), rip currents (Shepard and Inman, 1950; Sonu, 1973; Haller and Dalrymple, 2001; Yu and Slinn, 2003; Kennedy and Zhang, 2008), and the entire circulation cells (MacMahan et al., 2004). An understanding and predictive capability of the unsteady flow motions is very important and necessary for engineering practices. In order to study the instabilities of longshore currents and rip cu rrents, an unsteady circulation model extended from the steady version is deve loped and coupled with the wave model in the time domain, detailed descripti on given in Section 4.2.3. 4.1 Waveinduced Alongshore Currents over Barred Beach 4.1.1 Introduction The important role of alongshor e currents in the surf zone in coastal processes has long been recognized. Measurement and prediction of alongshore currents generated by obliquely PAGE 142 142 incident waves breaking near the shoreline have gained tremendous efforts by numerous scientists and coastal engineers. Models of alongshore currents are often based on the onedimensional alongshore momentum flux balance between wave, wind and tidal forcing, bottom stresses and lateral mixing (LonguetHiggins, 1970; Thornton and Guza, 1986). Breaking waves are the most important forcing mechanism in th e surf zone, although wind (Feddersen et al., 1998) and tidal (Houwman and Hoekstra, 1998) forcing can contribute significantly. Bottom friction also plays a very importa nt role in the waveinduced near shore circulation. In the surf zone, bottom stress could be the major force to balance the wave forcin g under certain conditions. Unfortunately, no rigorous formula accurately desc ribing the complicated near bottom processes is available. In practice, parameterized formul ations are usually used in modeling nearshore circulation. Lateral momentum mixing from wave breaking tur bulence is another factor in alongshore currents and more compli cated nearshore circulations. Other than using a constant eddy viscosity in and out the surf zone, the vari able parameterized expressions of lateral mixing have been proposed (e.g. zkanHaller and Kir by, 1999). The widely used formulations for bottom friction and lateral mixing, and their effects on the alongshore currents over a planar beach have been discussed in chapter 3. Barred beaches are another type of topography th at is often seen in the nature. For waveinduced longshore currents over barred beach es, laboratory experiments suggest that the maximum alongshore velocity is located on the bar (Reniers and Battjes, 1997). However, early numerical studies (e.g., Church and Thornton, 1993) predicted a strong, narrow jet on seaward side of the bar and nearzero flow in the bar tr ough, which differs from the laboratory results and field data. This failure may result from the alongshore variations in bathymetry and/or the neglect of wave surface ro llers in wave forcing. Variations of bathymetry in longshore direction PAGE 143 143 may alter the longshore current fiel d significantly, but which is not a subject of this study. The surface rollers carry the energy diss ipated from wave breaking and release it later as the rollers become smaller, thus shifting the maximum velocity shoreward and increasing the bar trough velocity as well. The roller model proposed by Stive and DeVriend (1994) is discussed in Section 4.1.2, and improvement in predictions of the measured alongshore currents at Duck, North Carolina on Oct, 14, 1994 by in cluding rollers in the wave forc ing, is also discussed in Section 4.1.2. Monochromatic wave drivers were traditi onally used in simulating waveinduced longshore currents. However, the monochromatic wave approximation is very crude and may overestimate or underestimate the radiation stre ss components significantly, which was discussed and quantitatively analyzed in Chapter 2. Due to the difference in wave radiation stresses, the resulting longshore currents are supposed to be different between fr om monochromatic waves and from irregular waves. However, the quantita tive change in longshore currents due to the presence of wave directional and frequency spreading is not well known. Using the current wave/circulation model, the effects of the wave frequency and directional spreading on the wavedriven alongshore currents are investigated, and the difference in al ongshore currents of DUCK94 owing to the presence of frequency/dir ectional bandwidth is demonstrated in Section 4.1.2. 4.1.2 Comparison to DUCK94 Field Experiment Data Models neglecting the wave rollers predicted current profiles with a maximum velocity seaward of the bar or even in the trough (e.g., Church and T hornton, 1993). Many researchers have proved that by including wave rollers can give better predic tion (e.g. Ruessink et al., 2001). The wave surface roller concept is introduced before presenting the simulation results. To PAGE 144 144 authors knowledge, narrowbanded waves have been usually assumed in current roller models, which means that the effects of frequency/directional spreading were not being accounted for. A roller model including the eff ects of wave frequency/directio nal spreading can give a more accurate representation. In the following, a brief liter ature review in wave roller models is given, followed by introducing the spectral version of the roller model. Surface Wave Roller As waves break, a volume of turbulent water is generated, which rushes down the front face of the waves. This body of aer ated water detached from the wave form is called the surface wave roller. The surface rollers carry the energy imparted from wave breaking shoreward. As a result, the wave momentum is imparted into th e water column and transferred shoreward, and consequently the surface rollers have an im portant impact on surf zone dynamics, including wave setup/setdown and waveinduced curren ts (Svendsen, 1984; Roelvink and Stive, 1989; Fredsoe and Diegaard, 1992; Smith et al., 1994). Svendsen proposed that the turbulent water body is pushed by the wave front with horizontal velocity equal to the wave phase speed, so the roller energy travels with wave phase speed. It also has been pointed out that as soon as the surface rollers are generated by wa ve breaking they are dissipated locally by the shear layer between the rollers and the wave motions. The onedimensional roller energy balance eq uation has been given by Stive and De Vriend (1994), and Reniers and Battjes (1997) independently. When extended to a twodimensional timedependent formulation, the energy balance for rolle r can be written as rr r b rEEDD t c+u (41) where rE is the roller energy density, derived from the kinetic energy density of the roller volume with unit crest le ngth (Svendsen, 1984a,b). c and u are the wave phase speed vector and PAGE 145 145 ambient current vector, respectively. The first term on the right hand side is the roller sink term, given by (Duncan, 1981; Dieggard, 1993): sinr rEg D c (42) where is the wavefront slope a nd usually is assumed to be 010 (Lippmann et al. 1996). brD stands for the wave breaking di ssipation and is the source term for the roller energy balance equation. Usually, an approximation of Equation (41), assuming that a wave spectrum is narrowbanded in frequency and direction (Lippmann et al., 1996; Ruessink et al., 2001), is used in practice, which is a gross approximation for br oad wave spectra. For a welldefined wave spectrum, for instance the Gaussianshaped wave sp ectrum, in fact, the formulation for the roller balance equation can be rigorously derived using the same strategy as to the wave action balance equation. The detailed algebras of derivation will not be reproduced here. In the surface roller energy balance equation, wave phase velocity c is weakly relevant to frequency since surface rollers only exist after wave breaking occurs, therefore, the equation only considering the directional spreading as following is used: 0,10,0 0,10,0//rrr r b r RIEcEEuEcEEvEDD tx y (43) 0,1RE and 0,1IE are the real and imagin ary part of wave moment0,1E as defined previously. 0,10,0/REE and 0,10,0/IEE contains not only information on th e mean wave direction but also the wave directional spreading. The Equation (43) is solved after solving the wave moments with the same numerical schemes. rE is set to be zero at the offshore boundary. PAGE 146 146 The formulation of roller radiation stress is similar to the wave energy stress. Since gcc in shallow water, the roller radiation stress is given by 22ijij ijrrkk SE k (44) The wave roller energy is coupled to the circulation model through the wave forcing in the same manner as the wave energy: ,,1 x yryyr r ySS x y (45) ModelData Comparison Alongshore currents driven by wave forcing over the barred beach are simulated and comparison to DUCK94 data is made. The bathymetry and experimental data used here were downloaded from the DUCK94 Experiment Data Se rver. The Duck94 data were collected from August to October 1994, during the near shore field experiments at the U.S. Army Corps of Engineers Field Research Facility (FRF) located in Duck, North Carolina. Directional wave statistics ( s igH peak period, and mean wave direction) were measured at 8m water depth. Pressure and velocity observa tions were obtained at more than 15 crossshore locations extending from near the shor eline to 8m depth. The bathymetry used in this study (as shown in Figure 41) was the bathymetry of daily minigrid 3D survey data at around 10:00 am, Oct ober 14, 1994. All the wave and current data are 3 houraveraged for a storm event. Th e 8m depth wave statistics are: m Hsig07.2 8.866pTs and 020 The pressure gage nearest to the offshore end of the measured bathymetry is the gage P19, located 347.4 m from shoreline. In order to take advantage of the available wave conditions at gage P19, the o ffshore boundary of the computational domain is PAGE 147 147 chosen to be at pressure sensor gage P19. The s igH at P19 is of 1.373 m, and 6.4pTs The mean wave direction at P19 is calculated using the Snells law. Alongshore uniformity in morphology, waves and currents is assumed through this numerical experiment. The infl uence of wavefront slope which determines how fast the roller energy dissipates (large for fast roller energy dissipa tion), is examined using the bathymetry and wave conditions of Duck94, Oct 14. The bed roughness ak is fixed to be 0.07 m and the mixing coefficient 0.6M for all numerical simulations. With no roller, the current jets are located on the seaward side of the sandbar and near the shoreline, with nearzero current in the trough area (50100 m). The rollers stored the en ergy from water breaking and released later, which results in shift of the maximum along cu rrent onshore and broadens the current jets by increasing the currents in the trough. With decreasing the maxv is shifted further onshore and velocity in the trough incr eases continually. For up to 0.2, the roller e ffect is not significant any more (Figure 42), consistent with the concept that for large dissipation rate (associated with large ) the roller energy is di ssipated quickly. The roller energy has a significant impact on the total radiation stress (sum of wave radiation stress and roller radiat ion stress) distributed along the crossshore distance from the shoreline, peak of ro ller energy located on where waves brea k strongly, i.e. where the trough of wave radiation stress located (Figure 43a). C onsequently, the crosss hore wave forcing peak shifts shoreward and the forcing distribu tion become broader. The magnitude of maxy is decreased considerably (Figure 43b). The wave model predicted the observed rmsH very well at all sensors. The rmsH errors for individual sensors vary between 0.01 m and 0.12 m, with an aver age of 0.04 m for all sensors. PAGE 148 148 The larger errors occur near the bar crest ar ea (Figure 44a). Figure 44b shows comparison of the predicted alongshore current (wit h and without roller model) with measured velocity profile. The roller model greatly improved the performa nce of the current model in simulating the alongshore current distribution by shif ting the current jets onshore and lifting th e trough velocity. Modeled alongshore velocity with 0.05 0.07akm and 0.6 M agree most well with measured data for this run. The 0.05 is consistent with values used by Walstra et al. (1996) and Ruessink et al. (2001), while, both0.6M and 0.07akm are equal to and slightly over the upper range suggested by zkanHaller and Kirby (1999). Terms in the ymomentum equation are plotted in Figure 45. Obviously, the wave forcing is balanced mostly by the bottom friction in this case. 4.1.3 Effects of Frequency/Directional Spreading For simplicity, the wavedriven nearshore circ ulation models (e.g., Church and Thornton, 1993; Ruessink et al., 2001) ofte n use the wave statistics (rmsH peak frequency and mean direction) to compute the radia tion stress and then the wave forcing. For waves with finite frequency/directional bandwidths, this approximation can result in significant errors (Battjes, 1972). Feddersen (2004) compared the exact radia tion stresses and showed that the component x yS could be overestimated up to 60% using the narrowbanded approximation, which is roughly consistent with the combined Duck94 and SandyDuck observations. As discussed by Battjes (1972), errors in radiation stresses will result in modeled current errors or, when fitting a model to observations, in biased model parameters such as bottom friction coefficient, eddy viscosity etc. The current wave model solves not only th e wave statistics but al so the frequency and directional bandwidths, and then the radiation stresses including the effects of frequency and directional spreadi ng can be obtained. PAGE 149 149 To demonstrate effects of frequency spread ing upon the modeled wave height, radiation stress, and finally the longshore velocity, the fi eld data of DUCK94, 10am, October 14 is resimulated assuming that the frequency spectrum is a Gaussian distributions with different values of frequency bandwidth S (/0,0.2,0.4mS are examined). As expected, no appreciable differences in predicted rmsH (Figure 46(a)) are observed between different frequency bandwidths because the wave breaki ng is so strong that other factor s are not very important. As discussed above, the radiation stresses are weakly dependent on frequency spreading in intermediate water, less than 10% variation even for a fairly large frequency bandwidth. In shallow water, the changes in radiation stress x yS with frequency spreadi ng are even less (Figure 46(b)). The small variation in x yS can only result in slight ch anges in alongshore velocity profiles (Figure 46(c)), the alongshore velocities decreasing s lightly with the increase of S Free parameters0.05 0.04akm and 0.5M were used for these experimental runs. Figure 47 demonstrates effect s of directional spreading upon the modeled wave height and alongshore velocity profile on DUCK94 bathymet ry. Due to the strong wave breaking, the rmsH is not obviously affected by dire ctional spreading (panel (a)). x yS changes considerably for 023 S (panel (b)) compared to the no directi onal spreading case (about 40% larger before waves break strongly), which is consistent with the analysis theoretical conclusion discussed above. The predicted maximum maxv for the no spreading case is 43% larger than with 023 S 4.2 Shear Stability of Alongshore Currents 4.2.1 Introduction In 1986, observations during the SUPERDUCK field experiment (OltmanShay, et al., 1989) revealed a surprising behavior of longs hore currents. The longs hore currents began to PAGE 150 150 oscillate with very low frequencies that typically fall into the lower end of the infragravity band 0.0010.01 f Hz However, the alongshore wavelengths of this motion are more than an order of magnitude less than the wavelengths predicted by edge wa ve or gravity wave theories, which suggested that these wave motions are not surface gravity waves (OltmanShay, et al., 1989). Bowen and Holman (1989), in a companion paper, interpreted the energetic, propagating perturbations in alongshore currents as shear instability of the mean currents us ing linear analysis, and their explanation has lead to these wavelike motions being referred as shear waves. This wave motion, occurring in the horizontal plane and moving with a speed less than the mean alongshore current, causes the alongshore current to move back and forth across the surf zone. To test the hypothesis that the newly observed wavelike motion is dynamically derived from the shear instability of mean longshore cu rrents, Bowen and Holman (1989) performed an analytical study of a simple, threepiece longshore current profile over a constant water depth. They showed that a shear inst ability leads to growing, nearly nondispersive shear waves propagating in the direction of th e alongshore current with a phase speed between and of the maximum value of the mean longshore current ()Vx, which was qualitatively in agreement with observations (OltmanShay, et al., 1989). The re storing mechanism of these shear waves is potential vorticity consisting of a relative vorticity and a background vorticity. The background vorticity is defined by the sh ear of the mean alongshore currents. Dodd, OltmanShay and Thornton (1992) extended the linear instability theory to include a linearized bottom friction and applied to more realistic alongshore currents and beach profile. They obtained good agreement of wavelengths and wave speeds from observations and from theore tical predictions based on the most unstable linear mode. PAGE 151 151 As shear waves reach to finiteamplitude, the assumptions of linearity would fail. In order to investigate the finiteamplitude behavior of shear waves, Feddersen (1998) carried out an analytical study on the weakly nonl inear shear waves using the st andard perturbation method. Using the same basestate longshore currents and ba thymetry as in the nu merical study of Allen et al. (1996), Feddersen obtained reasonably consistent results in the weakly nonlinear regime. However, this analytical model is mathema tically complex and becomes invalid when the instabilities are strongly nonlinear. Since numerical models based on the nonlinear shallow water equations can simulate shear waves from linear, weakly nonlinear to strongly nonlinear, numerical simulations of nonlinear instability of longshore currents ha ve received a great number of studies in literature. Falqus et.al (1994) performed numerical experiments of nonlinear shear wave s with the rigidlid assumption for a planar beach, including the bottom friction effects and tu rbulence mixing. They obtained some preliminary results on the finiteamplitude behavior of these shear instabilities. Using the rigidlid nonlinear shallow water equations on a planar beach with an idealized base mean longshore velocity profile and linearized bottom friction, Allen et al. (1996) conducted a set of numerical experiments and studied the effects of various bottom friction coefficients and varying longshore width of the modeling domain. They found that for flows from low to high nonlinearity (low nonlinearity defined by high fr iction coefficient with weakly unstable base longshore current), a wide range of behavior is observed. For weak nonlinearity, equilibrated steady finiteamplitude disturban ces are observed. As the nonlinear ity increases, the disturbances display modulated amplitudes and eventually irre gular eddies. They also found that for modeling domain substantially larger than the wavelengt h of the most unstable linear mode, unstable waves subsequently evolve into largerscale, nonlinear propaga ting wavelike disturbances. PAGE 152 152 Allen et al. (1996) also pointe d out that the measured mean current profiles in the ocean correspond to the final mean current in the presen ce of the fluctuations rather than fluctuationfree initial state. Slinn et al. (1998) extended the study of Allen (1996) over barred beach topography with a more realistic forced current computed usi ng the Thornton and Guza ( 1986) alongshore current model. They found that instabilities cause substa ntial lateral mixing of momentum in the surf zone and alter the initial current profile significantly. Soon afte r, zkanHaller and Kirby (1999) simulated the nonlinear instabilities of alongshore currents of the SUPERDUCK field experiment, including the bottom friction and lateral mixing effects. They obtained satisfactory results in the propagation velocities of the inst abilities as well as th e maximum mean alongshore currents including the shear instabilities. They ca rried out simulations fo r varying friction and mixing coefficients and pointed out that a smalle r friction coefficient leads to a stronger mean current and faster and more energe tic vortex structures and less en ergetic flow structures result from an increase in mixing coefficient. They al so found that shear instab ilities induce significant horizontal momentum mixing, as a result, the re sulting current profile in cluding effects of shear instabilities displays a smaller peak valu e and increased veloci ties in the trough. Numerical simulations of the finiteamplitude development of shear instabilities have gained great achievements and provided a better understanding of shear waves. The effects of friction and mixing on shear instab ilities have gained intensive studies by many investigators. Most of the studies, however, used the monochroma tic wave drivers. That how the instabilities of the mean alongshore currents driven by waves w ith frequency/directional spreading differ, has not been carefully studied to date. Since the wa ve forcing tends to be overestimated by using PAGE 153 153 monochromatic wave models (as discussed in previous section), bot h the resulting steady alongshore currents and shear inst abilities will alter consider ably for irregular waves. In this section, one of our goa ls is to investigate the effects of directional spreading of incident waves on the resulting sh ear instabilities of alongshore currents. Considering the fact that wave forcing simulated by the present re duced wave model for a narrow spectrum does not differ significantly from that of monochromatic waves, the effects of frequency spreading on shear instabilities are not very not iceable and will not be presented here. 4.2.2 Linear Stability Analysis Following Bowen and Holman (1989), the velocity field is assumed to consist of a steady longshore current () Vx and perturbations in alongshore curr ents. Therefore, the total velocity vector is (,,)(,,), ()(,,) x ytuxytVxvxyt u (46) where, (,,)vxytis the longshore pert urbation velocity and is assumed that vV. Substituting Equation (46) into the nonlinear shallow water equations and retaining only the linear terms, including bottom friction and lateral mixi ng, following equations can be obtained 2tyxuVuguhu hh (47) 2tyxyvVvVugvhv hh (48) where, is the free surface elevation, is the bottom friction coe fficient given by Equation (36). The above perturbation equa tions were derived under the ri gidlid approximation and have neglected the lateral mixing term, which have been proved reasonable for lin ear stability analysis (Allen et al., 1996). Under the assumption of nondivergence, a stream function (,,) x yt can be introduced, such that PAGE 154 154 /yuh /xvh (49) which allows the momentum equations to be combined and yields 2 2// /yy x xxyy x x yx xx x xxVhh tyhh Vh hhh (410) assuming a solution of the following form: ()(,,)Reikytxytxe (411) where, the simple alongshore harmonic is assumed. k is the longshore wavenumber and is a real number. But x and may be complex. Substituting Equation (411) into (410) gives the linear stability equation 2 24 2///// 2 / 0xxxx x x xx xxxx xx xxxxxx xVikhchhkhVhihk iih kkhk kk (412) Generally, the equation above need to be solv ed numerically. The fin ite difference method is utilized to solve the linear st ability problem and the following secondorder center difference scheme is used for spatial discretization: 11/2xii x (413) The following boundary conditions are impos ed at both ends of the domain: 0xx 0, x L (414) Discretizing Equation (412) a nd rearranging, we arrive at an algebraic eigenvalue problem c AB (415) where, only values of at 2,3,...,1 in are required to be solved. A and B are matrices. The above equation may be solved much more efficiently after the following operation: PAGE 155 1551c BA (416) A simple linear stability model programmed in Matlab is developed and applied to two cases. In order to verify the capability of the present m odel, the case of Bowen and Holman (1989) is reproduced as the first test. In the second test, a predicted steady alongshore current profile on a barred beach is studied including the bottom friction term with acrossshore variable This barred bathymetry is also used to study the nonlinear evolution of shear instability, which is presented in following Section 4.2.3. Case 1: Bowen and Holman (1989) A threepiece longshore current profile used by Bowen and Holman (1989) is shown in Figure 48. In this case, no bottom friction a nd no lateral mixing are assumed in order to compare with their an alytical results. Results of the linear stability calcu lations in terms of growth rate im for the most unstable mode as a function of alongshore wavenumber k are shown in Figure 49. The propagation velocity associated with the most unstable mode rc as a function of k is also shown. The alongshore wavenumber k of the most unstable linear mode is approximately 0.026m1 (and wavelength 02/242 km ), with a growth rate of about 0.00335 s1. The corresponding propagation velocity 0.335/rcms which gives a period of 02/()721rkcs These calculation results are in a very good agreem ent with Bowen and Holmans results. Figure 410 shows the spatial form of x y of the fastest growing waves. The domain in alongshore direction is chosen to be 484 m, which is two time s of the wavelength of the most unstable mode. As discussed by Bowen and Holm an (1989), progressive waves in the alongshore direction are observed. The pattern of total veloci ty, in which the fluctuat ion velocities are scaled PAGE 156 156 such that its peak magnitude equals the peak of the mean alongshore current for illustration purpose, is shown in panel (c). It is obvious that the alongshore currents meander along its peak position and secondary circ ulation cells can also be seen. Case 2: Longshore Currents on a Barred Beach In this particular test case and in next Section 4.2.3, a barred beach bathymetry approximately fitting topography measured at Duck, North Carolina, on October 11, 1990, given by Lippmann et al. (1999) is used 2 1111 12 1111tanh tanhexp5c cxx bxAbxbx hxA A AAAx (417) where 1tan(0.075) b 2tan(0.0064) b 112/ bb 80c x m is location of the sandbar. 12.97 A and 21.5 A This bathymetry is plotted in Figure 411. Monochromatic waves of 1.373 Hm 6.4 Ts and with the incident angle of 15 at offshore boundary are modeled ove r this barred bathymetry. The mean alongshore currents driven by wave forcing are simulated for vary ing bottom friction coefficients. Large bottom friction results in relatively weaker currents w ithout significantly chan ging the shape of the current profile (Figure 412(a)). The acrosss hore variations in botto m friction coefficients for 0.01,0.02,0.03akm are calculated and used during simula tions of the mean currents (Figure 412(b)), and they are also utilized in linear stability calculations. The resulting wavenumber of the most unstable linear mode for different bottom friction coefficients are about the same, 1 00.032 km which leads to a wavelength of 02/200 km (Figure 412(c)). However, the growth rate increases greatly as ak decreases. A growth time scale of 4.4 minutes is resulted for 0.03akm The corresponding propagation velocity rc increases considerably by decreasing ak PAGE 157 157 from 0.03 m to 0.01 m. As a result, the period 02/()rkc of shear waves is smaller for higher bottom friction (Figure 412(d)). For 0.03akm the period is about 4 minutes. 4.2.3 Numerical Simulation of Nonlinear Shear Instability 4.2.3.1 Unsteady Circulation Model The timeindependent 2D shallow water equati ons with the rigidlid approximation were used in the steady circulation model described in chapter 3. To simulate the time domain solution of nearshore currents, an unsteady version ex tended from the steady circulation model is developed based on the timedependent 2D shal low water equations. In the steady circulation model, the rigidlid approximation was used fo r simplicity. While, for unsteady flow problems treatment of the free surface is no t as difficult as for steady flow s using the pressurecorrection method. Therefore, the rigidlid approximation is not necessary in th e unsteady version of circulation model, and the effect of the free surface is included. The timedependent depthaveraged NSWE averaged over incident wave tim escales, including forcing and dissipation, can be written as 0 HuHv txy (418) x xxx txyHuHuuHuvgHF (419) yyyy txyHvHuvHvvgHF (420) where is the waveaveraged water surface elevation (setup), Hh is the total water depth. x y are bottom friction terms, given by Equation (36), x y are lateral mixing terms, given by Equation (314). In chapter 3, the SIMPLE algorithm used to solve steady shallow water equations have been described in great detail. Since the governin g equations for unsteady flows involve transient terms and the total water depth H is time varying too, the solution for unsteady shallow water PAGE 158 158 equations is required. The easie st and most direct way is to extend the SIMPLE algorithm to transient calculations. The discretized momentum equations now include transient terms, and the resulting discretized uequation at time level 1n outer iteration level m, becomes 1/tiiuu m PP unsteady steadyAAH and /tiin mm uu unsteady steadySSHu (421) where, iu P s teadyA and im u s teadyS are the main diagonal coeffici ents and source terms for steady flow problems respectively, more details referred to chapter 3. is the control volume, t being the time step. All other matrix coeffici ents remain the same. It is similar for the vequation and will not be presented here. For the continuity equation at time level 1 n outer iteration level m, the discretization equation is */t+0mn l lm (422) where l lm is the summation of mass fluxes computed using the intermediate velocities from the momentum equations, for details refer to chapter 3. 1 mm and / p g is the correction of the mean surface elevation from last outer iteration. Substituted into Equation (422), the main diagonal coefficients and the source terms for the pressure correction equations can be given by /gtpp PP unsteadysteadyAA and 1/tmmmn pp unsteadysteadySS (423) At each time level, the iterative procedure is applied until convergence is achieved. Considering that it is easy to get, only a li ttle effort required from the so lution method used in the present steady circulation model, and that it is very r obust, this transient solu tion method is used for PAGE 159 159 solving the transient equations although some other solution methods which could be more efficient, can be developed. 4.2.3.2 Simulation Results The same bathymetry described in previous section (see Equation (4 17) and Figure 411), is used in this section. Noslip boundary condi tions are incorporated at the offshore boundary and at the shore, and the offshore boundary is far enough from the surf zone so that the vortices can not reach the offshore boundary. The peri odic boundary condition is imposed in the alongshore direction. The linear stability analysis in previous section suggests the typical alongshore wavelength is about 200 m, although it may be affected slightly by using different bottom friction coefficients, mixing coefficients and with the presen ce of direction spreading. In this modeling effort, the computational domain is chosen to be 10 times of the predicted most unstable wavelength so that a sufficient number of waves can exist in the modeling domain and the flow is not influenced by the finite size in ydirection. The computation grids with uniform spacing (5dxdym) are used, and the implicit Euler method with 4dts is used. Sufficient iterations are used at every time level to en sure an accurate time evolution of the shear instabilities. Numerical simulations are carried out to investigate the effects of bottom friction, lateral mixing, especially the wave directional spread ing on the resulting mean alongshore currents and then their shear instabilities. Many different expe riment simulations are ca rried out and nine of them are presented here. Simulations are conduc ted for the directional spreading 0, 10, 20 while keeping bottom friction coefficient 0.03akm and the mixing coefficient 0.5 M Simulations are also carried out for varyi ng bottom friction coefficient with a fixed M and S PAGE 160 160 Simulation conditions and maximum of th e resulting mean alongshore current maxV are given in Table 41. All simulations have the initial condition cons isting of the mean base flow superimposed with a small fluctuating component (,,0)()*1() vxytVxSy (424) where, 0.001 is the perturbation amplitude and experiment results proved that the long time behavior of flows are inde pendent of the perturbation amplitude. The function() Sy is given by 1()cos2/J yj jSy jyL (425) where, 10 J is the number of times of the most uns table wavelength fits into the modeling domain. j are random phases. Initial values of crossshore velocity, surface elevation and wave forcing are not perturbed. This initial conditi on ensures that the mean alongshore currents are perturbed at the most unstable wa velength as well as longer wavele ngths can also coexist in the modeling domain. Different ways to set initial co nditions have been used in the literature of nonlinear shear instability simulation. For example, Slinn et al. (2000) started simulations with the fluid at rest but the wave fo rcing terms are perturbed. It is found that the flow behavior after long enough time run (4 hourt) is independent of the way to se t the initial condition. The way to set initial conditions used in this study requires less computer time for the instabilities to reach finite amplitudes. For incoming waves of 1.373rmsHm ,6.4 Ts 0 015 the predicted basestate steady alongshore velocity profiles for the tested experime ntal conditions (see Table 41) are shown in Figure 413. In general, large alongshore velocity results from smaller bottom friction coefficients and smaller mixing coefficients. Wave s with directional spreading lead to relatively PAGE 161 161 weak currents without altering the shape of the ve locity profile or shifting the peak location. Decreasing the mixing coefficient M will increase the maximum velocity over the bar, whereas decrease the velocity in the trough. As discussed above, the resulting steady mean currents are affected considerably by the existence of directional spreadi ng in wave directional spectra. Consequently, characters of shear instability of alongshore currents induced by monoc hromatic waves are supposed to be different from those of currents induced by waves with fini te wave directional bandwidth. For the purpose of investigating the effects of di rectional spreading on shear instabilities, simulation results with fixed ak and M but varying S are compared. Figure 414 shows the time series of velocities (,) uv at 420 x m 100 y m for different values of S while keeping 0.03akm and 0.5 M The time series show that fluctuations reach finite amplitude within about 90 min. For 020 S the instabilities reach finiteamplitude later than that for other two cases. This is because the basestate mean flow associated with 020 S is more stable. For all three cases, time series of acrossshore velocity show regular fluctuations similar to wavelike pattern with nearly constant amplitudes. For 00 S the amplitude of fluctuations in u is about 0.21 m and the period is about 5.5 minutes. For 010 S the amplitude is about 0.19 m a nd the period is about 6 minutes. For 020 S the amplitude decreases to 0.16 m and the period increases to 7 minutes. In summary, the amplitude is decreasing and the period is increasing as direction spreading S increases. This is consistent with the results of linear stability analysis. A longer period and smaller amplitude is resulted from a more stable base flow. Time series in alongshore velocity v at the same location indicate that the alongshore velo city begins to oscillate as the same time as fluctuations in u reaches to finite amplitude. Fluctuations in acrossshore velocity v are PAGE 162 162 irregular although oscillations in u are regular. This may be explai ned by the fact th at the kinetic momentum 221 2 uv at a fixed location is not constant over time due to the instability of the flow. Again, relatively weak oscillations are found for 020 S than those for smaller S For all three cases, the timeaveraged v is well less than the base mean flow velocity (neglecting shear instability) as expected. For analysis of the results of the numerical experiments, it is useful to introduce the vertical component of vorticity x yvu (426) and to define an alongshore average and a time average as follows 01 (,)(,,)yL yQxtQxytdy L (427) 1 (,)(,,)e st t esQxyQxytdt tt (428) where, (,,) Qxyt some variable. yL is the length of the computation domain in ydirection. Using these definitions, the velocities can be defined by (,),(,) uvuvuv (429) The spatial structure of the flow is best illu strated by snapshot of th e vorticity field. Figure 415 shows the contours of vorticit y at several times. It can be seen the disturbances are essentially of constant form propagating in the direction of steady mean alongshore current both for 00 S and 020 S cases. This flow character is similar to the flow described by Slinn et al. (1998) as equilibrated shear waves. The main differences between 00 S and 020 S are the magnitude of the vorticity and the propagation speed. For 00 S the time variability of the PAGE 163 163 flow is further demonstrated by contour plots of the alongshoreaverag ed alongshore velocity (,)vxt and the alongshoreaveraged pe rturbation kinetic energy density 221 (,) 2 uvxt shown in Figure 416. The approximately equilibrated na ture of the disturbances is indicated by the nonzero and nearly constant values of the v and 221 2 uv Convergence Test For numerical simulations, it is very important to verify that the computational results are independent to the grid size and the time step size. Because of limited computational resources and the tremendous time needed for each run, co nvergence tests were pe rformed for a single case E01. Spatial resolutions with 5.5dxdym, 5dxdym and 4.5dxdym were used to examine the effects of grid size on computationa l results. To test whether the computational results vary with varying time step size, simulations with time step sizes of 3.6,4,4.4dtsss were conducted. These tests had all parameters the same as the test case E01 other than the spatial and/or time resolution. Figure 417( a) shows comparison of timeand alongshoreaveraged alongshore velocity of simulation result s using different time step sizes, and panel (b) shows the comparison of perturbation kinetic dens ity. Similarly, Figure 418 shows the timeand alongshoreaveraged alongshore velocity profiles and perturbation kinetic density for different grid sizes. Both the mean alongshore velocity profile and the perturbation kinetic density distribution do not change much with varying the spatial resolu tion and time step size. In addition, fluctuating quantities were also compar ed between these simulations (not shown here) and they are very close to each other. Therefore, it further confirms that the simulation results are independent from grid size and time step size. Effects of Directional Spreading PAGE 164 164 The timeand alongshoreaveraged alongshore velocities for cases E01, E02, and E03 are shown in Figure 419(a). The alongshore averages are performed over the entire width of the domain, while the time averages are obtained neglecting the first 3 hours of the computed time series. The maximum values of the three velo city profiles are 1.134m/s, 1.0675m/s and 0.9m/s respectively, corresponding to maxV values of 1.350m/s, 1.258m/s and 1.026 m/s (see Table 41). For case E01 the maximum alongshore velocity d ecreases about 0.216m/s, while for E03 it only decreases 0.126m/s. This may suggest that more energy is dissipated due to the instabilities for more unstable basestate alongshore currents. The crossshore distributions of the perturbation kinetic density associated with these cases are shown in Figure 419(b). All three cases display similar shapes with two local maxima, shorewar d of the location of m ean current maximum and seaward of the trough of the veloc ity profiles. The perturbation kinetic energy decreases in the surf zone very quickly and comp letely disappears toward the offs hore direction. The differences in kinetic energy density for E01 and E03 ar e pronounced in the bar trough region. The maximum of kinetic energy density for E01 is as large as over 2 times of that for case E03. The timeand alongshoreaveraged acrossshore momentum flux 2hu for the three cases is plotted as a function of x coordinate in Figure 420. It is indicated that the acrossshore momentum flux decreases as S increases, with maxima located around the sandbar crest. In simulation E01, E02 and E03, the relatively large friction is used. As a result, all the computed flows display approximately steady wavelike motion and the presence of a finite directional spreading does not alter the flow regime although qua ntities describing the flows are certainly affected. Next, a smaller bottom friction coefficient ak is used in the second set of experiments, i.e. E04, E05 and E06 (for details see Table 41). PAGE 165 165 As the friction factor is decreased, the maximum velocity maxV of the resulting alongshore currents increases and the currents become more uns table. Time series of velocities for the three cases are shown in Figure 421. Compared to cases with 0.03akm the oscillations in acrossshore velocity u have became irregular (E04 and E05), and the amplitudes of oscillations in both acrossshore and alongshore veloc ity also increase. For case E06, however, approximately steady wavelike flow pattern characters of equilibrated shear waves are observed since the basestate flow is less unstable compared to E04 and E05, a nd the overall flow pattern is different too. The influences of directional spread ing can be further demonstrated by comparing the snapshots of vorticity field for case E04 (00 S ) and case (020 S ), see Figure 422. For case E04, the shear waves develop into longwavelength propagating nonlinear waves. The nonlinear disturbances form, interact and propagate in a complicated manner such that at different times several separate individual di sturbances with different alongs hore scales are observed. For example, at 9.00h t three disturbances with length scal es of approximately 200 m, 600 m, and 800 m can be seen (Figure 422(a)). The changes in the vorticity with time are also pronounced. For case E06 (Figure 422(b)), pr opagating regular disturbances are observed, which are similar to those found in cases with 0.03akm For case E04, the behavior of the flow is also illustrated in the (,)vxt and 221 (,) 2 uvxt contour plots in Figure 423. Small variations in v with time are observed in the bar trough area. Fluc tuations in kinetic energy density 221 2 uv on a very long time scale of about 23 hours are evident in Figure 423(b). The timeand alongshoreaveraged velocity profiles and the perturbation kinetic density for cases E04, E05 and E06 are plotted in Figure 424. The timea nd alongshoreaveraged momentum flux 2hu for 020 S PAGE 166 166 is much narrower and the magnitude is apparently smaller compared to the other two cases, see Figure 425. To further examine the effects of wave dire ctional spreading on the shear instabilities, another set of experiments (E11, E12 and E13) were performed. In these three experimental runs, a constant mixing coefficient 0.3M was used. Nonzero directional bandwidth 010 S was used in E12 and 020 S in E13. The bottom friction factor 0.0235akm for test E12 and 0.012akm for E13 were obtained by fitting the basestate velocity profiles to that of test case E11, in which 0.03akm (see Table 41). Comparison of snaps hots of the vorticity field for the three cases is shown in Figure 426. For all th ree cases, time varying nonlinear shear waves propagating in the alongshore dir ection are observed. However, shear waves with larger alongshore wavelengths from 400 m to 900 m are developed in case E12 and E13 compared to the case E11. The overall flow patterns of the th ree cases are apparently different from each other although the different model parameter set tings produce quite close velocity profiles. Strong energy dissipation due to shear instabilities is evident for all the three cases by comparing their timeand alongshoreavera ged alongshore velocities to the basestate longshore velocity profile (Figure 427a). While, the v profiles are close, which may suggest that the amount of energy dissipated is about the same although thei r flow patterns are different. The timeand alongshoreaveraged pertur bation kinetic energy 221 2 uv is shown in Figure 427b, and it also indicates that energy dissipation due to the in stabilities is close for the three cases. Effects of Lateral Mixing Simulation results of experimental runs E04, E07 and E08, in which fixing 0.01akm and 0 S are compared to investigate the e ffects of lateral mixing coefficient M on shear PAGE 167 167 instabilities. For these three cases the steady basestate alongshore current profiles (Figure 413) clearly show that the velocity gradient is smaller for larger M Simulations of nonlinear evolution of shear instab ilities of these basestate velocity profiles are carried out with the method discussed above. Time series of velocities for E04 (0.5M ), E04 (0.3M) and E04 (0.2 M ) are plotted in Figure 428. It can be seen that initial pe rturbations took longer to the reach finite amplitudes for larger M showing the damping nature of the lateral mixing terms. It is evident that the instabilities are more energetic for lower M This is partly because the basestate velocity profile is more unstable with lower M but also because the weaker eddy viscosity mixing effects reinforce the instabilities. For 0.3 M and 0.2 M the oscillations in velocities display underlying longer time s cales of about 45 minutes. For 0.3M this longer time scale flow behavior is even approxi mately stable. No obvious long tim escale oscillations are observed for 0.5 M Snapshots of the vorticity field for case E04, E07 and E08 at 10.0 th are shown in Figure 429. As M decreases from 0.5 to 0.2, the maxima of both the positive and negative vorticity increase. For 0.5 M the vorticity field is relatively dispersive, nearzero vorticity existing in most region of the domain. For 0.3M areas with strong positive vorticity are obviously narrower. For 0.2 M the features display a longe r alongshore scale and vortex shedding occurs. The timeand alongshoreaveraged perturbation kinetic energy density (Figure 430) also shows that the motions are more energetic over most of the computation domain for smaller M Figure 431 shows the acro ssshore momentum flux 2hu for the three cases. The acrossshore flux of the acrossshore momentum is not altered significantly in the region between bar crest and the shoreline by varying the mixing coefficient M while it increases dramatically as M decreases in the region from sandbar crest to the offshore boundary. In addition, local PAGE 168 168 maxima of 2hu are observed around 250 x m for 0.2,0.3M Finally, it is somewhat surprising to find that the timeand alongshoreaveraged velocity v is very close for varying M although the their basic state velo city profiles are quite different and the instability climate is also very different (see Figure 432). Especially, the v over sandbar and in the trough is almost the same for the three cases. The current V profile neglecting the instabilities is remarkably different, suggesting that th e instabilities induc e momentum mixing. Effects of Bottom Friction Fixing the mixing coefficient and the directi onal spreading, simulations with varying bottom friction factor ak are carried out and analyzed to study the effects of bottom friction on the instabilities. This set of expe riments consists of run E01 (0.03akm ), E04 (0.01akm ), E09 (0.02akm ) and E10 (0.005akm ). Again, the steady basestate alongshore current profiles corresponding to these four cas es can be found in Figure 413. Time series of velocities at the monitoring location are shown in Figure 433. Generally, as the bottom friction coefficient ak decreases, the amplitude of oscillations in the time series of velocities increases, and the fluctuati ons become more irregular. For 0.03akm the oscillations are in regular, wavelike motions and are approximately steady. For 0.01akm the motions are more energetic, and oscillations with th e longer time scales in time series of v can be roughly seen. At the very small friction factor 0.005akm oscillations are initia lly of frequencies of about the order of 5 min. These waves, however, break down eventually after about 2.5th and develop to motions with longer time scales of about 50 min. The alongshore velocity oscillations are very energetic and reach to about 0.5m/s at times with the undistur bed steady velocity of about 2.5m/s. The spatial flow structures are fu rther demonstrated by snapshots of the vorticity PAGE 169 169 field at several times (see Figure 434). Generall y, one largescale distur bance is observed at all five times. The disturbance is strongly nonlinear and propagates in the al ongshore direction, at the same time, it travels back and forth in the acrossshore direction irregularly with time. In addition, the offshore extension of the vorticity disturbance is greatly increased compared with higher ak reaching the offshore boundary. The timeand alongshoreaveraged velocity v profiles and steady current profiles V for the four cases are plotted in Figure 435(a). For all four cases, the mean velocities are considerably damped and the shape of veloci ty profiles are also altered somehow by the instabilities. The instabilities play the roles of momentum mi xing, stretching out the steady velocity profiles in the acrossshore directions. As ak decreases, the velocity profile is altered more significantly by the instabil ities. The time and alongshoreaveraged kinetic energy density (Figure 435(b)) suggests that strong oscillat ions occur at the location of about 420 m. For 0.005akm the kinetic energy is much larger around the bar crest area than that for higher friction, and there is no the sec ondary maximum in this case as in cases with higher frictions. Figure 433 indicates that th e acrossshore momentum flux increases substantially as ak decreases. For 0.005akm strong acrossshore momentum flux presents outside the surf zone, indicating that eddies break away from the mean currents towards the offshore direction. 4.3 Simulation of Rip Currents There have been a great number of efforts in simulating the rip currents using numerical models since the introduction of radiation stress by LonguetHiggins and Stewart in a series of papers (1960, 1961, 1962, 1964). The depthaverag ed twodimensional shallow water equations are usually used as the governing equations in rip current models, and the wave forcing described by the gradients in radiation stress is associated as an input to drive the circulation. Numerical PAGE 170 170 studies have advantages over fi eld or laboratory experiments in respect of allowing direct computation and thus rip currents can be investigated in more details. A variety of models in rip currents have been developed in past decades, early studies including Bowen (1969), Noda (1974), Liu and Mei (1976), and Ebersole and Dalrymple (1980), among others. These studies have generally assumed the currents are weak such that their effects on waves are negligible. Some of them also neglected the nonlinea r convective terms for simplicity. More recently, more s ophisticated modeling efforts have been reported. Haas et al. (1998) conducted numerical simulations of a labor atory rip current system They found that the offshore extension of rip currents can be signi ficantly reduced by including the wavecurrent interaction. Chen et al. (1999) ut ilized Boussinesq equations to simulate the same laboratory set up. In their model, the wave motion and induced currents were always solved simultaneously, therefore the mutual interaction is always include d. Therefore, it is not possible to compare the results with and without inter action. Yu and Slinn (2003) con ducted a numerical study of rip currents on a barred beach with gentle sinusoida l alongshore variations. Th ey investigated the effects of wavecurrent interac tion using the linear bottom friction model and ignoring the lateral mixing. Most the modeling efforts in rip currents to date have been for monochromatic (single frequency and unidirectional) waves. The wave forcing of irregular waves, as discussed in previous sections, can be differe nt to some degree. Therefore, the resulting waveinduced rip currents may be influenced somewhat by the presence of wave frequency and directional bandwidth. It is still uncertain that how the fre quency/directional spreadin g affects rip currents. Using the present wave model, the effects of fr equency/directional spreading on rip currents can be studied in detail, which is very useful. PAGE 171 171 Oscillations or pulsations in rip currents over periods of minutes or tens of minutes have been noticed and recorded sin ce very early in rip current st udies. Shepard and Inman (1950) observed fluctuations with peri ods varying from 1.7 minutes to 7.8 minutes. Many other field observations (Mackenzie, 1958; Sonu, 1973; Smith and Largier, 1995; Huntley et al., 1998; Brander and Short, 2001; MacMah an et al., 2004) have also ob served energetic rip current velocity pulsations. In addition, Haller and Dalyrmple (2001) docum ented oscillations of the rip current jetlike motions in a laboratory study. Many field experiments have overlooked the amount of energy associated with the rip current pulsations du e to not sampling long enough to resolve them. Ogston and Sternberg (2003) reana lyzed some measurements using longer records and found that there is a great am ount of energy associated with the low frequency oscillations. Fluctuations in rip current strength at a location can significantly increase the danger to swimmer. Therefore, it is very important to study the ri p current pulsations and their causes. Two causes have been recognized to be res ponsible for rip curren t pulsations: timevarying wave groups and flow instabilities. In present study, we will conc entrate on flow instabilit ies rather than wave groups. In this section, two main topics are investig ated: wavecurrent interaction on rip currents, and rip current instabilities. The objectives of this section include: improving the understanding of the physical processes of wave current interaction, and its e ffects on rip currents; searching the relation between frequency/di rectional spreading and rip curre nt features; and finally, rip current stabilities. 4.3.1 Steady Rip Currents It has been known that wavecu rrent interaction could change the structure of wave forcing dramatically and is very important in numerical simulation of rip curr ents (Haas et al., 1998; PAGE 172 172 Chen et al., 1999; Yu and Slinn, 2003). Numeri cal simulations including effects of bottom friction and lateral mixing were carried out in order to obtain a better understanding of the physical processes involved in the wavecurrent interaction as it affects waves and rip currents. For cases where currents are strong, convergence c ould be difficult to achieve using the twoway coupling of the wave model and the steady versio n of the circulation model as discussed in chapter 3. Here, we used the wave model coupled with the unsteady version of the circulation model described in Section 4.2.3. Specifically, both the wave equations and the shallow water equations were solved in the time domain. At each time step, the wave model computed the instantaneous wave forcing which was then used to drive the circulation model, and the wave field and velocity field computed we re then used as the initial values for the next time step. This time marching process was continued until a convergent solution had been reached both for the wave field and for the velocities. The velocity fields for cases without waveinteraction were computed by using the steady circulation model along with the wave model. The bathymetry from Kennedy and Zhang ( 2008) was used thorough this section. The bathymetry was parameterized as a planar beach with a superimposed bar. The bar is Gaussian in crosssection and itself gives the alongshore variations in bathymetry, i.e. 2 2 00,e x p / 2bb bhxyhxxhyxx (430) where 0h is the offshore water depth, and 0 x is its corresponding crossshore coordinate, is the slope of the planar beach, 80b x m is the crossshore location of the bar crest, and b is the standard derivation of the Gaussian bar, which determines the width of the bar in acrossshore direction. For all experiments, 50/(32)mb was used. The longshore variation of the bar height is given by bh PAGE 173 17321,be e bhyCnUMA (431) where Cn is a Jacobian elliptic function, and eU and e M are the elliptic argument and modulus respectively. e M is in the range of 0 to 1, and a larger value of e M gives a narrower rip channel. bA is the bar amplitude and equal to 1m for all cases. For most cases, 0.95eM (a narrow rip channel) was used. In this section, 0.95eM unless stated otherwise. The argument ()2(/1/2)eyeUyyLK where ()eeKM is the complete elliptic integral of the first kind. Figure 437 shows the rip channel bathymetry with 0.95eM The noslip wall boundary condition was app lied to both the offshore boundary and the shoreline boundary. In the alongs hore direction, periodic bou ndary conditions are imposed on ,, uv because of the topography and the flow pattern we are interested in. For incoming waves, an equivalent deep water wave height of 01mrmsH and thus the RMS wave height at the offshore computational boundary is computed according to 000cos/cosrmsrmsg gHHcc with angle 00 period 10mTs For all the experiments presented below, 5m xy and the lateral mixing coefficient 0.8M were used. For cases including the wavecurrent interaction, 0.2st was used (relatively small because of the stability of th e wave model). For most cases, the computational domain is 300xLm and 200yLm To conduct a comprehensive study of the rip current system, a number of cases were simulated and one case was selected as the base case for detailed study. Departures from this base case were then examined to evaluate the importance of various parameters. These variations PAGE 174 174 included the effects of bottom friction, wave heig ht, rip channel width an d frequency/directional spreading as well. Base Case In the base case, monochromatic waves with eq uivalent deep water wa ve height is 1 m and the bottom friction factor 0.08akm Figure 438 shows the comput ed wave field without (a) and with (b) the wavecurrent interaction for the base case. Figure 439(a) shows velocity field without taking into account the wavecurrent interaction. For th e same parameters, the results with the wavecurrent interaction are given in Figure 439(b). Some typical features on the influences of the wavecurrent interaction upon the wave field and the circulation system may be summarized as follows 1. Without the wc interaction, wave shoaling occu rs as the normally incident waves propagate toward the shore. After approach ing to the bar region, waves re fract from the rip channel and focus over the bar crest area, which leads to st ronger breaking over the bar than in the rip channel (Figure 438(a)). The longshore variation in wave he ight then resulted into a classic current pattern with circulation cells and a narrow offshoredirected jet. 2. With the interaction, noticeable changes in wa ve height due to the ambient currents are observed: in the region from 90 x m to 160 m around the rip channel, the offshoredirected currents shorten the waves and c onsequently the wave height is increased. On the other hand, waves are lengthened by the onshore current in the area from the bar to the shore, which decreases the intensity of wave breaking and cons equently wave height is slightly increased. Owing to the refraction effects by the circulatio n cells, dramatic changes in wave direction are observed. The changes in wave field due to the wc interaction results in a forcing effect PAGE 175 175 opposite to that due to the topography. As a resu lt, with the interacti on the rip currents are noticeably weaker and circulations are mostly within 150 m from the shoreline in this case. 3. Without the interaction, rip cu rrents reach to the offshore boundary even with the relative large parameter values of bottom friction factor ak and lateral mixing coefficient M With the interaction, the offshore extent of rip cu rrents is greatly redu ced. Figure 440(a) shows the comparison of the uvelocity along the center of rip channel with and without the interaction. Rip currents are also noticeably broader with the inter action and uvelocity comparison along the longshore transect at 100 x m is given in Figure 440(b). Figure 441 shows comparison of the wave setup along th e rip channel and the bar crest. It was observed that the wc interaction did not alter the wave setup significantly and the maximum difference is only about 1 cm. Bottom Friction In modeling nearshore currents, the bottom fr iction is usually a source of uncertainty because the bottom friction coefficient is usua lly unknown and also the formulations of bottom friction are roughly approximated. To test the effects of bottom fric tion, several different values of the bottom friction factor ak were used while other parameters were remained the same as in the base case. Figure 442 shows th e uvelocity profiles along the rip channel cente r with various values of ak The shapes of velocity profiles are quite similar for all the bottom frictions and have maxima (offshore velocities) located at 105 x m and minima (onshore velocities) located around50 x m The maximum offshore velocity is not sensitive to the bottom friction coefficient and changes only about 10% from a small friction coefficient case (0.05akm ) to an unrealistic large friction case (0.2akm ). As ak increases, the onshore velocities decrease apparently although the offshore velocities onl y decrease slightly. Figure 443 shows the PAGE 176 176 offshore extent and the maximum/minimum velo cities versus the bottom friction factor ak The offshore extents are given by the ac rossshore locations where the o ffshore velocities decrease to 0.01m/s. We can see that the offshore extent of the rip currents is w eakly dependent on the bottom friction. It is also show n that all the velocity extreme values are smaller for a larger bottom friction, which is certainly reasonable si nce more dissipation asso ciated with a larger friction. It should be also noted that bot h the velocities and the offshore extent only change slightly as the bottom friction varies in a very wide range, which definitely suggest that the bottom friction becomes much less influential to the prediction of the circulation system. This is not the case when the effects of wc interaction are not taken into account because the bottom friction is the only major dissipative effect on the flows. Re sults of the numerical experiments without the wc interaction also suggest th at the strong dependence of the predicted flow on bottom friction (not shown). With the interaction, decreasing th e bottom friction will result in a stronger current, which in turn leads to decreased forcing to the current. This selfgenerated dissipation is dominant in deep waters where bottom friction is relatively small. It also explains why the bottom friction has more impact on velocities in the region from the alongshore bar to the shoreline. Wave Height As wave heights increase, the breaking induced currents are expected to become stronger. However, wavecurrent interaction acts as an op posing factor to focus wave breaking in the rip channel and thus reduce forcing, and breaking will saturate offs hore of the bar for large wave heights, so the overall behavior may be very co mplex and is not well understood. To test the effects of the wave heights, another set of e xperiments was performed over equivalent deep PAGE 177 177 water wave heights of 0.6, 0.8, 1.0, 1.2 and 1.5 m. Again, normal incident waves are used and the bottom friction factor ak is fixed to be 0.08 m in these experiments. In Figure 444(a), the mean water level (wave setup/setdown) relative to the still water level along the line 100 y m for various wave heights is plotte d versus the acrossshore distance from the shoreline. The magnitude of wave se tup increases monotonically with the deep water wave height, and the setup at the shore line increases approximately linearly with 0H It is also noted that as 0H increases, the onset location of wave setdown moves seaw ard considerably. Figure 444(b) shows the uvelocity prof iles with different wave heights. For 00.8m H the maximum rip current velocity is located approximately at 105 x m For 00.6m H it moves shoreward slightly to 95 x m The maximum velocity increases significantly as 0H varies from 0.6 m to 0.8 m but only slightly as 0H increases from 0.8 m to 1.5 m. This is because for 00.8m H wave height distributions are similar, while for 00.6m H wave breaking occurs only in very shallow water and results in a different wave hei ght distribution pattern. It is also evident that the maximum rip current velocity increases with wave height, as shown in Figure 445. The onshore maximum velocity shows a little variati on as wave height incr eases and the maximum velocity of the feeder currents only increases slightly with wave height. The offshore extent of rip currents initially decreases quickly with 0H and then shows a slow decrease when 01.2m H Rip Channel Width Rip currents have been observed for a wide range of topographies. Almost any geometry that induces a significant longs hore variation in wave forcing will see rip currents when the incident waves are close to s horenormal. Kennedy et al. (2008) showed in a laboratory PAGE 178 178 experiment that rip currents with varying gap widths had similar peak velocities although the flows could be very different. In this set of experiments, the eff ects of relative rip channel width on the present system are investigated by varyin g the rip channel width. A relative rip channel width may be defined as the ratio of rip channe l width to rip current periodic length. For the bathymetry described above, this is controlled by th e elliptic modulus e M which for all previous tests was set to 0.95eM (relative rip channel width 0.336) Figure 446 shows the longshore variation of bar/rip channel with varying e M 0eM corresponds to a sinusoidal variation (relative width 0.5). The narrowest rip channel used had 0.99eM (relative width 0.2). The relative rip channel width is taken at the longsho re location where the bar height decreases to half its maximum height. In this set of experiments, wa ves are normally incident with 01 Hm and the bottom friction factor ek is fixed to be 0.08. Figure 447 shows the acrossshore uvelocity profiles for various 0,0.50,0.75,0.95,0.99eM Generally, the stronger curren t is observed for a narrower channel (larger e M ), and the velocity profile does not ch ange its acrossshore shape as the rip channel width changes. The longs hore uvelocity profile at 100 x m changes significantly as e M varies from 0 to 0.99, as shown in Figure 448. This implies that the rip current becomes narrower as the rip channel becomes narrower, whic h can be seen even more clearly in Figure 449. As rip channel becomes broader, it is also observed that the maximum velocity of the rip current decreases noticeably, while the maximu m feeder velocity increases slightly. The maximum onshore velocity does not change monotonously with rip channel width, initially decreasing then beginning to increas ing with rip channel width. Frequency/Directional Spreading PAGE 179 179 Monochromatic waves have been used in all th e previous tests. As discussed in previous sections, with certain frequency and directional bandwid th the resulting wave radiation stresses could be considerably different from that of monochromatic waves, especially for irregular waves with a relatively broad dir ectional spreading. Consequently, the wave forcing and then the flow pattern could be different. The effects of frequency/dire ctional bandwidth on waveinduced longshore currents and on the nonlinear instability of longshore currents are investigated in Section 4.2. To investigate the frequency/direct ional spreading on this waveinduced rip current system, a series of numerical experiments was conducted with different values of frequency bandwidth S and directional bandwidth S For all irregular wave cases, a Gaussianshaped spectrum was used with mean wave period of 10 s and mean incident angle of 00 (shorenormal). The equivalent deep water wave height 01mrmsH and the friction factor 0.08akm are used in these simulations. To examine the effects of frequency bandwidth, results from cases with relative frequency bandwidth /0.1,0.2,0.3mS are compared to those of the monochromatic wave case. Comparison of the uvelocity profiles of the re sulting rip currents in Figure 450 does not show appreciable differences between these cases. This is not a big surprise partly because as shown in previous sections the wave field simulated by the present wave model is not significantly affected by a relative small freque ncy bandwidth with the Gaussian spectrum assumption. It is also because the topography induced wave tran sformations are dominant as irregular waves propagate over such an irregular bathymetry. In addition, the assumption that an input spectrum conserves its spectral shape used in the wave model may not hold here because strong wave focusing and scattering have occurred in these cases. PAGE 180 180 To asses the effects of the directional bandw idth on this rip current system, simulations with different directional bandw idth were carried out. A numerical instability problem was encountered for larger directional bandwidth (015 S ), and it is believed because wave caustics may occur when the tested direct ional spectra are too broad. So the directional bandwidth used here is limited to be no larger than 10. Figure 451 shows the uvelocity profiles at 100 x m corresponding to the directional bandwidth 00000,4,8,10 S It can be seen that the uvelocity profile of the rip current only show s a little variation with directi onal bandwidth, with almost the same maximum velocity of the rip current but the maximum onshor e velocity decreasing noticeably with S These observations are confirmed in Fi gure 452, and it also shows that the maximum velocity of the feeder currents do not change much either. 4.3.2 Stability of Rip Currents Field measurements indicate that rip currents can have very long peri od oscillations (e.g. Shepard and Inman, 1950; Sonu 1973; Huntley et al. 1988; MacMahan et al., 2004, and many others). The lowfrequency motions of rip curren ts have also been r ecognized in laboratory experiments (e.g. Haller a nd Dalrymple, 2001; Kennedy and Thomas, 2004). Two mechanisms have been suggested to be res ponsible for the pulsati ons in rip currents: long period modulations in incoming wave heights i.e. wave grouping, and instability of the hydrodynamic flow itself (Haller and Dalrymple, 2001). In this study, we will restrict ourselves to fl ow instabilities. Haller and Dalrymple (2001) performed an analytical study to investigate the linear stability characteristics of rip currents in a laboratory experiment, and they modeled the vicinity of the rip by considering it as a slowly varying, long narrow jet and obtained good agreement with laboratory data in the vicinity of the rip neck. To obtain typical pictures of instability over the whole circulation domain, Kennedy and Zhang (2008) developed a numerical model to perform PAGE 181 181 the linear stability analysis to a waveinduced ri p current system by using the similar strategy as used to the longshore currents (see Section 4.2.2). The effects of wa vecurrent interaction, lateral mixing and bottom friction were included in their model. Numerical simulations of the lowfrequency mo tions of rip current and circulation cell have not received sufficient studies to date. Slin n et al. (2000) simulated the flow instabilities over rip channel topography. Howe ver, all computations were performed using a very large incident angle of 45 degree, which resulted in long shore currents rather than rip currents for all cases. Yu and Slinn (2003) modeled instabilities of rip currents induced by normal (or nearly) incident waves breaking over artificial rip channel bathymetries. They found that the instabilities of the cellular circulations are sensitive to the angle of wave incidence and rip channel spacing. They also observed that the inclusion the wave current interaction not only reduced significantly the mean forcing available to the rip current but also changed the instabilities of rip current: onset of instabilities occurs at the nearshore region, rather th an offshore at rip heads with nocurrent interaction. These findings are very useful and can provide us useful insights into some typical features of rip current in stabilities. While they used a gen tle rip channel bathymetry with the sinusoidal perturbation in the ydirection and small bar heights (0.1 m and 0.2 m used), and the resulting currents were rela tively weak, with the maximum rip current velocity of about 0.2m/s0.3m/s for most cases. For strong rip currents, for example in which the rip velocity is as large as 0.5m/s or larger, the st ability characteristics may be quite different from those for weak rip currents and were not suffici ently studied. In addi tion, the lateral mixing effects were not included in their model and they used the simple linear bottom friction model with a constant bottom friction coefficient. In this study, we speci fically simulate the inst abilities of rip currents PAGE 182 182 with a maximum rip velocity of about 0.5m/s, including the effects of lateral mixing and bottom friction with acrossshore variable coefficients. In the previous section, the simulation result s showed that steady strong rip currents with a maximum offshore directed velocity of around 0.5m/s or larger resulted with the fairly large bottom friction 0.05akm and lateral mixing 0.8 M for both cases includ ing wc interaction and neglecting the interaction. Using the same wa ve parameters and bathym etry, the rip currents will become unsteady and the instabilities begin to grow if the bottom friction coefficient and the lateral mixing keep decreasing. The wave field is expected to vary temporally due to the existence of the timedependent circulation fiel d, and the resulting wave forcing in turn will affect the current field. This feedback mechanism, i.e. wavecurrent interaction, can be simulated using a timedependent wave model along with an unsteady circulation model. While, numerical difficulties have been encountered in simulating the unsteady energetic rip current system using the present coupled wave/circulati on model. Possible reasons that ar e believed to be responsible for the failure include: (1) possible occurrence of wave ray crossing, caustics and/or wave blocking, which is beyond the capability of th e present wave model; (2) The explicit time differencing used in the wave model is not compatible very well with the implicit time differencing scheme in circulation model. Because of the difficulties, the instabilities of rip currents are simulated excluding the effects of wavecurrent interac tion in this study. Specifically, the steady wave forcing is used to drive the unsteady rip current syst em instead of the timedependent forcing. 4.3.2.1 Experimental Plan Although the feedback between waves and curr ents can not be fully simulated by the present wave/circulation model due to some numerical difficulties when the currents are PAGE 183 183 unsteady and energetic, the impact of wavecu rrent interaction on the resulting unsteady rip currents can be investigated to some degree through comparing the unsteady flows driven by wave forcing with and without th e wavecurrent interaction. Two sets of numerical experiments were performed: one completely neglecting the wavecurrent interaction; the other using timeinvariant wave forcing in which th e wavecurrent interaction is included. Specifically, in the first set of simulations, the wave forcing was calcula ted based on the steady wave field simulated by neglecting the effects of the ambient current on waves, and the wave forcing was assumed invariant as the resulting rip cu rrents become unsteady. In the second set, the wave forcing was from a steadystate solution of the wave and current system with sufficiently large bottom dissipation and lateral mixing coe fficient including the wavecurrent interaction (as did in the previous section). The wave forcing then was used to drive the circulation. Rip currents are known to be strongest with normally incident waves. However, both the mean properties and the stability of rip currents may change considerably with relative small changes in incident angle. In their study cases, Yu and Slinn (2003) found th at instabilities of the rip currents were sensitive to the angle of wave incidence, when the flow is dominated by cellular circulations a small longshore flow te nds to destabilize th e circulations, e.g., 01 and 03. As the longshore component becomes strong, it suppresses the circulation cells, and the flow starts behaving like a longshore currents. For relative strong rip currents, it is still uncertain whether the angle of wave incidence will affect th e development of rip current instability in the same way. Another set of experiments were perform ed to investigate the effects of wave incident angle on rip current instability. Expe riment setup and parameter settings are given in Table 42. In test cases T03 and T04, the wave forcing is computed based on th e wave field without including the wavecurrent in teraction and remains invarian t during the simulation of the PAGE 184 184 unsteady rip currents. In all other test cases, the wave forcing is from the simulation results of steady state wave and rip current system with th e wavecurrent interaction being accounted for, in which the bottom friction factor0.06akm and the lateral momentum mixing coefficient 0.7M were used. 4.3.2.2 Results To decrease the influence of the latera l boundaries on the unsteady flow motions, simulations are carried out for the bathymetry composed of three rip channels. For all simulations, the computation domain is: 300xLm and 600yLm Rip channels are located at 100y, 300 and 500 m. The rip channel bathymetry an d incident wave conditions are the same as used in the base steady case in previous sec tion except that the incide nt wave angle may be not normal for some cases. The noslip wall cond ition is applied to both the offshore and onshore boundary, and the lateral boundaries are periodic. For all cases, the grid size of (5dxdym) is used. Because the wave model does not need to be solved and the robust implicit discretization scheme is incorporated in the circulation model, a fairly large time step of 2 s is used for most cases. All computations cover durat ion of 15 hours. The time averag es are obtained using the last 10 hours of the time series. The lateral mixing coefficient 0.3 M was used for all the simulations. Convergence Convergence tests were performe d for test cases T01 and T08. For each test case, simulations with several grid sizes of 6,5,4 dxdymmm were conducted to investigate the effects spatial resolution upon the simulation results, and different time step sizes of 2.5,2.0,1.6 dtsss were used to examine the effects of the time step size on computational results. Figure 553 shows the comparison of the timeaveraged vorticity field between PAGE 185 185 simulations with different grid sizes and time step sizes for test case T01. For the different time steps, the vorticity fields are very close. Using different grid sizes, while, evident differences in vorticity field appear. For the finer grid resolu tions, the rip currents reach further offshore, especially, for the grid size of 4 dxdym Comparison of the acrossshore velocity profiles along the rip channel also suggest s that the time step sizes have been small enough, while the grid resolutions affect to some degree the mean ve locity (see, Figure 554). For the test case T08, comparison of the timeaveraged vorticity field is shown in Figure 555. It can be seen that simulations with these different time step sizes and grid sizes are very clos e in the mean vorticity field. Figure 556 shows the comparison of acros sshore velocity profiles along the rip channel for these simulations. The time series of velocities at several monitored locations were also compared and no significant difference was obs erved. Since it is difficult to show the comparison of time varying quantities, only the co mparison of mean quantities is presented. Generally, these convergence tests verify th at the time step of 2 s is small enough and using small time steps will not affect the comput ational results much. The grid resolution has more impact on the simulation results, and using finer resolutions than 5 dxdym could change the results to some degree. But the overall flow pattern will not alter significantly. As the result of negotiation between computational time and accuracy, the grid size of 5 dxdym and 2dts are used in all the following unsteady rip current simulations. WaveCurrent Interaction The effects of wavecurrent interaction on the time evolution of rip currents are examined by comparing test cases T01, T02, T03 and T04. The bottom friction factor 0.02akm is used in T01 and T03, and a smaller friction 0.01akm is used in T02 and T04. Figure 457 shows the time series of velocities (,) uv at the location of (,)(100,302.5) x ymm for the four cases. The PAGE 186 186 flows are time dependent for all the four cases. For 0.02akm (T01 and T03), fluctuations in velocities grow to be apparent after several hours of simulation (about 5 hours for T01 and 8 hours for T03). For the smaller friction 0.01akm cases (T02 and T04), the instabilities grow much more quickly. Regular oscillations with th e period of 12 minutes in velocities are observed for the simulation results of case T01 after the fl ow developing into an equilibrated state, and oscillations are somewhat irregula r for all other three cases. For wave forcing without the effects wavecurrent interaction (T03 a nd T04), fluctuations in velocities display longer time scales (20~50 minutes) and are more irregular compared to T01 and T02. Comparison of vorticity between case T02 and T 04 is shown in Figure 458. The snapshots are separated by a 12minute time interval. For test case T02, the vorticity is mainly concentrated between the shoreline and 150 m fr om the shoreline. For wave forcing without wavecurrent interaction (test case T04), however, the spa tial distribution of vorticity is all over the computation domain, positive and negative vorticity existing even in the adjacent region of the offshore boundary. In inspection of temporal varia tions of the vorticity for T02 shows that the very low frequency (VFL) motions in rip currents are more or less trapped within the offshore extent of 150 m and they do not propagate seaward away from the surf zone. In contrast, in case T04 the VFL motions move back and forth from the surf zone to the offshore boundary. The strength of the vorticity for the two cases is surprisingly close alt hough the shape and position are quite different from each other. It should be al so noted that the flow patterns at different rip channels are quite different and the flows are complex for both test cases. The crossshore mass flux and momentum flux associated with a rip current system may be as a measurement of the intensity of the rip current s, and may also be of interest in the study of sediment transport phenomena and then morphol ogical evolution in the nearshore areas. The PAGE 187 187 magnitude and spatial distribution of fluctuations of the rip current can be examined by the timeaveraged turbulence kinetic energy, 221 2 TKEuuvv Figure 459 shows the spatial distribution of the timeaver aged acrossshore mass flux hu timeaveraged acrossshore momentum flux huu and timeaveraged turbulence kinetic energy TKE for test cases T01, T02, T03 and T04. For cases with wave forcing neglecting the wavecurre nt interaction (T03 and T04), more water and more momentum are tr ansported offshore through the rip channels compared to case T01 and T02. Although energetic, irregular oscillations in velocities are observed in T01 and T02, the mean flow pattern (not shown) is similar to the steady state rip currents. It should be also noted that hu and huu are generally slightly smaller for the smaller bottom fiction cases, which may imply that more energy is dissipated by the instabilities. The turbulent kinetic energy shows that the fluc tuations are mainly lo calized in the feeder currents and rip currents, with lo cal maxima of TKE located in th e rip channels for case T01 and T02. For case T03, however, the TKE is associated with the offshore rip heads. The TKE for test case T04 suggests that oscillations occur over the entire circulations with local maxima associated with the feeder and the rip currents. It is also evident that a smaller bottom friction results into more TKE by comparing T01 and T02, and T03 and T04. Wave Incident Angle To examine the response of the resulting unstead y flow to small deviation in angle of wave incidence from the shorenormal di rection, computational experime nts were performed for wave incident angles of 1, 2, 3, 4, 6 and 8, with the same parameter settings in the wave model and circulation model (for details see Table 42). In these simula tions, the wave forcing used is PAGE 188 188 calculated based on the simulation results of stea dy wave fields with wa vecurrent interaction being accounted for, in which 0.06akm and 0.7 M Figure 460 shows the time series of velocities (,) uv for simulation cases with various wave incident angles at the location of 100 x m and 302.5 y m where is seaward of the transverse sandbar around the rip channel cente r. At the very small angle of incidence 01 oscillations are initiated intermediately and deve lop to an approximately constant pattern after about 4 hour simulation. The time series show that the oscillations have the period of O(10 min) with magnitude of about 0.2m/s. Oscillations in the acrossshore velo city also display the underlying longer time scales of about 80 minutes. When 02 significant differences in time series of velocities from the case 01 are observed: high frequency oscillations are observed at first. After about 3.5 hours, oscillations begin to display much longer time scales of about 53 minutes with smaller magnitudes; the mean value of the acrossshore velocity is also reduced significantly. At 03 oscillations in velociti es can still be roughly s een, but the magnitudes of fluctuations decrease further to very small values. As further increasing to 40, a very different flow pattern results and energetic fluctuations with time scales of about 135 minutes are observed. It appears that with in one cycle of the periodic case, the acrossshore velocity u is close to 0.1m/s most of the time and increases to over 0.5m/s rapidly in a short time period, meanwhile the longshore velocity v is also at its peak value. At a still larger angle of 06 evident longshore character in the mean flow is observed, but the fluctuations in velocities have smaller magnitudes, especially in the crossshore velocity, and shorter time scales of about 80 minutes. At the largest incident wave angle of 08 in this set of experiments, the longshore current character is more pronounced and the unsteady flow develops into a wave motion propagating in the longshore directio n. The time series of velocities suggest the os cillations are PAGE 189 189 harmonic with a period of about 26 minutes. It should be also noted that both the mean value and the oscillating magnitude of the acrossshore velocity are small for this case, whereas the mean longshore velocity is about 0.45 m/s. This also implies that the flow becomes longshore current dominated. From plots of time series of velocities uv, i.e. Figure 460 along with Figure 457 (00, case T01), it seems that the resulting flow is quite sensitive to the wave incident angle. For the case with shorenormal incidence, the fl ow instabilities take a few hours to develop and very strong oscillation in veloci ties are observed. A small incident wave angle tends to trigger the instabilities much quickly, however, the magnitudes of oscillations decreases, e.g. 01 02 and 03. Irregular oscillations in velocities with the time scales of from 1 to 2 hours are observed and longshore character begins to be evident for larger wave angles, i.e. 04, 06. As further increasing the angle of inciden ce to a certain value, i.e. 08 longshore currents become dominant and the flow acts as the typical shear waves. Snapshots of vorticity for cases with the angles of incidence 02 04and 08, each of which represents a typical flow pattern (see Fi gure 460), are compared at several times, as shown in Figure 461. For 02, the sequence of vorticity fiel ds roughly covers one cycle of periodic oscillations in Figure 458. Let us fo cus on the patches of vorticity around the north rip channel. During the first half of the cycle, the patch of the positive vorticity close to the shoreline stretches clockwise with decrease in strength and the negative patch below shrinks. In the second half of the cycle, this process reverses. Th e vorticity patches around other rip channels do not appear to change their shapes and position significantly. For 04 snapshots of vorticity only cover about half cycle of the periodic cas e in Figure 460. During the time from 14.0th to 14.6 th the positive patch close to th e shoreline around th e south rip channel stretches to the PAGE 190 190 north and also seaward and the strength is decr easing. In the meantime, the positive patch around the north channel undergoes an opposite proces s. In addition, the vorticity pattern around the center rip channel does not alter very much during this time period. At 15.0th when both the acrossshore and longshore veloci ties are at their peak values (see Figure 460), however, strong circulation cells are observed and the offshore patch the positive vorticity turns clockwise toward the offshore direction. For 08, the vorticity fields suggest the circulation cells have virtually disappeared and unsteady longshor e currents have been develope d both close to the shoreline and over the sandbar. Fluctuations in the vorticity fields appear to have the longshore wavelength of 200 m, which is also the l ongshore length scale of the rip ch annel bathymetry. This suggests that the topography has its signature on the flow instabilities. The timeaveraged vorticity field, acrossshore mass flux hu and turbulent kinetic energy 221 2 TKEuuvv for the set of experiments are shown in Figure 462. As the wave incident angle increases from 1 to 8, the rip cu rrents gradually ti lt away from the shorenormal direction and th e circulation cells develop to longshore currents. At 02 and 03, the similarity of timeaveraged mean flow patterns at different rip channels breaks down, and the cause for this is not clear. At other angles, th e similarity is observed. The timeaveraged acrossshore mass flux decreases considerably as increases from 1 to 8. The timeaveraged turbulent kinetic energy pr ovides a good measure of the in tensity and distribution the fluctuations duet to the flow instabilities. Str ongest turbulent kinetic energy is observed when 04 and 06, which is consistent with energetic osci llations in the time series of velocities shown in Figure 453. For 01, the TKE is mainly associated with the rip currents and shows similar pattern at different rip channels. For 02 and 03, the TKE is small and can be only PAGE 191 191 seen at one rip channel. At 08 the TKE is not energetic but clearly associated with the meandering longshore currents, nearly absent from the rip channels. 4.3.3 Summary The waveinduce rip currents ove r an artificial rip channel ba thymetry have been simulated using the coupled wave/circulation model in this section. Both the steady state rip currents and unsteady motions (instabilities) of rip currents we re examined. In the simulations of steady state rip currents, the wavecurrent interaction is fully represented by using the unsteady circulation model along with the wave model through the ra diation stress concept. However, the wave feedback was not included in simulating the unsteady rip currents due to some numerical difficulties, and the stationary wave forcing was used. The bottom friction model with acrossshore variable friction coefficient is used in these simulations. The lateral momentum mixing is also included. From simulation results for the steady state rip currents, several conclusions can be drawn. For 10 s waves with equivalent deep water wave height of 1 m, fairly strong rip currents with maximum offshore velocity over 0.5m/s result over this artificial rip channel bathymetry (see Equation (430) and Figure 437) even for mode rately large bottom friction and lateral mixing. Both the wave field and the current field are a ltered considerably by including the wavecurrent interaction. Generally, the wc interaction produced an opposite fo rcing effect to that due to topography. As a result, the rip cu rrents are weaker, slightly broade r and dramatically reduced in offshore extent with the wc interaction being acc ounted for. Second, with the interaction the rip currents are much less sensitive to changes of the bottom friction. This is because that decreasing the bottom friction will result in a stronger current, which in turn leads to decreased forcing to the current. This selfgenerated dissipation m echanism supersedes the effects of bottom friction PAGE 192 192 significantly. These two conclusions are in good agreement with th e previous simulation results by Yu and Slinn (2003). It is also found that larger incident wave s lead to a larger wave setup and wave setup at the shoreline increases appr oximately linearly with the wave height. For incident wave heights leading to a similar wavebreaking pattern, how ever, the resulting rip currents alter only slightly. The rip channel width appears to have a bigger impact on rip currents compared to the incident wave height. A narrow rip channel leads to a strong and narrow rip current without changing the shape of the acrossshore velocity profile along the rip channel center. Moreover, it is suggested that the presence of limite d wave frequency and direction spreading does not change the rip current apparently by the results of simulation experiments. For rip current instabilities, we find that the very low frequency motions in rip currents are confined in the surf zone by using the wave forc ing with wc interaction rather than moving back and forth from the surf zone to the offs hore boundary as predicte d without including wc interaction in the wave forcing. Finally, the instabilities are sensitive to the angle of wave incidence in terms of periods and magnitudes of the oscillations in velocities. The unsteady rip currents display a wide range of flow patterns as the wave inci dent angle varies from shorenormal to 8. Under the current experiment conditi ons, the strongest fluctua tions with very large time scales (from 1 to 2 hours) are resulted from the wave incident angle of 04 and 6. As increases from 0 to 3, the energy content of the instabilities decreases continually. At larger values 06 and 8, circulation cells begin to disapp ear and the flow acts like shear waves with the spatial length scale in longshore direction equa l to that of the bathymetry. PAGE 193 193 Table 41. Experimental conditions and paramete r settings for shear in stability simulations. ak is the bottom friction coefficient, M the lateral momentum mixing coefficient, S the directional bandwidth of incident waves, maxV the maximum longshore velocity of the base mean longshore currents. Run name ak (m) M S (0) maxV (m/s) E01 0.03 0.5 0 1.350 E02 0.03 0.5 10 1.258 E03 0.03 0.5 20 1.026 E04 0.01 0.5 0 1.815 E05 0.01 0.5 10 1.694 E06 0.01 0.5 20 1.385 E07 0.01 0.3 0 1.974 E08 0.01 0.2 0 2.088 E09 0.02 0.5 0 1.499 E10 0.005 0.5 0 2.192 E11 0.03 0.3 0 1.452 E12 0.0235 0.3 10 1.451 E13 0.012 0.3 20 1.440 Table 42. Experimental conditions and parameter settings for rip current instability simulations. WFinteraction stands for whether wave forcing is computed with the wavecurrent interaction being accounted for; 0 is the angle of wave incidence, ak the bottom friction coefficient, M the lateral momentum mixing coefficient. Run name WFinteraction 0 (0) ak (m) M T01 Yes 0 0.02 0.3 T02 Yes 0 0.01 0.3 T03 No 0 0.02 0.3 T04 No 0 0.01 0.3 T05 Yes 1 0.02 0.3 T06 Yes 2 0.02 0.3 T07 Yes 3 0.02 0.3 T08 Yes 4 0.02 0.3 T09 Yes 6 0.02 0.3 T10 Yes 8 0.02 0.3 PAGE 194 194 0 50 100 150 200 250 300 350 400 5 4 3 2 1 0 1 Distance from shoreline(m)h(m) Figure 41. Sea bottom elevation relative to m ean sea level versus crossshore distance on October 14, 1994 at DUCK, NC and selected crossshore locations of pressure sensors and/or current meters. 0 50 100 150 200 250 300 350 0 0.5 1 1.5 Longshore Current(m/s)Distance from shoreline(m) Figure 42. Predicted alongshore curre nt versus crossshore distance without roller model (points) and including the rollers with va rious wavefront slope values: 0.02 (dashdotted line), 0.05 (solid line), 0.1 (longdashed line), 0.2 (shortdashed line). 0 50 100 150 200 250 300 350 0 200 400 600 Radiation Stress(N/m)(a) 0 50 100 150 200 250 300 350 5 0 5 10 x 103 Distance form shoreline(m)Wave forcing(m/s2) (b) Figure 43. (a) Radiation stress x yS the wave radiation stress (shortdashed line), the roller radiation stress (longdashed line), and the to tal (solid line); (b) Wave forcing in crossshore direction, no roller (dashed line), including roller model with 0.05 (solid line). PAGE 195 195 0 50 100 150 200 250 300 350 0 0.5 1 1.5 Hrms(m)(a) 0 50 100 150 200 250 300 350 0 0.5 1 1.5 v(m/s) (b) 0 50 100 150 200 250 300 350 4 2 0 Distance from shoreline(m)h(m)(c) Figure 44. Modeldata comparison: (a) meas ured (circle) and modeled wave height rmsH ; (b) alongshore velocity, measured(circle), no roll er (dashed line), roller model with 0.05 (solid line); (c) bathymetry used. 0 50 100 150 200 250 300 350 6 4 2 0 2 4 6 x 103 Distance form shoreline(m)Momentum balance(m/s2) Figure 45. Alongshore momentum balance: wave forcing (solid line), bottom friction (longdashed line), lateral mixing (shortdash ed line), and residua l (dashdotted line). PAGE 196 196 0 50 100 150 200 250 300 350 0 0.5 1 1.5 Hrms(m) (a) 0 50 100 150 200 250 300 350 0 200 400 600 Sxy(N/m) (b) 0 50 100 150 200 250 300 350 0 0.5 1 v(m/s) Distance from shoreline(m) (c) Figure 46. Effects of frequenc y spreading on the modeled (a)rmsH (b) x yS and (c) alongshore velocity. Bathymetry, offshore rmsH mean period and direction are the same as in Figure 44. Single frequency (shortdashed line), /0.2mS (solid line), /0.4mS (longdashed line). 0 50 100 150 200 250 300 350 0 0.5 1 1.5 Hrms(m) (a) 0 50 100 150 200 250 300 350 0 200 400 600 Sxy(N/m) (b) 0 50 100 150 200 250 300 350 0 0.5 1 v(m/s) Distance from shoreline ( m ) (c) Figure 47. Effects of directi onal spreading on the modeled (a)rmsH (b) x yS and (c) alongshore velocity. Bathymetry, offshore rmsH mean period and direction are the same as in Figure 44. No directional spreading (shortdashed line), 00.211.5 Srad (solid line), 023 S (longdashed line). PAGE 197 197 0 50 100 150 200 250 300 1 0 1 2 x(m)v(m/s) Figure 48. The artifici al longshore current profile of the study case Bowen and Holman (1989), a constant water depth of 1 m is used. 0 0.01 0.02 0.03 0.04 0.05 0 1 2 3 4 x 103 im(s1) 0 0.01 0.02 0.03 0.04 0.05 0.3 0.32 0.34 0.36 0.38 cr(m/s)k(m1) Figure 49. Results of linear st ability calculations for the alongshore velocity profile (Figure 48) in terms of growth rate (top) and propag ation velocity (bottom) versus alongshore wave number. The thick dashed lines show wave number of the most unstable linear mode. PAGE 198 198 Figure 410. Spatial form of shear wave asso ciated with the most unstable linear mode (1 00.026 km ). (a) stream function pattern; (b) fl uctuation velocity field; (c) total velocity pattern, the shear wave has been scal ed so that its peak magnitude equals the peak of the mean alongshore current. 0 50 100 150 200 250 300 350 400 450 500 6 4 2 0 2 x(m)h(m) Figure 411. Bathymetry for a barred beach using the formulation given by Slinn et al. (2000), which is an approximate fit to topography measured at Duck, North Carolina, on October 11, 1990. PAGE 199 199 0 0.5 1 1.5 2 V(m/s) (a) 0 100 200 300 400 500 1 2 3 4 x 103 x(m) (m/s) (b) 0 2 4 6 x 103 im(s1) (c) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.5 1 1.5 k(m1)cr(m/s) (d) Figure 412. Linear stability calculations for waveinduced alongshore currents over the barred beach (Figure 411) for varying bed apparent roughness, 0.01akm (short dashed lines), 0.02akm (long dashed lines), 0.03akm (solid lines). (a) Mean velocity profiles, (b) bottom friction coefficient, (c) growth rate versus alongshore wavenumber, (d) propagation velocity versus alongshore wavenumber. 0 50 100 150 200 250 300 350 400 450 500 0 0.5 1 1.5 2 2.5 x(m)V(m/s) E01 E02 E03 E04 E05 E06 E07 E08 E09 E10 Figure 413. The predicted steady alongshore velocity () Vx for various experimental conditions (see Table 41). PAGE 200 200 0.2 0 0.2 u(m/s) 1 1.2 1.4 v(m/s) 0.2 0 0.2 u(m/s) 1 1.2 v(m/s) 0.2 0 0.2 u(m/s) 0 1 2 3 4 5 6 7 8 9 10 0.8 1 v(m/s)Time(hour) S=0oS=10oS=20o Figure 414. Time series of velocitiesu, v at ,(420,100) x ymm with 0.03akm and 0.5 M for different values of wave directional spreading S from top to bottom 0000,10,20 S (Run name E01, E02, E03). PAGE 201 201 Figure 415. Snapshots of the contour plots of vorticity at times separated by 0.25 h for 0.03akm and 0.5 M Top panels for 00 S bottom panels for 020 S Figure 416. The alongshoreaver aged alongshore velocity (,)vxt (m/s) (top) and the alongshoreaveraged perturbati on kinetic energy density 221 (,) 2 uvxt (m2s2) (bottom) as a function of the acrossshore coordinate x and time t for 0.03akm ,0.5 M and 00 S PAGE 202 202 0 100 200 300 400 500 0 0.5 1 PAGE 203 203 0 50 100 150 200 250 300 350 400 450 500 0 0.5 1 1.5 PAGE 204 204 0.4 0.2 0 0.2 0.4 u(m/s)S=0o 1 1.5 2 v(m/s) 0.4 0.2 0 0.2 0.4 u(m/s)S=10o 1 1.5 2 v(m/s) 0.4 0.2 0 0.2 0.4 u(m/s)S=20o 0 1 2 3 4 5 6 7 8 9 10 0.5 1 1.5 v(m/s)Time(hour) Figure 421. Time series of velocitiesu, v at ,(420,100) x ymm with 0.01akm and 0.5M for different values of wave directional spreading S from top to bottom 0000,10,20 S Figure 422. Snapshots of the contour plot of vorticity at times separated by 0.25 h for 0.01akm and 0.5 M (a) for 00 S (b) for 020 S PAGE 205 205 Figure 423. The alongshoreaver aged alongshore velocity (,)vxt (m/s) (top) and the alongshoreaveraged perturbatio n kinetic energy density 221 (,) 2 uvxt (m2s2) (bottom) as a function of the acrossshore coordinate x and time t for 0.01akm ,0.5 M and 00 S 0 100 200 300 400 500 0 0.5 1 1.5 PAGE 206 206 0 50 100 150 200 250 300 350 400 450 500 0 0.01 0.02 0.03 0.04 PAGE 207 207 0 100 200 300 400 500 0 0.5 1 1 5 PAGE 208 208 Figure 429. Snapshots of the contour plot of vorticity at 10ht for 0.01akm and 0 S and varying lateral mixing coefficient M 0 50 100 150 200 250 300 350 400 450 500 0 0.05 0.1 0.15 0.2 0.25 PAGE 209 209 0 50 100 150 200 250 300 350 400 450 500 0 0.01 0.02 0.03 0.04 PAGE 210 210 0.5 0 0.5 u(m/s)ka=0.03 1 1.5 v(m/s) 0.5 0 0.5 u(m/s)ka=0.01 1 1.5 2 v(m/s) 0.5 0 0.5 u(m/s)ka=0.005 0 1 2 3 4 5 6 7 8 9 10 0.5 1 1.5 2 2.5 v(m/s)Time(hour) 0.5 0 0.5 u(m/s)ka=0.02 1 1.5 v(m/s) Figure 433. Time series of velocitiesu, v at ,(420,100) x ymm with 0.5 M and 0oS for different values of bottom roughness from top to bottom 0.03,0.02,0.01,0.005akm (case E01, E09, E04 and E10 respectively). Figure 434. Contour plots of vorticity at times separated by 0.25 h for 0.005akm and 0.5 M (run E10). PAGE 211 211 0 100 200 300 400 500 0 0.5 1 1.5 2 2.5 PAGE 212 212 0 100 200 300 0 50 100 150 200 6 4 2 0 x(m) y(m) h(m) Figure 437. The rip channel bathymetry with 0.95eM Figure 438. Contour plot of the computed wave height and wave dire ction. (a) without the wavecurrent interaction; (b) with the wavecurrent interaction. Parameters are 00 and 01 Hm in deep water, 0.8M and 0.08akm PAGE 213 213 Figure 439. Computed velocity fields: (a) withou t the wavecurrent interaction; (b) with the wavecurrent interaction. Parameters are 00 and 01m H in deep water, 0.8M and 0.08akm 0 50 100 150 200 250 300 0.5 0 0.5 1 x(m)u(m/s) (a) 0 50 100 150 200 0.2 0 0.2 0.4 0.6 0.8 y(m)u(m/s) (b)without wc interaction with wc interaction Figure 440. (a) Acrossshore pr ofiles of uvelocity at 100 y m the center line of rip channel. (b) Longshore velocity pr ofiles of uvelocity at 100 x m 10 m seaward of the longshore bar. Parameters are 00 and 01 Hm in deep water, 0.8 M and 0.08akm PAGE 214 214 0 50 100 150 200 250 300 0.05 0 0.05 0.1 0.15 0.2 x(m) (m)y=100m,with interaction y=100m,no interaction y=0m, with interaction y=0m, no interaction Figure 441. Computed mean water level (wave setup) relative to the still water level. 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 0.6 x(m)u(m/s) ka=0.05m ka=0.08m ka=0.1m ka=0.2m Figure 442. The uvelocity profiles of rip cu rrent with different values of bottom roughness ak 0 0.05 0.1 0.15 0.2 0.25 0.3 0.4 0.2 0 0.2 0.4 0.6 velocity(m/s) 0 0.05 0.1 0.15 0.2 0.25 0.3 0 50 100 150 200 250 300 ka(m)offshore extent(m) umin umax vmax offshore extent Figure 443. Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent of rip current, with the bottom roughness ak PAGE 215 215 0.1 0 0.1 0.2 0.3 (m) (a) 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 0.6 x(m)u(m/s) (b) H0=0.6m H0=0.8m H0=1.0m H0=1.2m H0=1.5m H0=0.6m H0=0.8m H0=1.0m H0=1.2m H0=1.5m Figure 444. Effects of the incident deep wate r wave height on (a) Wa ve setup along the center line of rip channel; (b) Acrossshore profiles of the uvelocity at 100 y m 0 0.5 1 1.5 2 0.4 0.2 0 0.2 0.4 0.6 velocity(m/s)H0(m) 0 0.5 1 1.5 2 0 50 100 150 200 250 300 offshore extent(m) umin umax vmax offshore extent Figure 445. Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent of rip current, with incident wave heights. PAGE 216 216 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 y/LyRelative bar height Figure 446. Variation of bar/ri p channel in the longshore coordi nate with changing elliptic modulus e M Values are shown for 0,0.50,0.75,0.95,0.99eM with rip channels monotonically decreasing in width e M increases. The solid curve corresponds to 0.95eM which is the value for most of the test cases. 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 0.6 x(m)u(m/s) Me=0 Me=0.5 Me=0.75 Me=0.95 Me=0.99 Figure 447. The uvelocity profiles of rip current for different rip channel widths. 0 20 40 60 80 100 120 140 160 180 200 0.4 0.2 0 0.2 0.4 0.6 y(m)u(m/s) Me=0 Me=0.5 Me=0.75 Me=0.95 Me=0.99 Figure 448. The longshore pr ofiles of uvelocity at 100 x m different rip channel widths. PAGE 217 217 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.4 0.2 0 0.2 0.4 0.6 velocity(m/s)Relative ri p channel width 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0 50 100 150 200 250 300 Width(m), offshore extent(m) offshore extent width umin umax vmax Figure 449. Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv and the offshore extent and the width of rip current, with relati ve rip channel width. 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 0.6 x(m)u(m/s) S/ m=0 S/ m=0.1 S/ m=0.2 S/ m=0.3 Figure 450. The uvelocity profiles of rip curr ent for different values of the wave frequency bandwidth S 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 0.6 x(m)u(m/s) s=0o s=4o s=8o s=10o Figure 451. The uvelocity profiles of rip curren t for different values of the wave directional bandwidth S PAGE 218 218 0 2 4 6 8 10 12 0.4 0.2 0 0.2 0.4 0.6 S(o)velocity(m/s)umin umax vmax Figure 452. Variations of the maximum rip velocity maxu maximum onshore velocity minu and maximum velocity of the feeder current maxv with the wave directional bandwidth S Figure 453. Convergence tests for the test case E01. The timeaveraged vorticity fields for different grid sizes (a) and for di fferent time step sizes (b). PAGE 219 219 0 50 100 150 200 250 300 1 0.5 0 0.5 1 x(m)u(m/s) dt=2.5s, dx=5m dt=1.6s, dx=5m dx=6m, dt=2s dx=4m, dt=2s dx=5m, dt=2s Figure 454. Comparison of acrossshore veloci ty profiles along the rip channel center for different grid sizes and time step sizes for test case T01. Figure 455. Convergence tests for the test case E08. The timeaveraged vorticity fields for different grid sizes (a) and for different time step size (b). PAGE 220 220 0 50 100 150 200 250 300 0.4 0.2 0 0.2 0.4 x(m)u(m/s) dt=2.5s, dx=5m dt=1.6s, dx=5m dx=6m, dt=2s dx=4m, dt=2s dx=5m, dt=2s Figure 456. Comparison of acrossshore veloci ty profiles along the rip channel center for different grid sizes and time step sizes for test case T08. 0 0.5 1 0.5 0 0.5 (a) 0.5 0 0.5 1 1 0 1 (b) 0 0.5 1 0.5 0 0.5 (c) 0 5 10 15 0 0.5 1 Time(hour)u(m/s) 0 5 10 15 1 0 1 Time(hour)v(m/s)(d) Figure 457. Time series of velocities u (left) and v (right) at the location of ,100,302.5 x ymm (a) for T01, (b) for T02, (c) for T03. (d) for T04. PAGE 221 221 Figure 458. Snapshots of vorticit y contours at various times. (a) simulation results of T02; (b) simulation results of T04. Figure 459. Contour plots of tim eaveraged acrossshore mass flux hu timeaveraged acrossshore momentum flux huu and timeaveraged 221 2 TKEuuvv (a) for test case T01, (b) for T02, (c ) for T03, (d) for T04. PAGE 222 222 0 0.5 1 0.5 0 0.5 =1o 0 0.5 0.5 0 0.5 =2o 0 0.5 0.2 0 0.2 =3o 0 0.5 1 0 0.5 1 =4o 0 0.5 0 0.5 =6o 0 5 10 15 0.2 0 0.2 Time(hour)u ( m / s ) 0 5 10 15 0 0.5 Time(hour)v(m/s)=8o Figure 460. Time series of velocities u (left) and v (right) for simulation results with offshore incident angles 0000001,2,3,4,6,8 around the rip channel cen ter, at the monitored location of ,100,302.5 x ymm PAGE 223 223 Figure 461. Snapshots of contour plot of vorticit y at various times for the experiment cases of 02 04 and 08 PAGE 224 224 Figure 462. Timeaveraged vortici ty (top), acrossshore mass flux hu (middle), and timeaveraged turbulent kinetic energy 221 2 TKEuuvv (bottom) for simulation cases with various wave incident angles. PAGE 225 225 CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS 5.1 Conclusions This study sought to develop numerical mode ls to study water waves and waveinduced currents in nearshore area. A reduced wave spect ral model was developed based on the evolution equations of wave moments of wave energy density spectrum. This new momentsbased wave model is different from the traditional monochrom atic wave models or the full spectral wave models because the dependent variables to be solved are the moments of wave spectrum. The present wave model can serve as an intermediate ground in the field of nearshore wave modeling between the two traditional types of wave models. Generally, the wave model is computationally cheaper than full spectral wave models and mo re accurate than monochromaticbased wave models for wave spectra with finite wave fre quency and directional band widths, but relatively compact shapes with clear peak directions and frequencies. To study the waveinduced currents, a finite volume, steady twodimensional circ ulation model was first developed using the pressurecorrection solution method. Coupling of the wave model and the steady circulation model was utilized to simulate steady wave field and waveinduced current cases. The simulation results showed good agreement with the DUCK 94 fi eld data of waves and longshore currents. It has long been known unsteady flow motions in near shore currents can result from time invariant input waves. To be able to simulate the uns teady flow motions, a tim edependent circulation model extended from the steady solver was cr eated and coupled with the wave model through the radiation stress concept. Shear instabilitie s of the longshore curren ts over a barred beach (approximating the DUCK90 bathymetry) were si mulated and compared to the results of analytical linear stability analys is. Rip current instabilities were also studied using the coupled PAGE 226 226 wave/circulation model system. Overall, achi evements achieved through this study can be summarized as follows The monochromatic wave models including the mildslope wave models and the Boussinesq wave models, and the full spectral wave model (e.g. SWAN) were reviewed. The standard spectral wave action balance e quation was rederived using two different approaches. Starting with the standard wave action bala nce equation, the defi nition of the general wave moment (,)lEtx was introduced and the basic id ea to build a momentsbased wave model was presented. Based on the idea, th e evolution equation of the general wave moment using a weighting f unction with the form of 1 nime was derived. A simple fiveparameter wave model, which gives the RMS wave height, mean frequency, mean direction, and the frequency and directi onal bandwidths, was developed with some assumptions and approximations. One major assumption in this model is that the input wave spectra were assumed to be Gaussians haped distributions both in frequencies and directions and the input wave spectrum was assumed to reserve its spectral shape. Simulation of the evolution of an initial bum p in wave heights for the irregular wave conditions showed good agreement with analyt ical solutions for both the wave heights and the wave spectra when the wave spectra are relatively narrowbanded. As wave spectra become broader, the model still predicted the wave heights reasonably well, however, the mean frequency, and especially the frequency bandwidth were not predicted well because the shape of wave spectrum at some locations is far from a Gaussian distribution. Simulations of nonbreaking wave transformation on a sloping beach were conducted and compared with the analytical solutions. It was shown that the reduced model simulated not only the shoaling and refraction of irregular waves very well for this simple bathymetry, but also the limited evol ution of the frequency bandwidth and the strong evolution of the direc tional bandwidth. Numerical test s of the wave model on the relative complex bathymetry were performed, and reasonable results were obtained. The wave model was used to investigate wa ve transformation due to the presence of ambient currents. For monochromatic waves, currentinduced wave transformation was simulated very well by comparing with analytical solutions. The wave model was then used to investigate effects of wave fre quency and directional bandwidths upon evolution of waves with the existence of artificial ambient currents. Generally, the wave field is more sensitive to directional spreading th an frequency spreading. For waves with a considerable directional spread, the wave tr ansformation due to ambient currents can be very different from that of monochromatic waves. The radiation stresses com puted based on the monochromatic approximation are often used to drive the waveinduced currents in the nearshore circulation models. This approximation, however, may cause significant errors in the radiat ion stress components when the waves have finite frequency and directional bandwidths. Formulations of radiation stresses including th e effects of frequency and directional bandwidths were PAGE 227 227 derived, and a much more accurate estimate than the monochromatic approximation was provided. A 2D finite volume steady circulation m odel based on the nonlinear shallow water equations was developed using the pressure correction solution method. The model was built on the boundary fitted structured nonor thogonal computational grids and can be used to simulate flows in complex geomet ries. The models performance, convergence and accuracy were investigated using an arti ficial test with the analytical solutions available. It was proved that the model has 2nd order accuracy for both Cartesian and nonorthogonal computations grids. The model wa s further tested using the well known liddriven cavity flow case, and the computed results showed a good agreement with the benchmark solutions (Ghia et al., 1982). The wave model and the steady circulation m odel were coupled and used to simulate wave transformations and waveinduced longshore currents over artificial and field bathymetries. Simulation results of the coupl ed models were compared to the DUCK94 field data. Both the wave field and current field were well predicted by the models. A wave roller model including the effects of wave directional spreading was derived, simulation results showed that inclusion of surface rollers improved the prediction of the acrossshore velocity profile by shifting the maximum velocity shoreward and increasing the bar trough velocity over the barred beach Effects of wave frequency and directional spreading were studied by numerical experi ments on the DUCK94 bathymetry. Generally, the wave height distribution is relatively insensitive to wa ve frequency and directional spreading, however, the magnitude of the resulting longshore currents was considerably decreased by the presence of a finite directional bandwidth. The linear shear instabilities of longshore curr ents over a barred beach were studied by an analytical stability analysis model including the bottom friction. The analytical results suggested that the shear instabilities of alongshore currents give rise to alongshore propagating shear (vorticity) wave s, which is consistent with previous studies. However, studying the fully nonlinear dynamics of finiteamplitude shear waves requires numerical simulations. A timedependent circulation m odel extended from the steady version was created and coupled with the wa ve model. Using the coupled models, sets of simulations were performed to investigate the effects of bottom friction, lateral mixing, especially wave directional spreading on shear instabilities. Numerical results suggested that more unsteady, energetic instabilities are associated with a smaller bottom friction, a weaker lateral mixing and a stronger mean current Shear instabilities can cause significant damping of the mean current by inducing th e horizontal momentum mixing. A finite wave directional bandwidth leads to a less unstable flow, consequently, the shear instability character can be altered quantitiv ely and qualitatively due to the presence of the directional bandwidth. Steadystate waveinduced rip currents on an artificial rip channel bathymetry were simulated using the wave mode l along with the unsteady circulation model. Simulation results showed that the wavecurrent intera ction has significant impacts on both the wave field and the resulting circul ation system. Generally, the wc interaction resulted in a forcing effect opposite to that due to the t opography. As a result, the rip currents were PAGE 228 228 noticeably weaker and slightly broader and the offshoredirected rips were mostly restricted in the surf zone when including the wc interaction. In addition, the wavecurrent interaction is of such significance that the resulting rip currents are much less sensitive to a change of bottom friction coefficient. It was also observed that the strength and offshore extent of rip currents are subjec t to the incident wave height and the rip channel width. Limited spreadi ng in wave frequency and dire ction do not have significant effects on the resulting rip currents becau se topographyand currentinduced wave transformations are dominant factors in generating the rip currents. The instabilities of rip curre nts were investigated by nume rical experiments focusing on the following two aspects: effects of wavec urrent interaction on the development of rip current instabilities; and dependency of uns teady rip currents on the angle of wave incidence. The simulation results suggested that low frequency fluctuat ions in rip currents are restricted in the surf zone when usi ng the wave forcing with the wavecurrent interaction rather than moving back and fort h from the surf zone to offshore boundary as predicted without including wavecurrent interaction in the wave forcing. With no wavecurrent interaction, the timeaveraged acrossshore mass and momentum transported offshore through the rip channels may be overestimated, and offshore extent of transport processes could be overpredicted considerab ly too. It was also observed that the instabilities of rip currents ar e sensitive to the angle of wa ve incidence. The unsteady rip currents displayed a wide range of flow patte rns as the wave incident angle increases from shorenormal to 8. The energy content of the instabilities initially decreased as the angle continually from 0 to 3 degrees. Most en ergetic fluctuations with periods of about 2 hours were observed for wave in cident angle of 4. Circula tion cells began to disappear and the longshore current character became pr onounced as the angle further increased to 6 and 8. 5.2 Discussion and Recommendations for Future Work In this study, the reduced wave spectral model and the circulation model were developed and used to simulate the nearshore waves and cu rrents. Satisfactory simulation results have been achieved by using the newly developed models, especi ally in the aspect that the models allow for the direct investigation of the effects of th e wave frequency and di rectional spreading on nearshore waves and currents. However, there are uncertainties associated with various assumptions, approximations and simplifications in the models, especially in the wave model. There are also some approximations that have been made in the applications of the models. The following recommendations are made for further studies. Derivation of the evolution equations of wa ve moments, by integrating the wave action balance equation multiplied by some weighting f unctions over frequencies and directions, PAGE 229 229 does not need any prior assumptions. Howeve r, closure assumptions are necessary in order to evaluate the unsolved integrals in terms of wave moments. In this study, the input Gaussianshaped spectra have been a ssumed exclusively. More realistic wave frequency spectra such as JONSWAP (Hasselmann, 1973), TMA (Bonuws et al., 1985) and PM Spectrum (Pierson and Moskowitz, 1964) should be used in the future. One major difficulty in developing the mo mentsbased wave model is that more dependent wave moments will be involved th an the number of evolution equations. To define a close system, extra relations be tween different moments are required. For simplicity, the spectral shape of the input spectrum was assumed to reserve in this study. This, however, is not the case, especially fo r the directional spect rum when strong wave focusing and/or defocusing take place. One po ssible improvement of this is to use an external estimator for the extra wave moment s using the available wave moments, such as the maximum entropy technique. The effects of wave diffraction are not accoun ted for in the present wave model, like the typical full spectral wave model. To accomm odate the wave diffraction, one could use the phasedecoupled refractiondiffraction approximation initially proposed for the SWAN wave model by Holthuijsen et al. (2003) For a complete simulation of nearshore waves, wind input, energy dissipation due to whitecapping, and wavewave interactions should also be included. The present wave model is limited to wave spectra with single peak frequency and direction. One possible improvement is to de velop a more sophisticated model that can deal with wave spectra with multipeaks in both frequency and direction. In both the steady and unsteady circulation models, the SIMPLE algorithm was used as the solution method. One could use enhanced versions of SIMPLE to make the model more efficient and more robust. For the unst eady circulation model, a more sophisticated time differencing scheme rather than the implicit Euler method should be used for accuracy and stability. Although the circulation models were deve loped on the boundaryf itted nonorthogonal structured grid system, tests on realistic coastal topographies us ing nonorthogonal grids were not conducted. More verifications of the model on computation domain with complex external and/or internal boundaries s hould be performed before any application. In the present circulation mode ls, a stationary computation grid was used. To model the shoreline movement, capability of deali ng with the dynamic boundary should be built into the models. One simple way to do th at is following the dryingwetting scheme developed by Casulli and Cheng (1992). In simulating the nonlinear shea r instabilities of longshore cu rrents, the effects of wavecurrent interaction on the devel opment of were not included. It has been found that wavecurrent interaction can significantly alter th e finite amplitude behavior of the shear instabilities (zkanHaller a nd Li, 2003). Therefore, for a complete study, wavecurrent interaction should be accounted for. PAGE 230 230 In simulating the rip current instabilities, th e wavecurrent interac tion was only partially represented by using the time invariant wa ve forcing including the wavecurrent interaction. In reality temporal changes in the wave field will be induced due to the changes of the current field, and so does the wave forcing. Therefore, the full representation of wavecurrent interaction is necessary for a complete study. The present models have been mostly used in artificial or idealized systems. More simulations of field data should be conducted in future work. PAGE 231 231 APPENDIX A DERIVATION OF MOMENTS OF THE GAUSSIANSHAPED S PECTRA With the assumption of the Gaussianshaped wave spectra 22 0,01111 (,,;,) exp exp 22 22mmFxytE SS SS (A1) According to its definition, the general mo ment of the wave spectrum density (A1) ,nm E can be written as 2211 // 22 ,0,0 011 ee 22mmSS nim nmEE ded SS (A2) First, look at the integral 21 / 2 0emS nd Let mt S which gives 2 211 / 22 0/eem mSt n n m SdStSdt (A3) In the case /mS is not too small, the integral a bove can be approximated as follows 2 2211 1 / 22 2 0/eeem mStt nn n mm SdStSdtStSdt (A4) Because the exponential function 21 2et decreases quickly to 0 as t increases to about 3. With the narrowbanded assumption, /mS in fact, is always larger th an 3. Therefore, approximation from Equation (A3) to (A4) is reasonable. According to the binomial theorem, we have 0! !()!n n nkknk mm kn St St knk (A5) Substituted into Equation (A4) gives PAGE 232 232 2 21 1 / 1 2 2 0 0 1/2 1 0! ee !()! 2 1/2, ()=0,2,4... = !()! 0, mn St nn k k n k m k n nk nkk m kn dStd t knk n Sn kn k knk otherwise (A6) Where, 1/2 nk is the gamma function. Similar process is applied to the integral 21 / 2emSd Substituting mt S into the integral and assuming that both /mS and /mS are not small so that the approximation followed can be made 2 2211 1 / / () () 22 2 /em m mm mS Stt imSt imSt im Seed eeSdteSdt (A7) And, 2 221 ()/2 2e2cossinmt imSt mS mmedtemim (A8) So, 2 221 / /2 22cossinmS mS im mmeedSemim (A9) Substituting equations (A6) and (A9) into (A2), the moment ,nm E is finally given by 220,0 /2 1/2 1 0cossin* 2 2 1/2, ()=0,2,4... !()! 0, otherwisemS nmmm n nk nkk m kE Eemim n Sn kn k knk (A10) PAGE 233 233 As mentioned above, /mS /mS and /mS were assumed to be no less than 3 in order to take advantage of the defin ite integral containing e xponential functions To maintain the assumption valid, following limitations are forced: /0.25mS /5mS and /5mS (A11) After deriving the general form of wave moment nmE it is straightforward to give the moments appearing in th e current wave model: 1,00,0mEE (A12) 22 2,00,0 mEES (A13) 2/2 0,10,0cossinS mmEEei (A14) 2/2 0,10,0cossinS mmEEei (A15) 22 0,20,0cos2sin2S mmEEei (A16) 29/2 0,30,0cos3sin3S mmEEei (A17) PAGE 234 234 APPENDIX B EVALUATION OF THE INTEGRALS IN GOVERNING EQUATIONS This appendix gives the deriva tions of integrals in the equation of the general m oment nmE Integrals in the equations of the reduced spectral wave model are straightforward to obtain from the three general integrals. Recall Equation (247) 00 1 0 1 0nim nim nm nimEe F d d i m c e F d d t nceFdd gcU (B1) 1. First integral 0 0 0,0 0 (1)(1)(1)(1) 0,0 01 cos,sin ,( ) 22 1 ,( ) 2nim nim gn m iiii nim gn m ni m i mi m i m g nmIe F d d ce F d du v E eeee ceEMDdduvE i cMdeeieeEDduvE gcu 0,10,10,10,1 01 ,, 2n mmmmg nmEEiEEcMduvE (B2) Using the Taylor expansion approximation of g c we have: 0 2 2 2 00 0 1 000 2 21 2 2 0001 2 1 2 2m m m m m mn g gg nnn gmm g nnn g m g nnn mmcGd cc cGd Gd Gd c cGdGdGd c GdGdGd (B3) PAGE 235 235 Inserting the above equation back in to Equation (B2), we arrive at 2 2 ,1,1,1,1 2 2 1,11,1 1,11,1 2 2,12,111 1, 22 1 2 1 2m m m m mgg nmnmnmnmgm m gg nmnmnmnm m nmnmcc IEEiEEc cc EEiEE EEiE 2 2,12,1, 21 2mg nmnm nmc Eu v E (B4) 2. Second integral Similar to the first integral, the second integral can be given by 0 0,10,1 0,10,1 0 ,,2,2 ,2,2,,2,22 11 () 22 11 22 444 1 2mnim n gm m m m nmnmnm nmnm nmnmnm mIceFdd hh ccGdiEEEE hxy ui u vv EEE EEEEE yx yx cc cc h 2 2 ,1,1,1,1 2 2 1,11,1 1,1,1 2 2 2,1 21 2 1 2 1 4m m m m mmn m n m n m n m mn m n m n m n m nmncch h iEEEE xy cc cch h iEEEE hxy cch iEE hx 2,1 2,12,1 ,,2,2 ,2,2 ,,2,21 2 4 1 2 44mn mn mn mn mn m nmnm nmnmnmhu EEEEE yy iuvv EEEEE xyx (B5) PAGE 236 236 3. Third integral 1 0 0, 0,0,20,2 00 0,20,2 0,0,20,23 11 () ()2 24 2 11mmnim gg m nn mmm mmmmm gg mIceFdd cc E u hGdGdEEE hccx uv v iEEEEE yx y cc h hcc U U 2 2 2 22 1,2, 22 2 2 ,,2, 21 22 1 2 11 2 24m mm m mmmg mn m ggg mn mn m ggg m m nmnmnmc E c ccc EE ccc ccc u EEE cccx 2 2 ,2,2 ,,2,2 2 1,1,21,2 1,21,2 1,1,21,2 1 2 4 2mmgg nmnm nmnmnm m nmnmnm nmnm nmnmnmcc uvv iEEEEE yxycc uu v EEEiEE xy x v EEE y 2 22 2 2 2 2 2 2,22,2 2,2,22,21 2 8 2g nmnmnm nmnm nmnmnmc u EEE cx uvv iEEEEE yxy (B6) For convenience, let 2 2 21 1 2m m mgg gm mcc Ac 2 22m mgg mcc A 2 21 3 2mgc A 2 2 211 1 2m m mmmcc cc Bcc h 2 21 2m mmcc cc B h 2 21 3 2mcc B h PAGE 237 237 2 2 21 1 2mmmggg mmccc C ccc 2 22mmgg mcc C cc 2 21 3 2mgc C c ,, 1 1, 1 1 nm nmnm nmnmhh XiEEEE xy ,,2,2 ,2,2 ,,2,21 22 4nm nmnmnm nmnm nmnmnmuu vv YEEEiEEEEE xy xy ,,2,2 ,2,2 ,,2,211 22 444nm nmnmnm nmnm nmnmnmui u vv ZEEEEEEEE yx yx Then, we obtain much shorter expr essions for the three integrals ,1,1,1,1 1,11,11,11,1 2,12,1 2,12,1 ,12 1, 22 3 2nmnmnmnm nmnmnmnm nmnmnmnm nmAA IEEiEEEEiEE A EEiEEuvE (B7) ,1 ,2 ,,1 2123 2nm nmnmnm I BXBXBXZ (B8) ,,1 ,2 ,1 ,2 ,1 310 523 123nmnmnmnm nmnmnmIhCEECECE h CYCYCY U (B9) After lots of algebra, the governing equations for the reduced wave model are given by 2 2 0,0 0,10,1 0,0 2 2 0,0 2 0,00,0 21 ,, 2 11 0 22m m mmg g RI ggc EcSEEEuv t cc E ShEYh cchh UU (B10) 2 22 1,0 0,10,1 1,0 2 2 2 1,0 2 1,00,0 21 ,, 2 11 2 0 2m m m mm mgg gm m RI ggg m mcc EcSSEEEuv t ccc E S ShEYh ccchh UU (B11) PAGE 238 238 2 24 2,00,10,12,0 2 2 0,10,1 2,0 2 2 23 ,, 2 3 2, 2 9 32 2m m m mmmg g RI g m RI ggg mc EcSSEEEuv t c SEE hE h ccc S ccc U 2 2,00,01 0 hEYS h U (B12) 2 0,10,10,10,1 2 0, 0, 2 2 0, 2 0,0, 2 21 ,, 222 11 22 2 1 2m m m m mg mmmm mg m m mm gc EEEE EcSiuvE t E cc imccSXZh hh c c U 2 0,0, 21 0mg mmc ShEY ch U (B13) in which, 0,0 0,00,2 0,2 0,00,21 2RIIuu vv YEEEEE xy xy 0, 0,10,1 0,10,1 mmmmmhh XiEEEE xy ,,2,2 ,2,2 ,,2,21 22 4nm nmnmnm nmnm nmnmnmuu vv YEEEiEEEEE xy xy 0, 0,0,20,2 0,20,2 0,0,20,211 22 444mm m mm mm m mui u vv ZEEEEEEEE yx yx PAGE 239 239 APPENDIX C DISCRETIZATION OF TERMS IN TH E SHALLOW WATE R EQUATIONS Mass Flux The midpoint rule of approximation is used to compute the mass fl uxes. The method is demonstrated by applying to the east face of a t ypical 2D control volume (see figure 31), and it is the same to apply to other faces. As st ated previously, on the outer iteration level m all nonlinear terms are approximated by a product of o ld (from the preceding outer iteration) and a new value. Thus the mass flux through each CV face is evaluated using the latest available velocity field: 1em ee e SmhvndShvnS (C1) /eneseneseenyyixxjS (C2) Where en is the outward unit normal vector at the face e and eS is the surface area. Unless stated otherwise, all variables in this section are about the m th outer iteration. Convective Flux The convective flux of iu momentum equation through e face of CV is: 1 ,,eecm ie i i eie SSFhuvndShuvndSmu (C3) The CV face value of iu used in this expression is not necessary to be the same as that used to calculate the mass flux. Linear in terpolation is the easiest se cond order approximation, known as the central difference scheme (CDS). Some iteration methods fail to converge when applied to the algebr aic equation system derived from central difference approximations of convective fluxes since the matrices may not be diagonally dominant. For large Peclet numbe r flow case, central difference approximation may lead to physically impossible solutions. An approach called deferred correction is very PAGE 240 240 useful to ensure the convergence of iterative solver when high order approximation is used. Deferred correction approach combines an explicit high order approximati on (here is CDS) and an implicit lower order approximation (here is the Upwind Differencing Scheme, UDS) in the following way: 1 ,,,, m cUDS CDSUDS ieeieeieieFmumuu (C4) The term in parenthesis on the right hand side of the above equation is evaluated using values from previous iteration. As convergence appro aching, the UDS contributions cancel out, the solution is still the second order approximation (CDS ). So this procedure always makes iterative solver used more robust and efficient, meanwhile, keeps the high order accuracy. Diffusive Flux The lateral momentum mixing terms given by Equation (314) involve second derivatives. Applying the finite volume method to the integral form of the momentum equations, the lateral mixing terms become diffusive fluxes. For convenience, we use the general term ihu for demonstration, and it is straightforward to ob tain the specific expression for each second derivatives from the general form. The midpoint rule approximation applied to the integrated diffusive flux gives: ,ed ieiie e SFhundShunS (C5) The spatial gradient of ihuat the cell face center can be define d in terms of the derivatives with respect to either global Cartesian coordinates x y or local orthogonal coordinates, nt: iiii ihuhuhuhu hu x ynt i j nt (C6) PAGE 241 241 Where, n and t denote respectively the nor mal and tangential coordina te direction to the cell face. A better way to approximate the diffusive flux is that the gradient evaluated with respect to the local orthogonal coordinates, nt, since only the derivative in the ndirection contributes to the diffusive flux: d i ie e e ehu FS n (C7) For the nonorthogonal grid, the line norm to th e cell face e may not pass through the CV center nodes P and E (see figure 33). As a result, the derivative of ihu with respect to the direction normal to cell face at cell center e can not be approxi mated using the values at nodes P and E. Muzaferija (1994) s uggested using the deri vative with respect to the coordinate t, along the line connecting nodes P and E, as an implicit flux approximation, if the t direction is nearly orthogonal to th e cell face, Equation (C7) can be approximated as ii d i EP ie e e ee e PEhuhu hu FSS L t (C8) where P ELdenotes the distance from P to E. Muzaferija (1994) suggested a deferred correct ion part to the diffusive flux, and then Equation (C8) becomes 1 m d ii i ie e e ee ee ehu huhu FSS n tt (C9) in which, 1gradm i i e ehu hun n (C10) PAGE 242 242 The gradient at cell face center e is calculated by interpolating the gradients at CV centers P and E. It should also be noted that if n t the last term in Equation (C9) becomes zero. When the nonorthogonality is not severe, th e method discussed above is second order accurate. The approximations discussed above as sumed that the line connecting nodes P and E passes through the cell face center e. In irregular grid case, however, the line connecting P and E may not pass through e, and the approximations for i ehuand i ehu tused above are second order accurate at a point e' rather than at the face center e. If e' is close to the corners se or ne, the approximation reduces to first order accurate. To give an effective cure, Ferziger and Peri c (2002) suggest using values at auxiliary points P', E' and e'. P' and E' lie at the intersections of the cell face normal n and straight lines connecting cell face centers n and s of each CV. e' lies at the intersection of line connecting P and E and line connecting CV corners, see figure 33. Then the modified diffusive flux is written as 1 m d iii ieee ee eeehu huhu FSS n tt (C11) where, ii i EP e PEhuhu hu L t (C12) ii i EP e PEhuhu hu nL (C13) P EL stands for the distance between P' and E'. The values of i Ehu and i Phucan be calculated either by bilinear interpolation or by using the gradient at CV center: iiiPP PPPhuhugradhurr (C14) PAGE 243 243 iiiEE EEEhuhugradhu rr (C15) Bottom Friction Term In this circulation model, the bottom fric tion term has the following general form b xu b yv (C16) Where, is the bottom friction coefficient and the formulation is given by Equation (36). For the circulation mode l, the values of at every grid points are known as a priori. The bottom friction terms are linear expression s in terms of velocities, so they can be easily treat implicitly. In disc retized momentum equations, implicit treatment of bottom friction term increases the diagonal dominance of coefficient matrix and thus usually has a positive effect on many iterative solution methods. The midpoint rule approximation applied to the volume integral of the bottom friction term gives: ii Pudu (C17) where is the volume of the interested CV. Wave Forcing Term As for the circulation model, the wave forci ng is an input quantity as a surface force term in the momentum equations. This term only cont ributes to the source terms of the discretized momentum equations. The second order midpoint rule approximation is applied: In the xmomentum equation: 1xx xx y SSiSjndS (C18) In the ymomentum equation: 1yy xy y SSiSjndS (C19) PAGE 244 244Pressure Term The pressure term in momentum equations can be treated either as surface forces (conservative approach) or as body forces (nonconservative approach). In the first case (as surface forces): 11, ,,, mm P l S l lhh SPPindSPndSlewns (C20) In the second case (as body forces): 11 mm P d ii PhPhP SP d xx (C21) Equation (C21) is essentially equiva lent to (C20) if the derivative 1 m i P hP x is calculated using Gauss theorem. PAGE 245 245 APPENDIX D TWO NUMERICAL ISSUES FOR THE CIRCU LATION MODEL: UNDERR ELAXATION TECHNIQUE AND MOMENTUM INTERPOLATION METHOD UnderRelaxation Technique For the algebraic equa tion for the variable on outer iteration level m P mm PPll lAAQ (D1) in practice, only a fraction of th e wouldbe difference is allowed to be changed for stabilization consideration: 11 mmnewm, 01 (D2) where is called the underrelaxation factor. new is the solution of Equation (D1) and can be written as m P ll new l P PQA A (D3) Substituting Equation (D3) into (D2), it gives: *1 P1P PQ A mmm P PllPP lA AQA (D4) where P A and P Q are modified main diagonal elements and source terms respectively. When convergence approaches, the terms involving will cancel out and the solution of the original problem will be obtained. Compared to the explic it application of the underrelaxation technique according to Equation (D2), this treatment has been shown to have positive effects on many iterative solution methods since the diagonal domi nance of the coefficient matrix is increased (Patankar, 1980). PAGE 246 246 The optimum underrelaxation factors, which require least iterations to converge, may vary from case to case. A good strategy is to us e a small value at the be ginning of iterations and increase it towards unity as convergence is appr oached. For solving th e incompressible NavierStokes equations based on the pressurecorrec tion method, the two momentum equations usually use one underrelaxation factor u and the pressure correction e quation should be assigned its own factor P Pantanker (1980) suggested using 0.5u and 0.8P in the SIMPLE algorithm. Ferziger (2002) proposed that 1.1 P u gives the best results for the liddriven cavity flow problem. Momentum Interpolation Method The collocated grids have advantages over stag gered grids in respect of reducing required storage memory and computational costs, and in respect of dealing with complex geometries. Nonetheless, the main drawback of the collo cated grids is the occurrence of nonphysical oscillations of pressure and/or velocities, causing severe numerical instability in some cases. The spurious oscillations are mainly due to the d ecoupling between the pressu re and the velocities. To address the problem, Rhie and Chow (1983) first proposed to employ a special interpolation practice, the momentum interpolation method (MIM ), instead linear interpolation for calculating the velocities at the faces of CVs. The success of this method reported by Rhie and Chow (1983) has stimulated more research efforts on this su bject. Following the same basic idea, Peric(1985) and Majumdar (1986) have worked on subseque nt refinements. Later, Majumdar (1988) observed that the converged results are underre laxationdependent during trial runs. Such dependence on underrelaxation factor is certainly an undesirable f eature of the MIM. To remove the dependency, Majumdar (1988) analyzed th e reason for the dependency and proposed a modified momentum interpolati on method illustrated for onedimensional flow. In the following, PAGE 247 247 we will give a brief review of the momentum in terpolation method and then present the modified version for 2D steady flow based on the same principle of method proposed by Majumdar (1988). For simplicity, they will be demonstrated on Cartesian grids. Taking out the pressure term from the source term, the algebraic equations for discretized u momentum equation can be written as P Pllp lAuAuSuSu (D5) where pSu is the source term due to pressure gradient. And, P peewwh SupSpS (D6) From equations (D5)(D6), we can get ll l eeww P PP P PP PPSuAu pSpS h u AA (D7) The right hand side of Equation (D 7) consists of two parts, the influences of velocities at the surrounding CV centers denoted by the first term, and the influence of the pressure gradient represented by the second term Following the expression, the u at the CV center E can be written as ll l eeww EE E E PP EESuAu pSpS h u AA (D8) Mimicking Equation (D7) and (D8), we can expr ess the interface velocity at cell face e as ll l EP eee e PP eeSuAu p p hS u AA (D9) This equation gives the original mome ntum interpolation method (OMIM). PAGE 248 248 When the underrelaxation is in corporated for stability pur pose, Equation (D9) becomes to: 00' 11ul l l EP euee e ueEeP PP eeSuAu pp hS uf u f u AA (D10) where 0 P u and 0 Eu are values at the prev ious iteration level. /ePePE f LL is the linear interpolation factor. Su is the modified Su (see Equation (D6)). Equa tion (D10) is still the formulation of the OMIM, and the only differen ce from Equation (D9) is the incorporation of underrelaxation. Majumdar (1988) suggested that a slight change of Equation (D10) could eliminate the underrelaxation dependent 0' 1ul l l EP euee e ue PP eeSuAu pp hS uu AA (D11) This expression for cell face velocities can ach ieve a unique solution that is underrelaxation factor independent. For a better understanding of Equation (D11), we will take a close look at it. It is easy to show that: 0' 1ul l up l P uP PPSuAu Su uu AA (D12) Using Equation (D12), we can get 00' (1)1 (1) (1)ul l l e eEePueEeP P e eeww eeww EEPP ue e PP EPSuAu fufufufu A pSpS pSpS hh ff AA (D13) Substituting Equation (D13) into (D11), we finally arrive at: PAGE 249 249 0001 1 11ue ew wue ew w eE eP PP EP eeEeP u EP ee ueeEeP P epSpS pSpS fh fh AA ufufu pp hS ufufu A (D14) The expression of cell face velocities by Equatio n (D14) is composed of two parts. The part from linear interpolation is given by the firs term on the right of side of the equation. The part from momentum interpolation is denoted by all terms in the cu rly braces. The last term in the braces can be called MMIM correction term since it is intr oduced by using the modified momentum method. Without this te rm Equation (D14) becomes to be equivalent to the MIM. Obviously, the part from momentum interpolation can be regarded as a correction to the linear interpolation. The first three terms in the braces have the function of smooth ing the pressure field, and then the unrealistic pressure oscillation could be removed. As iterations approach to convergence, 0 P puu 0 EEuu and 0 eeuu Pulling out the MMIM correction term and reorganizing, Equation (D14) can be rewritten as a formulation th at every term is containing the multiplier u Divided by u it yields an expression of eu that is not a function of u Therefore, the dependency of underrelaxation fa ctor is certainly elim inated. Verification by numerical tests is shown in Section 3.5.2. PAGE 250 250 LIST OF REFERENCES Allen, J. S., P. A. Newbe rger, and R. A. Ho lman (1996), Nonlinear shear instabilities of alongshore currents on plane beaches, J. Fluid Mech., 310, 181213. Battjes, J. A. (1972), Setup due to irregular waves, in Proc. 13th Int. Coastal Eng. Conf., ASCE, 19932004. Battjes, J. A., and J. P. F. M. Janseen (1978) Energy loss and setup due to breaking of random waves, in Proc. 16th Int. C onf. Coastal Engineering, ASCE, 569587. Berkhoff, J. C. W. (1972), Computation of combined refraction diffraction, Proc., 22nd Int. Conf. on Coastal Engrg., Vancouver, Canada, 471490. Bouws, E., H. Gunther, W. Rosenthal, and C. L. Vincent (1985), Similarity of the wind wave spectrum in finite depth water, part I Spectral form, J. Geophys. Res., 90(C1), 975986. Booij, N. (1981), Gravity waves on water with nonuniform depth and current, Doctoral dissertation, Technical University of Delft, The Netherlands, 131 pp. Booij, N., R. C. Ris, L. H. Holthuijsen (1999), A thirdgeneration wave m odel for coastal regions, Part I, model description and validation, J. Geophys. Res. 104(C4), 76497666, doi: 10.1029/98JC02622. Bookm an, C. A., T. J. Culliton, and M. A. Warren (1999), Trends in U.S. coastal regions, 19701988, NOAA. Bowen, A. J. (1969), Rip currents: 1. Theoretical investigations, J. Geophys. Res., 74, 54675478. Bowen, A. J. and R. A. Holman (1989), Shear instabilities of the mean longshore current 1. Theory. J. Geophys. Res., 94, 18023. Brander, R. W., and A. D. Short (2001), Flow kinematics of low energy rip current systems, J. Coastal Res., 17(2), 468481. Casulli, V. and R. T. Cheng (1992), Semiimplicit finite difference methods for threedimensional shallow water flow, Int. J. for Numerical Methods in Fluids, 15, 629648. Chawla, A., H. T. zkan, and J. T. Kirby ( 1998), Spectral model for wave transformation and breaking over irregular bathymetry, J. Waterway Port Coastal Ocean Eng., 124(4), 189198. Chen, Q., J. T. Kirby, R. A. Dalrymple, A. B. Kennedy, and M. C. Haller (1999), Boussinesq modelling of a rip current system, J. Geophys. Res. 104(9), 2061720637. PAGE 251 251 Church, J. C. and E. B. Thornton (1993), Effects of breaking wave induced turbulence within a longshore current model, Coastal Eng., 20, 128. Dalrymple, R. A. (1975), A mechanism for rip currents generation on an open coast, J. Geophys. Res., 80, 34853487. Dean, R. G., and R. A. Dalrymple (1998), Water Wave Mechanics for Engineers and Scientists, World Scientific, N. J., 353 p. Dean, R. G., and R. A. Dalrymple (2002), Coastal Processes with Engineering Applications, Cambridge University Press, Cambridge UK, 475 p. Deigaard, R. (1993), A note on the three dimensi onal shear stress distribution in a surfzone, Coastal Eng., 20, 157171. Dodd, N., J., OltmanShay, and E. B. Thornton (1992), Shear instabiliti es in the longshore current: A comparison of observation and theory, J. Phys. Ocean., 22, 6282. Duncan, J.H. (1981), An Experimental Inves tigation of Breaking Waves Produced by a Towed Hydrofoil, Proc. R. Soc. Lond., A, 377, 331348. Ebersole, B. A., and R. A. Dalrymple (1980), Numerical modelling of nearshore circulation, paper presented at Intern ational Conference on Coasta l Engineering, ASCE, Sydney. Eldeberky, Y., and J. A. Battjes (1996), Spectra l modeling of wave breaking: Application to Boussinesq equations, J. Geophys. Res., 101(C1), 12531264 Falqus, A., V. Iranzo, and M. Caballeria (1994), Shear instability of longshore currents: effects of dissipation and nonlinearity, Proc. 24th Intl. Conf. Coastal Engng, ASCE, pp 19831997 Feddersen, F. (2004), Effect of wave directional spread on the radiation stress: Comparing theory and observations, Coastal Eng., 51, 473481. Feddersen, F., R.T. Guza, S. Elgar, and T.H.C. Herbers (1998), Alongshore m omentum balances in the nearshore, J. Geophys. Res. 103, 15667. Ferziger, J. H. and M. Peric (2002), Computational Methods for Fluid Dynamics, Springer, New York. Fletcher, C. A. J. (1991), Computational techniques for fluid dynamics. vol., 1. Springer, Berlin Fredsoe J., and R. Deigaard (1992), Mechanics of Coastal Sediment Transport. Advanced Series on Ocean Engineering, Vol. 3. World Scientific, 369 pp. PAGE 252 252 Ghia, U., K. N. Ghia, and C. T. Shin (1982), HighRe solutions for incompressible flow using the NavierStokes equations and a multigrid method, J. Comput. Phys., 48, 387411. Gobbi, M. F., J. T. Kirby, and G. Wei (2000), A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O(kh)4, J. Fluid Mech. 405, 181210. Goda, Y. (1985), Random seas and design of maritime structures, Univ. of Tokyo Press, Tokyo, Japan. Grassa, J. M. (1990), Directional random waves propagation on beaches, Proc., 22nd Int. Conf. on Coastal Eng., Delft, The Netherlands, 798811. Haas, K. A., and I. A. Svendsen (2000), Threedimensional modeling of rip current systems. Research Report CACR0006, Center for Appl ied Coastal Research, University of Delaware. Haas, K. A., I.A. Svendsen and M.C. Haller (199 8), Numerical modeling of nearshore circulation on a barred beach with rip channels. In Proc. 26th Int. Conf. on Coast. Eng., ASCE, 801814. Haas, K. A., I. A. Svendsen, M. C. Haller and Q. Zhao (2003), Quasithreedimensional modeling of rip current systems. J. Geophys. Res., 108 (C7), 3217, doi:10.1029/2002JC001355. Haller, M. C., and R. A. Dalrympl e (2001), Rip current instabilities, J. Fluid Mech., 433, 161192. Hasselmann, K., T.P. Barnett, E. Bouws, H. Ca rlson, D.E. Cartwright, K. Enke, J.A. Ewing, H. Gienapp, D.E. Hasselmann, P. Kruseman, A. Meerburg, P. Mller, D.J. Olbers, K. Richter, W. Sell and H. Walden (1973), M easurements of windwa ve growth and swell decay during the Joint North Sea Wave Pr oject (JONSWAP), Dtsch. Hydrogr. Z. Suppl., 12, A8. Hedges, T.S. (1987), An approximate model fo r nonlinear dispersion in monochromatic wave propagation models, Coastal Eng., 11, 8789. Holthuijsen, L. H. (2007), Waves in oceanic and coastal waters, Cambridge University Press, New York, 387 pp. Holthuijsen, L.H., N. Booij, and R.C. Ris (1993) A spectral wave m odel for the coastal zone, Proc. 2nd Int. Symposium on Ocean Wave Measurement and Analysis. ASCE., New York, pp. 630. Holthuijsen, L.H., A. Herman, and N. Booij N. (2003), Phasedecoupled refractiondif fraction for spectral wave models, Coastal Eng., 49, 291305. PAGE 253 253 Houwman, K. T., and P. Hoekstra (1998), Tidal ellipses in the nearshore zone (3 to 10 m); modelling and observations, In Proc. 26th Int. Conf. on Coast. Eng., pp. 773786, ASCE, New York. Huntley, D. A., M. D. Hendry, J. Haines, and B. Greenidge (1988), Waves and rip currents on a Caribbean pocket beach, Jamaica, J. coastal Res., 4, 6979. Izumiya, T., and K. Horikawa (1987), On th e transformation of directional waves under combined refraction and diffraction, Coast. Eng. in Japan, Tokyo, Japan, 30, 4965. Jonsson, I. G., O. Skovgaard, and T. S. Jac obsen (1974). Computation of longshore currents. Proc. 14th conf. on Coastal Eng., ASCE, 699714. Kennedy, A.B., and D. Thomas (2004), Drifter meas urements in a laboratory rip current, J. Geophys. Res., 109, C08005, doi:10.1029/2003JC001927. Kennedy, A. B., and Y. Zhang (2008), The stabi lity of wavedriven rip current circulation, J. Geophys. Res., 113, C03031, doi:10.1029/2006JC003814. Kennedy, A. B., Y. Zhang, and K. A. Haas (2008 ), Rip currents with varying gap widths, J. Waterway Port Coastal Ocean Eng., ASCE, 134(1), 6165. Kirby, J.T. and R.A. Dalrymple (1984), Verification of a parabolic equation for propagation of weaklynonlinear waves, Coastal Eng., 8, 219232. Kirby, J. T. (1986), Rational approximations in their parabolic equation method for water waves, Coastal Eng., 10, 355378. Kirby, J.T. and R.A. Dalrymple (1986), An a pproximate model for nonlinear dispersion in monochromatic wave propagation models, Coastal Eng., 9, 545561. Kirby, J. T. and R. A. Dalrymple (1994), Combined Refraction/Diffraction Model REF/DIF 1, Version 2.5. Documentation and User's Manual, Research Report No. CACR9422, Center for Applied Coastal Research, Depart ment of Civil Engineering, University of Delaware, Newark. Kirby, J. T., Wei, G., Chen, Q., Kennedy, A. B. and Dalrymple, R. A., FUNWAVE 1.0. Fully nonlinear Boussinesq wave model. Documentation and user's manual, Report CACR9806, Center for Applied Coastal Research, Departm ent of Civil and Environmental Engineering, University of Delaware. Lascody, R. L. (1998), East central Florida rip current program, Natl. Weather Dig., 22(2), 2530. Li, R., Yan, Y. and Cao, H. (2003), Nonlinear dispersion relation in wave transformation, China Ocean Eng., 17, 117122. PAGE 254 254 Lim, S. H. L., and E. S. Chan (2003), Surf Zone Wave Modeling Imp rovement of Dissipation by DepthInduced Wave Breaking in Swan, the Tenth OMISAR Workshop on Ocean Models, Hanoi, Vietnam. Lippmann, T. C., A. H. Brookins, and E. B. Thornton (1996), Wave energy transformation on natural profiles, Coastal Eng, 27, 120. Lippmann, T. C., T. H. C. Herbers, and E. B. Thornton (1999), Gravity and shear wave contributions to nearshore infragravity motions, J. Phys. Ocean., 29(2), 231239. Liu, P. L.F. and C. C. Mei. (1976), Water motion on a beach in the presence of a Breakwater, I and II. J. Geophys. Res., 81, 30793094. LonguetHiggins, M. S., and R. W. Stewart (1960), The changes in the form of short gravity waves on long waves and tidal currents, J. Fluid Mech., 8, 565583. LonguetHiggins, M. S., and R. W. Stewart (1961), The changes in the form of short gravity waves on steady nonuniform currents, J. Fluid Mech., 10, 529549. LonguetHiggins, M. S., and R. W. Stewart (1962 ), Radiation stress and mass transport in gravity waves, with application to surfbeats, J. Fluid Mech., 13, 481504. LonguetHiggins, M. S. and R.W. Stewart (1964), Radiation stresses in water waves: a physical discussion with applications, DeepSea Research, 11, pp. 529. LonguetHiggins, M.S. (1970), Longshore currents generated by obliquely incident sea waves. J. Geophys. Res., 75, pp. 6778. LonguetHiggins, M.S. (1975), On the joint distribution of the periods and amplitudes of sea waves, J. Geophys. Res., 80, pp. 2688. MacMahan, J. H., A. J. H. M. Reniers, E. B. Thornton, and T. P. Stanton (2004), Surf zone eddies coupled with rip current morphology, J. Geophys. Res., 109, C07004, doi:10.1029/2003JC002083. Madsen, P. A., R. Murray, and O. R. Srensen ( 1991), A new form of Boussinesq equations with improved linear dispersion characteristics, Coastal Eng., 15, 371388. Madsen, P. A., and O. R. Srensen (1992), A ne w form of Boussinesq equations with improved linear dispersion character istics, Part 2: A slow ly varying bathymetry, Coastal Eng., 18, 183204. Majumdar, S. (1988), Roles of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids, Numer. Heat Transfer, vol. 13, pp. 125132. PAGE 255 255 Majumdar, S. (1986), Development of a FiniteVolume Proc edure for Prediction Fluid Flow Problems with Complex irregular Boundaries, Rep. 210/T/29, SFB210, University of Karlsrule, Karlsrule, Germany. McKenzie, P. (1958), Rip current systems, J. Geol., 66, 103113. Mei C.C. (1983), The applied dynamics of ocean surface waves. In: Advanced Series on Coastal Engineering vol. 1, World Scientific, Singapore, 734 pp. Muzaferija, S. (1994), Adaptive finite volume method for flow predictions using unstructured meshes and multigrid approach, PhD Thesis, University of London. National Research Council (NRC), Meeting Research and Educat ional Needs in Coastal Engineering, Washington, DC, National Academy Press, 1999. Newberger, P. A., and J. S. Allen (2007), Forcing a threedimensional, hydrostatic, primitiveequation model for application in the surf zone: 1. Formulation, J. Geophys. Res., 112, C08018, doi:10.1029/2006JC003472. Noda, E. K. (1974), Waveinduc ed nearshore circulation, J. Geophys. Res., 79(27), 40974106. Nwogu, O. (1993), An alternative form of the Boussinesq equations for nearshore wave propagation, J. Waterway Port Coastal Ocean Eng., 119(6), 618638. Ogston, A.S., and R.W. Sternberg (2003), Nearshore sediment suspension and flux in response to lowinfragravity frequency and episod ic processes: exploratory analysis, Coastal Sediments Proceedings. OltmanShay, J., P.A. Howd, and W.A. Birkemei er (1989), Shear instab ilities of the mean longshore current 2. field observations. J. Geophys. Res., 94(C12), doi:10.1029/89JC01392, 1989 zkanHaller, H. T., and J. T. Kirby (1999), Nonl inear evolution of shear instabilities of the longshore current: A comparison of observations and computations, J. Geophys. Res., 104(C11), 2595325984. zkanHaller, H.T., and Y. Li (2003), Wavecurre nt interaction effects on shear instabilities, J. Geophys.Res., 108(C5), 3139, doi: 10.1029/2001JC001287. Pangchang, V. G., G. Wei, B. R. Pearce, and M. J. Briggs (1990), Nu merical simulation of irregular wave propa gation over shoal, J. Waterway Port Coastal Ocean Eng., 116(3), 324340. Patankar, S. V. (1980), Numerical heat transfer and fluid flow. McGrawHill, New York. PAGE 256 256 Patankar, S. V., and D. B. Spalding (1972) A calculation procedure for heat mass and momentum transfer in threedimensional parabolic flows, Int. J. Heat Mass Transfer, 15, 17871806 Peregrine, D. H. (1967), Long waves on a beach, J. Fluid Mech., 27, 815827. Peric, M. (1985), A finite volume method for the predicti on of three dimensional fluid flow in complex ducts, PhD Thesis, University of London. Pierson, W. J., and L. Moskowitz (1964), A pr oposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. J. Geophys. Res., 69, 51815190. Putrevu, U. and I.A. Svendsen (1999), Threedimensional dispersion of momentum in waveinduced nearshore currents, Eur. J. Mech B. Fluid., pp 409. Radder, A. C. (1979), On the parabolic equation method for waterwave propagation, J. Fluid Mech., 95, 159176. Rattanapitikon W., and T. Shibayama (1998), Energy dissipation model for regular and irregular breaking waves, Coastal Engineering Journal, 40(4), 327346. Reniers, A.J.H.M., and J. A. Battjes (1997), A laboratory study of longshore currents over barred and nonbarred beaches, Coastal Eng., 30, 122. OReilly, W.C. and R.T. Guza (1991), Comparis on of spectral refraction and refractiondiffraction wave models, J. Waterway Port Coastal Ocean Eng., 117, 199215. Rhie, C. M. and W. L. Chow (1983), A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation, AIAA Journal, vol. 21, pp 15251532. Ris, R.C., L.H. Holthuijsen and N. Booij (1994), A spectral model for waves in the near shore zone, Proc. 24th Int. Conf. Coastal Eng., Kobe, Japan, pp. 6878. Ris, R.C., N. Booij and L.H. Holthuijsen ( 1999), A thirdgeneration wave model for coastal regions, Part II: Verification, J. Geophys. Res., 104(C4), 76677681. Roelvink, J.A., and M. J. F. Stive (1989), Barg enerating crossshore flow mechanisms on a beach, J. Geophys. Res., 94, pp 4785. Ruessink, B.G., J. R. Miles, F. Feddersen, R. T. Guza and S. Elgar (2001), Modeling the alongshore current on barred beaches. J. Geophys. Res., 106, 2245122464. Sancho, F. (1997), Unsteady nearshore currents on longshore varying topographies, Ph.D. thesis, University of Delware, 1997. PAGE 257 257 Shepard, F. P., K. O. Emery, and E. C. LaF ond (1941), Rip currents: A process of geological importance, J. Geol., 49(4), 337369. Shepard, F. P., and D. L. Inman (1950), Nears hore circulation related to bottom topography and wave refraction, Tran. Amer. Geophys. Union, 31(4), 555565. Sleath, J. F. A. (1984), Sea Bed Mechanics. John Wiley and Sons, New York. Slinn, D. N., J. S. Allen, and R. A. Holman (2000), Alongshore currents over variable beach topography, J. Geophys. Res., 105(C7), 16,97116,998. Slinn, D. N., J. S. Allen, P. A. Newberger, and R. A. Holman (1998), Nonlinear shear instabilities of alongshore curr ents over barred beaches. J. Geophys. Res., 103, 18,35718,379. Smith, J. A., and J. L. Largier (1995), Observa tions of nearshore circ ulation: Rip currents, J. Geophys. Res., 100, 1096710975. Smith, J.M., M. Larson, and N. C. Kraus (1994 ), Longshore current on a barred beach: field measurements and calculation. J. Geophys. Res., 98(C12), 22,717,731. Sonu, C. J. (1973), Field observation of nearshore circulation and meandering currents, J. Geophys. Res., 77(18), 32323247. Srensen, O. R., H. A. Schffer, and P. A. Madsen (1998), Surf zone dynamics simulated by a Boussinesq type model, III. Waveinduced horizontal nearshore circulations, Coastal Eng., 33 (23), 155176. Stive, M. J. F., and H. J. De Vriend (1994), Shear stress and mean flow in shoaling and breaking waves, in Proceedings 24th international Coastal Engineering Conference, pp. 594608, ASCE, New York. Stone, H. L. (1968), Iterative solution of im plicit approximations of multidimensional partial differential equations. SIAM J. Numer. Anal., 5, 530558. Svendsen, I.A. (1984a), Wave hei ghts and setup in a surf zone. Coastal Eng., 8, pp 303. Svendsen, I.A. (1984b), Mass flux and undertow in a surf zone. Coastal Eng., 8, pp 347. Svendsen, I. A. and U. Putrevu (1994), Nearshore mixing and dispersion, Proc. Roy. Soc. Lond. A, 445, 561576. Thornton, E. B. (1970), Variation of L ongshore Current Across the Surf Zone, Proceedings of the 12th Coastal Engineering Conference, pp 291308. PAGE 258 258 Thornton, E. B., and R. T. Guza (1983), Transf ormations of wave height distribution, J. Geophys. Res., 88, 59255938. Thornton, E.B., and R. T. Guza (1986), Surfzone longshore currents and random waves: field data and models. J. Phys. Ocean., 16, 7, 11651178. Tolman, H.L. (1991), A thirdgeneration model for wind waves on slowly varying, unsteady and inhomogeneous depths and currents, J. Phys. Ocean., 21, 6, 782797. Tolman, H. L. (1999), User manual and system documentation of WAVEWATCHIII version 1.18. Technical report, NCAR Ocean Modeling Branch. 110 p. Tolman, H. L., and D. Chalikov (1996), Source terms in a thirdgeneration windwave model, J. Phys. Ocean., 26, 24972518. Van Dongeren, A., F. Sancho, I. A. Svendsen, and U. Putrevu (1994), SHORECIRC: A quasi 3D nearshore model In Proceedings 24th Coastal Engineering Conference, pp 27412754. Vincent, C. L., and M. J. Briggs (1989), Re fractiondiffraction of i rregular waves over a mound, J. Waterway Port Coastal Ocean Eng., 115(2), 269284. Walstra, D. J. R., G. P. Mocke, and F. Smit (1996), Roller contributions as inferred from inverse modeling techniques, in Proceedings 25th international Coastal Engineering Conference, pp. 12051218, ASCE, New York. WAMDI group (1988), The WAM model a third ge neration ocean wave prediction model, J. Phys. Ocean., 18, 17751810. Wei, G., J. T. Kirby, S. T. Grilli, and R. Subramanya (1995), A fully nonlinear Boussinesq model for surface waves, 1. Highly nonlinear unsteady waves, J. Fluid Mech., 294, 7192. Whitford, D. J. (1988), Wind and wave forcing of longshore currents across a barred beach. PhD dissertation, Nav. Postgrad. School, Monterey, California. Whitham, G. B. (1965), Linear and nonlinear waves. Wiley, New York, 636p. Willebrand, J. (1975), Energy transport in a nonl inear and inhomogeneous random gravity wave field, J. Fluid Mech., 70, 113126. Yu, J., and D. N. Slinn (2003), Effects of wavecurrent interaction on rip currents, J. Geophys. Res., 108(C3), 3088, doi: 10.1029/2001JC001105. Zang, Y., R. L. Street, and J. R. Koseff (1994), A nonstaggered gr id, fractional step method for timedependent incompressible Navierstokes equations in curvil inear coordinates. J. Comput. Phys., 114, 1833. PAGE 259 259 BIOGRAPHICAL SKETCH Yang Zhang was born as the first child of three of Zhijun Zhang and Mengshu Huang in 1979 in Yanshi, China. He obtaine d his bachelors degree in ci vil engineering in 2001. Shortly after, he was admitted to the physical oceanogr aphy program at Hohai University, where he obtained his Master of Science degree in 2004. U pon graduation, he joined the Second Institute of Oceanography, National Oceanic Ad ministration of China as an assistant researcher. In June 2004, he was admitted to the University of Florida and he decided to pursue his Ph.D. degree in coastal and oceanographic engin eering. His research interests involve ocean wave modeling, waveinduced currents, wavecur rent interaction, and coas tal sediment transport and morphological evolution. 