1 APPLICABILITY OF SUSPENDED SEDIMENT CONCENTRATION FORMULAE TO LARGE-SCALE BEACH MORPHOLOGICAL CHANGES

Study of beach morphological changes under storm conditions and its prediction capability are of paramount importance in coastal zone management. Seabed sediment is picked up violently in and outside the surf zone due to suspension mechanisms, therefore a considerable amount of sand is transported in coastal waters due to such mechanisms. For the construction of an accurate beach morphological model, it is necessary to elucidate the sediment suspension and to introduce it properly into the modelling of sediment transport. Jayaratne and Shibayama (2007) developed a complete set of explicit theoretical formulae to predict the time-averaged concentration on sandy beaches due to three suspension mechanisms: a) vortical motion over wave-generated sand ripples, b) from sheet flow, and c) turbulent motion under breaking waves. The present paper focuses on the development of a quasi-3D beach deformation model using the sediment concentration models of Jayaratne and Shibayama (2007), the bed load model of Watanabe (1982), the wave propagation model of Onaka et al. (1988), the nearshore current model of Philips (1977) and the undertow model of Okayasu et al. (1990) to predict the large-scale morphodynamics of sandy beaches. The predicted beach profiles and total sediment transport rates were compared with two sets of large-scale laboratory experimental data [Kajima et al. (1983); Kraus and Larson (1988)] and Seisho beach at Kanagawa Prefecture, Japan. It can be concluded that the present numerical model is capable of predicting sediment transport direction, onoffshore sand bar formation and the general trend of the beach profiles of large-scale erosiveand accretive-type sandy beaches to a satisfactory level.


INTRODUCTION
Utilisation of the coastal environment has steadily been increasing due to human activities such as transportation, tourism, leisure and habitation over the recent past.Protection of coastal areas from wind-induced high waves and extreme events exacerbated by global warming has become a task of paramount importance today.Coastal zone management requires a quantitative and qualitative predictive capability that represents beach profile changes due to such events.To deal with this problem, numerical tools have gradually evolved as a powerful technique to predict total cross-shore sediment transport rates (bed load and suspended load) and beach profile evolution.But due to the complexity of the sediment transport processes, a full description of these mechanisms (e.g.suspension) has not reached a satisfactory level.Even at the current stage of research, the research is still far from complete and there is no thorough understanding of all phenomena induced by coastal waves, currents and tides on coastlines or on coastal structures.
Sediment suspension usually occurs in and outside the surf zone due to vortices generated by sand ripples, moving of the bottom layer with high bed shear stresses and turbulence generated by wave breaking (Jayaratne and Shibayama 2007).When surface waves propagate into the shallow water region, sand ripples are formed at the bottom of the sea.Vortex ripples enhance the separation of the flow and production of turbulence, which induce dramatic changes of the velocity distribution in the bottom boundary layer together with a considerable amount of energy dissipation and sand agitation (Sato 1987).Sediment suspension over vortex ripples thus plays a paramount role in various coastal processes such as sediment transport, wave attenuation and mass transport due to waves and currents.The characteristics of waves, currents and sediment particle movement in the sheet flow regime are evidently completely different over the ripple-covered bed.The rippled seabed gives way to the sheet flow region in the breaker zone, in particular under extreme storm conditions.Under such conditions, when large wave heights cause large near-bed velocities, such ripples are washed away and the sand is suspended over the thin layer close to the bed.Sand movement in the surf zone is more difficult to describe as it is influenced by the turbulence and the organised motion generated by the breaking waves as well as by the asymmetric oscillatory motion under shallow-water waves.The mechanism of sediment suspension at the breaking point is particularly important, where the turbulence and large-scale vortices that are generated agitate a large amount of sediments to form a suspended sand cloud with extremely high concentration (Sato et al. 1990).
Therefore, developing an accurate complete model to obtain large-scale beach morphological changes, dominated by sediment suspension is rather a challenging task.The present paper discusses the elements of the proposed quasi-3D beach deformation model including the time-averaged suspended sediment concentration models of Jayaratne and Shibayama (2007).Finally, the applicability and robustness of this model is examined by large-scale experimental and natural beach data.

