ON THE CHOICE OF RANDOM WAVE SIMULATION IN THE SURF ZONE PROCESSES

In this paper, the two common approaches to account for wave randomness, the spectral approach and the wave-bywave approach, are compared through numerical experiments conducted with the coupling of a surf zone hydrodynamic model and a bedload sediment transport model. Special attention is paid to the wave nonlinearity and net cross-shore bedload transport predictions. The two approaches are found to have negligible difference in their predictions of certain average hydrodynamics, such as wave heights, set-up and undertow. However, the wave-bywave approach outperforms the spectral approach in the wave nonlinearity prediction, and the two approaches differ significantly in their predictions of wave-induced net cross-shore bedload transport which strongly depends on wave nonlinearity. This suggests the necessity of using the wave-by-wave approach. The computational efficiency of the wave-by-wave approach is also discussed.


INTRODUCTION
In the surf zone, the cross-shore sediment transport is of crucial importance to the prediction of beach bathymetry evolution.The leading driving force is the random waves which have intense interaction with the sea bottom as they shoal, break and travel across the surf zone as broken waves.Since the net cross-shore sediment transport over a wave cycle is the delicate balance between the onshore and offshore components, to successfully model the propagation of the sea waves, especially to account for their randomness, is one of the most important tasks facing any surf zone process model.Two main approaches, the spectral approach and the wave-by-wave approach, are commonly used to represent wave randomness.The spectral approach is based on a single characteristic wave, while the wave-by-wave approach represents the random waves as a series of individual periodic waves.In this paper, the surf zone hydrodynamic model by Tajima and Madsen (2006), and the bedload transport model by Gonzalez-Rodriguez and Madsen (2007), are coupled using both aforementioned approaches.By comparing the simulation results, especially the wave nonlinearity characteristics and net bedload transport rate predictions, the practicality of the two approaches is discussed.Since this work is primarily aimed at modeling the gravity waves and their associated net cross-shore bedload transport, some equally important components of surf zone processes, such as those in the swash zone and contributions of infragravity waves, are not included at the present time.

APPROACHES TO REPRESENT WAVE RANDOMNESS
The spectral approach usually starts with certain assumptions on the shape of the wave spectrum, e.g. the narrow-banded spectrum and Rayleigh distributed wave heights.With these assumptions, the governing equation may be integrated analytically over frequency space to obtain predictions of a single representative mean wave such as a wave with the significant or root-mean-square wave height and the average or peak period.All the subsequent simulations, including the wave set-up, undertow, wave nonlinearities and sediment transport, are consequently based on this single representative wave.Therefore, the response time of the surf zone system is essentially assumed to be infinite when the spectral approach is used to simulate random waves.Typical examples are Thornton and Guza (1983) and Tajima and Madsen (2006).These models have been shown to be reasonably accurate in their hydrodynamic predictions.
Alternatively, the wave-by-wave approach is based on the joint wave height and period probability distribution of random waves.The general concept is to select a number of individual periodic waves, each of them with a certain probability, to represent the whole distribution.In terms of numerical modeling, the selected waves are simulated with periodic wave models and their contributions are averaged over their probabilities and periods.The basic assumptions inherent in this approach are that the wave-wave interaction is negligible and the scale of the surf zone system response time is the same as that of the individual wave period.Most models using the wave-by-wave approach have shown good agreement with observations (Mase and Iwagaki, 1982;Mizuguchi, 1982;Dally and Dean, 1986;Dally, 1992;Van Rijn and Wijnberg, 1996;Grasmeijer and Ruessink, 2003).Nevertheless, some researchers also found that the wave-by-wave approach does not outperform the spectral approach in predicting some hydrodynamic quantities such as the wave heights and currents (Grasmeijer and Ruessink, 2003).Considering the low computational efficiency, the wave-by-wave approach seems less attractive.However, some other quantities, especially those which have a strong nonlinear dependency on the wave heights, like the wave asymmetry and skewness and the sediment transport rate, are suspected to be poorly predicted by the spectral approach with a single mean wave.On the other hand, since the wave-by-wave approach is able to account for the variability among individual waves, this approach is believed to be more realistic to model these highly nonlinear quantities.Tajima and Madsen (2006) developed a theoretical model (TM06 hereafter) to predict the surf zone hydrodynamic processes driven by either obliquely-incident or normally-incident waves along a long straight beach.The model consists of three parts: the wave model, the surface roller model and the current model.The wave and surface roller models are based on the linear wave energy balance equation.In the current model, the momentum equation is integrated above the wave trough level and over the entire water column, respectively, to determine the mean shear stresses at the wave trough level and the sea bottom.The vertical current profile is analytically calculated by introducing a simple turbulent eddy viscosity model.The wave and current interaction in the bottom friction formulation is considered by using Madsen's (1994) modified wave-current boundary layer model.The TM06 model is extended to the random wave case using the spectral approach with the assumptions that the wave spectrum is narrow-banded (periods are all equal to the peak period) and with wave heights following the Rayleigh distribution, and the wave energy balance equation is then analytically modified to predict the root-mean-square wave.All the subsequent simulations, including wave nonlinearity and waveinduced current, are based on this single representative wave.

