~ ) Pergamon Computers & Geosciences Vol. 20, No. 6, pp. 905-917, 1994 Copyright © 1994 Elsevier Science Ltd O098-3004(94)EOOI3-J Printed in Great Britain. All rights reserved 0098-3004/94 $7.00 + 0.00 COMPUTATION OF LONGSHORE CURRENTS AND SEDIMENT TRANSPORT FERNANDO J. CAVIGLIA* Comisi6n de Investigaciones Cientificas de la Provincia de Buenos Aires y Facultad de Ciencias Naturales y Museo, UNLP, Casilla de Correo 45 (1900) La Plata, Argentina (Received 26 July 1993; accepted 28 December 1993) Abstract--A numerical model for the steady-state profile of the longshore current induced by regular, obliquely incident, breaking waves is presented. The wave parameters must be given at an arbitrary depth. A rapid convergent numerical algorithm is described for the solution of the governing equation. The model is solved using a nonlinear bottom friction law in which the friction coefficients are a function of the bottom roughness which is computed at each point using an empirical formula. The predicted current profiles are combined with some of the known formulae for sediment transport computations. Key Words: Longshore currents, Sediment transport, Numerical model, Forward sweep implicit technique. INTRODUCTION If breaking waves approach a straight coastline at an oblique angle, they induce a longshore current in the surf zone. The current acts somewhat analogous to a river, transporting sediment mobilized by the break- ing waves, and is confined to a zone with a width of two to three times the width of the surf zone. The prediction of such longshore currents and the associ- ated sediment transport is of prime importance for coastal engineers. The purpose of this paper is to present a numerical model of longshore current in combination with some of the known formulae for predicting sediment trans- port which can be used for engineering purposes. The model is solved on a desktop computer and requires only data available in any engineering project. LONGSHORE CURRENT MODEL According to the generally employed assumptions in longshore current modeling, namely (i) longshore uniformity in waves and bottom topography, (ii) negligible bottom friction in the cross-shore direction, and (iii) applicability of linear-wave theory, the verti- cally integrated, time-averaged momentum equations for nearshore water motion (Mei, 1983) become pgd ~ dS,.x = d x ' ( 1 ) d F.4h d d V ] _ zb__Ly = 1 dS.,. ~x L -~xJ p p d x ' (2) *Address for correspondence: Laboratorio de Oceangrafia Costera y Estuarios, Casilla de Correo 45 (1900) La Plata, Argentina. for the cross-shore x component and for the long- shore y component respectively (see Fig. 1), where p is the water density, g the acceleration resulting from gravity, d = h + r/ the mean water depth, h the still water depth, t/the elevation of the free surface above the still water level, Ah the eddy viscosity coefficient, V the longshore current velocity and %y the bottom frictional stress. The radiation stress components Sxx and Sxy are given by Sxx = E[n (cos 2 c~ + 1) - ½], (3) sin Sxy = T E~, (4) where Fx = - E C n cos c~ is the wave energy flux towards the shoreline, C the phase velocity of waves, (kd/Sh2kd), ct the angle that the wave front n = ~ + makes with the shoreline, E =~pgH 2 the energy density of the waves per unit area, H the wave height, k = 21t/L the wave number and L the wavelength. The phase velocity is calculated with the approxi- mated expression introduced by Visser (1984) 1 - as --~< 0.36 C = L° d , (5) / ° - = Co as 0.36 L 2 n L0 > where T is the wave period and the subscript 0 denotes a value on deep water. Wave set-up Equation (1) governs the mean water surface displacement r/. According to the specific assumptions of this paper namely, (iv) constant beach slope, (v) 905 Y F. J. CAVIGLIA proportionali ty between the wave height and the mean water depth in the surf zone, H = yd, (vi) negligible wave set-down offshore of the plunge line, xo, and (vii) applicability of the shallow water approximation, n = 1, cos g ~ 1 and C = x / ~ , the equation for the mean water surface displacement can be integrated to obtain = K s ( x p - x ) , (6) where K = [1 + 8/372] l, y the breaker index and s the beach slope. If Equation (6) is evaluated at x = x~, the following relation can be obtained x , = K x r, (7) where xs is the line of maximum wave set-up. The plunge line is calculated with the following expression (Svendsen, 1987; Larson and Kraus, 1989) xp = xb - 3Hb, (8) where Hb is the breaking wave height, Xb = (db/s + xs) wave crest >X shore s t i l l -wa te r line line PllUnneg e b r ea ker line Z 906 XS r/ xp xb h beach Figure 1. Definition sketch for coordinate system and nearshore region, where rt is wave set-up, h still-water depth, xb width of breaker zone, x, distance from shoreline to still-water line, xp distance from shoreline to plunge line, and a angle of incidence of waves. Longshore currents and sediment transport 907 the breaker zone width and db the depth at breaking (db= Hb/7 ). Substitution of (8) into (7), taking into account the expression for Xb, gives 3~'2 [db -- 3nb]. (9) x ~ = 8 L s From (6) and (7), the mean water depth can be written as = { (i-@)sg~bbxb as g ~ < P , (~ - ~ ) s ~ as g > P (10) where the following dimensionless variables were introduced Xp X s 3 = d 5 x P - - - , ~ (11) X b Xb Xb Driving force The longshore driving force resulting from obliquely incident wave is 0Sxy = sin cti 0F~ (12) dx C O x ' where the subscript i denotes a value at an arbitrary depth. Substituting the expressions of the wave energy flux towards the shoreline and the energy density of the waves per unit area, in combination with the assump- tion (v) and the shallow water approximation, into (12) gives I I-'q ~S~y sin 0t i 5 pg3/272 X 3/2 as x ~< Xb = ] 0 LXb -1 ~x C~ as x > Xb L (13) Visser (1984) introduced the following assump- tions: (viii) the dissipation of wave energy takes place shoreward of the plunge line, (ix) dFx/dx is pro- portional to x 3a shoreward of the plunge line and vanishes offshore of this line, and (x) the transport of wave energy towards the shore is given by the transport of wave energy predicted by linear wave theory in the cross-shore x-direction. Under these assumptions follows from Equation (13) that • f I-£ -13/: dS:y s,ncti~--Dbk-fi] as g < ~ P (14) ~-x = c"(O as ~ > P ' where 5 1 D b = ~ : - - E~ Cin, cos oc i. rxb Bottom friction term Under the combined action of waves and currents the time-averaged bottom friction in longshore direc- tion can be expressed as TbY-- 1FI;oT p 8 [ V2 + 2~Ur. V sin ~ cos cot -q- (~Um) 2 COS 2 (.0/] 1/2 * (V + GUm sin ~ cos cot) dt, (15) where co = 2rr/T is the angular frequency, u m = n / T H / S h k d the maximum orbital velocity near the bottom and ¢ = 2(Fw/F) l/: with Fw being the Jonsson friction coefficient and F the Darcy-Weisbach friction coefficient. Because the elliptic integral (15) must be solved frequently in the numerical longshore current program in this paper the square wave approximation introduced by Nishimura (1988) is used• Then, the time-averaged bottom friction in longshore direction becomes ~ ( z2sin2~z' / Zby= F Z + ~ I V , (16) P Z = ½[(V 2 + z 2 + 2zV sin ~)2,,2 + (V 2 + z: - 2zV sin ct)l/2]. and z = (2/n)¢Um. The most popular formulae for predicting the friction coefficients are I [ / r \ore-1 r exp -5.977 + 5.213/--/ I as -- ~< 0.63 F w = \ ab] J ab r | 0 . 3 as - ->0 .63 L ab (17) where ab=um/co is the horizontal water particle excursion at the bottom and r the equivalent diameter of the roughness elements. Equation (17) was suggested by Jonsson (1980) and Equation (18) is based on Nikuradse's experiments. According to Nielsen (1985) the total roughness, r, is calculated from the following expression r = 190Dx/O' - 0.05 + 8 =-, (19) Ar where q~ is the ripple height, 2r the ripple length, D the mean sediment grain diameter, O' = s F w ~ " the so-called skin friction Shields parameter, F;, the Jonsson friction coefficient based on r = 2.5D, ~P = u~/[(s - l)gD] the mobility parameter and s the relative sediment density (s = Ps/P, where Ps is the sediment density). If 0 ' ~< 0.05 then r = 2.5D, which is the roughness for a flatbed of fixed sand 908 F.J. CAV1GLIA grains. Nielsen (1981) recommended the following formulae for prediction of tipple geometry r/_~ = 0.342 - 0.34,ff-07, (20) 2r q~= ~0 .275- 0.022x/~ as ~ ~< 10 (21) ab ~, 21~P-LSs as ~ >I 10" This author recommended to base the tipple calculations on the significant wave height for field situations. The proposed expressions for ripple geometry calculation are valid for nonbreaking wave conditions. Ripples are washed out when the mobility parameter, ~, is larger than about 200-250 (see Dingier and Inman, 1976). Van Rijn (1989) assumed that the tipple existence is limited to values of 7 t smaller than 250. Then, from Equation (19) r = 190Dx/~ 7 - 0.05, (22) which is the formula recommended by Nielsen (1985) for the prediction of the roughness of flat, mobile sand beds. Under breaking wave conditions the mobility parameter, in general, will be larger than 250. Further- more, it is assumed that the longshore current has no influence upon the tipple geometry. Lateral mixing The first term of Equation (2) is an empirical closure relation describing lateral mixing resulting from Reynolds stresses. Visser (1984) expressed the eddy viscosity coefficient as A h = Mqd, (23) where M is a constant of order 3 and q the character- istic turbulent velocity. The nondimensionai form of A h is Ah "e4h ~ " Z ' (24) where A0 can be defined in arbitrary form. The eddy viscosity coefficient resulting from this parametetization decreases rapidly offshore of the plunge line. This result seems to be in agreement with the laboratory measurements of Deguchi, Sawaragi, and Ono (1992), who determined that Ah becomes maximum near d/db = 0.7-0.8 and decreases rapidly toward offshore. NUMERICAL SOLUTION Substitution of (14), (16), and (24) into (2) yields d27 d 7 l f o g 3 ( £ ) as . ~ < P d~ ~ t - g ~ ( ~ ) ~ - g : ( ~ ) 7 = in which 7 = V~ Vo, 1 d-4h~ g ' ( if) = -4h2 d.f ' as . f > P ' (25) (26) with F f Z z2 sin 2 ct'~ g2( i ) = Cl 1 F . f l 3/2 (27) (28) 1 x~, Dbx~, sin ct i C l = - - - and V 0 = - - - - 8 Aod b pAodb Ci The differential equation to be solved is nonlinear because g2($) contains V through Z. The differential equation is solved as a linear equation. This is done by solving (25) and (27) with an iteration procedure. The canonical form of the differential equation for (25) for calculation of the grid cell number is g4, Vi - 1 + gs, Vi + g6, Vi + 1 = gT~, (29) g4, = 1 l - ( 2 + g2i A.~ 2) - - ~gl , A.~, gsi = g6, = 1 + lgli A.f and og3, AX 2 as £ ~< P gT,= as ~ > P" Equation (29) is solved for each calculation cell, with i extending from 2 to N - 1, for a grid encom- passing N cells. Boundary conditions for 7 must be provided at cells 1 and N. The boundary conditions are 7 = 0 as £ = 0 or i = l 7 = 0 as £ ~ or i ~ with ~ = (i - 1)A.f. Because for £ > 2 the longshore current velocity 7 becomes small (Visser, 1984) it is assumed that 7 = 0 as £ = 3 or i = 2 4 1 if A£ = 1/80. For i = 2-240, Equation (29) is a tridiagonal set of linear equations in the unknowns V = [ F z , 73 . . . . . 7N_ l] and can be written in matrix form as MV = G where G = [gT:, g73 . . . . . g7u ,]. The matrix M is tridiagonal with principal diagonal elements gs,, the elements to left of the diagonal are g4i and to the tight are gr,. All the other elements are zero. This system is solved easily by forward sweep and a back- ward substitution. The algorithm is simply Gaussian elimination. The backward substitution solves for V, and it is given by Vi = Ui - I4I,. 7i+! for i = 2-240, (30) where with g6i w~= gs, - Wi_ lg4, g7i - - U i - lg4, W gsi - i- lg4~ 724 ~ = W 1 = U~ = O. and U, Longshore currents and sediment transport 909 Table 1. Limiting wave height (from Van Rijn, 1989) Breaker index Wave steepness s = 1/100 s = 1/20 s = 1/10 H~o/L o = 0.002 0.60 0.70 0.80 0.004 0.60 0.70 0.80 0.006 0.60 0.70 0.80 0.010 0.60 0.70 0.80 0.020 0.50 0.65 0.80 0.040 0.40 0.55 0.70 0.060 0.35 0.50 0.70 The forward sweep reduces ~ to an upper-triangu- lar matrix with principal diagonal elements unity and elements to the right of the diagonal given by W i. The current velocity calculated with (30), previously multiplied by V 0, is replaced into (27) and the procedure described here is repeated again. This pro- cedure is carried out until the absolute value of the difference between V~ and the velocity calculated in the preceding iteration step be less than E = 0.01. In gen- eral, the solution converges after no more than eight iterations. FIELD SITUATIONS A question which arises in comparing monochro- matic wave models with field random waves is the determination of an appropriate wave height parameter for model calculations, Wu, Thornton, and Guza (1985) gave a detailed explanation for selecting the root-mean square wave height, H~ ...... as the ap- propriate one. Another related problem is the appropriate breaker index for random waves. Thornton, Wu, and Guza (1984) presented a criteria for random breaking waves based on field and laboratory data. From their results, Van Rijn (1989) derived the values of Vs(?s = Hbs/db , where Hbs is the significant breaking wave height) shown in Table 1. Multiple linear regression gives ?~ = 0.603 + 2.147s + (30.685s - 5.192)[~0~ ] . (31) where the correlation coefficient is 0.99. For a Rayleigh wave height probability distribution, the relationship between significant wave hei~_ht and root-mean square wave height is H s ~ x / 2 H ...... • Thornton, Wu, and Guza (1984) showed that this relation could be used with an error of about + 5%. From this relation, it is assumed that 7 ...... = ?s/X/~, (32) where 7 ...... is the breaker index based on Hr.~.~. When the input conditions are those at the breaker line a simple predictive formula is used to estimate the deep water wave height prior to determining the breaker index (Larson and Kraus, 1989) I-H0 -]-0.2, Hbs = 0.53H0~L-~0 j . (33) USING THE MODEL The program has been written in Microsoft Quick- BASIC 4.50. The user is requested to input the follow- ing parameters: the significant wave height (Hs in meters) or the root-mean square wave height (Ha ..... in meters), the angle that the wave front makes with the shoreline (~ in degrees), the wave period (Tin seconds), the depth in which wave measurements were made (meters) if the wave conditions are not at the breaker line, the beach slope, the mean sediment diameter (D in millimeters), the particle diameter 090 (10% by weight exceeded in size, millimeters) which is not strictly necessary for the program to be run, and the sediment density (Ps in kilogram per cubic meter). With these parameters the program computes the longshore cur- rent and sediment transport distributions, the mean longshore current in m/s and the longshore sediment transport rate in ma/s and m3/day. The program includes a routine which allows a screen view of the longshore current and sediment transport distributions. This routine requires a VGA or M C G A display adapter. As an example of the program capability a print screen of the resulting longshore current distribution for the Thorn ton and Guza (1986, 1989) 4 February data are shown in Figure 2. Wave conditions at 9. l-m depth were: Hr.~.s.=0.52m, T = 14.2sec and = 18.4 °, with s = 0.038, D = 0.23 mm (Gable, 1981) and Ps = 2650 kg/m 3. The data are superimposed on the graph. DISCUSSION The specialized problem of studying longshore current generation in the surf zone is solved by the use of a forward sweep implicit technique. The program presented here is short, fast, easy to use, and versatile. The numerical model can be adapted easily 1.0 Longshore current d is t r ibut ion V 0 = 1.14 m/s X b = 45.0 m H b = 0.76 m o~ b = 7.8 ° V T =14.2 s V-- 0.5 -- S = 1:26.3 O f m .m • m ~'~m m'-~,,,~% ~ 0 i 2 x/X b Figure 2. Print screen of computed longshore current distri- bution ( ) for Thornton and Guza (1986, 1989) February 4 data (n) . 910 F.J . CAVIGLIA to o ther eddy viscosity coefficients, wave energy dissipat ion models, b o t t o m shear stress theories, etc. The present model for longshore current computa t ions was combined with any of the known formulae for predict ing sediment t r anspor t distri- but ion. Some of them are described briefly in the Appendix. REFERENCES Bijker, E. W., 1967, Some considerations about scales for coastal models with movable bed: Publ. No. 50, Delft Hydraulics Laboratory, Delft, The Netherlands, 147 p. Deguchi, I., Sawaragi, T., and Ono, M., 1992, Longshore current and lateral mixing in the surf zone (abst.): Proc. 23rd Coastal Engrg. Conf., A.S.C.E. (Venice, Italy), Book of Abstracts, Paper No. 172, p. 353-354. Dingier, J. R., and Inman, B. L., 1976, Wave-formated ripples in nearshore sands: Proc. 15th Coastal Engrg. Conf., A.S.C.E. (Honolulu, Hawaii), p. 2109-2126. Gable, C. G., ed., 1981, Report on data from the Nearshore Sediment Transport Study experiment on Leadbetter Beach, Santa Barbara, California, January-February 1980: IMR Ref. 80-5, Univ. California, Inst. Marine Resources, La Jolla, California, 314 p. Jonsson, I. G., 1980, A new approach to oscillatory, rough turbulent boundary layers: Ocean Engrg., v. 7, no. 1, p. 109-152. Larson, M., and Kraus, N. C., 1989, SBEACH: numerical model for simulating storm-induced beach change: Rept. 1, U.S. Army Coastal Eng. Res. Center, Washington, D.C. Tech. Rept. CERC-89-9, 255 p. Mei, C. C., 1983, The applied dynamics of ocean surface waves: John Wiley & Sons, Inc. New York, Ch. 10, p. 451-503. Nielsen, P., 1981, Dynamics and geometry of wave gener- ated ripples: Jour. Geophys. Res., v. 86, no. C7, p. 6467 6472. Nielsen, P., 1985, A short manual of coastal bottom bound- ary layers and sediment transport: Public Works Dept., Tech. Memo. 85/I, N.S.W., 56 p. Nishimura, H., 1988, Computation of nearshore current, in Horikawa, K., ed., Nearshore dynamics and coastal processes: Univ. Tokyo Press, Tokyo, Japan, p. 271-291. Svendsen, I. A., 1987, Analysis of surf zone turbulence: Jour. Geophys. Res., v. 92, no. C5, p. 5115-5124. Swart, D. H., 1976, Predictive equations regarding coastal transport: Proc. Coastal Engrg. Conf., A.S.C.E. (Hon- olulu, Hawaii), p. 1113-1132. Thornton, E. B., and Guza, R. T., 1986, Surf zone longshore currents and random waves: Field data and models: Jour. Phys. Oceanography, v. 16, no. 7, p. 1165-1178. Thornton, E. B., and Guza, R. T., 1989, Models for surf zone dynamics, in Seymour, R. J., ed., Nearshore sedi- ment transport: Plenum Press, New York, p. 337-369. Thornton, E. B., Wu, C. S., and Guza, R. T., 1984, Breaking design criteria: Proc. 19th Coastal Engrg. Conf., A.S.C.E. (Houston, Texas), p. 31 41. U.S. Army Corps of Engineers, 1984, Shore Protection Manual: Wtrwy. Expr. Station, Vicksburg, Mississippi, v. 1, Ch. 4, p. 89-99. Van Rijn, L. C., 1989, Sediment transport by currents and waves: Handbook, Rept. H 461, Delft Hydraulics, Delft, The Netherlands, Ch. 2, p. 15-19 and Ch. 6, p. 15-21. Visser, P. J., 1984, A mathematical model of uniform longshore currents and the comparison with laboratory data: Rept. No. 84-2, Dept. Civ. Engrg., Delft Univ. Tech., Delft, The Netherlands, 151 p. Watanabe, A., 1992, Total rate and cross-shore distribution of longshore sand transport (abst.): Proc. 23rd Coastal Engrg. Conf., A.S.C.E. (Venice, Italy), Book of Ab- stracts, Paper No. 297, p. 619-620. Wu, C. S., Thornton, E. B., and Guza, R. T., 1985, Waves and longshore currents: Comparison of a numerical model with field data: Jour. Geophys. Res., v. 90, no. C3, p. 4951-4958. APPENDIX Sediment Transport Formulae Because details on the applied formulae are given in the original papers, the formulae are described only briefly in this Appendix. Adapted Engelund-Hansen Formula The formula for the total sediment transport qT of the adapted Engelund-Hansen method is 0.05chd qT = Vp2gS/2 A2 ~ , (AI) where Zr=~pFV2[ 1+1{ ~ ~ - ) J ' Um~2] C~=8g and A =s--1 . F Bijker Formula Bijker (1967) suggested to use the Kalinsk~Frijlink for- mula for bedload computation C. L /aT, j where/~ = (Ch/Cgo) 15, C90 = 18 log (12d)/Dgo. Bijker (1967) used the Einstein expression for suspended load. Swart (1976) proposed a modification of the procedure for calcu- lations of suspended load. He used the following expression (z)b c = c, a ' (A3) where c is the sediment concentration at height z, c,, the concentration at a reference height a and 09 a 0.013z with Z = wJO.4V,, r, V,, r =(r r /p) I/2 and w, the sediment fall velocity. The bedload is assumed to take place in a layer with thickness a = r above the bed. From Bijker (1967) qb C , - 6.35 V, r ' (A5) so the suspended load becomes q' ] q s - 6.35(]-- - 1 V ' (A6) where V, = [1/8 F] I/2 V. The total load is qv = qb + q~" Watanabe Formula The sediment transport formula proposed by Watanabe (1992) is q T : A c [ ( z ~ V - ] , (A7, where A¢ is a dimensionless coefficient of order 2, r the maximum value of bottom friction under the combined action of waves and longshore current with the roughness length equal to the grain diameter, Tcr : ( P s - - p ) g D O c , the critical shear stress for the onset of general sand movement and 0or the Shields parameter. Here r, is used instead of r. Longshore currents and sediment transport 911 CERC Formula The CERC formula is an integral method which relates the total sediment transport across the surf zone with the longshore component of the wave energy flux at the breaker line (U.S. Army Corps of Engineers, 1984) e,b Qr = K ( p x _ p ) g a " (A8) where K is a dimensionless constant, a ' = 0.65 the sand porosity and Pzb = EbCb sin ct b cos ct b the longshore com- ponent of the wave energy flux. The subscript b indicates that all evaluations are to be made at the breaker line. A value of K = 0.77 must be used if the computation of Pub is based on the H, ..... wave height. For reducing the wave height to the breaker line is used the usual breaking criterion for spilling breaking waves, V = 0.78. 10 Program Listing A QuickBASIC 4.50 program to compute longshore currents and the associated sediment transport S -Beach slope (non-dimensional); HOI -Wave height to an arbitrary depth (meters); HB -Wave height at the breaker line (meters); TA -Angle that wave front makes with the shore llne (degree); T -Wave period (seconds); DO1 -Depth of the wave measurements (meters); D50 -Mean sediment diameter (millimeters); Dg0 -Particle diameter such that 10% by weight exceeded in size (millimeters); RSED-Sediment density (kilogram per cubic meter). V(I)-Nondimenslonal current velocity CLS LOCATE 12, 15: INPUT •Do you want to see the instructions-(Y or N)"; P$ IF P$ = •N • OR P$ = •n • GOTO 5 CLS : LOCATE 4, 25: PRINT •Input parameters: • LOCATE 8, i0 PRINT •A-Significant wave height or root-mean square wave height (m). • LOCATE 9, I0 PRINT =B-Angle that wave front makes with the shore line (degree)." LOCATE I0, i0: PRINT "C-Wave period (s). • LOCATE ii, I0: PRINT =D-Beach slope. • LOCATE 12, I0: PRINT "E-Mean sediment diameter, DS0 (mm)." LOCATE 13, i0 PRINT •F-Particle diameter such that 10% by weight exceeded" LOCATE 14, i0: PRINT " in size, Dg0 (mm). It is not strictly necessary for the" LOCATE 15, i0: PRINT = program to be run." LOCATE 16, I0: PRINT =G-Sediment density (kg/m3)." LOCATE 17, i0 PRINT "H-Depth of wave measurements (m). It is necessary only when the" LOCATE 18, i0: PRINT " input conditions are not at the breaker line. • LOCATE 25, I0: INPUT "Press ENTER to continue"; P$ DIM H(241), ALF(241), FC(241), FW(241), UM(241), D(241), AH(241) DIM W(241), U(241), V(241), VI(241), QB(241), QW(241), QH(241), R(241) CLS : LOCATE 8, 15: PRINT "Wave height parameter" LOCATE 11, 10: PRINT •Significant wave height ........ (1)" LOCATE 12, I0: PRINT •Root-mean square wave height ...(2) • LOCATE 13, i0: INPUT • •; D$: CLS LOCATE 12, I0: INPUT •Are conditions at the breaker line-(Y or N)•; C$ CLS : LOCATE 12, 15: INPUT "Beach slope'; S CLS : LOCATE 12, iS: INPUT "Wave height=; HOI CLS : LOCATE 12, 15: INPUT =Wave angle=; TA CLS : LOCATE 12, 15: INPUT •Wave period•; T CLS : LOCATE 12, 15: INPUT "Mean sediment diameter•; D50 CLS : LOCATE 12, 12 INPUT "Is the sediment diameter Dg0 available-(Y or N)=; PI$ CLS : IF PI$ = •YW OR PI$ = =y= THEN LOCATE 12, 20: INPUT =Dg0•; Dg0 CLS : LOCATE 12, 15: INPUT "Sediment density•; RSED IF C$ = my. OR C$ = "y" GOTO 15 CLS : LOCATE 12, 15: INPUT •Depth of wave measurements•; DO1 912 F. J. CAVIGLIA 15 IF D$ = "i" THEN H01 = H01 / (2 ^ .5) D50 = D50 / 1000z Dg0 = D90 / 1000 T01 = TA * ATN(1) / 45 IF C$ = "Y" OR C$ = "y" THEN HB = HOI: TB = T01 LO = 9.8 * (T ^ 2) / (8 * ATN(1)): CO = L0 / T AAI = .603 + 2.147 * S BBI = 30.685 * S - 5.192 IF C$ = "N m OR C$ = =n e GOTO 20 HO = ((2 ^ .5) * HB * (LO ^ -.24) / .53) ^ (1 / .76) MU = AAI + BBI * (HO / L0): MU = MU / (2 ^ .5): HO = HO / (2 ^ .5) PB = HB / MU CB = ((9.8 * PB) ^ .5) * (I - PB / LO) GOTO 30 Deep water wave condition computation 20 IF DO1 >= LO / 2 THEN TO = T01: HO = HOlt GOT0 25 CO1 = ((9.8 * D01) ^ .5) * (i - D01 / LO) IF DO1 / LO > .36 THEN CO1 = CO LOI = C01 * T ALF5 = SIN(T01) * CO / CO1 TO = ATN(ALF5 / ((i - ALF5 ^ 2) ^ .5)) KK = 8 * ATN(1) * DO1 / L01 SHI = 2.7182818# ^ (2 * KK) : SH2 = 2.7182818# ^ -(2 * KK) SH = (SHI - SH2) / 2: GX = 1 + 2 * KK / SHt CGOI = .5 * CO1 * GX HO = H01 * ((2 * CGOI * COS(T01) / (CO * COS(T0))) ^ .5) Breaking wave computation 25 MU -- AAI / (2 ^ .5) + BBI * (HO / LO) HB = HO DO HBI = HB PB = HB / MU CB = ((9.8 * PB) ^ .5) * (1 - PB / LO) ALF7 = SIN(T0) * CB / CO T B = A T N ( A L F 7 / ( ( I - A L F 7 ^ 2 ) ^ . 5 ) ) HB = HO * ((CO * COS(T0) / (2 * CB * COS(TB))) ^ .5) LOOP WHILE ABS (HBI - HB) > .001 30 LP = 3 * HB XS = 3 * (MU ^ 2) * (PB / S - 3 * HB) / 8 XB = PB / S + XS XP = XB - LP XS -- XS / XB. P = XP / XB Plunge line search DX= 1/ 80 J-- 0 PP = P DO PP = PP - DXt J = J + 1 LOOP UNTIL PP <= DX P = J* DX IP = J + 1 J = 0: PP = 0 CGI = CGOI: TI = T01: CI = COlt HI = HOI" EI = 9.8 * (HOI ^ 2) / 8 EB = 9.8 * (HB ^ 2) / 8: EO = 9.8 * (HO ^ 2) / 8 IF C$ = .yu OR C$ = uy. THEN EI = EB: TI = TB. CGI = CB IF C$ = .ym OR C$ = my. THEN HI = HB." CI = CB IF HOI = H0 THEN EI = E0" TI -- T0: CGI -- CO / 2: CI = CO" HI = H0 DB = 2.5 * EI * CGI * COS(TI) / (P * XB) M = 3 Bl = 3 / (M * (((i - XS / P) * S) ^ 2)): B2 = 3 / ((S ^ 2) * M) R1 = -.5 + (.25 + BI) ^ .5:R2 = -.5 - (.25 + B2) ^ .5 AI2 = R2 - R1 * (i - XS / P) A1 = ((i - XS / P) * 1.5 - R2) * (P ^ -RI) / AI2 A2 = ((P- XS) ^ (i- R2)) * ((1.5- RI) / P) / AI2 Longshore currents and sediment transport 913 AHB = CB * PB DELTA = (ESED - 1025) / 1025 FOR I = 2 TO 240 X = (I - I) * DX Non-dimensional mean water depth computation IF X <= P THEN D(I) = (1 - XS / P) * S * X * XB / PB: GOTO 35 D(I) = (X - XS) * S * XB / PB Phase velocity computation 35 C = ((9.8 * PB * D(I)) ^ .5) * (i - D(I) * PB / L0) IF D(I) * PB / LO > .36 THEN C = CO Incidence's angle computation ALF = SIN(TI) * C / CI ALF(I) = ATN(ALF / ((i - ALF ^ 2) ^ .5)) Wave height computation L011 = C * T. KK = 8 * ATN(1) * PB * D(I) / L011 IF I <= 81 THEN H(I) = HB * D(I): GOTO 40 SHI = 2.7182818# ^ (2 * KK). SH2 = 2.7182818# ^ -(2 * KK) SH = .5 * (SHI - SH2) : GX = 1 + 2 * KK / SH: CG = .5 * C * GX H(I) = HOI * ((CGI * COS(TI) / (CG * COS(ALF(I)))) ^ .5) Maximun orbital velocity computation at the bottom 40 SHI = 2.7182818# ^ KK: SH2 = 2.7182818# * -KK SH = .5 * (SHI - SH2) US(I) = 4 * ATN(1) * H(I) / (T * SH) Non-dimensional eddy viscosity coefficient computation 45 IF I > IP THEN QCIN -- (A2 * ((X - XS) ^ R2)) ^ (1 / 3): GOTO 45 QCIN = (AI * (X ^ El) ÷ (X / P) ^ 1.5) ^ (I / 3) AH(I) = M * QCIN * D(I) * PB * (DR A (i / 3)) / AHB Ripple and roughness computations 50 UMB = US(X) * (2 ^ .5) AB = UMB * T / (8 * ATH(1)) F1 = (UMB ^ 2) / (DELTA * 9.8 * D50) FWI = 2.7182818# A (5.213 * ((2.5 * D50 / AB) ^ .194) - 5.977) T1 = .5 * F1 * FWI IF F1 >= 250 THEN ET = 0: GOTO 50 IF F1 < I0 THEN ET = AB * (.275 - .022 * (FI ^ .5)) : GOTO 50 ET = 21 * AB * (FI ^ -1.85) EL = .342 - .34 * (TI ^ .25) IF T1 < .05 THEN R(I) = 2.5 * DS0." GOT0 55 R(I) = 8 * EL * ET + 190 * D50 * ((TI - .05) ~ .5) Friction coefficients computations 55 FC(I) = 8 * ((2.5 * LOG(12 * PB * D(I) / R(I))) A -2) ABC = E(I) / AB IF ABC > .63 THEN FW(I) = .3: GOTO 60 FW(I) = 2.7182818# A (-5.977 + 5.213 * (ABC ^ .194)) 60 NEXT I VO = DB * (XB ^ 2) * SIN(TI) / (AHB * PB * CI) C1 = .125 * (XB A 2) / (AHB * PB) v(1) = o W(1) = 0, U(1) = V(1) DO • Computation of the elements of the tridiagonal and • upper-trlangular matrices 914 F.J . CAVIGLIA 65 J = 0 FOR I = 2 TO 240 X = (I - i) * DX EP = 2 * ((FW(I) / FC(I)) ^ .5) FI = EP * UM(I) ZIX = FI / (2 * ATH(1)) Zl = (VO * V(I)) ^ 2 + 2 * ZIX * SIN(ALF(I)) * VO * V(I) + ZIX ^ 2 Z2 = (VO * V(1)) ^ 2 - 2 * ZlX * SIN(ALF(I)) * VO * V(I) + ZIX ^ 2 Z = .5 * (Zl ^ .5 + Z2 ^ .5) FT = Z + ((ZIX * SIN(ALF(I))) ^ 2) / Z G2 = C1 * FC(I) * FT / (AH(I) * D(I)) GII -- AH(I + I) * D(I + i) - AH(I - i) * D(I - i) GI2 = 2 * DX * AH(I) * D(I) G1 = GII / GI2 IF I > IP THEN G3 = 0: GOTO 65 G3-- ((X / P) ^ 1.5) / (AH(I) * D(I)) G4 = 1 - .5 * G1 * DX G5 -- -(2 + G2 * (DX ^ 2)) G6 = 1 + .5 * G1 * DX G7 = -G3 * (DX ^ 2) W(I) = G6 / (G5 - G4 * W(I - i)) U(I) = (G7 - G4 * U(I - i)) / (G5 - G4 * W(I - I)) NEXT I Backward substitution for computing the non-dimensional longshore current velocity V(241) = 0 FOR I = 240 TO 2 STEP -i v(I) = u(x) - w(x) * v(x + 1) IF ABS(V(I) - VI(I)) > .01 THEN J = J + 1 Vl(X) = v(I) NEXT I LOOP WHILE J <> 0 I = 0: MM = 0 70 DO MM = MM + V(I) * VO I = I + 1 IF I <= 81 GOTO 70 LOOP WHILE (V(I) * VO) > .01 FIN -- I - 1 VM = MM / FIN Fall sediment velocity computation VISC = .0000014 B = 9.8 * DELTA * (D50 ^ 3) / (VISC ^ 2) IF B < 39 THEN WS = 9.8 * DELTA * (D50 ^ 2) / (18 * VISC): GOTO 75 IF B > i0000 THEN WS = (9.8 * DELTA * DS0 / .91) ^ .5: GOTO 75 WS = ((9.8 * DELTA) ^ .7) * (D50 ^ 1.1) / (6 * (VISC ^ .4)) 75 BB = 0: HH = 0: WW = 0 80 Computation of the critical shear stress for the onset of sediment movement DAST = D50 * ((9.8 * DELTA / (VISC ^ 2)) ^ (1 / 3)) IF DAST <= 4 THEN TICR = .24 / DASTt GOTO 80 IF DAST <- i0 THEN TICR = .14 * (DAST ^ -.64): GOTO 80 IF DAST <= 20 THEN TICR = .04 * (DAST ^ -.i): GOTO 80 IF DAST <= 150 THEN TICR = .013 * (DAST ^ .29): GOTO 80 TICR = .055 TBCR = DELTA * 9.8 * D50 * TICR Sediment transport computation FOR I = 2 TO FIN EP = 2 * ((~W(I) / FC(1)) ^ .5) TC = .125 * FC(I) * ((V(1) * VO) ^ 2) TR = TC * (I + .5 * ((EP * UM(I} / (V(I) * VO)) ^ 2)) VR = TR ^ .5: VC = TC ^ .5 85 Longshore currents and sediment transport BiJker formula CHE = (8 * 9.8 / FC(I)) ~ .5 IF PI$ = "N" OR PI$ = "n" GOTO 85 IF TR <= TBCR GOTO 85 Bll = WS / (.4 * VR) B1--1.05 * (BII ^ .96) * ((R(I) / (D(I) * PB)) A (.013 * BII)) BI2 = 1 - B1 Z1 = ((D(I) * PB / R(I)) ^ BI2) - 1 FC90 = 8 * ((2.5 * LOG(12 * D(I) * PB / D90)) ^ -2) CHg0 = (8 * 9.8 / Fcg0) ^ .5 NU = (CHE / CH90) ^ 1.5 QBB = 5 * DS0 * VC QB = QBB * (2.7182818# ~ -(.27 * 9.8 * DELTA * DS0 / (NU * TR))) QS = QB * Z1 * V(I) * VO / (6.35 * BI2 * VC) QB(I) = QB + QS IF QB(I) > QB(I - i) THEN QBMAX = QB(I) BB = BB + QB(I) 915 Watanabe formula 90 IF TR <= TBCR GOTO 90 QT = 2 * (TR - TBCR) * V(I) * V0 / 9.8 QW(I) = QT IF QW(I) > QW(I - i) THEN WAMAX = QW(I) ww = ww + Qw(1) Adapted Engelund-Hansen formula QHI = ((9.8 * DELTA) ^ 2) * D50 * (9.8 ^ .5) QH(I) = .05 * CHE * (TR ^ 2) * V(I) * VO / QHI IF QH(I) > QH(I - I) THEN QEMAX -- QH(I) HH = HH + QH(1) NEXT I CERC Formula 95 CB1 = (9.8 * HB / .78) ^ .5 IF C$ = "Y" OR C$ = "y" THEN EBI = EB: TBI = TB: GOTO 95 HB1 = (((.78 / 9.8) ^ .5) * (HI ^ 2) * CGI * COS(TI)) ^ .4 CB1 = (9.8 * HB1 / .78) A .5 EBI = .125 * 9.8 * (HB1 ^ 2) ALF8 = SIN(TI) * CBI / CI TBI = ATH(ALF8 / ((1 - ALF8 ^ 2) ^ .5)) PLB = EBI * CBI * SIN(TBI) * COS(TBI) CERC = .77 * PLB / (DELTA * 9.8 * .65) QH = HI{ * DX * XB QW -- WW * DX * XB QB = BB * DX * XB CLS QBY = INT(QB * 86400). QHY = INT(QH * 86400) QWY = INT(QW * 86400): CERCY = INT(CERC * 86400) LOCATE 5, 5 PRINT "Root-mean square wave height at the breaker line =" LOCATE 5, 55. PRINT USING "##.##'; HB LOCATE 6, 5: PRINT "Angle of incidence of wave at the breaker line =" LOCATE 6, 54: PRINT USING "##.#'; TB * 45 / ATN(1) LOCATE 7, 5: PRINT "Period of wave = " LOCATE 7, 22: PRINT USING "##.#'; T LOCATE 8, 5" PRINT "Breaker zone width =" LOCATE 8, 26: PRINT USING "###.#'; XB LOCATE 9, 5: PRINT "Breaker index =" LOCATE 9, 21: PRINT USING "#.##-; MU LOCATE i0, 5: PRINT "Mean longshore current velocity =" LOCATE I0, 39: PRINT USING "#.##'; VM LOCATE 12, 15: PRINT "Sediment transport" LOCATE 13, 14: PRINT • (m3/seg and m3/day) " IF PI$ = "N l OR PI$ = "n" GOTO i00 LOCATE 15, 5: PRINT "BiJker formula =- LOCATE 15, 40. PRINT QB: LOCATE 15, 55: PRINT QBY 916 F.J. CAVIGLIA I00 LOCATE 16, 5: PRINT "Adapted Engelund-Hansen formula =" LOCATE 16, 40: PRINT QH: LOCATE 16, 55: PRINT QHY LOCATE 17, 5: PRINT "Watanabe formula =" LOCATE 17, 40: PRINT QW: LOCATE 17, 55: PRINT QWY LOCATE 18, 5z PRINT "CERC formula =" LOCATE 18, 40: PRINT CERC: LOCATE 18, 55: PRINT CERCY LOCATE 25, 5: INPUT "Paper output-(Y or N)"; YY$ TB = TB * 45 / ATN(1) IF YY$ = "N" OR YY$ = "n" GOTO ii0 LPRINT TAB(10); "Breaking wave height = "; USING "#.##"; HB LPRINT TAB(10); "Angle of incidence at breaking = "; USING "##.#"; TB LPRINT TAB(10); "Period of wave = "; USING "##.#"; T LPRINT TAB(10); "Breaker zone width = "; USING "###.#'; XB LPRINT TAB(10); "Breaker index = "; USING "#.##"; MU LPRINT TAB(10); "Mean longshore current velocity = "; USING "#.##"; VM LPRINT LPRINT TAB(30); "Sediment transport" LPRINT TAB(29); "(m3/seg and m3/day)" LPRINT IF PI$ = "N" OR P15 = "n" GOTO 105 LPRINT TAB(10); "Bijker formula="; TAB(45); QB; TAB(60); QBY 105 LPRINT TAB(10); "Engelund-Hansen formula="; TAB(45); QH; TAB(60); QHY LPRINT TAB(10); "Watanabe formula="; TAB(45); QW; TAB(60); QWY LPRINT TAB(10); "CERC formula="; TAB(45); CERC; TAB(60); CERCY ii0 CLS LOCATE 12, 15: PRINT "Do you want to see the longshore current" LOCATE 13, 15 INPUT "and sediment transport distributions-(Y or N)"; P$ IF P$ = "N" OR P$ = "n" GOTO 150 115 SCREEN ii CLS : LOCATE 12, 8: PRINT "Longshore current distribution ....... (i)" LOCATE 13, 8: PRINT "Sediment transport distribution ...... (2)" LOCATE 14, 44: INPUT ; A2 CLS : IF A2 = 1 GOTO 125 120 IF P15 = "N" OR PI$ = "n" GOTO 120 LOCATE i0, 15: PRINT "BiJker formula ................... (I)" LOCATE ii, 15: PRINT "Watanabe formula ................. (2)" LOCATE 12, 15~ PRINT "Adapted Engelund-Hansen formula..(3)" LOCATE 13, 47: INPUT ; A3 125 CLS LINE (100, 350)-(500, 350): LINE (100, 350)-(100, 360) LINE (500, 350)-(500, 360): LINE (300, 350)-(300, 360) LINE (200, 350)-(200, 355): LINE (400, 350)-(400, 350) LOCATE 24, 12.5: PRINT "0.0": LOCATE 24, 62: PRINT "2.0" LOCATE 24, 37.25: PRINT "I.0": LOCATE 25.4, 36.5: PRINT "X/ Xb" IF A2 = 1 THEN X$ = "V": Y$ = "---": Z$ = "Vo" IF A2 = 1 THEN LOCATE 9, 53: PRINT "Vo = "; USING "#.##'; VO IF A2 = 1 THEN LOCATE 9, 65: PRINT "m/s': GOTO 130 LOCATE 4, 20: PRINT "Sediment transport distribution" X$ = "q": Y$ = "---": Z$ = "qo" IF A3 = 1 THEN Q01 = QBMAX IF A3 = 2 THEN Q01 = WAMAX IF A3 = 3 THEN Q01 = QEMAX LOCATE 9, 53: PRINT "qo = "; USING "#.#####'; Q01 LOCATE 9, 65: PRINT "m": LOCATE 8.5, 65.5: PRINT 2 LOCATE 9, 67: PRINT "/s" 130 IF A2 = 1 THEN LOCATE 4, 20: PRINT "Longshore current distribution" IF A2 = 1 GOTO 135 IF A3 = 1 THEN LOCATE 5, 27: PRINT UBiJker formula" IF A3 = 2 THEN LOCATE 5, 26: PRINT "Watanabe formula" IF A3 = 3 THEN LOCATE 5, 20: PRINT "Adapted Engelund-Hansen formula" Longshore currents and sediment transport 135 LOCATE I0, 53: PRINT "Xb ="; USING "###.#"; XB LOCATE i0, 65: PRINT "m" LOCATE 14.5, 4: PRINT X$: LOCATE 14.501, 3: PRINT Y$ LOCATE 15.5, 3.2: PRINT Z$ LINE (i00, 350)-(100, I00): LINE (100, 350)-(90, 350) LINE (i00, I00)-(90, I00): LINE (i00, 225)-(90, 225) LINE (I00, 287.5)-(95, 287.5): LINE (I00, 162.5)-(95, 162.5) LOCATE 22.51, 8.5: PRINT "0.0": LOCATE 6.51, 8.5: PRINT "i.0" LOCATE 14.51, 8.5: PRINT "0.5" LOCATE II, 53: PRINT "Hb ="; USING "##.##"; HB LOCATE ii, 65: PRINT "m" LOCATE 12, 53: PRINT "ub = "; USING "##.#"; TB LOCATE 12, 65: PRINT .o. LOCATE 13, 53: PRINT "T = "; USING "##.#"; T LOCATE 13, 65: PRINT "s" LOCATE 14, 53: PRINT "S = I:"; USING "###.#'; 1 / S IF A2 = 1 THEN LINE (405, 111)-(548, 111): LINE (405, 111)-(405, 235) IF A2 = 1 THEN LINE (548, 111)-(548, 235) IF A2 = 1 THEN LINE (405, 235)-(548, 235): GOTO 140 LINE (405, 108)-(560, 108): LINE (405, 108)-(405, 235) LINE (560, 108)-(560, 235): LINE (405, 235)-(560, 235) 140 V1(1) = 350 FOR I = 2 TO FIN IF A2 = 1 THEN A1 = V(I) IF A3 = 1 THEN A1 = QB(I) / QBMAX IF A3 = 2 THEN A1 = QW(I) / WAMAX IF A3 = 3 THEN A1 = QH(I) / QEMAX IF I > 161 GOTO 145 X = i00 + 200 * (I - 1) / 80 Xl = 100 + 200 * (I - 2) / 80 VI(I) = 350 - 250 * A1 LINE (Xl, Vl(Z - 1))-iX, Vl(I)) NEXT I 145 LOCATE 30, i0: INPUT "More graphics-(Y or N)"; P$ IF P$ = my. OR P$ = "y" THEN A2 = 0: A3 = 0: GOTO 115 150 A2 = 0 : A3 -- 0 CLS : SCREEN 0 LOCATE 12, 15 : INPUT "More calculations- (Y or N) "; XX$ IF XX$ = "Y" OR XX$ -- "y" THEN CLEAR: GOTO I0 CLS : END 917