MODEL FORMULATION
The proposed quasi-3D beach deformation model primarily consists of a five-tier hierarchy of modules: a) wave propagation, b) nearshore current field, c) vertical distribution of undertow, d) sediment transport models (bed load and suspended load), and e) bottom topography change.The following sub sections explain each module and its formulations in detail.

a) Wave Propagation Model
Evaluation of wave transformation in the shallow water region has been researched using numerical models by many researchers.Berkhoff (1972) used the mild-slope equation to study the combined effect of wave refraction and diffraction on gradually varying bottom topography.Watanabe and Maruyama (1986) proposed a time-dependant version of the mild-slope equation which can deal with wave breaking and attenuation in the surf zone.None of the above work included the wave-current interaction effect.Booij (1981), Liu (1983) and Kirby (1984), to name a few, proposed wave field with varying current models.However, these models consisted of elliptic-type differential equations which are in general difficult to solve numerically.
In order to overcome this issue, Ohnaka et al. (1988) proposed a new set of time-dependant mildslope equation by separating the surface elevation (ζ) and the flow rate vector (Q) from the elliptical model of Kirby (1984).The wave propagation model in the present beach deformation model is based on the solutions of Ohnaka et al. (1988).Equations ( 1)-( 3) give the surface elevation (ζ) and x and y directional flow rates (Q x and Q y ).These three equations were solved by using a finite difference scheme.The study area was divided into grid cells and a staggered mesh scheme was employed.To avoid unphysical reflection at the offshore boundary, it was set in terms of the surface elevation as described by Ohnaka et al. (1988).
where ζ=surface elevation, t=time, x=cross-shore direction, y=alongshore direction, U=steady current in x direction, V= steady current in y direction, Q x =flow rate in x direction, Q y =flow rate in y direction, ω=apparent angular frequency, ζ=intrinsic angular frequency, c=wave celerity, The nearshore current system includes longshore currents, shoreward-directed currents and rip currents.By assuming the nearshore current field varies slowly in space and in time on a scale that was large compared to a typical wavelength and wave period of sea and swell waves, Phillips (1977) obtained derivations for mean water level ) ( and mean velocities in x and y directions (U and V) from the mass conservation and the Navier-Stokes equations.Equations ( 4)-(6) show ,  U and V respectively.These three equations were also solved by using a finite difference scheme.
where h=water depth, F x , F y =friction terms in x and y direction, M x , M y =diffusion terms in x and y direction and R x , R y =radiation stress terms in x and y direction.

c) Vertical Distribution of Undertow Model
It is indispensable to estimate the distribution of the undertow in order to predict the sediment transport and the material diffusion in the surf zone.Many of the present undertow models have been developed from local properties of waves and turbulence in the surf zone.Okayasu et al. (1990) developed a cross-shore vertically 2D distribution model ) , ( v u by considering the accurate description of wave attenuation, energy distribution of turbulence including energy of the organized large vortices, based on the actual wave breaking phenomena.By considering mean shear stress (η), eddy viscosity (ν e ), molecular viscosity (ν) and steady currents (U and V) in x and y directions, the following set of equations [Eqs.( 7)-( 10)] was derived by Okayasu et al. (1990).
vertical distributions of undertow in x and y directions, a ηx , b ηx , a ηy , b ηy =coefficient of shear stresses in x and y directions, a ν =eddy viscosity, ν=molecular viscosity, z / =vertical elevation measured from the seabed, C 1 , C 2 =integral constants which were obtained in terms of the mass transport by waves and d t =water depth at wave trough.

d) Sediment Transport Models
In order to evaluate quantitatively the bottom topography change caused by waves and currents, total sediment transport rate must be calculated.Due to the complexity of the sediment transport problem, governing equations of sediment movement based on fundamental physics have not been established.For simplicity, the sediment transport formulae for cross-shore and longshore transport have been developed based on a macroscopic approach to oscillatory and steady fluid motion.The following two sub sections discuss i) a bed load due to mean currents and wave action, and ii) suspended sediment transport rate models used in the present study.

