UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
NUMERICAL MODELING OF WAVEINDUCED CURRENTS USING A PARABOLIC WAVE EQUATION By HARLEY STANFORD WINER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 I BNPlsTY OF FLORIDI IBRARLES ACKNOWLEDGEMENTS The writing of a dissertation is a long process involving interaction with and support from many people. The author would like to thank all those who in one way or another offered encouragement, support, advice and helped in the problem solving process. The contributions of those who are not mentioned by name are in no way diminished by the lack of individual acknowledgement. First, the author wishes to express deep appreciation to the members of his committee: the chairman, Dr. Hsiang Wang, who enticed him to come to Florida and enroll in the doctoral program, and the cochairman, Dr. James Kirby, who shared his knowledge of wave equations and numerical modeling, along with Dr. Robert Dean, Dr. Bent Christensen and Dr. Joseph Hammack. Each of these men has been both a teacher and a friend. The author wishes to also thank Dr. Peter Sheng who at one time served on the supervisory committee but due to a prior commitment was unable to be at the defense and officially be on the committee. The author would like to thank the staff of the Coastal and Oceanographic Engineering Laboratory at the University of Florida and the staff of the department of Coastal and Oceanographic Engineering who helped and aided the author in innumerable ways. Thanks are also extended to Ms. Lillean Pieter who helped prepare the figures and Dr. Claudio Neves, a fellow student, who manipulated the ITE3X software so as to produce a document in thesis mode. A special acknowledgement is given to the author's parents who encouraged and pushed him to complete this doctorate. The author would also like to thank his teachers, friends and fellow students who offered encouragement and interesting discussions. Lastly the author thanks his soulmate and partner, Esther DeJong, who managed to be patient with each extension of the final date. Baby cakes, it's Phinally Done. I love you. TABLE OF CONTENTS ACKNOWLEDGEMENTS ................................. ii LIST OF FIGURES ..................................... vi ABSTRACT ................. ...................... ix CHAPTERS 1 INTRODUCTION ..................................... 1 1.1 Literature .................. .................. 4 1.2 General Description of the Model ........................ 8 2 CIRCULATION MODEL ............................... 16 2.1 Governing Equations .............................. 16 2.2 Radiation Stresses ................... .............. 24 3 WAVE MODEL ..................... ............... 31 3.1 Governing Wave Equation .............................. 31 3.2 A Parabolic Approximation .... ............. ........... .. 36 3.3 Energy Dissipation Due to Wave Breaking ..................... 45 4 FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDITIONS ..... 49 4.1 Numerical Solution of the Governing Equations ................ 49 4.2 Finite Differencing of the Parabolic Wave Equation .............. 52 4.3 Boundary Conditions in the Circulation Model ................ 56 4.4 Boundary Conditions in the Wave Model .................. 59 5 RESULTS ............... .......... ................ 63 5.1 Wave Setup ................... .... ............ 64 5.2 Longshore Currents .................. ............... 68 5.3 Waves On a RipChannel ............................. 70 5.4 ShorePerpendicular Breakwater .................. ..... 76 5.5 Shore Parallel Emerged Breakwater .................. .... 88 5.6 Comparison with Experimental Results of Gourlay ............. 92 6 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY ... 103 6.1 Improved Wave Formulations ........................... 104 6.2 Including Wave Reflection ............................. 105 6.3 Extending the Model to Include Sediment Transport ............. 105 APPENDICES A DERIVATION OF THE DEPTH AVERAGED EQUATIONS .......... 107 B DERIVATION OF THE RADIATION STRESS TERMS .............. 113 C LAGRANGIAN DERIVATION OF THE GOVERNING WAVE EQUATION 117 D OBTAINING EQUATION (3.79) .......................... 120 BIBLIOGRAPHY ................... ................. 124 BIOGRAPHICAL SKETCH ................... ........... 129 LIST OF FIGURES 1.1 Computational domain of the model ...................... 8 1.2 Flow chart of computer program ...................... 11 1.3 Flow chart of circulation model ....................... 13 1.4 Flow chart of wave model .............................. 14 2.1 Longshore current profiles for different friction formulations for an 11 second, 1.028 meter wave at 17.28 degrees to a 1:25 plane beach. ..... 22 2.2 Experimental results for wave setup and setdown. Reprinted from Bowen, Inman, and Simmons (1968), Journal of Geophysical Research, vol. 73(8) page 2573, copyright by the American Geophysical Union. 30 4.1 Definition sketch of grid nomenclature and coordinate system. ...... 56 4.2 Startup function according to Eq. 4.43 with C1=30 and C2=3. ..... 60 4.3 Waves passing through a lateral boundary unaffected by the boundary. 61 5.1 Comparison of the numerical solution of Vemulakonda et al. (1985) for setup with experimental data of Bowen (1968). Source: Vemulakonda et al. (1985) . . . ... 65 5.2 Comparison of numerical results from the model with experimental re sults of Hansen and Svendsen (1979). Wave setup and wave heights for a 2 second wave on a beach slope of 1:34 and incident wave height of 7 centimeters. ............... .. .............. 67 5.3 Nondimensional longshore velocities vresus the nondimensional dis tance offshore from the analytical solution of LonguetHiggins. Reprinted from LonguetHiggins (1970b), Journal of Geophysical Research, vol. 75(33) page 6793, copyright by the American Geophysical Union. .... 68 5.4 Comparison of experimental results of Visser (1982) and numerical result for longshore currents due to a 1.1 second wave on a slope of 1:20 for the deep water conditions Ho = .065m and 0. = 20 degrees. ........ 69 5.5 Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as shown. Bottom: Depth contours as used in the present model. Grid spacing is 5 meters. .................. ...... ... 71 5.6 Velocity vectors from the results of Ebersole and Dalrymple for a wave angle of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec. Source: Ebersole and Dalrymple (1979). ................ 72 5.7 Velocity vectors obtained from the present model using the same input conditions as for Ebersole and Dalrymple. ............... 73 5.8 Mean water surface elevation contours showing that the water surface is lower above the trough. Values indicated are in millimeters. ........ 74 5.9 Depth contours in meters for the case of a bar with a gap. ........ 75 5.10 Velocity vectors for the case of waves over a bar with a gap. Dashed lines indicate depth contours. Grid spacing equal to 10 meters. ..... 75 5.11 Depth contours for the case of normally incident waves upon a shoal. Grid spacing is .05 meters. Contour interval is .02 meters. ....... 77 5.12 Velocity vectors for the case of waves breaking upon a shoal. ....... 78 5.13 Maximum velocity across the top of a shoal versus incident wave height for waves over a shoal on a plane beach. Incident wave height less than 9.1 centimeters are nonbreaking waves. Incident wave heights greater than 9.1 centimeters produced waves breaking over the shoal. ....... 79 5.14 Position of the groin in relation to the grid rows and columns. ...... 80 5.15 Results for the numerical model of Liu and Mei (1975). Deepwater wave height equal to 1 meter, deepwater wave angle equal to 45 degrees. Source: Liu and Mei (1975). .................. .... 83 5.16 Results from the present numerical model using the same conditions as were used by Liu and Mei. Offshore boundary condition of constant wa ter surface level which precludes tangential flow at the offshore boundary. 84 5.17 Results from the present numerical model using the same conditions as were used by Liu and Mei. Offshore boundary condition of no onoffshore flow. Grid spacing equal to 20 meters. ................. 85 5.18 Experimental setup for the tests conducted by Hales (1980). Source: Hales (1980) . . . . .. 86 5.19 Experimental measurements of the middepth averaged velocities obtain by Hales. Measurements at one foot intervals. Source: Hales (1980). 87 5.20 Ripcurrent velocities adjacent to the down wave side of the jetty. Com parison of experimental and numerical results ............ 88 5.21 Velocity vectors of the currents obtained numerically for the same con ditions as for the Hales experiment. Grid spacing equal to one foot. 89 5.22 Definition sketch of grids for the case of a shore parallel breakwater. .. 90 5.23 Wave crests in the lee of a semiinfinite breakwater in a region of constant depth. . . . . ... 91 5.24 Velocity vectors obtained from the model for a 100 meter long breakwater 100 meters offshore for a 1 meter, 6 second wave for increasing angles of incidence . . . . 93 5.25 Stream function contour lines obtained from the Liu model for a 10 second 1 meter wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope. Source: Liu and Mei (1975). 94 5.26 Flow lines obtained from the present model for a 10 second 1 meter wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope. .................. ... .. 95 5.27 Velocity vectors obtained from the present model for a 10 second 1 meter wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope. .................. ... .. 95 5.28 Setup contour lines obtained from the present model for a 10 second 1 meter wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope. .... .......... 96 5.29 Physical layout for the experiment of Gourlay. Source: Gourlay (1974), Proceedings 14th Coastal Engineering Conference, copyright by the Amer ican Society of Civil Engineers, reprinted with permission. ....... 97 5.30 Computational domain for the experiment of Gourlay. Grid spacing equal to .1 meter ................... . 98 5.31 Comparison of experimental and numerical results for the experiment of Gourlay: Wave Setup Contours. Top figure is from Gourlay (1974) Pro ceedings 14th Coastal Engineering Conference, copyright by the Amer ican Society of Civil Engineers, reprinted with permission. The bottom figure is the computer results. .............. . 100 5.32 Results from the experiment of Gourlay: Contours of current veloci ties and streamlines. Source: Gourlay (1974), Proceedings 14th Coastal Engineering Conference, copyright by the American Society of Civil En gineers, reprinted with permission. ............... 101 5.33 Results from the numerical model: Streamlines of the flow. ........ 101 5.34 Results from the numerical model: Contour lines of the current veloci ties . . . . 102 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL MODELING OF WAVEINDUCED CURRENTS USING A PARABOLIC WAVE EQUATION By HARLEY STANFORD WINER August 1988 Chairman: Dr. Hsiang Wang Major Department: Aerospace Engineering, Mechanics and Engineering Science Prediction of waves and currents in the nearshore region in the presence of proposed breakwaters is needed in the design process in order to optimize the placement of the structures. A numerical model is developed which solves for the wave field, the depth averaged mean currents and the mean water surface elevation within a given computa tional domain with the presence of breakwaters. The model iterates between a wave model and a circulation model. The wave model uses a parabolic approximation to a combined refractiondiffraction mild slope equation which includes currents. The circulation model uses an alternating direction implicit method solution to the depth averaged equations of momentum using the radiation stress terms obtained through solution of the wave model as the forcing terms. Each of the terms in the depth averaged equations of momentum involve assumptions that are discussed in the report. The model is constructed in modular form so that improvements can easily be incorporated into the model. The model was employed for several test cases such as rip channels, bar gaps, and both shore perpendicular and shore parallel breakwaters. Good agreement was obtained when the model was used to simulate a laboratory experiment of Michael R. Gourlay ("Wave setup and wave generated currents in the lee of a breakwater or headland," Proc. 14th Coastal Eng. Conf., Copenhagen, 1974, pp. 19761995). Suggestions are made for improving the model and for including sediment transport in the model. CHAPTER 1 INTRODUCTION As waves propagate they are transformed by variations in currents and depths, along with interaction with obstructions. The major processes of wave transformation include shoaling, reflection, refraction, diffraction and dissipation. Shoaling is the change of wave energy, propagation speed and wave height due to changes in the water depth and/or current field. Refraction is the change in the direction of wave propagation due to the spatial gradient of the bottom contours and/or current field. Diffraction is the phenomenon of wave energy being transmitted laterally along a wave crest. Dissipation is the process of energy loss due to a multitude of possible causes such as bottom percolation, damping due to viscous effects, turbulence, and wave breaking. As waves approach a shoreline and break, momentum exchanges can produce currents and mean water level variations. These currents can be in the form of longshore currents, rip currents, or circulation cells in the lee of obstructions. They in turn modulate the incident wave field. The process can be quite complicated depending upon the complexity of the bathymetry, the ambient current field, and the presence of obstructions. Coastal engineers and planners have many reasons for wanting to obtain accurate pre dictions of waves and currents and water elevations in the coastal xooe. For instance, pre dictions of extreme conditions are needed for site planning and erosion mitigation studies. Harbor authorities need wave predictions for structure design and for operational decisions. This dissertation develops a numerical model for obtaining predictions of waves, currents and mean surface elevations in a given domain with the offshore wave conditions and the bathymetry as the input. The model uses a linear wave equation and depth integrated equations of momentum and mass conservation. Comparison with experimental results 2 indicates that the model produces useful and valid predictions of the waves and currents in the nearshore region in the presence of obstructions such as groins and breakwaters. The predictions obtained through use of the computer model herein reported are not intended to be the sole resource for the engineer. A good engineer or planner will seek answers from multiple sources and cross check them, recognizing the limitations of each method used to obtain the answer. Predictions of waves and currents can be obtained for different cases analytically, through experiments or by numerical solutions. Analytical solutions are obtained only for simple cases involving idealized beaches using simplified assumptions which are not al ways satisfied in nature. LonguetHiggins and Stewart (1963) give an analytic solution for wave setup for steady state conditions on a plane beach with straight and parallel con tours, assuming linear waves approaching the beach normally, and spilling breakers where the wave height shoreward of the breaker line is a constant times the total depth of the water column. LonguetHiggins (1970a) gives a cross shore profile of the longshore current for a plane beach with straight and parallel contours, for linear waves approaching the beach at an angle, spilling breakers, and linearized bottom friction, while neglecting the influence wave setup has upon the total depth of the water column. In a companion paper, Longuet Higgins (1970b) extended the solution to include lateral diffusion of the current momentum. There are other analytical solutions for waves and currents but all are limited to very simple cases. The literature is replete with laboratory investigations of waves and currents. Physical models have advantages and disadvantages. The advantages of experimentation are that results can be obtained for complex situations where there is uncertainty about the exact form of the governing equations or about the precise form of any of the component terms of the equation. Any reasonable amount of detail can be included. Visualization of the physical processes provides an opportunity to identify problems which otherwise might not be evident. 3 However there are limitations to the use of laboratory experiments. There is the problem of scale effects which means that all scaling requirements cannot be satisfied simultaneously; thus, all features may not be interpreted properly to the prototype. For example, the time scale dictated by Froude modelling laws in which gravitational forces dominate is different from that obtained from Reynolds modelling laws in which viscous forces are most important. Another disadvantage of physical models is the cost of labor to construct and run large scale experiments. Also, because of the cost of construction, it is not always very easy to change details in the model. Thus, it is not always possible to isolate one parameter for study. For example, an investigation of the influence of the size of the gap between segmented breakwaters might require many reconstructions of the rubble mound structure. The large size of a model means that it is not easily stored once the initial experiment is completed. The options are to either dismantle the model upon completion of testing or to incur some maintenance costs. The significance of this is that the model cannot always be easily retested if, for example, additional input information becomes available or man or nature makes a substantial alteration of the prototype after the model has been dismantled. Also the results of a model for one prototype are not always transferable to other prototype situations. Often the available space for the physical model will require the use of a distorted model in which the vertical scale is different from the horizontal scale. This presents the additional problem that for a distorted model, when modelling a situation where both refraction and diffraction are equally important, the laws of scale modelling dictate different time scales for refraction and diffraction. For refraction the time ratio is T, = 1/TL, where the subscript r is for the ratio of the value in the model to the value in the prototype and T and Lv are time and vertical length respectively. For diffraction T, = iVf in shallow water and T, = v\/L in deep water, where LH, represents the horizontal length scale. Numerical models do not have many of the disadvantages found with physical models. Any mechanism for which the governing equations are known can be modeled numerically. 4 There are no limits in terms of scale effects. Numerical models, once developed, can be adapted easily to many different prototype situations. There is no problem in the storage or duplication of the model. Today computer facilities are far more available than laboratory facilities. Computer models easily allow for sensitivity testing of the importance of an individual variable and any specific feature of the model can be monitored. Numerical models also have their disadvantages. In fact some of their advantages are also their disadvantages. For example, being able to monitor any feature also means that the computer will deliver only the information that is explicitly specified. Thus an investigation of one aspect will not indicate problems with other aspects. Though any mechanism for which the governing equation is known can be modeled, the model is only as good as the assumptions involved in the equations. For example, the forces exerted on the water column by the complex interaction between the fluid motion and the bottom are not totally understood and are often greatly simplified. Another example is the use of depth integrated twodimensional models of fluid flows when the reality is threedimensional. Both of the above are incorporated in the present study. Knowing the correct equation does not always make for successful computer modelling since the complexity of the equations may make for difficulties (for example solving for the complete threedimensional flow field using the NavierStokes equations). Computer time and cost for large programs can often be significant. 1.1 Literature This section is a brief history and literature review of the attempts of the coastal engineering profession to obtain predictions of waves and currents in the nearshore region given the offshore wave conditions. The extension of Snell's law which governs optical wave refraction to the analogous problem of water waves over straight and parallel contours has been known for some time. Johnson, O'Brien, and Isaacs (1948) reported graphical methods for constructing wave refraction diagrams. Two methods were reported, one producing wave fronts and the second 5 producing wave rays. Pierson (1951) showed that the method using wave fronts introduced more error than the wave ray method. These methods proceeded upon the assumption that the contours are locally straight. The wave rays obtained were important in that the wave height was determined to vary inversely as the square root of the spacing between the wave rays. Using these ray construction techniques O'Brien (1950) was able to demonstrate that the April 1930 destruction of the Long Beach Harbor breakwater on a seemingly calm day most probably resulted from the focusing of long wave energy from a South Pacific storm. Munk and Arthur (1952) showed a method to obtain the wave height directly from the curvature of the ray and the rate of change of the depth contours. The wave ray tracing techniques did not include the effects of currents and had the disadvantage (from a numerical modeling point of view) that wave heights and directions were obtained at positions along a ray and not at prespecified grid positions. Noda et al. (1974) solved this problem by devising a numerical scheme that solved a wave energy equation and a statement for the irrotationality of the wave number to obtain wave height and direction at grid points. This scheme also included the effects of current refraction but neglected diffraction. Numerical models of wave induced currents were developed by Noda et al. (1974) in which time dependent and advective acceleration terms were ignored. Birkemeier and Dal rymple (1976) wrote a timemarching explicit model which neglected lateral mixing and advective acceleration. Ebersole and Dalrymple (1979) added lateral mixing and the non linear advective acceleration terms. Vemulakonda (1984) wrote an implicit model which included lateral mixing and advective acceleration terms. All of the above used the numer ical scheme of Noda et al. (1974) to obtain wave heights and wave angles. Each of these models neglected diffraction. Diffraction of water waves represents some very difficult mathematical problems. His torically most solutions were based upon the assumption of a flat bottom, since to quote Newman (1972), The restriction to constant depth h is important.... Problems involving wave propagation in a fluid of variable depth are of great interest, and the resulting 6 diffraction effects are among the most important of all ... but these problems are less tractable than those involving a fluid of constant depth. (page 2) Penney and Price (1944) showed that Sommerfeld's solution of the optical diffraction prob lem as described by Bateman (1944) is also a solution of the constant depth water wave diffraction problem. Penney and Price (1952) solved the constant depth diffraction problem for breakwaters assuming a semiinfinitely long, infinitesimally thin barrier. The solution is in terms of Fresnel integrals which need to be evaluated numerically. The Shore Protec tion Manual (SPM,1977) produced many graphical solutions based upon Penney and Price (1952) that give the wave height and direction in the lee of breakwaters. The problem of combined refraction and diffraction according to the SPM could be handled by assuming a constant depth for a few wavelengths so as to obtain the diffraction effects and then continue with the refraction effects. Perlin and Dean (1983) used essentially the same procedure in a numerical model which calculated the longshore thrust in terms of the breaking wave height and angle. Liu and Mei (1975) solved for the wave field around shore parallel and shore perpendicu lar breakwaters by using a curvilinear coordinate system where one coordinate corresponded to the wave ray and the other to a curve of constant phase. Their solution was for "an am plitude modulation factor D(f,f) which accounts for diffraction" (page 72), where C and f are the two coordinates. An equation for D(e, f) was obtained and a boundary value prob lem was established by matching conditions in the incident zone and the shadow zones. For constant depth the solution for D(C, f) is in terms of Fresnel integrals. Upon obtaining the wave field the radiations stresses were obtained and used to force a circulation model that solved for the stream function and mean water surface level while ignoring the nonlinear advective acceleration terms and the lateral diffusion of current momentum. Although Liu and Mei produced results for the wave field, currents, and setup resulting from the diffrac tion of waves by breakwaters, their study is not a wavecurrent interaction model since there are no current terms in the wave equations and no feedback of current terms from the circulation model. 7 Berkhoff (1972) derived an equation governing combined refraction and diffraction with the assumption of a mild slope. Radder (1979) applied the splitting method of Tappert (1977) to obtain a parabolic equation governing the propagation of the forward propagat ing component of the wave field. Booij (1981) extended the work of Berkhoff to include the effects of currents and also obtained a parabolic approximation. Kirby (1984) made corrections to the equation of Booij. Parabolic wave equations have recently been used in waveinduced circulation models. Liu and Mei (1975) solved for the wave field by matching conditions between the shadow zones and the illuminated zone in the lee of a structure. The wave information was then used to force a circulation model which formulated the governing depth integrated time averaged momentum equations in terms of a stream function. However this model did not have current feedback in the wave equation. Yan (1987) and the present study both include a current feedback for the parabolic wave equation and solve directly for the mean currents. Prior to the use of parabolic wave equations to obtain the wave forcing in circulation models, refraction wave models using a numerical scheme derived by Noda et al. (1974) were used. Birkemeier and Dalrymple (1976) produced an explicit numerical solution of the depth averaged equations of momentum and continuity using the wave model of Noda et al. (1974). Their model did not include the nonlinear advective acceleration terms or any mixing. Ebersole and Dalrymple (1980) included the advective acceleration terms and lateral mixing in an explicit model which also used the refraction scheme of Noda et al. Vemulakonda (1984) Produced an implicit model which also was based upon the wave refraction scheme of Nods et al. Yan (1987) used a parabolic approximation to the mild slope equation which accounts for combined refraction and diffraction. The present model Improves upon Yan's model in numerical efficiency, by introducing structures and most importantly by obtaining the radiation stress terms directly from the complex wave amplitude. J=N J=N1 J=N2 LL C 0LLO T J=3 J=2 J=l _____LflTERflL BOUNDPRT_____/   __ __ / ^ _ /  / __ /  ___ /  ___ /  /* __ /  _ /  _ /  _ / ___ /  ___ /  ___ /  ___ /  ^ __ /  ~ ~ / _ / i 1 LATERAL n 7 BOUNDARY Y x x .  Figure 1.1: Computational domain of the model. 1.2 General Description of the Model This section is a description of the numerical model that has been constructed to calcu late the wave and current field within a given domain. This domain usually is bordered by a shoreline along one boundary and an offshore region along the opposite boundary. The other boundaries are sufficiently far from the region of interest (i.e. a structure or shoal, etc.) so as to not affect the region of interest. This will be discussed further in the section on boundary conditions. The computational domain for the purposes of this investigation consists of a rectilinear coordinate system. Figure 1.1 shows the computational domain. The model consists of two parts, a circulation model and a wave model. The circulation model computes the mean water surface level and the depth averaged mean currents using the radiation stresses from the waves as the forcing terras and includes bottom friction, advective acceleration terms and the lateral diffusion of current momentum. The wave model determines the radiation stress terms by first solving for the complex amplitude as a 9 function of the total depth and the depth averaged mean currents, and then computing the radiation stresses from the complex amplitudes. The model starts with initial conditions of zero mean currents and zero mean water surface elevation. The offshore incident wave height is initially zero and then is smoothly increased up to its given value using a hyperbolic tangent curve. This is done so as to minimize any transients due to a shock startup. The model iterates between the wave model and the circulation model until steady state convergence is obtained. Each of the two parts of the model will be described separately. The circulation model uses an alternating direction implicit (ADI) method of solution. Three equations are solved for three unknowns. The unknowns are q7 the mean water surface level, and U and V, the two components of the depth averaged mean current. The equations are the continuity equation and the x and ydirection momentum equations. Each iteration of the circulation model solves for the three unknowns on the next time level. In other words each iteration advances the values of i, U and V by a time increment, At. For the purposes of this investigation the model solves for steady state conditions using the time marching character of the circulation model. The model can be used to solve for time dependent wave fields, but this would require the substitution of a time dependent wave equation in the wave model. The ADI method first solves, grid row by grid row, for U at the next time level and an intermediate value of fi, by solving the x direction equation of motion and the continuity equation using an implicit numerical scheme. This solution is in terms of the known values of U,V,and f7 at the present time level. Once all the grid rows are solved in the x direction, the y equation of motion and the continuity equation are solved, again by an implicit scheme, grid column by grid column, for V and q7 at the next time level in terms of U at the new time level, qf at the intermediate level and V at the present level. For each of the implicit schemes in the two directions, the governing equations yield a tridiagonal matrix. The tridiagonal matrix and the method of solution will be discussed later. The wave model also yields a tridiagonal matrix equation. The governing equation is a parabolic equation for the complex amplitude A. The complex amplitude contains phase 10 information meaning among other things that y variations of the wave field are contained within the complex amplitude. The term "parabolic equation' indicates that there are first order derivatives in the x direction and second order derivatives in the y direction. The method of solution is to solve for the amplitude on the 1+1 grid row in terms of the known amplitude values on the I grid row. The solution scheme starts on the first grid row using the prescribed offshore incident wave and then advances grid row by grid row to the last grid row. To advance to the next grid row the amplitude values on the 1+1 grid row are implicitly formulated using a CrankNicolson scheme. This will be further discussed in section 4.2 which describes the finitedifferencing scheme for the parabolic wave equation. Once the complex amplitudes are known the radiation stress terms are computed in terms of the absolute value of the amplitude and the gradients of the amplitude. This is described in section 2.2. The circulation model and the wave model are the heart of the program; however there is much more. In all there are some 27 subroutines in the program. These subroutines cover everything from input and output to convergence checks and flooding of the shoreline. The main program is compartmentalized into many subroutines so as to facilitate any changes or future upgrading. For example, the bottom shear stress, a term of importance in the equations of motion, currently uses a crude approximation. If a better approximation is developed that is also computationally efficient, it will be an easy task to replace the present subroutine that computes the bottom shear with the newer subroutine. There are many places where the model can be improved with better assumptions and better approxima tions. These will be discussed in detail in the final chapter outlining recommendations for further study. The modular construction of the model also allows for ease in substituting one wave model for another. Figure 1.2 shows a flow chart of the computer program. The subroutine INPUT is called to initialize the values of the unknowns and to input the water depth, the offshore wave height and direction and the position of any structures. The three subroutines WAVE, Figure 1.2: Flow chart of computer program. 12 SHEAR, and LATMIX compute the radiation stresses, the bottom shear stress, and the lateral mixing coefficients. These terms are the forcing terms in the equations of motion and are used in the subroutine CIRC, which is the circulation model, to solve for q, U and V. There is a decision tree which directs the program to the three subroutines WAVE, SHEAR, and LATMIX at selected intervals. This is purely arbitrary and is done in the interest of computer speed. It is assumed that the forcing terms do not change at a great rate and can thus be updated periodically. This updating is more frequent during the startup time while the offshore incident wave height is being increased. The CHECK subroutine is to check the model for convergence. Since steady state conditions are assumed, this subroutine indicates convergence when the absolute value of the percentage difference between the updated values of U, V, and qi and the previous values is less than an arbitrarily small specified value. The UPDATE subroutine simply updates the values of the unknowns for the next iteration. The FLOOD subroutine allows for the addition or elimination of grid rows at the shoreline to allow for beach flooding due to setup. Figure 1.3 is a flow chart of the circulation model. The subroutines XCOEF and YCOEF determine the coefficients of the tridiagonal matrix equation for the general cases in the x and y directions. The subroutine TRIDA solves the tridiagonal matrix equation using the doubleelimination scheme described by Carnahan, Luther, and Wilkes (1969). Subroutines XCOJ1 and XCOJN are adaptations of XCOEF to determine the coefficients of the tridiagonal matrix for the special case along the lateral boundary where J=1 or J=N. Likewise YCOIM is an adaptation of YCOEF for the I=M grid row which borders the shoreline. These adaptations are necessary because boundary conditions necessitate different numerical formulation of some of the derivative terms. Similarly with the presence of an emerged impermeable barrier bordering the shore side of the JJJ grid row, special subroutines YCJJJ and YCJJJ1 are used on the grid rows I=JJJ and I=JJJ+l so as to properly include the boundary conditions for these grid rows, and for shore perpendicular groins subroutines XCJGR and XCJGR1 incorporate the boundary conditions for a groin Figure 1.3: Flow chart of circulation model. Figure 1.4: Flow chart of wave model. located between the J = JGR and the J = JGR+1 grid columns. Figure 1.4 is a flow chart of the wave model. Subroutine WAVNUM determines the wave number k, the group velocity C, and the intrinsic frequency a at each of the grid centers. These terms are used in the governing parabolic equation. Subroutine AMPCO determines the coefficients of the tridiagonal matrix and CTRIDA solves the matrix equation for the complex amplitude on the next grid row. CTRIDA is essentially the same as TRIDA except that it has been modified to handle complex numbers. The subroutine WW is called by the AMPCO subroutine in order to determine the dissipation coefficient following the work of Dally (1980,1987). Chapter 2 describes the circulation model in greater detail. Individual sections will cover the governing equations and the formulation of each of the component terms. Chap ter 3 describes the wave model, the derivation of the linear wave equation, the parabolic 15 approximation, the manner of including wave energy dissipation due to wave breaking and the method of obtaining the radiation stress terms directly from the complex amplitude. The fourth chapter gives details of the numerical methods used to obtain solution of the governing equations and the boundary conditions employed. The fifth chapter presents the results of applications of the model and comparisons with other studies both numerical and experimental, as well as some analytical solutions to some idealized cases. The final chapter summarizes the study and makes recommendations for future research. CHAPTER 2 CIRCULATION MODEL 2.1 Governing Equations The governing equations for the circulation model are the depth integrated time aver aged horizontal equations of momentum au au au a, j 1 1 S+U a +V +g + I 7,I ot ox 9y ax pD pD 1 (as,, as,, la, + a + +  = (2.1) pD a x ay pay av uav + av an 1 1 t _x ax ay Ty pD pD I a + a + =0 (2.2) pD ax ay pax and the continuity equation a# a a a + (UD)+ (VD) O 0 (2.3) where U = x component of the mean current V = y component of the mean current q = mean water surface elevation p = water density D = the total water depth = h. + 4 h. = local stillwater depth rT = the lateral shear stress due to turbulent mixing rt. = x component of the bottom shear stress 17 r&y = y component of the bottom shear stress r,. = x component of the surface shear stress r,. = y component of the surface shear stress and Ss,, Sy, and S. are the radiation stress components which arise from the excess momentum flux due to waves. It should be noted that U and V include the wave induced mass transport. These equations are obtained by integrating the local x and y momentum equations and the continuity equation over the depth of the water column and then time averaging the results. This is demonstrated in Appendix A, which follows closely the work of Ebersole and Dalrymple (1979). For greater details of the derivation the reader is referred to Ebersole and Dalrymple (1979) or Phillips (1969). These are the same equations that have been used by other modelers such as Birke meier and Dalrymple (1976), Ebersole and Dalrymple (1979), Vemulakonda (1984), and Yan (1987). Birkemeier and Dalrymple used an explicit numerical scheme and neglected the advective acceleration terms and the lateral mixing. Ebersole and Dalrymple added the advective acceleration and lateral mixing terms but still used an explicit numerical scheme. Vemulakonda introduced an implicit numerical formulation for these equations. Each of these investigators used a refraction wave model based upon the work of Noda et al. (1974) which solves for the wave height and wave angle at given grid points. Yan (1987) used a combined refractiondiffraction wave equation to determine the complex wave amplitude (magnitude and phase) and then solved for the wave angle in terms of the slope of the wa ter surface. Each of these studies obtained the radiation stress components in terms of the wave height and wave angle. The present study differs from the previous ones in that the radiation stress terms are obtained directly from the complex amplitude and its gradients. This will be discussed in section 2.2. The use of Eqs. (2.12.3) involves many assumptions. The primary assumption is that the flow field can be represented in two dimensions using depth and time averaged 18 equations. This assumption implies that vertical fluctuations of the flow which are zero in the time average can be ignored. Many details of the flow may be missed by this assumption of purely horizontal flow. Of course reducing the problem to a twodimensional flow has the distinct advantage of reducing an essentially intractable problem to a simpler tractable form. Some of the terms in Eqs. (2.12.3) also involve certain assumptions. These terms are the bottom friction, the surface wind stress and the lateral mixing. These are all phenomenon that are not completely understood so that assumptions are made concerning the complex forces that are represented. These assumptions usually involve the use of empirical coefficients. In this study the wind forces are ignored so the assumptions involved with the surface stress will not be discussed in detail here. It would be a relatively simple task, though, to add a wind surface stress to the model, using a quadratic shear stress as formulated by Birkemeier and Dalrymple (1976) and also used by Ebersole and Dalrymple (1979) following the work of Van Dorn (1953) in determining the value of the empirical coefficient. The interaction of an oscillatory fluid flow with a real bottom is very complex. Many physical phenomena may be at work. With sand bottoms there can be percolation into the sand and frictional effects resulting in wave energy dissipation. With mud bottoms where the mud is essentially a dense fluid layer there can be damping and internal waves along the mudwater interface. Wave energy losses due to bottom effects are not included in the model but can easily be incorporated into the wave dissipation subroutine. Resuspension and settling of particles also constitutes a complex phenomenon with density changes and momentum transfers. All of these effects are glossed over in the model by assuming a rigid impermeable bottom. As with all of the models discussed above the bottom friction term is represented in the quadratic form. Noda et al. (1974), Birkemeier and Dalrymple (1976) and Vemulakonda (1984) all used the weak current formulation as proposed by LonguetHiggins (1970a) in which rbs = pF luor IU (2.4) rT = pF us V (2.5) where F is a drag coefficient (of the order of .01) and Iuorbl is the time average over a wave period of the absolute value of the wave orbital velocity at the bottom. From linear wave theory it is found that 2H (2.6) Uo,.b[ = Tsinh kh where H is the wave height and T is the wave period, k is the wavenumber, and h = D is the total water depth. Ebersole and Dalrymple used a more complete quadratic formulation rb = pFti~\t (2.7) where the overbar represents the time average over a wave period and tri is comprised of both the mean current and the wave orbital velocities. The total velocity vector ti is expressed as ut = (U + acos 8)I+ (V + asin )}j (2.8) so that its magnitude is given by \t\ = /U2 + V2 + a2 + 2Ucos 0 + 2Vasin 9 (2.9) The wave orbital velocity a is expressed as u= tm cos at (2.10) where Um is the maximum wave orbital velocity at the bottom which is found to be iH =T sinh kh (2.11) The x and y components of the time averaged bottom friction are thus rb = pF" (U + +mcosocosart). I uld(t) (2.12) pF +c( 2b = r (V + "msin cosat) Ilutd(at) (2.13) 2w  where tiutl is written as Itui = /U2 +V + 6f2C08a + 2Ucos2 co't cos9 + 2Vmcos tsi9 (2.14) The integrals in Eqs. (2.12) and (2.13) are then evaluated by a 16 point summation using Simpson's Rule. In the present study several different formulations are used and compared. The first is the formulation given by Longuet Higgins (1970a) and used by several previous investiga tors. The second is a modification of the LonguetHiggins formulation with partial inclusion of current strength in the friction formula. rb, = pF uo,,bl U + pFI7IU (2.15) rbV = pF I[u,b V + pF1f7\V (2.16) It is thought that this formulation has some value since in the lee of structures where there is relatively weak wave activity, there may still be relatively strong currents generated by the gradients of the setup outside of the sheltered region. The third formulation is similar to that of Ebersole and Dalrymple in that the full quadratic expression is used. However the wave orbital velocities at the bottom are expressed as gradients of the complex amplitude and the resulting integral is determined using an eight point Gaussian integration. It is found that an eight point Gaussian integration of a cosine squared integration over an interval of 2r is accurate to within 1 percent. The quadratic formulation of the bottom friction is used rb = pFuIUi (2.17) where tr is composed of both the mean current and the wave orbital velocities. The wave orbital velocities are expressed as the real part of the gradient of the complex velocity potential $ which is expressed as S= igA(z,y)cosh k(h + z) ,,(2.18) o cosh kh 21 where A is the complex amplitude, w is the absolute frequency, a is the intrinsic or relative frequency following the current, and ho, h and k are as previously explained. Equations governing P and a parabolic approximation solution for A are presented in the next chapter. The total velocity vector Ut is expressed as t L=[U I+R(!) T+[v+R(l)]; (2.19) where !9 indicates the real part, so that its magnitude is given by 2+V2+ ("iZ )+ + [W ,)]2I+[ (2.20) The x and y components of the time averaged bottom friction are thus /T /av pF,: [U () k Ide (2.21) [V= +F It[("y)].I ii^dt (2.22) As in Ebersole and Dalrymple, the assumption is made that only the gradients of A are retained and the gradients of the depth and spatially varying quantities such as k and a are ignored. This is essentially the assumption of a fiat bottom and very slowly varying currents. Figure 2.1 plots the magnitude of the longshore current as a function of the distance offshore for the various friction formulations discussed above. These velocity profiles were obtained numerically for an eleven second wave with a deep water wave height of 1.028 meters (1 meter at offshore grid of the computational domain, using a AX of 10 meters and 30 grid rows) and a deep water wave angle of 17.28 degrees (10 degrees at the model boundary) approaching a 1:25 plane beach. A friction factor, F, of .01 was used for each case. Using a friction formulation that ignores the contribution of the waves results in a very weak friction and thus a large longshore current. Using the Longuet Higgins formulation gives a friction slightly less than the more proper time averaged formulation, while a combination of the LonguetHiggins formulation and the strong current model (no wave contribution) yields a slightly stronger friction. For practical numerical considerations the latter formulation is U' ,  r = pFlU\ ud a. *** *** r = pPFlu"b. ., ....................... r= f for lt l dt U.. /I \  r lpF\\u\+  LUi / W r = pF u u.b +iY] a ,r. / '/* . POINT I )S o 50. 100. 150. 200. DISTANCE FROM SHORE (M) Figure 2.1: Longshore current profiles for different friction formulations for an 11 second, 1.028 meter wave at 17.28 degrees to a 1:25 plane beach. most often used since it requires less computation time and more importantly the Gaussian integration involving the real part of the gradient of complex amplitudes at times produces stability problems. The lateral mixing term in Eqs. (2.1) and (2.2) also involve several simplifying assump tions. The first assumption is found in the derivation of the two depth integrated equations of motion where it is assumed that the turbulent shear stress term pu' is independent of depth. The other assumption is that the process of momentum transfer due to turbulent fluctuations can be represented by a product of a mixing length parameter and derivatives of the mean current. This is what is done in the present study and also in previous studies. Ebersole and Dalrymple (1979), Vemulakonda (1984), and Yan (1987) all used the formu lation given by LonguetHiggins (1970b) in which he developed an analytical expression for the cross shore distribution of the longshore current velocities using a linearized friction formulation as described above, linear waves at an angle to the normal of a plane beach, 23 and a spilling breaker assumption. LonguetHiggins wrote = ( aU aV) (2.23) in which e, and ey are mixing length coefficients. LonguetHiggins reasoned that the mixing length parameter e, must tend to zero as the shoreline was approached. He assumed that e, is proportional to the distance from the shoreline, X, multiplied by a velocity scale which he chose as the shallow water celerity fgTi. Thus e. is written as ez = NX Vg (2.24) where N is a dimensionless constant for which LonguetHiggins chose the limits as 0 < N < 0.016 (2.25) Obviously, for physical reasons there needs to be some limit on the scale of these eddies. Ebersole and Dalrymple arbitrarily chose the value of er near the breaker line to be the maximum value and for e. to be constant offshore from the breaker line. The same limit upon ex is used in this study. The coefficient ey is assumed to be a constant and it is arbitrarily assigned the maximum value of e,. In the present study the further assumption is made in computing rT that the cross derivatives are negligible (a common assumption in fluid mechanics). Therefore, the lateral mixing terms become 1 (a 2 aU) (2.26) in the x direction and a(e av) (2.27) p az aZ \ aZ in the y direction. This means that the lateral mixing term is represented by a purely lateral diffusion of momentum. That the lateral mixing term is a lateral diffusion term has some significance in that the diffusion term is expressed explicitly and thus there results a time step limitation to an implicit model. This will be discussed in section 4.1 dealing with the finite differencing of the equations governing the circulation model. 24 2.2 Radiation Stresses In a series of papers in the early 1960s LonguetHiggins and Stewart (1962, 1963, 1964) introduced and developed the mathematics of the radiation stresses. These stresses which can be conceptually thought of as the "flux of excess momentum" can be obtained through integration of the momentum equations over the water column as is shown in Appendix A. Using linear wave theory the radiation stress components are found to be S,. = E[n(cos'2+l)1] (2.28) S = E n(sin29+1) 1 (2.29) SV = Encosgsine (2.30) where E is the wave energy equal to IpgH2, n is the ratio of the group velocity to the wave celerity and 6 is the angle the wave rays make with the x axis. If the wave height and wave angle are known at each grid point, the radiation stresses and their gradients are readily determined for use in the governing equations. To determine the wave height and angle at each grid point Birkemeier and Dalrymple (1976), Ebersole and Dalrymple (1979) and Vemulakonda (1984) each used the following scheme which is based upon the irrotationality of the wave number and an equation for the conservation of wave energy. The irrotationality of the wave number determines changes in the wave angle due to refraction while the energy equation determines changes in the wave height due to shoaling and refraction. This scheme was first introduced by Noda et al. (1974). The irrotationality of the wave number V x k=O (2.31) where I= Ik[cos 01+ IkI sin 89 (2.32) is expanded in Cartesian coordinates to obtain a# I ae 1 alkl alk coas + sin O+ sin  cos 0 (2.33) 9x ay \k\ jx IkI ay 25 For waves on a current the dispersion relation is w = o + = UJglkl tanh klh + Ikl coso + Ikl sin 0 (2.34) Equation (2.34) is differentiated to eliminate MI 1 and m from Eq. (2.33). Using a forward difference derivative in x and a backward difference derivative in y, Oj is obtained as a function of cos 9, sin 9, U, V, k,h and their derivatives as well as 0 at the adjoining grids. The sin 8 and cos 9 terms are expressed in terms of the surrounding four grids using a 2nd order Taylor expansion. The solution is an iterative process alternating solution of the 0 equation with solution of the dispersion relationship (2.34) for k = \~I which is accomplished by a NewtonRaphson method until acceptable convergence is obtained for k and 9. Next an energy equation is solved for the wave height. A conservation of energy equation is given by Phillips (1969 ed. Eq. 3.6.21) as aE Ba auS + {E(U,+C,) + SI =0 (2.35) Assuming steady state conditions and expanding in Cartesian coordinates gives (U + C,coso) 1 + (V + C, sin ) + (U + Ccos ) aU CU 09V av + (V + Csin90)+ + ,Y +f~~,~~~Y = 0 (2.36) where &ff = Es,, = n(cos2 + 1)  SE = 1 1 vyy = ISvy = n(sin0+1 ) &V = rZ= z= = Y = ncosossinV Defining QiJ= ~{(+ + C, cos ) + (+Csin9)+ (2.37) where ?=a Vz + fzV lay.  26 and noting that U,V,C,, and a are known quantities; taking forward difference derivatives in x and backward difference derivatives in y and solving for the wave height, Hij, at the i,j grid yields (V +C, rin)H.,i_i (_ +C co. )i.+lj AY AX AX2. The computation is a row by row iteration proceeding from large i (offshore in the coordinate system of the above cited authors) to small i. During the computation of each Hij the breaking height is also calculated and if Hij exceeds this value, then the breaking wave height is used. Several breaking height indices are possible. Miche (1944) gave as the theoretical limiting wave steepness the condition Hb Hb= .142 tanh(kbhb) (2.40) where the b subscript indicates breaking conditions. Le MWhauth and Koh (1967) indicated from experimental data that the limiting steepness is Hb T = .12tanh(kbhb) (2.41) or Hb = (~2'.12tanh(kbhb) (2.42) This is the breaking height criterion used by Noda et al. (1974). Perhaps the simplest breaking height criterion is the spilling breaker assumption Hb = .78hb (2.43) which is used by Vemulakonda (1984) and Vemulakonda et al. (1985). The scheme described above of determining the wave angle 8 from the irrotationality of the wave number (2.33) and the dispersion relation (2.34) and the wave height from the energy equation (2.35) has its limitations. Principally this scheme does not take into account the effects of diffraction. This limitation is surmounted by use of a parabolic wave equation. However this presents another problem in that the wave angle is not determined 27 in a parabolic equation solution. Only the complex amplitude, which contains the mag nitude and the phase of the wave, is obtained. Yan's (1987) solution to this problem is to obtain 0 directly from the complex amplitude by assuming that the slope of the water surface resulting from the waves indicates the direction d wave propagation. This value for 0 and twice the absolute value of the complex amplitude for the wave height are then used in Eqs. (2.282.30) to obtain the radiation stresses. The value thus obtained for 0 is then also used in the dispersion relationship. For plane waves this works fine; however, in diffraction zones where shortcrested wave are encountered problems may arise, since the instantaneous slope of the wave at a particular location can be other than in the direction of wave propagation. In this report the radiation stresses are obtained directly from the complex amplitude and its gradients using equations developed by Mei (1972) and repeated in his book (1982). The basis of the derivation is to use the definition of the velocity potential to substitute for the wave induced velocities in the definition for the radiation stresses. For an inviscid incompressible fluid the instantaneous equations of motion are (_L+ T) = Vp pV, (2.44) V. = 0 (2.45) where q represents the velocity components, p is the pressure and i, is the z or vertical unit vector. Denoting the horizontal components of qas ua, a = 1, 2 and the vertical component of (as w, the boundary conditions are p = 0 (2.46) + usp = (2.47) on z = ), where q is the instantaneous free surface, and the bottom boundary condition up W (2.48) q is assumed to be composed of a time averaged part and a fluctuating part q3; i.e. = ( zy, ) + (z, y, z, t) (2.49) 28 Integrating a horizontal component of (2.44) vertically over the depth of the water column from z = h to z = q, using the boundary conditions (2.462.48) and then taking time averages (the same procedure as in appendix A) the following definition of the radiation stresses is obtained: Sap =p t up'dz + 6 \ pdz .u(f + h)' + (7 T (2.50) Jh afh 2 1 where the overbar is used to indicate the time averaged values and 6ap is the Kronicker delta function 1 for a ofl (2.51) 1ap  0 for a 6 Mei (1972) was able to express Eq. (2.50) in terms of the general complex amplitude function B by relating it to the velocity potential and the free surface height by = _iBCOsh k(ho + z) ,t (2.52) o cosh kh S= Beit (2.53) An expression for p is obtained by integrating the vertical component of Eq. (2.44) aw aw aw ap P + pusPW + PW z P9 (2.54) Integrating from any depth z to the surface, using Leibnitz's rule for interchanging the order of differentiation and integration, employing the two surface boundary conditions, Eqs. (2.46) and (2.47) and the continuity Eq. (2.45), Eq. (2.54) yields: p(z,y,Z,t) = pg(,r z) + p7f ,wdz+p uwz (2.55) Taking time averages gives )= z)+ d p(.)2 (2.56) Substitution of Eqs. (2.52),(2.53) and (2.56) into (2.50) and making the assumption that horizontal derivatives of k,h and a can be neglected (this is essentially a mild slope as sumption), the following is obtained after straightforward but tedious calculations which are shown in Appendix B: s 09 RP aBoB \ 1 + 2kh + 2 (.( ac_ 6{ (k ct2 sinh 2kh B [iIs2h + (IVB12 B12) khct kh 1) (2.57) where B* denotes the complex conjugate of B. In obtaining Eq. (2.57) use was made of the flat bottom Helmholtz equation VB +k B= (2.58) The general complex amplitude function B is related to the complex amplitude A which is solved for in the wave model by B = AeI(f l Idz) (2.59) There are some problems with Eq. (2.57). One may question the assumption of a fiat bottom; however the inclusion of the slope and gradients of the wave number and the intrinsic frequency leads to intractable mathematics. The use of the localized flat bottom is also implicit in the derivation of LonguetHiggins and Stewart to obtain Eqs. (2.28 2.30) since in carrying out the integration involved it is assumed that horizontal derivatives of the velocity potential are limited to derivatives of the phase function. For plane waves using a spilling breaker assumption, results using Eq. (2.57) differ little from those using Eqs. (2.282.30); however in diffraction zones Eq. (2.57) is superior for reasons mentioned previously. The major problem with the use of Eq. (2.57) is also a major problem associated with the use of LonguetHiggins's formulation. This problem is that in both formulations the initiation of wave breaking implies the dissipation of energy and thus the forcing of set up and longshore currents starting at the breaker line. The physical reality is somewhat different in that in the initial stages of breaking, energy is not immediately dissipated but rather is transformed from organized wave energy to turbulent energy and to a surface roller (Svendsen 1984). There is a space lag between the initiation of breaking and the dissipation 30 2.3 BEACH. 20 01 MEAN WATEI LEVEL, V 00 0 1 05 0 EXPERIMENT 05 / 05 BREAK / POINT I PLUNGE POINT ENVELOPE OF WAVE HEIGHT 6 BECH WAVE CREST,... BE 4 400 300 200 100 0 DISTANCE FROM STILL VWTER LINE ON BEACH, x (cms) Figure 2.2: Experimental results for wave setup and setdown. Reprinted from Bowen, In man, and Simmons (1968), Journal of Geophysical Research, vol. 73(8) page 2573, copyright by the American Geophysical Union. of energy. This can be seen in figure 2.2 which gives the results for the experiments of Bowen et al. (1968). Analytical and numerical solutions for the mean water surface elevation will give a setup starting at the break point while the solid line in the lower graph shows that setup measured in the experiment begins shoreward of the break point. This will be discussed further in Chapter 5. CHAPTER 3 WAVE MODEL 3.1 Governing Wave Equation Several methods have been used to derive the governing equation for the propagation of waves over varying topography and varying currents. Among the methods used are per turbation techniques, multiple scales, a Lagrangian and the use of Green's second identity. With the assumption of a mild slope and weak current conditions each of these methods should yield the same equation. For illustrative purposes the method of using Green's sec ond identity will be used here to derive a linear equation. A Lagrangian derivation of the same equation is presented in Appendix C. The linear equation is used since second order waves exhibit singularities as the water depth goes to zero. A velocity potential 4 is assumed such that the water particle velocities are given by V0, where V is the threedimensional gradient operator v +a a,+ C1 (3.1) v= oy aZ and ,3j, and k are the unit vectors in the xy, and z directions, respectively. The velocities are a superposition of steady currents and wave induced motion. The velocity potential is given by = =o + efi (3.2) where 0o is defined such that its gradient yields U, V, and W which are the steady current terms, and f4 is the wave component of the velocity potential. An undefined scaling e is used so as to separate the wave and steady current part of the velocity potential. The wave part of the velocity potential is composed of two parts: f, which contains the depth dependency and 4, which is a function of the horizontal coordinates and time. The depth dependency f is given by linear theory as f = cosh k(h + ) (3.) cosh kh where ho is the depth of the bottom referenced to the still water line, and h is the total water depth and k is the wave number vector. It should be noted that f is also a function of the horizontal coordinates since k and ho vary with horizontal position. For incompressible, irrotational flow the governing equation for the velocity potential is the Laplace equation V20 = 0 ho < z _< q (3.4) where rl is the instantaneous water surface elevation. Using a horizontal gradient operator Vh defined by a a. Vh TS + (3.5) the Laplace equation is written as V10 + 0 = 0 ho _< z << r (3.6) The governing equation is complemented by the following boundary conditions. On the free surface the kinematic free surface boundary condition is given by + Vn.V1 = 0Z=^ (3.7) The dynamic free surface boundary condition with zero atmospheric pressure O, + 2(V)2 + gZ = 0 z = z 7 (3.8) and the bottom boundary condition Vjh" Vhho = 0. z = ho (3.9) The derivation technique is to write Green's second identity for 4 and f as defined by Eqs. (3.2) and (3.3) respectively. fh f dz of, dz = ff fL. (3.10) 33 Then Eq. (3.6) is used to substitute for 0,, in the first term of Eq. (3.10) and the boundary conditions are used to substitute for 4, at the mean free surface and at the bottom in the third term of Green's identity. Equation (3.2) is used and terms of order e only are retained in order to obtain the governing equation for 4. In order to obtain 4, at q it is assumed that F1 = q + eq. (3.11) The kinematic free surface boundary condition, Eq. (3.7), is then expanded in a Taylor series about z = j, a[L"~ ] [+ + h v17 0] = (3.12) solving for f, q using Eqs. 3.2 and 3.11 while retaining only terms of order e yields Sq + U" V + V1, Vi W, (3.13) A total horizontal derivative D a Dt = t + V" V (3.14) is introduced so that the first two terms on the right hand side of Eq. (3.13) can be written as the total horizontal derivative of qi. The final manipulation of Eq. (3.13) is to recognize that a at the surface is the horizontal divergence of the horizontal current (= Vh U). Thus Dt .i<= + Vh. Vhq + ijVh (3.15) ij is eliminated from Eq. (3.15) by expanding the dynamic free surface boundary condition, Eq. (3.8), in a Taylor series about z = fj. [+ 2(V)2+z] = gq + [+ 0 (V)2] 1. ~Z S=F J =i 2 J= 19fC 2 Tza" J.1 Using Eq. 3.2 to substitute for 4 the order e terms yield a solution for i I= pt + *vh 1 = D 9 g Dt = 0 (3.16) (3.17) Substituting into Eq. (3.15) gives ID2 + VI. VI (3.18) g1 DO g Dt The first term of the Green's identity using Eq. (3.6) becomes flk f.dz= f fV24dz= f fV, (VIh)dz (3.19) Now Vho = I + ef Vh% + V f (3.20) and Vh (Vh,) = Vh. U + eVh (f Vh;) + teVh. Vhf + eCV^f (3.21) Retaining only terms of order e and dropping the last term since it is second order in the slope the integral is equal to f fVh. (fV )dz l fVhVhf dz (3.22) These two terms are handled using Leibnitz's rule of integration. Vh 1f(f VA)dz 4 Vhf (fVh=)dz + 7f Vh(fVh;)dz fho *~ho fho +Vhqf2 In Vh; + Vhhof2 h, VAh (3.23) Using this relationship and recognizing that f q is equal to one and the integral over the depth of f2 is CC,/g, where C is the wave celerity, C. is the wave group velocity and g is the gravitational constant, the first term is equal to Vh (CCOV ) + Vhq Vth + Vhhf2 Iho Vh (3.24) 9 where the e has been dropped. The second term of Eq. (3.10) is evaluated quite easily. fho fdz f f z e fcf..dz (3.25) 35 The first integral above is dropped since there are only steady current terms. The remaining term, dropping the e, is = f a f.dz j= k f2dz = k (3.26) $ho ho 8 The third term in Eq. (3.10) is f1. 0 h.= f I 14 I,o (3.27) The first term on the right hand side is equal to 4, s,1 since f lq= 1 and has already been determined. The second term on the right hand side of Eq. (3.27) is easily determined using the bottom boundary condition Eq. (3.9) 1'. Iho= fVh. Vhho Ih, (3.28) Using Eq. (3.20) this becomes f, Ih.= f 0. Vhho I_h, +ef2 Ih. Vh Vhho + efV hf V*hho It. (3.29) The first term is dropped since it not order e and the third term is clearly higher order in the slope. Thus, again dropping the e / Ilh.= f2 Iho V'. Vhho (3.30) This term cancels the bottom boundary term resulting from the Leibnitz integration above. The fourth term of Eq. (3.10) is also easily determined. sinh k(h0 + z) o"2 Ofo i= ' k Ih= k4 k tanh kh I=  (3.31) Evaluating 4 at the mean water surface and retaining only terms of order e yields h=  (3.32) Collecting terms from (3.24),(3.26),(3.18),(3.30) and(3.32) yields the following linear governing equation for 4. D2 + V ) ( + (Vh+ (2 k CC,)j = 0 (3.33) This equation was first derived by Kirby (1984). 36 3.2 A Parabolic Approximation Equation (3.33) is a hyperbolic wave equation. It can be altered to an elliptic equation by assuming that any time dependency is solely in the wave phase. Before doing this an energy dissipation term DD is introduced as follows. It is assumed that Eq. (3.33) holds for the case of no energy dissipation; then for the case with dissipation the left hand side is balanced by any energy that is dissipated. The governing equation is then written as D2 D;_ "t2 + (V d) v ,(cc,v5 ) + (o' kCcc,) = DD (3.34) In the derivation of the governing hyperbolic equation using a Lagrangian as presented in Appendix C the DD term can be viewed as a virtual work term. Following Kirby (1983) the DD term is assumed to be of the form DD = w (3.35) Dt where w is an undefined coefficient indicating the strength of the dissipation. Eventually the w term will be related to the energy dissipation due to wave breaking following the work of Daily, Dean and Dalrymple (1984). Using Eq. (3.35) in Eq. (3.34) and making the assumption that the only time dependency is in the phase, i.e. S= (3.36) reduces the hyperbolic equation (3.34) to the elliptic form 2iwU VA, + 7. V1(O? V^h) + (Vh, 1)(17 V1;) V1. (CC,VAh) +({ CCo iW(Vh* 1)}d = iow$ (3.37) where only the phase contribution to the horizontal derivative of 4 is retained in obtaining the term on the right hand side of (3.37). Before deriving a parabolic approximation to Eq. (3.37), the need for such an approxi mation should be discussed. There are two major computational drawbacks to numerically 37 solving an elliptic equation. The first is that solution is required simultaneously for each grid. This means that for a large computational domain, say m by n grids, it is necessary to invert an m x n by m x n matrix. If m and n are sufficiently large, problems may be encountered in terms of the internal memory storage, and even if the available storage is not a problem, the number of computations required to invert the matrix may make for an exceedingly slow solution. In some situations this may be acceptable; however for the purposes of a wave induced circulation model in which repeated solution for the wave field is necessary the program would take a very long time to run. The other disadvantage in volved in the numerical solution of an elliptic equation is that boundary conditions must be specified at all of the boundaries. In many applications the only known boundary condition is the offshore, or incident wave amplitude. For example, the domain of interest may be the wave field in the lee of an obstruction where the only down wave condition is that the waves freely pass through the down wave boundary. Or in the present study in which the wave amplitude is known to go to zero at the shoreline a problem ensues in that the shoreline is (for purposes of the circulation model) at the grid edge and not at the grid point. A parabolic equation eliminates both of these problems. Since the solution proceeds grid row by grid row where each solution uses the results of the previous grid row, no down wave information is required. Since only one grid row is solved at a time, the solution requires only that a tridiagonal matrix equation be solved to obtain values for the grid row. The only required boundary conditions are the conditions on the first grid row (usually the offshore boundary) and lateral boundary conditions which are usually a derivative condition that allows for the wave to pass through the boundary without any reflection at the boundary. There are several ways to transform the elliptic Eq. (3.37) to a parabolic equation. The most direct is to make certain assumptions concerning the length scale of variations in the x direction. At the heart of the parabolic approximation is the assumption that the direction of wave propagation is essentially along the x axis. For waves propagating at an angle to the x axis the proper form of ; is = _iA(x, V)ei(flfkOd,+fk wcdl) (3.38) and the proper form of the dispersion relationship is w = a + kcosoU + ksin9V (3.39) where 0 is the angle of the wave propagation relative to the x axis. The essence of a parabolic approximation is to assume that the angle that the wave rays make with the x axis is small so that the ksin 0 term can be neglected in Eqs. (3.38) and (3.39) and cos 0 is assumed to be unity. The manner of doing this varies with different authors. Yan (1987) retains the complete dispersion relation (Eq. (3.39)) but deletes the y dependency in Eq. (3.38) in the derivation of his governing wave equation. He solves for 0 by determining the slope of the instantaneous water surface. For long crested waves this method is satisfactory although requiring additional computations. However, this method may produce erroneous results in diffraction zones where the presence of short crested waves will produce a horizontal com ponent of the slope normal vector that is different from the direction of wave propagation. Kirby's (1983) method is to make the assumption that = ig A( 'cl k z, (3.40) so that the if kin Idy part of the phase function is absorbed into the amplitude function, A. Using Eq. (3.40) directly in Eq. (3.37) and further assuming that derivatives of A, are small compared to derivatives of the phase function (i.e. that << ikA, ) will yield a parabolic equation. For illustration the case of linear diffraction in the absence of currents and in a region of constant depth is examined. Substitution of Eq. (3.2) into the Laplace equation (3.6) yields the Helmholtz equation (3.41) 39 This equation is also obtained directly from Eq. (3.37) for the stated assumptions of no currents, no dissipation, and constant depth. Substitution of (3.40) into the Helmholtz equation gives A,, + 2ikA, + A,, = 0 (3.42) and employing the above mentioned assumption concerning derivatives of A, yields the linear parabolic equation 2ikA, + Ayy = 0 (3.43) Another method of obtaining a parabolic equation is to use a splitting method in which it is assumed that 4 is composed of forward moving and backward moving components that are uncoupled. 4 = + + (3.44) Uncoupled equations for 4+ and ; are developed and since the equations are uncoupled only the equation for i+ is retained. In the following the work of Kirby (1983) and Booij (1981) is closely followed. The second order elliptic equation :a (1fla ) + = o (3.45) can be split exactly into equations governing forward and backward disturbances 0+ an 0 which are given by S= +i+ (3.46) ax Q = i0 (3.47) (3.48) Booij chose the relation (3.49) 40 where f is an as yet unknown operator on 4. Substituting into Eq. (3.45) and neglecting derivatives of f alone leads to the result +( 1 ( ) + =0 (3.50) r / \ + 0 The elliptic equation (3.37) is now put into the form of Eq. (3.50). Letting p = CC, for notational convenience and expanding some of the terms of Eq. (3.37) in their spatial coordinates gives (;4.), + (Py)t + k2p; + (w2 o2 + i*(VI E7)) + io'W (VU,), (v2,), (UVW ), (UV,), + 2iwU Vh, = 0 (3.51) Booij chose to neglect terms involving the square of the mean current components assuming that they are smaller than terms involving p and then arranged the above equation as S+ P'(p + 2iwU) +k' (1+ )p = 0 (3.52) where M4 = (w2 ar2 + iw(Vh j) + iaw)$ + 2iwV>y + (pY)Y (3.53) Kirby points out that this choice by Booij leads to his inability to derive the correct wave action equation. Kirby instead does not drop any of the squares of the mean current components until after the splitting is complete. By assuming that the waves are oriented in the x direction the dispersion relation is written as o = w kU (3.54) leading to the result w, au = 2wkU (kU)2 (3.55) Using this relationship Kirby then arranged Eq. (3.51) as ez+ (p U2)I(p u 2), + (+ k i+ ,M U2) =o0 (3.56) where Mq = (2wkU + iw(Vi,. ) + iow} (UV4) (UV ), +{(p V2),}, + 2iw. .Vh (3.57) Comparison with (3.50) gives 72 = k(1+ k2(p U2) (3.58) t = (p U2) (3.59) or solving for y and t 7 = k(1+k2p M U2))2 (3.60) S= ki(pk( )+ M (3.61) Booij points out that the splitting matrix approach employed by Radder (1979) is equivalent to the choice yR = k l+2k:pUz) fR = k2(pU2)i Using the results (3.60) and (3.61) in Eq. (3.46) along with the definition (3.49) gives the parabolic equation i( (p1+ U ) += ikk*( (1+ (pU2) s k2( U2) By expanding the pseudooperators in the form (+ k2(p U2) (1+4k(p U2)) + (3.63) k2(p_U2) j k2+p l+ p ) (3.64 the result of both expansions may be summarized in the general equation Ia { (+ PM + P2M3 10 k1(p U2)2}=;(L M')  ik21 (3.65) T1 k2( U2) 1=kkrU) k2(p U2)I~ where D to Radder Bboij (3.66) and where P2 = P1 + j for both approximations. Kirby chose the approximation with P1 = 0 since for his derivation which was for weakly nonlinear waves it should correspond to the information contained in the nonlinear Schr&dinger equation. With P1 = 0 Eq. (3.65) becomes k(p U2)p + ~ U{2) + M (3.67) where Mr+ is defined by Eq. (3.57). An equation for the complex amplitude A is then obtained by making the substitution = (A)e,f (3.68) which after dropping squares of the components of the mean current results in (C,+U)A+ + C+U A+()VA, + ( A 01 )z or 2 or +i[ o r(+)] + A=o (3.69) Kirby then further modified the equation by a shift to a reference phase function using the relation A = A'ei(Ltf*kd) (3.70) where k can be taken as an averaged wave number, resulting in (C,+ U)A + i( k)(C, + U)A' + CI U) ^^^^/'h[CC# (A')&,] +(' VA+ + (V A' ,, + A'=0 (3.71) The shift to a reference phase function is intended to incorporate the phase information in the complex amplitude so that the instantaneous water surface can be easily recovered in the form rl(z,y)= R {A'e'f!}Ed (3.72) 43 where 9R designates the real part. This is an important consideration especially when ob taining the radiation stresses. It is important that in obtaining Eq. (3.69) that squares of the products of the current components be retained until the last step. The reason for this is seen in that terms such as kUV and kU2 are combined with terms of wV and wU using the dispersion relationship to give terms of aU and aV. Equation (3.71) has been successfully used by Kirby in the form given and with modifi cations. The modifications are the addition of a nonlinear term which makes the equation applicable for weakly nonlinear waves. However for present purposes Eq. (3.71) presents a certain problem which will shortly be corrected. The problem is that for waves propagating at an angle to the x axis it is not possible to obtain a correct form of the conservation of wave action. The case of no current on a beach with straight and parallel contours is given as an illustration. The equation for the conservation of wave action reduces to the conservation of wave energy flux. In mathematical terms this is [C cos a2] =0 (3.73) where a is the wave amplitude. For waves approaching the beach in a direction normal to the beach Eq. (3.71) reduces to C,A= + J(C,).A' = 0 (3.74) Multiplying Eq. (3.74) by the complex conjugate of A', and add to it the complex conjugate equation of (3.74) multiplied by A', yields [c,AI'] =0 (3.75) which is Eq. (3.73) for the special case of cos = 1. However it is not possible to obtain Eq. (3.73) from Eq. (3.71) for the general case of waves at angle 0 on a beach with straight and parallel contours. The importance of this inability to give the correct form of Eq. (3.73) is that if energy flux outside the breaker line is not conserved then gradients of the S,y radiation stress 44 will generate currents offshore of the breaker line in the circulation model. Although such currents in general will be of a small magnitude, they can cause some numerical computational problems. The problem though is easily corrected by using as the dispersion relationship w = o + kcos0U (3.76) and the substitution $= _A(Z, V) j(f koe0d ) (3.77) Rather than using a splitting method, the procedure will be to go directly from Eq. (3.51) using Eq. (3.77) and writing w2 _2 = 2wk cos 0 U (k cos 9 U)2 (3.78) in place of Eq. (3.55). Dropping the ()xx term and all terms containing products of the current components that cannot be reduced as explained above results in the following parabolic equation. (C, cos9 + U)A, + o os A+]VA + (V A tkC,(1 cos2 O)A +CC(a + = 0 (3.79) 2 [C? IVJ A  Now adding a phase shift A = AY'( f'oldz f tkcodz) (3.80) so that the free surface can be recovered by rl(X,Y) = A'e(f cidz) (3.81) where 0 is the wave angle obtained by Snell's law in the absence of any obstructions or y variation of the depth or current. The following governing equation is thus obtained. (c, cos + U)A +i(cos cos)(C cos + U)A' + a C cos + U +VA' + A' kC,( 2 $)A' iC,( + A'=O (3.82) (V) A! Lk,( WA (382 45 Details of the algebra and assumptions involved in obtaining Eq. (3.79) are given in Ap pendix D. For computational purposes it will be assumed that cos = cos 9. From Eq. (3.82) the correct form of the conservation of wave action can be obtained. For the case of no y variation of the wave number Eq. (3.82) reduces to (C, cos + U)X. + COS+U X+VA+ () XA rC^^I V a^^, 9)/ jkC,(1 cos')A' [c,() + A = 0 (3.83) Multiplying Eq. (3.83) by the complex conjugate of A and adding to it A times the complex conjugate equation of (3.83) yields [(C.cosg+ 12U) ] + [V4IA1] =0 (3.84) which is a closer approximation to the exact expression which should have a (C, sin 0 + V) term instead of V in the y derivative term. 3.3 Energy Dissipation Due to Wave Breaking Kirby (1983) found a method of relating the dissipation coefficient w to the work of Dally (1980) who proposed that the decay of energy flux in the surf zone was proportional to the amount of excess energy flux. The excess energy flux is the difference between the actual flux and a stable energy flux. In mathematical terms this relationship is expressed as (EC,) = [Ec, (EC,).] (3.85) where EC, = the energy flux E = wave energy = 1pgH2 C, = the rate of energy propagation which is " ok K = constant h = the water depth 46 (EC,). = the 'stable' energy flux. H = the local wave height The stable energy flux is obtained by assuming it to be a function of the height of a wave propagating over a flat bottom after breaking ceases, and that this height is a constant times the water depth, (0) = '7. The constants K and were determined by analysis of experimental data obtained from Horikawa and Kuo (1966). Kirby related Daily's work to the dissipation term w by deriving a conservation equation for the wave action E. Garrett (1967) and Bretherton and Garrett (1968) showed that for waves propagating in a moving and variable medium, the quantity that is conserved is the wave action (f), which is the wave energy divided by the intrinsic frequency. Kirby derived a wave action conservation equation from the linear hyperbolic wave Eq. (3.34) using the definition of W as given by Eq. (3.35): D2, + (Vh U)D# Vh(CC,Vh) + (o2 kACC,) = i;.w (3.86) Kirby follows the approach of Jonsson (1981) assuming that ; can be represented as i4 = Reie (3.87) where R and e are real. From linear theory it is obvious that R= IAIj. (3.88) To lowest order the derivatives of E are approximated as Ot = w (3.89) Vhe = lc (3.90) Substituting into (3.86) gives the complex equation D2R DR Dt2 + (Vh U) D Vh(CC,VhR) i(atR + 2 i {v, (fo})R + 2(oU) VhR + V, (kCC,)R + 2kCC, VAR} iorwo = 0 (3.91) 47 Since R and ( are real, the real and imaginary parts of (3.91) must individually equal zero. Taking the imaginary part and multiplying by R gives (oR'2)t af VhR' V.* (orf)R' CCk VhR2 Vh (kCC,)R owR2 = 0 (3.92) which using the substitution kCC, = aCg is written as (oR'), + V (oR2(0 + d ,)) = owR2 (3.93) Using (3.88) R2 may be replaced by R2= gA1' 2g (E) (3.94) 0' o \or} which yields the wave action conservation equation (E t + V. (E (0+d'I =W(E) (3.95) Assuming a time steady wave field and neglecting currents, (3.95) becomes (EC,), = wE (3.96) Comparing (3.96) with (3.85) and noting that (C,), = C., w is written as K E \ KC( (H (3.97) where H is the local wave height given by 21AI. In obtaining Eq. (3.97) Kirby neglected the effects of the currents. However, the same expression (3.97) is obtained when currents are included in the derivation. Dally (1987) extended his work to investigate wave breaking on currents and proposed the following equation: [(C, + U) ()] = [EC, (EC,).] (3.98) For the case of a steady state wave field propagating in the x direction the conservation of wave action Eq. (3.95) reduces to [(C + ,) ] (E) (3.99) 48 Setting the right hand sides of Eqs. (3.98) and (3.99) equal to each other results in Eq. (3.97). At first glance it may appear unusual that the same expression for the dissipation coefficient results for the case of currents as for the case of no currents. This can be explained in part by the fact that Daily's work was with waves where breaking was depth limited. For the case of current limited breaking (for example a stopping current), another expression is necessary. In this study Eq. (3.97) will be used, since in general the current field will be moderately weak. Such an assumption of weak currents was invoked in the derivation of the governing parabolic equation. CHAPTER 4 FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDITIONS 4.1 Numerical Solution of the Governing Equations The finite differencing of the governing equations in the circulation model is derived by a matrix analysis as suggested by Sheng and Butler (1982). A rational analysis of the complete Eqs. (2.12.3) including all the nonlinear terms is not possible. However guidance can be obtained through an analysis of the linearized equations. For this purpose the linear equations for the special case of a flat bottom and without any of the forcing terms will be examined. aq au av + = ++ = 0 (4.1) au +g = (4.2) av +9 =0 (4.3) at Ty Let W be a vector of q,U and V W= U (4.4) V Using Eq. (4.4), the governing equations are written as Wt + A W, + BW, = 0 (4.5) where A and B are the following matrices 0 D 0] O 0 D] A= g 0 0 and B= 0 0 0 (4.6) 0 0 0 g 0 0 To find an implicit finite difference scheme for the governing equations, Eq. (4.5) is finite difference as S + AWn+1 + BW+'1 = 0 (4.7) AT 50 This equation is solved for W"+1 in terms of W" (1 + 2A, + 2AX)Wn+l = W" (4.8) where A=, 1AAT6 and Ay = B AT (4.9) 2 AX 2 AY 4AA,W"n+l is added to the left hand side of the equation and 2AYWn is added and sub tracted to the right hand side of Eq. (4.8) so that the equation can be factored to get (1 + 2A,)(1 + 2AX)W"+1 = (1 2Ay)Wn + 2AyWn (4.10) An intermediate value W* is introduced so that the above equation can be written as two equations (1 + 2A.)W* = (1 2AX)Wn (4.11) (1 + 2AJ)W"+' = W* + 2AXWn (4.12) Note that by eliminating W* from the two Eqs. (4.11) and (4.12), Eq. (4.8) is recovered with the addition of an error of 4AAyAWn+" which has been introduced. Writing out the component equations of Eq. (4.11) gives A T AT ( r* + DAS6zU =qn DT Vn (4.13) A~X AY 6T U* +g6S t U= (4.14) V* = Vn gAT6,n" (4.15) (4.16) While writing out the component equations of Eq. (4.12) gives AT A^T qn+l + DA yVn+l q* + D A6vn (4.17) AY AY Un+l = U* (4.18) v+ + Tg + = V* +g,+1 (4.19) 51 Using Eq. (4.15) to eliminate V* in Eq. (4.19) and (4.18) to eliminate U* in Eq. (4.14) the following are obtained. For the x direction AT AT ** + D SU+1 = n D TSyV V (4.20) Un+1+ gAT S ^* = U" (4.21) and for the y direction + AT Vn+l AT @A+Y + DD6vV+ =Y (4.22) AT Vn+I +g ATV  (4.23) The above equations are numerically stable for any combination of AT, AX, and AY. However as a circulation model they are incomplete in that they do not contain any of the forcing terms that drive the circulation. In particular in this model wave forcing represented by gradients in the radiation stresses and bottom friction are of concern. Bottom friction is especially important because at steady state conditions gradients in the radiation stresses must be balanced by either gradients in the water surface elevation or by bottom friction forces. In the absence of friction forces the model could yield continuously accelerating currents. The method of adding the forcing terms is somewhat arbitrary in that the rational derivation of the finite differencing scheme of Sheng and Butler cannot be expanded to include these essentially non linear terms. The addition of these terms to the model also limits the latitude of the ratio of AT in a manner that is not analyzable. That both the friction and the radiation stresses are nonlinear should be recognized. The nonlinearity of the velocity friction term is quite evident since it is represented by the quadratic term pFIlIUa. The radiation stress terms, which represent excess momentum flux induced by waves, are functions of the wave height squared, based on linear wave theory. In shallow water the wave height H is a function of the total water depth, h. + q so that H2 is a function of q2. All the forcing terms and the nonlinear acceleration terms are added at the present or nth time level with the exception of the bottom friction term. The quadratic friction term 52 is finite difference as pFl nfUn+l and the friction term as suggested by LonguetHiggins is difference as pF ju.bl U"n+, where juorbl is in terms of values at the n time level. However, when using the integral quadratic expression, it is evaluated entirely in terms of the present nth or known values. The nonlinear advective acceleration terms are added at the nth time level using central difference derivatives. In the x direction momentum equation these terms are Ui U+j U j + (j + Vj + Vj+ + Vj+) U+l U (4.24) 13 2AX 2AY and in the y direction equation they are I (Uij + Uiij + Llij+1 + Uilj+ll+lj Vi,1 + VijVij+1 V141 (4.25) U, 2AX 2AY (4.25) In the model of Vemulakonda (1984) which is based upon the WIFM (Waterways Ex periment Station Implicit flooding Model) (Butler 1978) model, a three time level leapfrog scheme is used. This gives the advantage that all of the forcing terms (and the nonlinear acceleration terms) are added at the time centered level. However this has some disad vantages in that additional computer memory and computational time are required and leapfrog schemes can exhibit instabilities. In the present study a two time level scheme is used and the wave forcing terms as well as the nonlinear acceleration terms are added at the present time level. Since, for all of the applications in this report, the model is run until steady state conditions are obtained, it is reasoned that there is no great necessity to properly time center all of the terms. If the use of the model is extended to time dependent applications, this may need to be corrected. 4.2 Finite Differencing of the Parabolic Wave Equation Dividing Eq. (3.82) by (C, cos 0 + U) yields [CiCosG+U1 A. +i(kcosik osA)A'+ CAo+U + V A' 2(Ccos# + U) (CCos9 + U) ______ (V A .kC,(1co8'9) , 2(C.cos +U) c( ), + 2 (C, co +U) ^2(C,+cos+ U) [+UA'= 0 (4.26) 2(C co 0 IT CCR(Ao') V 2(C,; Cos 0 + Uj 53 The following notation is used to describe the finite differencing of Eq. (4.26): F' represents the value of F(i,j) where F is any term and i and j designate the grid row and grid column. Because of the staggered grid system used for the circulation model as previously explained, in which the velocity components are specified at the grid edges, the velocity terms U and V are averaged across the grid so as to produce a grid centered term for use in the wave model. Thus U' in the wave model is actually (U + U1+I) and similarly for V and other grid designations. A CrankNicholson finite differencing scheme is used for Eq. (4.26). All terms not involving an x derivative are an average of the value of the term on the i grid row and on the i +1 grid row. Terms having y derivatives are centered at i, j on the i row and i +1, j (A'+IA') on the i + 1 grid row. The one term with an x derivative is written as 'AX 1. Thus all terms are centered at i + 2,j. The finite differencing of Eq. (4.26) is written as A'+I A' AX +' [(,+' cos 8'+ k+' cos ,+')A, + ( cos 1' kj cos #)A;] +2 1{ [(C, cos 9 + U)/]j1+' [(C, cos + U)/o]' ( + A.) +2AX [ (C,cose + U)/O]1 + [(CcoBs + U)/VO J Vau \i [Ai+Ij A (tl Va ______ A _ +2 (Ccos+U)) 2AY \(Ccos9+U)) 2AY J +co + u tV a [(V/)i+ (V/a)1+ Ai] I { [_ ___o)_k k, [ (cos )1 K 1s ] 2 (C (\ A)]+1 ( + U' 24 + +1 4 2(C, cos ( + + MY + [1 ( wuoA) 1 (A +I 1 A+ +i + 0 (4.27) (c Cos e+ )'1 (C, cos0+ U). 54 where p = CC, for notational convenience and the y derivative of p (A) is written as f( [(A) ()' + [(A) ( A) P (A)] + + e [ 2(AYap (4.28) In order to evaluate w.1+ by the method to be explained in section 3.3 an intermediate value of A'+1 is obtained by the explicit formulation A'' + i(k cos k f cos P))A; + V [A +l AAI] + Co [(C,o cos 9 + )U) S=(Ccos+C )/2J. 2AY + (2 c) [ ) l 1. i [A(s.cos'e)lA, S2(C cos + 2AY 2 (,cos + U) J C o+)+ w =0 (4.29) 2 (Cpcosf + U)i 2(Ccos + U)l When Eq. (4.27) is written out using Eq. (4.28) and wi+1 obtained through use of (4.29) there results for any particular values of i and j an equation for the three unknown values Aj, A"1, and A'.' The equation is put into the form ai AS+1 + bi A1 + c At' = di (4.30) where the coefficients aj, bi, and ci are known quantities involving the properties of the environment (depth, current, celerity, etc.) and the di involves the environmental properties and the known values of the amplitude on the i grid row. When Eq. (4.30) is written for all the j values, j equal to 1 through N, for a specified i value, the ensemble of N equations takes the form of a tridiagonal matrix equation. For illustrative purposes let N be equal to four. The resulting matrix equation is a a b1 c A d (4.31) a2 b cs2 A.+' b C3 A',+' d 2 (4 .31) a4 b4 C4 A'+1 d4 Ai+1 55 Equation (4.31) contains four simultaneous equations with six unknowns. Two bound ary conditions should be specified. At each boundary a functional relationship is established between the grid inside the domain and the grid outside of the domain. For the wave prop agation model there are two types of boundary conditions that are readily available for use in a computer model. These are a totally reflective boundary condition which would be used to model a laboratory experiment in a wave basin, and nonreflecting noninterfering boundary which will allow waves to freely pass through the boundary. For the reflective boundary, the mathematical statement of the boundary condition is 8A = 0. At the lower boundary this is finite difference as A A = 0 which gives the functional relationship Al = Ao. (4.32) At the upper boundary the functional relationship for a reflective boundary is AN+1 = AN. (4.33) For the nonreflective boundary the mathematical statement is 1A = imA where m is the longshore component of the wavenumber. Since the lateral boundaries are sufficiently far from any diffraction effects, only refraction and shoaling effects are present. Thus using Snell's Law m = ksin e = ko sin e0. The non reflective boundary condition is finite difference at the lower boundary as S A [A, + Ao] (4.34) which gives (1+ inAY) Similarly at the upper boundary the functional relationship used to eliminate AN+1 is 1 + mAY) AN+1 = 2) AN (4.36) The use of either pair of boundary conditions, (4.32) and (4.33) or (4.35) and (4.36) to eliminate Ao and AN+I will then yield an N by N tridiagonal matrix equation comprising N  J.. 4 J=j+1 ~4 vij y J= j1 I=il I=i I=i+l Figure 4.1: Definition sketch of grid nomenclature and coordinate system. equations for N unknowns. Returning to the illustrative example where N is equal to four, the matrix Eq. (4.31) becomes bl + 7la, cl Ai+l di a2 b2 C2 A2+1 d2 (4.37) as b3 C3 A' (437) a4 4+ 7NC4 A+1 d 4[ where 7\ and IN represent the functional relationships as described above on the j=1l and j=N boundaries. Equation (4.37) is solved easily for the unknown A'+1 using a double sweep method. In the computer program the subroutine CTRIDA uses an algorithm of Carnahan, Luther, and Wilkes (1969) to efficiently accomplish the double sweep solution of eq. (4.37). 4.3 Boundary Conditions in the Circulation Model As shown in figure 4.3 a staggered grid system is used so that velocities are expressed at the grid edge and qy and the wave properties are defined at the grid center. So that the wave model will have full grids it is desirable to place the boundaries of the model at the grid edges on all sides, with the exception of the offshore grid row where the wave conditions are an input or boundary condition for the wave model. Figure 1.1 can also be referred to. The offshore boundary condition that is used in the circulation model is that the mean 57 water surface level is not changed by any of the circulation within the model, mathematically ETA(1,J) = 0 J = 1,N (4.38) This condition necessitates that the offshore boundary be sufficiently far from any setup or setdown. Other possible offshore conditions are a noflow condition or a radiation condition. The no flow condition would conserve water within the computational domain so that the volume of water involved in set up would be compensated by an equal volume of setdown offshore. In an open ocean environment the open ocean can be conceived of as an infinite reservoir. However in modeling laboratory experiments a no flow condition as the offshore boundary condition would be the most appropriate. A radiation boundary condition as proposed by Butler and Sheng (1984) imposes the condition (+C = 0 (4.39) at ax This condition allows long waves to freely propagate out of the computational domain without any reflection. This tends to eliminate any seiche within the computational domain which is induced by the startup transients. This condition was applied at one time in the present model but it did not preserve the offshore zero setup condition and was changed in favor of the more rigid zero water surface elevation condition. If the model is adapted for use with time varying wave trains (an adaptation that the model can easily accommodate) then the radiation condition will become imperative. At the shoreline the boundary condition is that there is no flow into or out of the beach. Since the model allows for flooding of the beach the index of the shoreline changes. The model allows for differential flooding along the shoreline, for example, due to differences of wave intensity in the lee of structures. For each grid column in the y direction the last grid is the IWET(J) grid, where IWET(J) is allowed to change with changes of the mean water level. The numerical statement of the no flow boundary condition at the shoreline is J=1,N (4.40) U(IWET(J)+IJ) = 0 58 This boundary condition is used in other models, though not all models allow for differential flooding of the shoreline. Yan (1987), whose model was used to investigate waves and currents around a tidal or fluvial discharge, used the no flow condition away from the discharge and defined U(IWET(J)+I,J) as equal to the discharge per grid divided by the water depth at the grids where discharge occurred. Other conditions are also imposed at the shoreline boundary. These are that the wave height and the longshore velocity are both zero. These conditions are imposed at the shoreward grid edge of the IWET(J) grids even though the wave height and longshore velocities are specified a half grid space form the shoreline. This is done by writing x derivatives of the radiation stress terms and of the V velocities that need to be positioned at the center of the IWET(J) grid in terms of three points. The three points are at the IWET(J)1 and IWET(J) grids and the zero value at the shoreline at the edge of the IWET(J) grid. For the lateral boundary conditions several alternatives are available. Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) both used a periodic or repeating beach. This meant that computations were performed between the second and Nth grid, and values outside of the computational domain are determined according to Q(I, 1) = Q(I, N) (4.41) Q(I, N + 1) = Q(I,2) for all I (4.42) for any quantity Q. At one point this was used in the present model until the problem of what should have been transient edge waves generated by wave obstructions during the startup were found to be trapped in the model. The periodic boundary condition allowed for disturbances propagating out of the domain at one lateral boundary to enter the domain at the other boundary. For straight and parallel contours, the periodic beach condition worked well in the present model but with the presence of breakwaters another condition had to be used. For the lateral boundary conditions there are two choices, either closed or open bound 59 aries. Closed boundaries are easily applied by the condition of no flow across the boundary. Such a choice of boundary conditions would accurately duplicate the physical conditions ex isting in a laboratory basin (or for prototype conditions within a harbor). Open boundaries are obtained by specifying derivative conditions. In the present study an open boundary is obtained by specifying that there are no gradients of the depth and the wave proper ties across the lateral boundaries. This necessitates that the lateral boundaries are placed sufficiently far from any obstructions or any y gradients caused by the obstructions. An alternative open boundary condition would be that there is no gradient of the flow. This would, however, put the position of the boundary condition at the center of the bordering grid rather than at the edge of the grid and could possibly introduce error at the boundary. 4.4 Boundary Conditions in the Wave Model In the wave model only lateral boundary conditions are needed for computational pur poses. The offshore condition is merely the specified wave amplitude and wave angle. No condition is needed at the shoreline since the parabolic wave model marches grid row by grid grow and needs only the lateral conditions. The offshore condition is that each time the wave model is called the wave amplitude on the first grid row is specified as A(1,J) = C(t) Ao(J) J1,N (4.43) where Ao(J) is defined in the input subroutine according to Ao(J) = eikosineoJAY (4.44) where ko is the wave number along the first grid row and 0o is the angle the wave ray makes with the x axis. This offshore condition necessitates that the water depth be constant along the offshore grid row. This is another reason for choosing the offshore condition in the circulation model to be a constant zero. With the presence of wave obstructions and current deflectors such as jetties it is possible that with other offshore conditions the mean water surface elevation could vary along the offshore grid row. The use of one condition in 60 50. 100. 150. 200. 250. TIME (SEC.) Figure 4.2: Startup function according to Eq. 4.43 with C1=30 and C2=3. the wave model which is contrary to the condition in the circulation model can introduce instabilities in the overall model. The interaction of the two models is more sensitive than either of the models individually. The C(t) term in equation 4.43 is a startup function which smoothly goes from zero to unity so as to bring the incident wave height from zero to Ho and minimize any shock loading of the model. Mathematically it is given by C(t) = .5 [I + tanh(Z C,)] (4.45) where Ci and C2 are arbitrary constants to adjust the slope and offset of the curve. Fig ure 4.2 shows equation 4.45 for C1 and C2 equal to 30 and 3 respectively. The lateral boundary conditions are either open or reflective. The reflective condition is expressed as aA S=0 (4.46) while the open condition is expressed as S= imA (4.47) where m is the longshore component of the wavenumber which in the absence of diffraction effects remains constant m = kc sin$o = ksin 9 (4.48) Direction of 0 Wave Propagation Figure 4.3: Waves passing through a lateral boundary unaffected by the boundary. according to Snell's law. Equation 4.46 or 4.47 are used to eliminate the amplitude outside of the domain from the tridiagonal matrix equation. To illustrate equation 4.47 is finite difference at the lower lateral boundary as A(I,1) A(I,0) A(I,1) + A(I,0) AY 2 (4.49) AY 2 Solving for A(I,0) 1 im*AY A(I,0) = m A(I,1) (4.50) 1 + im S Thus A(I,0) is eliminated in terms of A(I,1). similarly at the upper boundary A(I,N+1) is eliminated in terms of A(I,N). The same procedure is used with equation 4.46 which yields the simple relationships A(I,0) = A(I,1) (4.51) and A(I,N+1) = A(I,N) (4.52) 62 It is important that the boundary conditions 4.46 and 4.47 be used as prescribed above. The parabolic equation gives N equations containing the N+2 unknowns A(I+1,J) for J=O,N+1. The two boundary conditions are thus used to reduce the number of un knowns to a number equal to the number of equations. If equation 4.47 is written in terms of A(I+1,1) and A(I+1,2) at the lower boundary and in terms of A(I+I,N1) and A(I+1,N) at the upper boundary and the governing wave equation is written for J=2,N1, then errors may be introduced at the boundaries. This is especially true for the case of a reflective boundary. Figure 4.3 shows waves freely passing through the lateral boundaries without any noticeable reflection for the case of a plane wave with a period of 10 seconds propagating over constant depth at an angle of 10 degrees to the x axis. CHAPTER 5 RESULTS The numerical model combining the iterative interaction of the circulation model and the wave model as described in the previous chapters was tested for several cases in which either analytical or experimental results are available. The model was also run for cases where previous numerical results are available. The two cases involving analytical solutions are for setup induced by normally incident waves on a plane beach and longshore currents generated by incident plane waves at an angle on a plane beach. The first section of this chapter covers setup, comparing the results with the analytical solutions of LonguetHiggins and Stewart (1964) and the experimental results of Bowen et al. (1968) and Hansen and Svendsen (1979). The second section covers longshore currents and compares results with the analytical solution of LonguetHiggins (1970a,b) and experimental results of Visser (1982). The model was run for the case of a ripchannel on a plane beach. Comparisons with the numerical results of Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) are given in the third section. A rip current was also generated for the case of normal waves breaking on a shore parallel bar with a gap. No comparisons are available for this case. The model was run for the cases of shore perpendicular and shore parallel breakwaters represented in the model as infinitesimally thin impermeable barriers. For the shore per pendicular barrier, the results are compared with the experimental results obtained at the Waterways Experiment Station as reported by Hales (1980). For the shore parallel barrier the results are compared with the numerical results of Liu and Mei (1975) who did not include a wave current feedback mechanism. In the final section the model was used to simulate the laboratory experiments of 64 Gourlay (1974) in which the shoreline in the lee of a breakwater was made to curve and meet the breakwater as occurs with tombolo formation. The results of the model were in close agreement with Gourlay's experimental values. 5.1 Wave Setup LonguetHiggins and Stewart (1964) gave a simple but elegant derivation for the wave induced setup on a plane beach assuming linear shallow water theory. For the steady state case of normally incident waves on a plane beach without any y variations and without any currents the depth integrated equations of motions reduce to 9a + I aSz =0 (5.1) ax pD ax In shallow water for normally incident waves the radiation stress term from linear wave theory, Eq. (2.28), becomes s,, = E(n + 1) = = pgH2 (5.2) A spilling breaker assumption is made so that the wave height H is a constant ratio x to the water depth D = ho + f, H = i(ho + ) (5.3) Substituting equations (5.2) and (5.3) into Eq. (5.1) gives a 3 aK2(ho + 0)2 3 ,2 aho + ) . ax 16(ho + ) ax a8 8x 1x) .) Solving for ^ gives the result 9 Ji Oho (5.5) ax I 3,#2 ax that the slope of the water surface inside the breaker line is a constant times the slope of the beach. Previous models are able to produce the results of the analytical expression but to date none, including the present, have matched the experimental results. In figure 5.1 it is seen that the slope of the mean water surface produced by the wave setup is nearly 3.0 NUMERICAL SOLUTION 2.0 0 EXPERIMENTAL DATA OF BOWEN, ET AL. / 1.5 Break 1.0 Point /o 0.5 a __ ____. SWL / 0.5 400 300 200 100 0 DISTANCE FROM STILLWATER LINE IN BEACH, x CMS Figure 5.1: Comparison of the numerical solution of Vemulakonda et al. (1985) for setup with experimental data of Bowen (1968). Source: Vemulakonda et al. (1985) a constant. For the numerical results the setup starts at the breaking point, whereas in the experimental results the water surface elevation is nearly constant for some distance shoreward of the breaking point and then there is a change in the slope of the water surface elevation. The numerical results follow the analytical solution. This is an area that needs further investigation. The reason for the difference between the analytical and experimental results is that the analytical formulation is predicated upon the radiation stress term being a function of the wave height, so that as the wave height is reduced after the breaking point so is the radiation stress term. As pointed out by Svendsen (1984) the radiation stress term remains nearly constant after breaking for a distance roughly equal to five or six times the breaking depth. This distance is termed by Svendsen as the outer or transition region in which there is a rapid change of wave shape. In the inner region the wave shape according to Svendsen remains nearly constant and is characterized by a surface roller or bore. By incorporating the surface roller Svendsen is able to obtain close agreement for wave height and setup to the experimental results of Hansen and Svendsen (1979). However this close 66 agreement is just for the inner region. For the outer region there is no viable solution other than that the radiation stress term and the water surface elevation remain nearly constant. Further there is no analytical separation of the outer and inner regions other than the a poster definition of where there is a sharp change in the slope of the water surface elevation. This definition is not useful in a model that seeks to determine the water surface elevations. Svendsen's work shows promise and should be pursued; it is not however at this point viable for inclusion in a numerical circulation model. For one it needs to be expanded for the case of waves at a angle (this need will be discussed further in the next section). Another limitation is that Svendsen's model is for monotonically decreasing depths, which is a severe restriction (for example excluding barred profiles). Svendsen conceptually addresses the question concerning the spacial difference between the initiation of the setup in the analytical and experimental results. The analytical model assumes that once breaking commences energy dissipation also begins. Svendsen states that the organized potential energy of the wave prior to breaking is transformed into forward momentum in the outer region. In keeping with the definition of the radiation stresses as the flux of excess momentum, i.e. momentum that would not be present in the absence of the wave, it is easy to conceptualize that the radiation stress term should remain nearly constant just after breaking. The problem is in assuming that the wave form in the breaking zone is characterized the same as outside the breaking zone. The author, at present, offers no solution other than to recognize the existence of the problem. In the present study results for the setup on a plane beach with normally incident waves differ somewhat slightly from the analytical solution as given by LonguetHiggins and Stewart. This is because the present model does not use the spilling breaker assumption. Instead wave energy is dissipated according to a model suggested by Dally (1980) which is explained in section 3.3. Thus the slope of the mean water surface is not a constant times the bottom slope. A result from the experiments of Hansen and Svendsen (1979) show this clearly in figure 5.2. The figure shows a slight curvature in the wave height after breaking mm 140 120 100 80Model 40 Experiment 20 0 5 10 15 205 20 Figure 5.2: Comparison of numerical results from the model with experimental results of Hansen and Svendsen (1979). Wave setup and wave heights for a 2 second wave on a beach slope of 1:34 and incident wave height of 7 centimeters. and in the slope of the setup. Although the problem of where the setup begins is still present, the curvature of the mean water surface is the same for both cases. It is also quite clear that the diminution of the wave height is not constant as is assumed by the spilling breaker assumption. In figure 5.2 the breaking wave height in the experiment is much greater than in the model. This is because linear wave theory cannot adequately predict wave shoaling. Ven gayil and Kirby (1986) were able to duplicate the shoaling recorded in the experiments of Hansen and Svendsen (1979) by numerically solving coupled equations for the many har monic components. Their results showed that in the shoaling process prior to wave breaking, wave energy was transmitted to the higher frequency components. Their work however is not adaptable for use in the current model since it is not known how to dissipate energy 0.0 V_ I I : t I I 0.0 0.5 1.0 1.5 X 2.0 Figure 5.3: Nondimensional longshore velocities vresus the nondimensional distance off shore from the analytical solution of LonguetHiggins. Reprinted from LonguetHiggins (1970b), Journal of Geophysical Research, vol. 75(33) page 6793, copyright by the Ameri can Geophysical Union. from the various frequency components once breaking is initiated. 5.2 Longshore Currents LonguetHiggins (1970b) gives an analytical solution for the longshore current generated by obliquely incident waves on a plane beach. The solution involves a spilling breaker assumption and assumptions concerning the form of the bottom friction and the lateral mixing. For steady state conditions the y equation of motion reduces to a balance between the gradient of the Sxy term and the bottom friction with the lateral mixing serving to distribute the current. Details of the solution can be found in the cited paper. Figure 5.3 gives the results of the analytical solution of LonguetHiggins. In the figure the ordinate V is the nondimensional velocity v/Vo, where v is the actual velocity and Vo is the velocity at the breaking point obtained for the case of no lateral mixing. = (gh,)I n (5.6) The abscissa X is the nondimensional distance z/zb, where zb is the distance from the V(m.sec) V I I i' ":" x(m) 1 2 3 4 5 6 Wave SetUp Line Breaker Line Figure 5.4: Comparison of experimental results of Visser (1982) and numerical result for longshore currents due to a 1.1 second wave on a slope of 1:20 for the deep water conditions H, = .065m and 8o = 20 degrees. shore to the breaker line. Previous models have reproduced the results of the Longuet Higgins solution by incor porating all of the same assumptions (i.e. spilling breakers and bottom friction). Since the present model does not contain the same assumptions comparison will not be made with the analytical results of LonguetHiggins. Instead comparison is made with the experimental results of Visser (1982). For both the analytical and the numerical results it should be noted that the same inadequacy is present as was discussed in the previous section. In both cases it is implied that wave energy dissipation commences as soon as the wave begins breaking. Thus both analytical and numerical solutions initiate a longshore current starting at the breaker line. It is reasonable to expect that the same holds for the y component of the onshore flux of excess momentum as for the x component. That is, that the S,, term also remains nearly constant over a certain unspecified distance as the organized momentum flux of the non breaking wave becomes turbulent momentum upon breaking. This is illustrated in figure 5.4 which shows a comparison of the experimental results of Visser (1982) and the model results for the same test conditions. The maximum current is nearly the same in both the model and the experiment but the maximum is further offshore in the model. This can be attributed to the initiation of a longshore current in the model at the breaker line. 70 5.3 Waves On a RipChannel Both Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) examined the case of waves incident upon a shore with the presence of a ripchannel. The bottom .topography was given according to a formula first used by Noda et al. (1974) h(z, y) = mx {1 + Ae~i)' sin10 (y ztan) (5.7) where 8 is the angle the ripchannel makes with the x axis which is, in their coordinate system, directed offshore from the still water line, m is the beach slope, A is an amplitude of the bottom variation and A is the length of the repeating beach. The origin of the coordinates is at the intersection of the still water line and the trough as is shown in figure 5.5 which gives the depth contours for the case where the slope is .025 and 8 is equal to 30 degrees. Due to the restriction of the present model that there be no y gradients at the lateral boundaries the above formulation is modified to be h(x,y = mX(l + AeS()'sino (yxtanl)} I(y ztanl)l<  (5.8) mx 1(y xtan#)l > r Also shown in figure 5.5 are the contours obtained through use of the modified Eq. 5.8 for p equal to 30 degrees and a .025 slope. The two formulas differ only in that the repeating rip channels in the Noda formulation are removed. For application in the present model the presence of the trough in the farthest offshore grid row is removed so that the depth in the offshore grid row remains constant in compliance with the requirements of the boundary conditions as discussed in sections 4.3 and 4.4. Figure 5.6 shows the results obtained by Ebersole and Dalrymple for the case of a wave angle of 30 degrees a wave height of 1.0 meters and a period of 4 seconds. Figure 5.7 gives the results for the present model. Even though the scales of the two figures are different it is evident that there is agreement between the two. The length scale of the eddy is about 30 meters for both cases. In the caption to the figure of his results, Ebersole indicated that the lateral mixing coefficient may have been large. The present author does not believe Figure 5.5: Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as shown. Bottom: Depth contours as used in the present model. Grid spacing is 5 meters. 72 . .a S 0 0 0 0 0 5 S 00 40 d 0 S 0 S *. 0 S 0 % 0 0 S 0 m s O % S *a 0 0 0 60 *0 0 0 g s 0 0 0 0 0 0 . 0 0 . of 30 de res *Scl fr l oc i *v is I ic = 2 S e E l D y m 0 5 94 .( 9 7 ) Sm 24m 24 of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec. Source: Ebersole and Dalrymple (1979). ..  . . .. .. . .. I % Y ~ ~ e~ c4~.lcc~   LC ~4e~ e .. . . . . . .. .. . .. .. .. . . . . .. .. .. .. .. .. .. .. .. .      .            I25 m I Figure 5.7: Velocity vectors obtained from the present model using the same input condi tions as for Ebersole and Dalrymple. 0.635m/s Maximum Vector  ..... ...... .. ....... .....o o.o.. o.... ................ S....... ....... ...... .................................................. '. .... .... ................................. .......... .. .. I100mI 1 0 0 m Figure 5.8: Mean water surface elevation contours showing that the water surface is lower above the trough. Values indicated are in millimeters. this to be the case. In the present study the lateral mixing coefficient was reduced to as low a value as 2 meter2/second without substantially changing the magnitude of the eddy adjacent to the rip trough. However Ebersole did find a significant difference between the case of no lateral mixing and with lateral mixing. Results for the present model for the case of no lateral mixing are not available since the model sometimes exhibits an instability without the stabilizing influence of the lateral mixing terms. Figure 5.8 shows the mean water surface elevation contours obtained from the model which show that the water mean surface is lower above the trough. The agreement between the results of the Ebersole model and the present model for the case of waves incident upon a ripchannel should not be surprising since diffraction in this case is small compared to current and depth refraction. Another test result is the case of normally incident waves breaking on a bar in which 0 .0 0.0 0.0 .... .......... ........ ........................................................................... .. 1 1 1 1 !i 1 t I I I I I I I I I 1 Figure 5.9: Depth contours in meters for the case of a bar with a gap. 0.622m/s . .. Maximum Vector aaaaaaaaaa a   L l... I>.. ... . ....... .... S I I,. Figure 5.10: Velocity vectors for the case of waves over a bar with a gap. Dashed lines indicate depth contours. Grid spacing equal to 10 meters. indicate depth contours. Grid spacing equal to 10 meters. 76 there is a gap in the bar. Figure 5.9 gives the contour lines for this test case. The topography was created by adding an inverted parabola on to a plane beach with a slope of 1:40. The gap was created by suppressing the bar with a symmetric smooth tangent function. Figure 5.10 shows the velocity vectors resulting from a 6 second wave with an incident wave height of 1 meter. The contours of the bottom topography are represented by dashed lines. The driving force of the circulation is the return flow through the gap in the bar of the water that is carried over the bar by the breaking process. These results may not be totally correct given the arguments of section 5.1. If there were no gap in the bar the waves breaking on the bar would push water shoreward until sufficient setup was produced. However with a gap in the bar lateral flow is allowed and setup cannot be maintained. Kirby (conversation with the author) has remarked that he has observed the initiation of a very strong current in the direction of the waves produced by waves breaking on a sand mound in a small wave basin. This was tested in the model using normally incident waves on a slope of .02 with a shoal described by the following equation which is the shoal used in the Delft experiments of Berkhoff et al. (1982). .45 z < 5.82m h(z, y)= .45 .02(5.82 + z) (gs)2 + (,)2 > 1.0 (5.9) .45 .02(5.82+) .51 + x.3 (.)2+(Y)71..0 Figure 5.11 shows the depth contours for this test case. Figure 5.12 shows the velocity vectors obtained for the case of a barely breaking wave. For the case of a non breaking wave the velocity vectors plotted to the same scale are not discernable. Figure 5.13 plots the max velocity obtained in the shoal region versus the incident wave height. At a wave height of about 10 centimeters there is a sharp change in the slope of the curve indicating increased velocities once breaking commences. 5.4 ShorePerpendicular Breakwater The model is tested here for the case of waves interacting with a shoreperpendicular breakwater. The breakwater is represented as an infinitesimally thin impermeable barrier located in between the J = JBRK and J = JBRK+I grid columns and extending from Figure 5.11: Depth contours for the case of normally incident waves upon a shoal. Grid spacing is .05 meters. Contour interval is .02 meters. 0.162m/s Maximum Vector Figure 5.12: Velocity vectors for the case of waves breaking upon a shoal. L.) >6 CC CL) C, o .00oo 5.00 10.00 15.00 WRVE HEIGHT (CM) Figure 5.13: Maximum velocity across the top of a shoal versus incident wave height for waves over a shoal on a plane beach. Incident wave height less than 9.1 centimeters are nonbreaking waves. Incident wave heights greater than 9.1 centimeters produced waves breaking over the shoal. J = JBRK Y EHL EE E11EE11TI ~i~1 Figure 5.14: Position of the groin in relation to the grid rows and columns. the shore boundary to the I = IBRK grid row. Figure 5.4 illustrates the position of the breakwater. Mathematically, this thin impermeable barrier behaves as an internal solid boundary that prevents any activity on one side of the barrier from being transmitted to the other side. This requires modifications to both wave and current models. In the wave model an internal boundary condition is imposed so that at the breakwa ter face wave energy is totally reflected. This is accomplished by imposing the reflection condition 6A =0 (5.10) on both sides of the breakwater from grid row I = IBRK to the shoreline, or A(I,JBRK) = A(I,JBRK+I) for all I > IBRK When the wave equation is expressed in finite difference for J = JBRK, A(I, JBRK+I), which is on the other side of the barrier, is eliminated by virtue Eq. 5.11. Likewise when the wave equation is applied to the other side of the barrier A(I, JERK) is also eliminated. Further (5.11) E H 1\ 1 81 reflective conditions are imposed at the barrier by assuming that the derivatives of any property perpendicular to the breakwater has a zero value at the breakwater. In the circulation model the internal boundary condition of impermeability at the bar rier is imposed by setting the velocity component normal to the barrier equal to zero. V(I,JBRK+I) = 0 for all I > IBRK (5.12) To insure that currentinduced momentum is not diffused across the barrier, subroutines were introduced in the computation for the x sweep for the two grid columns adjacent to the barrier. Two sample cases were run for comparison purposes. The first case is a breakwater extending 400 meters from the shoreline on a beach with a composite slope of 1 on 10 for the first one hundred meters and a slope of 1 on 100 beyond that point. This is a geometry selected by Liu and Mei (1975). A deep water wave height of 1 meter and a deep water wave angle of 45 degrees was used by Liu and Mei. Using linear shoaling and refraction this gives a wave height of .84 meters and a wave angle of 31.15 degrees as the wave input conditions at the offshore grid row 800 meters from shore in a depth of 16.95 meters. The solution method of Liu and Mei was to solve for the wave field in the absence of currents and mean water level changes and then use the wave field as a constant forcing in the solution of the circulation. By ignoring the lateral mixing and advective acceleration terms and assuming steady state conditions, and by expressing the two horizontal components of the depth integrated flow as derivatives of a stream function V, the two horizontal equations of motion were reduced to one Poisson equation v2, f(z, ) (5.13) where the f term contains the forcing due to water level gradients, friction and wave forcing. This equation represents a boundary value problem with two unknowns 0 and f. Upon specifying boundary values an initial guess for the water surface !I is made and then Eq. 5.13 is solved for V'. With this new solution for 0 the value of o7 is updated by solving steady 82 state equations of motion. This process is repeated until convergence. Figures 5.15 and 5.16 show the results for Liu and Mei and the present study respec tively. Figure 5.15 shows the contours of the stream function whereas figure 5.16 shows flow lines which give the direction of the flow but not the intensity of the flow. The present model shows a separation eddy at the tip of the breakwater which is absent in Liu's results. Both models show a series of rip cells in the up wave side of the jetty which result from the interaction of the incident and reflected waves. There are not as many rip cells in the present model as the weaker ones farther upstream are buried by the stronger longshore current. Neither the Liu model nor the present one show a rip current at the downwave side of the jetty. In figure 5.16 the flow lines in the offshore region are directed in the onshore offshore direction which is different from the stream function contour lines in figure 5.15. This is because of a difference in the offshore boundary condition in the two models. In the Liu model the offshore condition is a constant 0 which means that the flow must be tangential to the offshore boundary. In the example used in the present model the offshore condition is a constant water surface (q = 0) which forces the tangential velocity to be zero at the offshore boundary. An attempt to change the offshore condition to a zero normal flow condition did not accurately conserve water within the computational domain, even though convergence was achieved. Figure 5.17 shows the flow lines resulting from using no onshore flow at the offshore boundary. The second case is to simulate the experiments conducted at the Waterways Experiment Station and reported by Hales (1980). A jetty perpendicular to the shoreline extended 15 feet out from the still water line on a plane beach with a slope of 1 on 20. The plane beach extended beyond the jetty tip until a depth of one foot was reached and further offshore the laboratory basin was a constant depth of one foot. Waves were sent at an angle by use of a moveable wave paddle and wave guides which were set up to follow the wave rays at the ends of the wave paddle. Figure 5.18 shows the experimental setup for these tests. These wave guides in essence made the experiments to be in a closed basin and greatly *401 p 4p ^ 4A 0 Cp 0 (p 00 0 0 0 0 I I I I i i A W W ft N  O V O wP O II OW ; 0 0 00 0080*0 0010 0 Figure 5.15: Results for the numerical model of Liu and Mei (1975). Deepwater wave height equal to 1 meter, deepwater wave angle equal to 45 degrees. Source: Liu and Mei (1975). 0 a 00 01111111fIWIIU Figure 5.16: Results from the present numerical model using the same conditions as were used by Liu and Mei. Offshore boundary condition of constant water surface level which precludes tangential flow at the offshore boundary. 450 ,jjei' I aol o" I E equal to 20 meters. Figure 5.17: Results from the present numerical model using the same conditions as were used by ]Liu and Mci. Offshore boundary condition of no onoffshore flow. Grid spacing equal to 20 meters. Figure 5.18: Experimental setup for the tests conducted by Hales (1980). Source: Hales (1980). I el f __ 2 N > a a W p A A 4 A a 4 A 4 4 4 4 M A 1 A A A q VERGEE MIDDEPTH VELOCITIES COINED REFRACTION OND DIFFRQCTION 1.0 FT/SEC ONGLE OF INCIDENCE = 30 DEGREES WOVE PERIOD a 1.0 SECONDS WOVE: HEIGHT NEAR GENERATOR 0.139 FEFT Figure 5.19: Experimental measurements of the middepth averaged velocities obtain by Hales. Measurements at one foot intervals. Source: Hales (1980). affected some of the circulation patterns away from the jetty. Figure 5.19 shows the average middepth velocities recorded on the down wave side of the jetty. Note the strong offshore current next to the boundary. Since the boundaries are less than two jetty lengths from the jetty, one must ask whether the circulation near the jetty was also affected. The numerical simulation of the Hales experiments was performed using open lateral boundaries, since to fully represent the experimental conditions in which the lateral boundaries are neither straight (i.e. curving wave rays) nor perpendicular to the shoreline boundary would require mapping techniques that are not yet included in the capability of the model. Figure 5.20 shows a comparison of the experimental and numerical results obtained for the magnitude of the rip current adjacent to the down wave side of the jetty for the case of a 1 second wave with an angle of incidence of 30 degrees and a one foot depth wave height & Middepth 0.3o Bottom s ... Model Results ) 0.2 0 .J wuj Shoreline w 0.1 / X > 0 2 4 68 10 12 4 DISTANCE FROM SHORELINE (ft) End of Breakwater (15 ft from shoreline) Figure 5.20: Ripcurrent velocities adjacent to the down wave side of the jetty. Comparison of experimental and numerical results. of .139 ft. One would expect the depth averaged velocity to be somewhere between the middepth velocity and the bottom velocity. The model appears to slightly underpredict the rip current. In figure 5.19 a circulation cell in the lee of the groin is barely discernable. In figure 5.21 the circulation cell is very obvious. 5.5 Shore Parallel Emerged Breakwater In order to solve for the wave and current field resulting from the interaction of waves and an emerged breakwater several boundary conditions at the breakwater are incorporated into the model in addition to the boundary conditions already discussed in sections 4.3 and 4.4. The breakwater is represented as being infinitesimally thin and impermeable. Figure 5.22 shows the position of the breakwater on the shoreward side of the I = IBRK grid row. The impermeability condition implies that neither waves nor currents pass through the breakwater. The no current condition is easily imposed in the circulation model by making U(IBRK+I,J) equal to zero for the entire length of the breakwater. In addition the circulation model is modified so that momentum diffusion and derivatives do not take place _______ __..... .... a...s.aa. aI I IS1iS . . .. ... I l I l l l i .... . . .. I . . .. . I . . ..... 11 1 i 1 1 . . . .. f i . I .. . .l i a.. . . .i l . ... ..... .............. ............... ... ..... ......... .I lll 1 ................. ............. tttttl'II I1ll 1 ..................................... .i I .............................. 1.. .1 ...... ........... ........ .............. .1 . . .. . . . . . ... . 0.4 J*I I I l I I I .... . l S. . . . . .. . * . . . . o ... .. ........ ... 8 8 8 8. . .* II . . . ......................... . .........%tl..l.... .. . . .......... ..... . . . .. . .......... ....... ...... .. .. ... .. ........ l l 1 ....~*&Jil  ' ',..... r r( '.'' .' **^'fftfffffff'~'''' ^ ^ f f f f f f f ^ ^ ^ f ^ f f f * ' ' . .. .***** <''*'^.._________>11 > < I f I < t f * Figure 5.21: Velocity vectors of the currents obtained numerically for the same conditions as for the Hales experiment. Grid spacing equal to one foot. x II II Figure 5.22: Definition sketch of grids for the case of a shore parallel breakwater. across the barrier. In the wave model an adaptation is made so as to prescribe the condition that the wave height is zero on the lee side of the breakwater. This condition on the lee side of the breakwater is physically not realistic since common observation reveals that waves actually propagate along the lee side of a breakwater. The parabolic model however marches waves in but one direction. This assumption however will give good results several grids away from the breakwater. Figure 5.23 shows the crests of the diffracted waves in the lee of a semiinfinite breakwater for a region of constant depth. This figure shows that the zero wave height condition on the lee side of the breakwater does not adequately represent the waves along the breakwater or allow wave energy to properly propagate in the lee of the breakwater; however, important information can be obtained through use of the marching parabolic wave equation method. To set the wave height to zero on the lee side of the breakwater the ACOEF subroutine is modified so as to advance the wave one half delta x from the center of the I = BRK grid row to the position of the breakwater at the edge of i\ 