MODEL DESCRIPTION
As the wave propagates into the surf zone, its nonlinear features start to increase so that the shape of the wave no longer is symmetric, but it becomes skewed (peaked crest and flat trough) and asymmetric (forward-leaning crest).This is also observed in the near-bottom wave orbital velocity profile.Gonzalez-Rodriguez and Madsen (2007) characterized the nonlinear near-bottom wave orbital velocity variation by the parameters: T c , T t , T cn , T cp , T tn , T tp , U c and U as defined in Fig. 1.Based on this, two parameters are used to quantify the two nonlinear features: skewness, Skw=2U c /U-1, and asymmetry, Asy=1-2T c /T.Notice that for a sinusoidal wave, these two parameters are both zero.TM06 predicts these two parameters based on linear wave characteristic and topography information with a set of empirical transform relationships which were developed from numerical experiments based on the Boussinesq type equations.The two nonlinearity features will yield an imbalance in the bottom shear stress and consequently induce a net bedload transport.
Gonzalez-Rodriguez and Madsen (2007) developed a conceptual bedload transport model (GRM07 hereafter) to predict the bedload transport which includes the effect of these two nonlinearities.It uses Madsen's (1991) bedload formula to calculate the instantaneous bedload transport rate: where q SB (t) is the volume of sediment transported per unit time and width in the direction of the instantaneous bottom shear stress, β is the bottom slope in the direction of transport (upslope is positive), s is the ratio of sediment and water densities, τ cr,β is the critical Shields parameter for the initiation of motion, α β is given by: tan tan tan tan and 50 s   and 30 m   are the values of the angles of static and moving friction recommended by Madsen (1991).The instantaneous shear stress is given by: 1 where u b (t) is the instantaneous near-bottom velocity, l τ (t) is the time lag between τ b (t) and u b (t), and f w is the friction factor.Both l τ and f w are obtained by applying Madsen's (1994) modified wave boundary layer model.The near-bottom velocity time series u b (t) is reconstructed with nonlinear wave orbital velocity amplitude, skewness and asymmetry predicted by the TM06 model.Obviously, with a constant friction factor, a purely asymmetric but non-skewed wave cannot induce net bedload transport, which contradicts many experiments.To fix this problem, the model uses time varying f w and phase lag l τ .With reference to Fig. 1 for notation, it first calculates the friction factor at the crest of the wave based on an equivalent sinusoidal wave with the velocity amplitude equal to U c and the period equal to T cp .Then the friction factor at the wave trough is calculated analogously with velocity amplitude U t =U-U c and period T tn .The time-varying friction factor is finally obtained by linearly interpolating the two friction factors over the period.Same method is used to calculate the phase lag.The net bedload transport rate is obtained by numerically averaging the instantaneous transport rate over the entire wave period.The model is tested against many measurements and has been shown to have reasonable accuracy.
In this paper, we couple the TM06 model and the GRM07 model using both the wave-by-wave approach and the spectral approach.The TM06 model does not have a swash zone process model at the shore-side boundary where the low frequency motion becomes significant.Also since we are primarily interested in gravity wave induced processes and the associated net cross-shore bedload transport, especially the effect of different approaches to account for wave randomness, the long period motions, e.g.surf beat and swash zone processes, and the processes which are strongly related to long period motions such as the suspended sediment transport, are not included in this work but are left for future study.