i) Bed load models due to mean current and wave action
Adaptation of an appropriate bed load model into a large-scale beach deformation model depends on the parameters which decide the direction of sediment transport, applicability of ranges of wave, flow and sediment conditions, and limitations associated with them.By considering the bed load models of Sunamura and Horikawa (1974), Shibayama and Horikawa (1982) and Watanabe (1982), the Watanabe (1982) model was selected as the most suitable in the present study.The sensitivity study related to the selection of bed load model is given in Takayama et al. (2012).Watanabe (1982) generalised sediment transport rate formulae and separated them into transport rate due to current (q c ) and due to direct action of waves (q w ) in order to utilise past studies and to simplify the resulting formulae.It may seem unreasonable to treat the transport rates under a combined wave-current action as a linear summation of those due to waves and currents respectively.However, it is much more practical to separate sediment transport formulae under direct action of waves and mean steady current where the resulting formulae will be compatible with each other and consistent with previous sediment transport studies.
The transport rates in x and y directions, q cx and q cy are estimated by Eqs. ( 11) and ( 12) by modifying the longshore drift distribution formula proposed by Komar (1977).These models are applicable for general nearshore current fields and accountable for the critical shear stress effect.
The bed load transport rates due to the direct action of waves in a wave-current coexistent field (q w ) is more difficult to estimate than that of mean current because the net transport rate due to both forward and backward transport rates should be treated.This makes more difficult when the effects of wave breaking and superposition of various wave trains due to reflection and diffraction are taken into consideration.According to Watanabe (1982), net bed load transport rate under a single progressive wave train is estimated as; where  13) and using the concept of power model, Eq. ( 14) can be established for bed load sediment transport due to waves by Watanabe et al. (1986).
where A w =dimensionless coefficient, assigned a different value at each calculation grid and u b =amplitude of near-bottom orbital velocity.
In order to smooth an unrealistic discontinuity at the null point in the distribution of sediment transport rate, a transport direction function, F D was introduced by Watanabe et al. (1986).It is noticed that П c is a crucial parameter determining for a given test case whether model predicts erosive-, accretive-or combination of both types of sandy beaches therefore in the present study П c is determined case by case from the published large-scale experimental data.This is discussed in the Model Verification and Discussion section and further information can be found in Takayama et al. (2012).Finally, combining Eqs. ( 14) and ( 15), bed load transport rates due to wave action in x and y directions (q wx , q wy ) can be given by Eqs. ( 16) and (17 where α=wave direction angle.

ii) Suspended sediment load models
In the Jayaratne and Shibayama (2007) concentration models, bed reference concentration (c a ) and diffusion coefficient (ε s ) were evaluated in terms of wave, sediment and flow properties using dimensional analysis and a best-fit technique.Bed reference levels (z 0 ) were taken where concentration could be measured without affecting the bed.Concentration profiles [c(z)] were derived from the steady diffusion equation and efficiency factors were calibrated with the help of a large amount of published data.Concentration profiles over rippled bed and lower sheet flow layer (pick-up layer) tend to have an exponential form with water depth while the concentration profiles over upper sheet flow layer (suspension layer) and turbulence included by wave breaking behave in a power relationship with water depth.Equations ( 18)-( 29) show the bed reference concentration, diffusion coefficients and concentration profiles developed under three suspension mechanisms.The bed reference level for rippled bed was taken at z=η/2 while for pick-up layer in sheet flow, it was considered to be at z=25d.On the other hand, the bed reference level for suspension layer in sheet flow was taken at z=100d and the same reference level was considered for breaking agitation.Figures 1-3 illustrate measured and computed time-averaged suspended sediment concentration profiles over water depth under large-scale rippled bed, sheet flow and breaking agitation respectively. Suspension under rippled bed: For suspension layer : where M s =parameter related to diffusion of the sediment. Suspension under breaking agitation: where û=wave orbital velocity at the location (grid point) to be considered, û b =parameter to represent the intensity of near-bottom flow at the wave breaking point (Sato et al. 1990), D B =average rate of energy dissipation due to wave breaking, computed from empirical relationships in Rattanapitikon and Shibayama (1998), k 1 =numerical constant depends on the wave period (T), k 2 , k 4 =numerical constants, k 3 =numerical constant and value in bracket varies from 0.3 at breaking point to 1.0 at transition point (Rattanapitikon and Shibayama 1996), x=position in the cross-shore direction and subscript b and t denote the distances from breaking and transition points respectively.

Suspended sediment concentration profiles [c(z)] evaluated above and vertical distributions of undertow velocities
) , ( v u obtained by Okayasu et al. (1990) are integrated with respect to z in order to obtain the suspended loads in x and y directions, q sx and q sy [Eqs.( 30) and ( 31 where δ s =elevation below the bottom surface where there is no effective movement of sand.

e) Bottom Topography Change
The change of local bottom elevation (z b ) or water depth (h) can be readily determined if the spatial distribution of sediment transport rates is given.The bottom topography change can be evaluated by solving the conservation of sediment mass equation.
where q tx and q ty =total sediment transport rates in x and y directions respectively and determined by the bed load models of Watanabe et al. (1986) under wave-current coexistent field and the suspended sediment concentration models of Jayaratne and Shibayama (2007) under rippled bed, sheet flow and breaking agitation.
The change of bottom topography in turn affects the wave and current fields.It is impracticable to handle a fully unsteady and interactive model of wave-current coexistent field and beach deformation therefore it is assumed in the present numerical model wave-current field remains unchanged for a certain time duration even if the local water depth changes.
Although the wave-current field varies with the beach deformation, this is not the sole mechanism to achieve equilibrium status of beach profiles.In nature, if the local slope becomes steep, sediment particles move downwards due to gravity.Therefore, in order to incorporate the local slope effect, Eq. ( 32) is modified by an additional term as follows (Watanabe et al. 1986).

MODEL VERIFICATION AND DISCUSSION
The simulation results of cross-shore beach profiles and distribution of on-offshore sediment transport rates were verified with the help of large-scale laboratory experimental data [Kajima et al. 1983;Kraus and Larson 1988].The natural beach data at Seisho beach in Kanagawa Prefecture, Japan was also used to obtain the alongshore distribution of bottom topography change under a storm event in 2007 (Typhoon Fitow).a) Large-scale laboratory beach profiles Kajima et al. (1983) measured beach profiles in a large-scale wave flume (205.0×3.4×6.0 m) at the Central Research Institute of Electrical Power Industry (CRIEPI) of Japan in 1981.The total number of test cases performed by Kajima et al. (1983) is 24, and all cases started with a uniform initial beach slope, ranging from 1/50 to 1/10.Kraus and Larson (1988) reported beach profile changes in a large wave flume (193.5×4.6×6.1 m) at the Coastal Engineering Research Centre (CERC) of the US Army Corps of Engineers Waterways Experiment Station, which were measured in 1956-1957 and 1962.The total number of test cases carried out by Kraus and Larson (1988) is 18 and all tests cases started with an initial beach slope of 1/15.Two tests were performed with an irregular initial beach slope.Tables 1  and 2 show the detail test conditions used for the calibration of bed load models and comparison of predicted cross-shore beach profiles.The bed load formulations given in Eqs. ( 11)-( 17) consist of various empirical coefficients.Among them, obtaining suitable values for the limiting parameter П c is a challenging task.By considering initial wave, sediment and beach conditions (H, T, d and m) given in Tables 1 and 2, П c was computed for each test case and plotted against the Sunamura and Horikawa (1974) parameter C. Figure 5 illustrates the categorisation of П c into four different zones: 1) Erosion, 2) Accumulation, 3) and 4) Combination of both.The limit parameter varies from smaller values in zone 3 to larger values in zone 4. According to the analysed data, in some test cases, the sediment transport direction is unexpectedly reversed, especially in the upper section of zone 4.This phenomenon was observed and supported by the laboratory experiments conducted by Shibayama and Horikawa (1982).Figures 6-7 illustrate the predicted large-scale beach profiles and total sediment transport rates superimposed with measured quantities.In Case 1.1 a larger on-shore directed sand bar is formed compared to the measured bar at t=9.0 hr while in Case 1.6 both predicted and measured sand bars at t=12.0 hr match quite well.Predicted total sediment transport rate distributions in Case 1.1 and 1.6 also match well with the measured transport rates, apart from some over-prediction at the bar crest levels.These two test cases were performed under moderate wave conditions.In Case 3.1 shoreward slope of the predicted bar is steeper than that of measured at t=6.5 hr but the offshore slopes of both bars match well.The sand bar in Case 4.4 is offshore-directed compared to the measured bar at t=10.0 hr and this could be due to the change of breaker line.Large deviation in sediment transport rates at the bar crest level can be seen in Case 4.4.The reason for this discrepancy is identified as the test Case 3.1 and 4.4 were carried out under large wave height and long period wave conditions.In general, predicted and measured beach profile shape of four tests cases of Kajima et al. (1983) follow the same trend.The causes can be due to the long period wave condition of Case 201 and large wave condition of Case 400.Due to the limitation of insertion of swash dynamics at the present level of study, a large discrepancy in beach slope and sediment transport rates in swash zone of Case 400 is noticed.Good agreement of beach profiles in Case 801 at t=5.0 hr can be noticed as this test case was carried out with moderate wave conditions.However, the predicted sand bar of this test case is offshore directed due to the change in breaker line.An over-prediction of offshore bar at t=9.5 hr can be observed in Case 901 and the reason can be due to the presence of large wave height and long period laboratory conditions.Apart from some localised under-prediction and over-prediction, the general trend of the measured beach profiles of Kraus and Larson (1988) are simulated satisfactorily well by the proposed numerical model.