REPRESENTATIVE INCIDENT WAVE SET FOR WAVE-BY-WAVE APPROACH
We apply the models to Test 3 of the moveable bed experiments which were conducted in the Large Scale Facility for Sediment Transport research (LSTF) (Wang, et al., 2002).In this test, the barred bottom profile is constructed from quartz sand of median diameter D 50 = 0.15 mm.The obliquely-incident waves follow the TMA spectrum with the significant wave height and the peak period 0.23 m and 3 sec, respectively.The incident angle is fixed at 10 degree.A plunging type breaker is mostly observed during the test.A long-shore circulation system is deployed to achieve long-shore uniformity.During Case 2 of this test which lasted 100 min, the time series of the surface level variations and the flow velocity are measured at ten locations along a cross-shore profile.For each measured time series, we perform a spectrum analysis and cut the parts outside the frequency range 0.2 Hz to 4 Hz to exclude the signal noise and the low frequency motions which are not considered in this work.
At the incident boundary, from a zero-down-crossing analysis of the measured surface level fluctuation, we can indentify roughly 3,000 individual waves.We first define the RMS-wave height H rms and average period T ave as: Where H i and T i are the wave height and period of the i-th observed wave, N is the total number of waves.After non-dimensionalizing the wave heights and the periods by H rms and T ave , respectively, a non-dimensional joint wave height, ξ, and period, τ, distribution is obtained, as shown in the left panel of Fig. 2. To select the representative incident wave set which is used in the wave-by-wave approach to represent this joint distribution, we first discretize the ξ-τ plane with rectangular grids of a certain resolution as shown in the right panel of Fig. 2. For each grid, we can choose a representative wave for which the wave height ξ ij , period τ ij and associated probability p ij are defined as: where N ij is the number of waves within the ij-th grid.Notice that, with this definition, the selected representative incident wave set will give exactly the same rms-wave height and average period as the observed waves at the incident boundary.With 0.1-by-0.1 resolution, a representative incident wave set with 247 waves is determined.To quantify the goodness of it, we can simply look at the standard deviations of the wave heights and periods which are defined as: A good representation should give values of these two parameters close to those based on the original observed waves which are defined as: Now with 247 waves, the relative errors between Eq. 6 and Eq. 7 are 0.33% and 0.30% for wave height and period, respectively.We will first take this error as insignificant and will return to this later.

HYDRODYNAMIC PREDICTIONS
We first compare the hydrodynamic predictions given by the two approaches.As shown in Fig. 3, both approaches reasonably predict the observed root-mean-square wave heights H rms .The relative error is roughly around 1-6% except for the two most shoreward points where both approaches underestimate the wave height by 20-50%.The reason for this is that these points are close to or within the region where swash zone processes are significant.This region is, as aforementioned, not included in the TM06 model.The difference between the two approaches is negligible, which has also been concluded in some previous work (Grasmeijer and Ruessink, 2003).
Similar to the H rms prediction, both approaches reasonably well predict the RMS near-bottom wave orbital velocity height U brms .The wave-by-wave approach is better than the spectral approach on the upslope side of the bar by reducing the relative error from roughly 7% to 1%.This region is where the strongest breaking occurs, so the good predictive skill gives confidence in the following current and bedload prediction.Unlike the H rms prediction, both approaches surprisingly give good predictions of U brms near the shoreline where the model is not totally applicable.The reason is unclear so we cannot comment on this at this point.In the surf zone, some hydrodynamic quantities such as the wave height and the wave orbital velocity have immediate response to the variation of individual waves.The concept of the wave-bywave approach, i.e. individual waves are simulated as if they are alone, is consistent with the nature of these quantities.Some other quantities, including the set-up and undertow, are believed to have slower response time, so that applicability of the wave-by-wave approach to their prediction is questionable.Similarly, the spectral approach which assumes the system response time to be infinite by using a single representative wave also has problems in its applicability.Nevertheless, since the two approaches are at the two limits of system response time, their predictions can quantify two limits within which the actual prediction should fall.The gap between the two limits also indicates how strong the effect of system response time is.
Fig. 4 shows the predictions of average wave set-up and undertow which, as mentioned before, have slow response time.Both approaches give good predictions of the wave set-up with the wave-by-wave approach slightly better.Their difference is generally minimal.Similar to the set-up, the undertow predictions based on the two approaches are also virtually identical.The models generally successfully predict the vertical undertow profiles, except at the most shoreward and the most seaward locations.The large overestimation at the most shoreward location is due to the swash zone processes which are not considered in this work.Meanwhile, the large error at the most seaward location could be caused by the downward penetration of the plunging breaker which will "push" the velocity profile shoreward at middle depth.
The closeness of predictions obtained by the two approaches, as shown in Fig. 3 and Fig. 4, generally indicates that the effect of system response time on the average hydrodynamic prediction is negligible.Considering computational efficiency, using the spectral approach is recommended for the purpose of predicting these quantities.
However, the difference between the predictions of average wave skewness and asymmetry given by the two approaches is large, as shown in Fig. 5.For asymmetry prediction, the spectral approach always overestimates it with a relative error that keeps growing shoreward to over a factor of 2. The wave-by-wave approach starts with a slight underestimation but the error is mostly within 20%.For skewness prediction, both approaches mostly overestimate it by over a factor of 2. However, unlike the spectral approach which has huge errors right from the incident boundary, the wave-by-wave approach reasonably predicts skewness until the crest of the sand bar (the third point from the left) where the most severe breaking occurs.Since the empirical formula to calculate the skewness is obtained through numerical experiments in which the waves are non-broken, the effect of wave breaking on the skewness is missing from the model.That could be a reason for the poor predictive skills of skewness.We will discuss the effect of overestimating skewness later.In general, we can conclude that the wave-by-wave approach has better ability to predict the average wave nonlinearities than the spectral approach.