b) Natural beach topography and longshore distribution
The proposed quasi-3D model is applied for a 2.5 km long coastal stretch of Seisho beach (40 km SW of Yokohama) located in Kanagawa Prefecture, Japan.This beach was heavily eroded and collapse of the by-pass was caused by the typhoon No 9 called 'Fitow' in September 2007. Figure 8(a) shows the location map of the coast, and 8(b) illustrates the aerial view of the beach.The study area covers a bridge called 'Kinpa', located in further west of the beach.Another feature in this beach area is that it consists of a large sand cliff at the offshore.Two field monitoring programmes were carried out by the Kanagawa Prefectural government in March and October 2007 in this beach and since there was no typhoon attack except due to typhoon Fitow, it was assumed that the large deformation to the beach profiles was due to the typhoon Fitow event.The offshore wave climate in this area was measured by the Hiratsuka Observation Tower, which is located 1.0 km from onshore to the Seisho beach at 20.0 m water depth.Table 3 shows the wave, sediment and computed parameters used in the numerical simulation.The initial topography data (water levels) just before the typhoon Fitow, which is the main input to the numerical model, was not available therefore it was made at 10×10 m grid points from the field measurements performed in March and October by using the Spline Interpolation Method.To determine sediment transport rates and bottom topography change, wave height and nearshore distributions should be predicted accurately.Figure 10 shows the model simulation results of (a) wave height distribution, (b) current distribution, and (c) aerial view of the entire study area of 2.5×1.0 km.Large wave heights and nearshore currents generated close to the Kinpa bridge and this could be due to the sudden change of water depth.Simulated results of combined effect of bed load and suspended loads (Fig. 11b) show better results than the bed load effect itself (Fig. 11c), in comparison to the measured beach topography at t=6.0 hr (Fig. 11a).A rather steep slope close to the shoreline with sand bar formation is observed in