NET CROSS-SHORE BEDLOAD TRANSPORT PREDICTION
Net cross-shore bedload transport is calculated with GRM07 model.The bottom shear stress required by the Madsen's (1991) bedload transport formula, as shown in Eq.1, should involve the combination of waves and currents.Therefore, we first decompose the total bottom shear stress into a wave part and a current part: The current shear stress is obtained from the TM06 model.The wave shear stress is obtained as in the original GRM07 model, except that the current shear stress, which is specified by τ cb , is added in Madsen's (1994) boundary layer model analysis.We also assume that the bottom condition is within the sheet flow regime and use Herrmann and Madsen's (2007) sheet flow roughness model, where  and cr  are the Shields parameter and its critical value for initiation of motion, and D n is the nominal diameter (D n ≈1.1D 50 ).The comparison between the two approaches is shown in Fig. 6.The wave-by-wave approach gives significantly larger net onshore transport rate than the spectral approach.This seems to contradict Fig. 5 where we see that the wave-by-wave approach gives smaller average wave nonlinearity and consequently the net bedload transport should be smaller.The reason is that the wave-by-wave approach can account for the variability of the wave nonlinearities among individual waves.Namely, larger waves can have larger wave nonlinearities, unlike the spectral approach in which every wave has the same wave nonlinearity as the single representative wave.Since most of the onshore sediment transport is caused by a small number of the very large waves, the wave-by-wave approach can yield a larger long-term average net onshore transport rate.This finding indicates that the net bedload transport has strong nonlinear dependency and the variability among individual waves is important.Because most natural sands including the sands used in the LSTF experiments can react immediately to the flow motion, the bedload transport is believed to have immediate response to variation of individual wave's nonlinear characteristics.This nature is consistent with the concept of the wave-by-wave approach.Meanwhile, the GRM07 model has been shown to have good accuracy (Gonzalez-Rodriguez and Madsen 2007Madsen , 2010)).Therefore, although we do not have measurements of the bedload transport, we are confident in concluding that the wave-by-wave approach is the more accurate.When using the pure wave-by-wave approach, the current shear stress in Eq. 8 also varies wave-bywave, so an immediate response is assumed for the current.However, as aforementioned, the current actually has slower response time.To check the influence of this, we modify the pure wave-by-wave approach by using a fixed average current shear stress for all the waves, which corresponds to assuming infinite response time for the current.Fig. 7 shows the modified bedload transport rate.The difference is negligible since the wave part dominates the bottom shear.This indicates that the response time effect of the wave-induced current on bedload transport prediction is generally insignificant.
As shown in Fig. 5, we over-predict the skewness by nearly a factor of 2, which may have significant influence on our bedload transport prediction.To quantify this, we artificially half the predicted skewness and re-calculate the bedload transport rate with the pure wave-by-wave approach.Fig. 8 shows the modified result.We see that the onshore transport is significantly reduced by nearly a factor of 2. This indicates that further improvement of skewness prediction is essential.