CONCLUDING REMARKS
The following conclusions and recommendations can be made from the study.1.A quasi-3D large-scale beach morphological model is developed using the suspended sediment concentration models of Jayaratne and Shibayama (2007) and applied to two large-scale laboratory test series and a natural beach in Japan subjected to a typhoon attack. 0 gravity of sand, L 0 =deepwater wavelength, П c =critical value of П (Limiting parameter) at the null point where the net sediment transport rate is zero and κ d =coefficient of order of unity which controls the degree of the change in transport rate around the null point.
20) where ν= kinematic viscosity, η=ripple height, λ=ripple length, u *wc =shear velocity under wave-current coexistent field, A b =orbital amplitude of fluid just above the boundary layer

Figure 1 .
Figure 1.Example of measured and predicted suspended sediment concentration profiles over large-scale rippled bed (Jayaratne and Shibayama 2007). Suspension from sheet flow: For pick-up layer : ) 100 ( d z 

Figure 2 .
Figure 2. Example of measured and predicted suspended sediment concentration profiles over large-scale sheet flow (Jayaratne and Shibayama 2007).The blue upwards arrow in the third panel indicates the pick-up layer while blue downward arrow in the same panel indicates the suspension layer.

Figure 4
Figure 4 depicts the flow chart of a quasi-3D large-scale beach deformation model developed from five-tier hierarchy of modules.The initial water depth (h), incident wave height (H), wave period (T) and average sediment diameter (d) are the main input parameters for the model.
and 400, sand bars at t=5.0 hr and 4.0 hr are under-predicted by the numerical model.

Figure 6 .
Figure 6.Measured and predicted cross-shore beach profiles and sediment transport rates (Kajima et al. 1983).

Figure 7 .
Figure 7. Measured and predicted cross-shore beach profiles and sediment transport rates (Kraus and Larson 1988).

Figure 8
Figure 8(a) Location map showing Seisho beach (blue circle) and Yokohama City (red circle), (b) Aerial view of eroded beach and partly damaged road embankment (Source: Yokohama National Highway Work Office, Ministry of Land, Infrastructure, Transport and Tourism, Kanto Regional Development Bereau).

Figure 9
Figure 9(a) Measured incoming wave during typhoon Fitow (1-10 September 2007 at Haratsuka Observation Tower), (b) Calibration of Limit parameter Пc against Sunamura and Horikawa (1974) parameter C for Seisho beach data, shown by red dot in region 1, and (c) Initial beach topography of the study area, t=0 hr.

Figure 10 (
Figure 10(a) Predicted wave height distribution, (b) Predicted nearshore current field, and (c) Aerial view of study area (purple) and severely damaged section (red).
Fig 11a-11c.The reason for this feature is due to the shifting of breaker line.The circles and arrows in the figure indicate the area of complete damage, which is near the Kinpa bridge.
Figure 11(a) Measured data at t=6.0 hr generated from field measurements in March and October 2007 (b) Predicted topography from combined effect of bed load and suspended loads at t=6.0 hr, (c) Predicted topography from the bed load effect only at t=6.0 hr, and (d) Sediment transport direction during typhoon event in 2007.

Figure 12
Figure 12(a) Predicted topography from combined effect of bed load and suspended loads at t=5.0 hr in x,y,z directions, (b) Predicted topography from the bed load effect only at t=5.0 hr in x,y,z directions.
2. The limiting parameter ofWatanabe et al. (1986)  П c is calibrated in each test case by considering theSunamura and Horikawa (1974) parameter C. 3. The numerical model is capable of predicting sediment transport direction, on-offshore sand bar formation and general trend of the beach profile of erosive-and accretive-type sandy beaches to a satisfactory level.4. It could handle variety of input wave, flow, tidal and sediment data at a reasonable computing cost. 5.The model will be further verified with the natural beach data of different coastal environments for wider applicability.