NUMBER OF WAVES IN WAVE-BY-WAVE SIMULATIONS
The computational time of the wave-by-wave approach is directly proportional to the number of the waves chosen to represent the random waves.With the 247 waves for 0.1-by-0.1 resolution, the total computational time is roughly 2 hours with a personal computer.Thus, it is desirable to have as few waves as possible.The easiest way is to reduce the resolution when discretizing the ξ-τ plane.Since most of the waves in the surf zone can be treated as long waves, the periods are less important than wave heights.Therefore, it is suggested to reduce the resolution for the period more.In Table 1, three cases for which the resolutions are 0.1-by-0.1 (the one initially used), 0.2-by-0.4 and 0.5-by-0.75 are presented.The number of waves can be significantly reduced even down to the order of 10.However, the relative error of the standard deviations of the wave heights and periods, defined by Eq. 6 and Eq. 7, also significantly increase from 0.3% to over 10%, which indicates the quality of the representative incident wave set dramatically decreases.
The upper panel of Fig. 9 shows the comparison of the predicted RMS near-bottom wave orbital velocity heights U brms for the three cases in waves has negligible effect on U brms .Even the 13-wave case only differs from the 247-wave case by a relative error of less than 10%.Similar conclusion can be made on the prediction of other hydrodynamic quantities, such as skewness and asymmetry.The lower panel of Fig. 9 shows the predicted net bedload transport for the two reduced cases and their comparison to the original 247-wave case.Negligible difference is observed between the 41-wave and the 247-wave cases.This indicates that by using 41 waves we have already achieved certain "convergence".Or in terms of the relative error of the two standard deviations in Table 1, an error less than 5% is sufficiently small to assure acceptable results.However, with 13 waves, there are several sharp peaks appearing along the profile, which indicates that the behavior of the wave set becomes akin to that of periodic waves.Although the magnitude of the net transport rate is still reasonably predicted, its cross-shore gradient, due to these peaks, is dramatically changed.Consequently, the rate of the change of the bottom morphology is significantly changed.Since our ultimate goal is to estimate the evolution of beach profiles, the discrepancy in the net bedload transport rate gradient is not negligible.Therefore, the necessary number of waves to reasonably predict the bedload transport is more than that required to predict the hydrodynamics.Otherwise, to utilize the bedload prediction of the wave-by-wave approach with very few (e.g.order of 10) waves, we need to apply some smoothing algorithm to the resulting net bedload transport rate.This issue will be pursued in future work.

CONCLUSION
Two main approaches to account for wave randomness, the spectral approach and the wave-bywave approach, are compared through numerical experiments with the coupling of a surf zone hydrodynamic model and a bedload transport model.The average predictions of hydrodynamics including wave height, near-bottom wave orbital velocity height, set-up and undertow are found to be virtually identical.This indicates that the system response time has negligible effect on most of the average hydrodynamics.However, the predictions of the wave nonlinearities, specified by skewness and asymmetry, are significantly different.The wave-by-wave approach gives much better prediction of asymmetry than the spectral approach, but both approaches dramatically overestimate the skewness by over a factor 2 with the wave-by-wave approach being slightly better.
The net bedload transport predictions given by the two approaches are also significantly different.Since the wave-by-wave approach affords better prediction of wave nonlinearity and can account for the variability of individual waves more realistically, we conclude that it gives more accurate net bedload transport prediction than the spectral approach and recommend using it for bedload calculations.The influence of the current response time is found to be negligible.This verifies the general applicability of the wave-by-wave approach.The effect of the overestimation of skewness is shown to be pronounced.In the future, improvement of the skewness prediction by the TM06 model will be pursued.
The necessary number of waves used in the wave-by-wave approach is briefly discussed.The hydrodynamic prediction is not very sensitive to the number of waves.The number can be reduced to the order of 10 without serious loss of accuracy.However, the net bedload transport prediction, which depends strongly on nonlinear wave characteristics, is found to be much more sensitive to the approach used to simulate the random waves.When the number of waves is reduced to the order of 10, periodicwave-like behavior is observed and the cross-shore gradient of the transport rate is changed to the extent that it cannot be used for the prediction of bathymetric changes unless it is either smoothed or based on a wave set larger than of the order 10.

Figure 1 .
Figure 1.Conceptual nonlinear near-bottom wave orbital velocity profile (onshore direction is positive)

Figure 2 .
Figure 2. Joint non-dimensional wave height and period distribution of Test 3 (the left panel shows the original distribution, the right panel shows the discretization with 0.1-by-0.1 resolution and the selected representative incident wave set)

Figure 3 .
Figure3.RMS wave height and RMS near-bottom wave orbital velocity heights predictions (dash lines: spectral approach; solid lines: wave-by-wave approach; dots: measurements)

Figure 7 .Figure 8 .
Figure 7.The effect of the current response time on net bedload transport prediction (dot line: wave-by-wave approach with fixed current shear stress; solid line: pure wave-by-wave approach)

Figure 9 .
Figure 9. RMS near-bottom wave orbital velocity height and net bedload transport prediction with reduced wave set (dash lines: 13-wave case; dot lines: 41-wave case; solid lines: 247-wave set)

Table 1 . Reducing the number of waves in the wave-by-wave approach
Table.1.It is observed that the reduction of the number of