Development, calibration and validation of a phase-averaged model for cross-shore sediment transport and morphodynamics on a barred beach

Simulating cross-shore sediment transport and associated sandbar migration is still a challenging task for phase-averaged coastal morphological models. Numerical studies have mostly relied on beach morphology prediction for calibration and validation, without examining in much detail the underlying hydrodynamics, sediment concentrations and transport rates. This paper reports on a new three-dimensional coastal morphodynamic model based on the hydrodynamic model of Zheng et al. (2017), combined with an advection-diffusion type suspended sediment transport model and the extended SANTOSS near-bed sediment transport formula of Van der A et al. (2013), to represent the key cross-shore transport mechanisms. The model is are calibrated based on comprehensive measurements from a large-scale laboratory experiment involving regular plunging breaking waves over an evolving sandbar, covering detailed comparisons on hydrodynamics, sediment suspension, transport rates, and bed level evolution. Separate validation using large scale wave flume experiments were also conducted to confirm the model ’ s performance on different conditions. Good agreements are obtained between measurements and model results, which demonstrates the model ’ s ability to reproduce cross-shore sediment transport processes under breaking waves correctly, given that the appropriate parameterizations for intra-wave processes are included. Model results also reveal the onshore near-bed transport is related to wave-induced near-bed streaming, wave skewness and asymmetry, and bed slope effects at different locations across the beach surface. Wave breaking-induced turbulence enhances the near-bed transport within the bed boundary layer which needs to be taken into account in order to achieve good model prediction skills.


Introduction
Sandbars can be found in many natural beaches around the world. They act as coastal barriers causing wave breaking and complex wavecurrent-sediment interactions. As important morphological feature and natural defence mechanism, sandbars have been of special interests to coastal engineers, especially in association with the more frequent extreme storm events and accelerated sea level rise in recent years. However, accurately predicting sandbar evolution and migration across the beach is challenging due to the variety of cross-shore sediment transport mechanisms at play, which may be onshore or offshore directed, and the balance of which determines the overall sediment transport rate direction and magnitude. This physical complexity is reflected in field observations that often see periodical motion of sandbars in the cross-shore direction, e.g. onshore migration under calm weather and offshore migration under storm conditions. Large-scale practical morphodynamic models, such as Delft3D, MIKE, TELEMAC and others, that describe physical processes involved (Van Rijn et al., 2013) based on phase-averaged quantities, have been widely employed in coastal engineering practice. These models work well for beach erosion and offshore directed sandbar migration under large storm events due to the wave breaking induced offshore sediment transport (Thornton et al., 1996;Hoefel and Elgar, 2003). In contrast, the onshore sediment transport which controls the beach recovery and sandbar shore-wards motion, is not well represented, see for example  and Van Rijn et al. (2013). In particular, two primary factors have been highlighted in more recent experimental and numerical studies that potentially contribute to the onshore transport: firstly, non-linear orbital velocity effects due to wave skewness (significant under 2nd Stokes waves in the shoaling zone) and asymmetry (significant under sawtooth-shaped waves in the surf zone), as well as related near-bed boundary layer streaming and Stokes drift; and secondly, breaking waves induced near-bed sediment transport (Christensen et al., 2019). Many of these intra-wave mechanism can be directly resolved in a phase resolving model. The major challenge in the widely used phase-averaged morphodynamic models, however, is to represent these processes through effective parameterisation, which directly affects the overall model accuracy and skills in modelling sandbar migration and the resultant beach profile changes.
Previous efforts focused on representing non-linearity of the nearbed orbital motion and the resultant transport mechanisms, and can be broadly classified into energetics-type transport formula approaches and bed shear stress based transport formulations. In the energetics transport approach, an extra term is often used to represent wave nonlinearity driven transport based on parameterisation of wave asymmetry/skewness, see for example Roelvink and Stive (1989), Elgar et al. (2001) and Hoefel and Elgar (2003), Hsu et al. (2006), Dubarbier et al. (2015) and Ferna'ndez- Mora et al. (2015). In the bed shear stress approach, the wave nonlinearity effect on sediment transport can be included through an additional parameterised transport term similar to the energetics approach, or by direct integration of the time-series of wave orbital velocity. For example, Henderson et al. (2004) used a one-dimensional phase-resolving model to represent the progressive wave-induced flow velocity and sediment transport within the wave bottom boundary layer. To include velocity asymmetry, wave boundary layer streaming and Stokes drift, the time series of the flow velocity measured at the top of the near-bed boundary layer was used to drive the simulations. Hsu et al. (2006) also used a one-dimensional phase resolving boundary layer model in combination with a quasi-steady bed shear stress parameterisation, in order to represent various cross-shore transport processes. Similar to Henderson et al. (2004), field measured velocity at 0.5m above bed was used to drive the onshore sandbar migration in the Duck94 experiment due to the wave velocity skewness and asymmetry. Ruessink et al. (2007) and Ruessink and Kuriyama (2008), on the other hand, used higher order nonlinear wave theory to estimate the near-bed intra-wave oscillatory flow velocities. A more recent study of Rafati et al. (2021) introduced an additional advection term depending on the orbital velocity skewness and asymmetry, to account for the wave nonlinearity driven suspended load transport in XBeach-Surfbeat, whereby the skewness and asymmetry was obtained from the parameterisation of Ruessink et al. (2012). Li et al. (2021) used a one-dimensional phase-resolving model to compute the instantaneous wave induced oscillatory flow and sediment suspension, with the wave asymmetry and skewness also given by Ruessink et al. (2012). To be readily implemented in a phase-averaged morphological model,  proposed the semi-unsteady formula, SANTOSS, which uses the "half-cycle" concept to describe onshore transport related mechanisms, e.g. the phase-averaged total transport rate is computed as summation of a positive transport rate during onshore period and a negative transport rate during offshore period. This approach has the advantage of being simple to implement for any intra-wave variations by compute separately the transport due to wave onshore motion and offshore motion, but with less computational demands than a traditional intra-wave simulation method. For example, Kalra et al. (2019) showed that cross-shore transport rate computed by SANTOSS formula in combination with the model COAWST due to compute the wave non-linearity is comparable with that from a CFD-DEM model. Rafati et al. (2021) also suggested the range of transport mechanism included in SANTOSS can be beneficial for the modelling of onshore sandbar migration. More recent study of Shafiei et al. (2022) demonstrated the usage of SANTOSS formula in simulating cross-shore sediment transport in CROCO morphological model for LIP experiments.
Apart from these wave non-linearity effects and wave-boundary layer interactions, many recent studies have shown that surface wavebreaking induced turbulence, as another important factor, has noticeable effects on not only sediment in suspension, but also on the near-bed (bedload) transport layer, see for example Scott et al. (2009), Ting and Nelson (2011), Yoon and Cox (2012), Sumer et al. (2013), Zhou et al. (2016), Christensen et al. (2019), and Van der Zanden et al. (2017a). However, there are still few methods to represent the breaking wave turbulence effects on the near-bed transport. Hsu and Liu (2004) added the near-bed turbulence kinetic energy under breaking waves to the Shields parameter to represent the enhanced sediment pick-up rate. Van der Zanden et al. (2017c) followed a similar approach as in Hsu and Liu (2004), but included the turbulence kinetic energy into the friction velocity. Rafati et al. (2021) used a depth-averaged advection-diffusion equation to compute the suspended sediment concentration. Similar to Van der Zanden et al. (2017c), the wave-breaking induced turbulence injection effects on the near-bed sediment transport are accounted by adding breaking turbulence term in the near-bed shear velocity, which enhances the sediment entrainment.
It is noted that so far the calibration and validation of these existing model concepts mostly rely on overall sandbar migration and model skills for beach profile changes only, due to the fact that comprehensive laboratory measurements are scarce (Dubarbier et al., 2015). The effectiveness of these parameterisations on wave non-linearity and breaking turbulence have not been systematically tested on all cross-shore processes, especially for the combined effectson the sediment transport in different modes. However, as shown in previous studies, the morphological changes are the result of a combination of many physical processes. Consequently, the modelling studies that are validated solely in terms of sandbar migration may not accurately reproduce all underlying physical processes that drive cross-shore sand transport (e.g. wave-induced orbital velocity skewness and asymmetry, boundary layer streaming, and wave breaking effects). Hsu et al. (2006) demonstrated that discrepancies in the computed transport from different processes may compensate each other and still lead to similar model skill in terms of beach profile changes. Results from Dubarbier et al. (2015) andFerna'ndez-Mora et al. (2015) support these findings and show that the best-fitted model parameters can vary largely from site to site, due to the lack of detailed representation of relevant physical processes. It is therefore important to calibrate model parameters and validate the model's prediction for a full range of physical quantities, including flow velocities, turbulence characteristics, sediment concentrations, and especially different modes of transport, so that the uncertainties in parameterisations of different mechanisms can be identified. More importantly, the detailed calibration and validation can also shed lights on the direct contribution from these processes towards the total sand transport rate at different locations along the beach, which has not been systematically studied so far.
Recently, Van der Zanden et al. (2016) and Van der Zanden et al.
(2017a) conducted extensive measurements of the sand transport processes under breaking waves over a moving bar. These data provides an opportunity to calibrate model parameters for fluid hydrodynamics, sediment transport and corresponding beach evolution simultaneously. Therefore, the present study aims to achieve the following objectives: firstly, implement the SANTOSS formula into the phase-averaged morphological model system of Zheng et al. (2017) to develop a practical tool for simulation of cross-shore sand transport under breaking waves. Secondly, conduct calibration based on Van der Zanden et al. (2016) and Van der Zanden et al. (2017a) measurements for cross-shore sediment transport processes, with particular focus on the wave non-linearity effects and wave-breaking turbulence induced sediment suspension, on-offshore transport rates and overall morphology changes. Results from these comparisons can quantify the effectiveness and validation of the parameterisations in representing cross-shore processes in the model. The model can then be validated against near full scale laboratory experiments on sandy bar migration under erosive and accretive conditions, e.g. the LIP11D data-set (Roelvink and Reniers, 1995), to demonstrate its model skills for the cross-shore transport and beach evolution. Thirdly, perform model tests on various key processes affecting on-/offshore transport rates, focusing on the near-bed region, in order to reveal their contribution to the overall cross-shore transport and the resultant beach profile changes. The remaining of the paper is organised as follows. The model system is summarised in Section 2 and Section 3. The detailed model calibration and validation are presented in Section 4 and Section 5 respectively. As application of the model, tests for different processes and simulation of an evolving sandbar are presented in Section 6. Discussions and conclusions are given in Section 7.

The hydrodynamic model
The three-dimensional nearshore hydrodynamic model of Zheng et al. (2017) is based on the ocean circulation model FVCOM (Chen et al., 2003), coupled to the unstructured-grid version of the third generation spectral wave model Un-SWAN (Booij et al., 1999), to enable the full representation of wave-current interactions. This model system has been validated against theoretical cases, laboratory scale and field experiments of oblique incident waves on a natural, sandy barred beach, illustrating the robustness and efficiency of the model for very different spatial scales and hydrodynamic conditions. For completeness, the hydrodynamic model system is briefly summarised in this section. Further details of the model system can be found in Zheng et al. (2017).
The phase-averaged flow velocities and water levels are resolved from the equations for conservation of momentum and mass. The vortex force concept of McWillams et al. (2004) and Uchiyama et al. (2010) is used to represent the wave-current interactions. The hydrodynamic model equations, including the vortex force formalism and the wave-current interaction (at the right-hand side of the equation), are given by:

∂V ∂t
in which the boldface typesets are used for horizontal vectors, while the vertical components are represented by a normal typeset; (V, w) and (V st , wst) are the Eulerian mean and Stokes velocities, respectively; f is the Coriolis parameter; φ is the dynamic pressure (normalised by the density ρ 0 ); F represents the non-wave non-conservative forces; F w represents the non-conservative wave-induced acceleration forces; (J, K) is the vortex force and K is the lower order Bernoulli head; ρ and ρ 0 are total and reference densities of sea water respectively; g is the acceleration due to gravity; K m and ν are the vertical eddy viscosity and molecular kinematic viscosity, respectively. The generic length scale (GLS) turbulence model of Umlauf and Burchard (2003) is modified to reproduce wave-breaking enhanced turbulence as follows: ∂k ∂t where k is the turbulent kinetic energy (TKE); ψ = (c 0 ∂z represent the turbulence production rates by shear and buoyancy, respectively; ε = (c 0 μ ) 3 k 3/2 l − 1 is the turbulence dissipation rate; σ k andσ ψ are the turbulent Schmidt numbers for k and ψ, respectively; F wall is a wall function; l = (c 0 μ ) 3 k 3/2 ε − 1 is the length scale, K h is the horizontal eddy viscosity, and C 1 , C 2 and C 3 are the coefficients that can be found in Warner et al. (2005). Following Craig and Banner (1994) and Feddersen (2012a,b), the wave breaking induced turbulence injection into the wave column is represented by a flux boundary condition at the surface: where ε b , ε r and ε wcap are energy dissipation rates from the depthinduced wave breaking, wave roller, and white-capping, respectively; α r denotes the fraction of wave dissipation converted into wave rollers, ζ c is the water surface level and D w is an empirical constant. Neumanntype surface boundary conditions for k and ψ are applied following Zheng et al. (2017). The involved surface roughness height z 0s in the surface boundary conditions directly affects the vertical distribution of TKE in the upper portion of the water column as described in Moghimi et al. (2016). For breaking wave conditions in particular, z 0s is connected to the length scale of injected turbulence, which is determined uniquely by the spectral properties of turbulence at the source. It is described by z 0s = H w H in this study, where H w is kept as a non-dimensional tuning parameter, H represents the (significant) wave height.
To represent the regular waves in the present study, the Un-SWAN model is specified with a narrow band wave spectrum with a single frequency and direction as described in Zheng et al. (2017). At the offshore boundary, the mean water level, wave height and period are specified. A non-gradient boundary condition is used for the wave-induced current and TKE. The approach of Soulsby and Davies (1995) is used to calculate the bed shear stresses due to wave-current interaction processes in the near-bed boundary layer, in which the bed roughness accounts for the bedform roughness following Grant and Madsen (1982), Nielsen (1986) and Li and Amos (2001). Further details of the calculation procedure can be found in Zheng et al. (2017).

Sediment transport model
In the present model, the total transport rate is computed as the sum of suspended load, q s , and near-bed transport rate, q nb : q tot = ( q tot,x , q tot,y ) = q s + q nb (6) q s is obtained by integration of the sediment flux from the top of the wave boundary layer (δ wbl ) towards the water surface as shown in Fig. 1.
where h is the still water depth; V, V st and V l are the wave-averaged Eulerian, Stokes and Lagrangian velocities respectively; C is the suspended sediment concentration. The sediment transport within the wave boundary layer is regarded as the near-bed transport rate q nb , which is calculated from the improved SANTOSS near-bed sand transport formula  as detailed in the following section.

Suspended sediment transport module
The suspended sediment concentration is computed through the 3D advection-diffusion equation. To account for non-uniform sediments, the advection-diffusion equation is calculated for several different grain size classes, based on the given sediment grain size distribution. The advection-diffusion equation for grain size class j is given as: where C j is the concentration of the sediment grain size class j; (V l , w l ) = (V, w) + (V st , wst) are the fluid Lagrangian velocities; K h is the vertical sediment diffusivity, which in the present study is equal to the turbulent eddy viscosity from the GLS turbulence model; F C j represents the horizontal diffusion of the sediment concentration; and w j s is the settling velocity (positive upwards) of sediment class j based on the settling velocity formula proposed by Soulsby (1997).
The process of sediment resuspension from the bed is represented by a pick-up function at the bottom boundary: where h is the local water depth and φ j p (θ) is a pickup function as suggested by Van Rijn (1984): where D * s is the non-dimensional grain size, θ cr is the critical Shields number determined by Soulsby (1997). At the water surface boundary, a no-flux boundary condition is applied for the sediment concentration.
Since the near-bed sediment pickup rate can be enhanced by wave breaking induced turbulence (Hsu and Liu, 2004), the method proposed by Van der Zanden et al. (2017c) is also included to represent this effect, e.g. the Shields parameter is given as: where e k is a turning parameter, k nbp represents the total near-bed TKE which contains both the breaking wave induced turbulence and TKE contributions from other sources (e.g. bottom friction, vertical shear velocity).
At this point we reiterate that the suspended load transport is computed by integrating the sediment flux over the region of the water column above the wave boundary layer. The current-and wave-related suspended transport taking place within the wave boundary layer, or near-bed layer, is accounted for through the SANTOSS formula as described in the following section.

SANTOSS formula, wave streaming and wave nonliearity
In the present study, the sediment transport within the wave boundary layer, q nb , is calculated based on SANTOSS formula ( , which accounts for the net wave-induced transport as well as the transport by mean currents within the wave bottom boundary layer. Following Dibajnia and Watanabe (1992), SANTOSS formula uses "half-cycle" concept to compute the phase-averaged total transport rate, in which the transport during the positive "crest" and the negative "trough" half-cycle are distinguished as follows: where s = (ρ s -ρ 0 )/ρ 0 is the relative density in which ρ s is the sand density, g is the gravitational acceleration; θ is the non-dimensional bed shear stress (Shields parameter), with subscripts "c" and "t" referring to the "crest" and "trough" half cycles respectively; T is wave period; T c is the duration of the crest (positive) half cycle and T cu is the duration of accelerating flow within the crest half cycle (Fig. 2); similarly, T t is the duration of the trough (negative) half cycle and T tu is the duration of accelerating flow within the trough half cycle. The variables Ω cc , Ω tc , Ω tt , Ω ct represent the sand load mobilised during the wave crest half cycle ("c") or wave trough half cycle ("t") and then transported during the wave crest half cycle or wave trough half cycle respectively, which are given as follows:  The parameters T c (T ori c ) and T t (T ori t ) represent the positive/crest and negative/trough flow duration defined in this study (in the original definition of ); similarly, T cu and T tu represent the duration of flow acceleration in positive and negative x directions. u δ is the current velocity at the top of wave boundary layer, u r,c (u r,t ) is the representative crest (trough) wave orbital velocity magnitude.
where m st = 11 and n st = 1.2 are two calibration coefficients; θ sf,i is the Shields parameter during the wave crest or trough half-cycle, where the subscript "i" is either "c" for crest or "t" for trough.
In the SANTOSS formula, the wave asymmetry effects are represented through three mechanisms: firstly, a higher/lower bed shear stress during the half cycle with higher/lower wave acceleration (see Eq. (21) in ; secondly, a greater travel distance for Ω tc compared to Ω ct to account for the fact that maximum velocity magnitudes are reached sooner after the 'offshore-onshore' flow reversal during the crest period, compared to the 'onshore-offshore' flow reversal in the trough period (see Eq. (12)); and thirdly, the grain size effects (phase-lag parameters) are affected with increased settling time during the crest half cycle and decreased settling time during the trough half cycle due to wave asymmetry (see the Eqs. (27) and (28) in . The bottom boundary layer streaming, on the other hand, is considered through an additional positive wave Reynolds stress term in the wave propagation direction at the edge of the wave boundary layer (see Eq. (22) in .

Wave breaking effects
The SANTOSS formula of  is originally developed to predict net near-bed sand transport rates under non-breaking waves. In order to apply it to breaking wave conditions in this study, the effect of breaking turbulence on near bed sediment transport needs to be incorporated. Similar to Reniers et al. (2004), the representative half-cycle orbital velocity used in the SANTOSS formula is adapted as following: where u r,i is the representative orbital velocity during the wave crest and trough half-cycle; γ k is a numerical calibration factor; k b is the near-bed wave-breaking induced extra turbulence kinetic energy. This effect only needs to be applied in the region where significant wave breaking occurs. In addition, the half-period concept in SANTOSS allows the added turbulence in Eq. (14) to be applied to the wave crest half-cycle and leads through enhanced shear stresses to an increased mobilised sediment load Ω i in the breaking region.

Bed slope effects
When the sea floor is sloping, the gravity contribution to the net force on the sand grain can increase or decrease the critical shear stress for sand mobilisation. This is particularly important to the transport of sand on the surface of onshore slope and offshore slope of the bar. According to Soulsby (1997), the critical bed shear stress (τ cr ) for the mobilisation of sand is modified according to the local bed slope. The corresponding critical shear parameter θ cr in Eq. (13) is altered at the back of onshore and offshore slope.

Modification for wave-current condition
As shown in Eq. (12) and Fig. 2, the SANTOSS formula calculates the near-bed transport rate based on the duration of each representative half cycle, T c and T t . The corresponding Shields parameters during each onshore and offshore period depend on the representative peak velocities u δ + u r , in which u δ is the mean current velocity at the top of boundary layer and u r is the wave orbital velocity. This approach works well for wave dominated conditions. However, in situations where the magnitude of the current velocity, u δ is close to, or exceeds, the freestream wave orbital velocity maxima u r , the combined velocity at the top of the wave bottom boundary layer remains negative or positive throughout the entire wave period. In these cases, the crest period T c (or trough period T t ) equals (approximately) zero based on the above definition, which means special treatment is required to prevent numerical instability (dividing by zeros) and to achieve a smooth transition between grid cells with unidirectional and bidirectional flows, see Veen (2014). Furthermore, experimental and theoretical evidences support the hypothesis that the near bed sediment transport within the boundary layer is largely driven by the on-offshore wave motion, even though the combined free-stream velocity remains onshore or offshore throughout the wave cycle, e.g. Davies et al. (1988).
In the present study, the representative half cycle periods is therefore defined based only on the wave orbital velocities as shown in Fig. 2, instead of using the combined wave and mean current velocity at the top of the wave boundary layer. Further discussion is presented in Appendix A. Such an approach removes many difficulties in defining various time period, e.g. T cu , T tu and similar terms during the calculation, see Veen (2014), and also reflects the typical near-bed oscillatory flow pattern under such extreme conditions as shown in Van der Zanden et al. (2016) (see their Fig. 6). Results from the present study suggest this approach works well as shown in later sections.

Input parameters
The SANTOSS formula requires several key parameters as input, including the wave boundary layer thickness (δ wbl ), mean current velocity (V δ ) and wave orbital velocity at the top of the wave boundary layer. Following the same approach as in , δ wbl is estimated using the expression as proposed by Sleath (1987). The velocity at z = δ wbl , V δ , is obtained by interpolating between the bed level and the model predicted mean flow velocity at the first model grid cell above the bed, assuming a logarithmic near-bed velocity profile. The analytical expression of Abreu et al. (2010) is used to prescribe the non-linear wave orbital velocity time series, whereby the parametrisation of Ruessink et al. (2012) is incorporated to determine the degree of skewness and asymmetry based on the local wave height, period and water depth (see details in Appendix B).

Morphological evolution module
The morphological evolution module computes the morphological changes by solving the sediment continuity equation, given by: where z bed is the bed level elevation measured from the deepest water depth in the study area (positive upwards); ε 0 is the sediment porosity; and q tot,x and q tot,y are the total volumetric sediment transport rates in x and y directions, respectively. The morphological changes can have significant influence on hydrodynamic processes and sediment transport when they are larger than a few percent of the water depth. Following Warner et al. (2008), the bottom boundary condition of the vertical velocity is equated to the rate of change of the bed elevation in order to account for the effect of morphological changes. The σ-coordinate is used in the present model so that the changes in bed level are reflected in the dynamic water depth and hence the computational mesh in the vertical directions. The interval of morphological updates is taken the same as the hydrodynamic time step.

Model calibration
The model system is calibrated against the comprehensive measurements from the large-scale mobile bed wave flume experiment of Van der Zanden et al. (2016) and Van der Zanden et al. (2017a). In these experiments, measurements of the hydrodynamics and sand transport processes were obtained with a high spatial and temporal resolution, especially around the breaker bar and in the near-bed region. In addition, the bed level evolution was recorded at regular intervals. With these detailed measurements, the performance of every module in the present model system (i.e. hydrodynamic module, wave module, sediment transport module and morphological evolution module) can be explicitly evaluated in detail.
Following a similar procedure as described in Zheng et al. (2017), a series of numerical tests were conducted with varying values of model parameters. In particular, the parameters for the hydrodynamic model, wave skewness and asymmetry, and SANTOSS formula have been tested previously in Zheng et al. (2017), Ruessink et al. (2012) and  and the same values as in these studies were kept in the present study. The calibration therefore focused on model parameters for representing effects of wave breaking and wave nonlinearity to the near bed sediment resuspenson and transport as main objective, including the wave roller parameter α r which controls the fraction of breaking waves that turn into wave rollers before dissipating; turbulent coefficient C 1 that used in the generic length scale equation of the GLS turbulence closure model; parameter D w in Eq. (5) which determines the proportion of wave dissipation energy transferred directly into TKE; the surface roughness height (z 0s = H w H s ) which controls the vertical distribution of TKE in the upper portion of the water column; the wave breaking parameter γ k which determines the breaking turbulence contribution to the near bed representative orbital velocity in Eq. (14); and e k in Eq. (11).

Experiments and model setup
The experiment was conducted in the large-scale CIEM wave flume in Barcelona (Van der Zanden et al., 2017b). Fig. 3 shows the experimental setup and the corresponding 12 measuring positions. The initial bed profile consisted of a bar-trough bed configuration which can be roughly divided into four sections: an offshore slope of the breaker sandbar (x = 35.0 − 54.8 m); a steeper shoreward bar slope (x = 54.8 − 57.5 m); a mildly sloping bed shoreward from the bar trough (x = 57.5 − 68.0 m); and a non-mobile sloping beach (x = 68.0 − 80.0 m). The coordinate system used in the following sections is defined as follows: the horizontal x-coordinate is positive towards the beach with its origin at the toe of the wave paddle in its rest position; the vertical zcoordinate is positive upward with z = 0 at the still water level. The water depth is h 0 = 2.55 m in the horizontal part of the flume, and the wave condition consists of monochromatic waves with wave period T = 4.0nulls and wave height H 0 = 0.85 m at the paddle (Table 1). The mobile bed consisted of medium sand (median diameter D 50 = 0.24 mm) with a measured settling velocity w s = 0.034nullm/s. Starting from the initial bed profile, the experimental duration was 90 minutes, comprising six 15-min runs, during which the bed evolved further. Outer flow velocities (i.e. higher than 10 cm above the bed) were measured by array of Nortek Vectrino acoustic Doppler velocimeters (ADV). In the near-bed region (i.e. below 10 cm above the bed), velocities were measured by a high-resolution Acoustic Concentration and Velocity Profiler (ACVP) (Fromant et al., 2019). Higher above bed, the suspended sediment concentrations were obtained with a six-nozzle Transverse Suction System (TSS). The bed profile is measured at the start of the experiment and after 30, 60 and 90nullmin. The wave broke at 53m, then plunged at 55.5m, and finally splashed around 58.5m. In the following discussion, the region seawards of breaking point is noted as shoaling zone, between breaking point and splash point as outer surf zone and shore-wards of splash point as inner surf zone. More details of this experiment can be found in Van der Zanden et al. (2016, 2017a. The model was setup with a domain covering a 70 m long (x = 10 − 80 m; cross-shore) by 0.5 m wide (alongshore) rectangular area, which is discretised using isosceles right triangles with grid size of 0.1 m in the horizontal direction and 31 uniform vertical sigma layers, resulting in a total of 4326 nodes and 7000 elements. Model tests were carried out with a time step of Δt = 0.01s, starting from initial still water until results were in hydrodynamic equilibrium, and then the results were compared with measurements. The measurements from 30 to 45nullmin was selected for the model calibration, during which the recorded bed evolution was small, and the bed is therefore assumed to be fixed as the measured profile at t = 30nullmin. This approximation enables the calibration to be focused on the time-averaged velocities, sediment concentrations and transport rates.
10 at x < 53 m and C 1 = 1.15 at x > 53 m. b Cross-shore varying value as shown in Fig. 5. c Linearly increases from 0 to 10 along the shoreward slope.
P. Zheng et al. error (RE), and normalised-root-mean-square error (NRMS) are calculated at each cross-shore measurement location, as well as the correlation coefficient (CCF) and model skill (Skill) which give an indication of the overall model accuracy. The absolute error and relative error are given as where M n and C n are the measured and computed (model) results, respectively, at n discrete points. The NRMS indicates the average deviation between the model predictions and the measurements and is defined as The correlation coefficient (CCF) and model skill (Skill) evaluate the coherence between the model results and the observations; a CCF or Skill value of 1 indicates a perfect agreement between model and observations whereas a value of 0 indicates complete disagreement. The CCF is given by where σ C and σ M are the standard deviations of the computations and measurements, respectively; the overbar represents the mean value. The Skill formulation is given as follows: Fig. 4 shows the typical results on the sensitivity tests of phaseaveraged mean flow (undertow) with various values of turbulence model parameter C 1 . Similar results were also obtained for testing of the other parameters. Through these sensitivity tests, the optimal values for the key model parameters are identified and presented in this section.

Optimal parameters
The optimal result found for D w in the present study is approximately 0.3, which is in the range for depth-induced breaking as shown in Govender et al. (2004); Huang et al. (2009);Feddersen and Trowbridge (2005); Feddersen (2012a,b). The value for H w and C 1 , however, were found varying across the beach. A slight increase in C 1 after the breaking point provided the best agreement with the measurements. The final value is therefore taken as C 1 = 1.10 at x < 53m and C 1 = 1.15 at x > 53m. These values are also within the range of variations between 0.9 and 1.4 as suggested in Warner et al. (2005). Similarly, varying H w across the surf zone as shown in Fig. 5 is required to achieve better agreement with the measured data on the wave driven mean current (undertow) and turbulence kinetic energy as shown in later sections. At the offshore side, the H w is found to be 0.2 prior x = 55m, and increases to 0.6 in the outer surf zone till x = 57m, following with a reduction to 0.4 in the inner surf zone. Most of the previous studies proposed 0 < H w < 1.0 (Terray et al., 1999;Burchard, 2001;Umlauf and Burchard, 2003;Moghimi et al., 2013), whist there are still some studies proposed H w > 1.0 (Kantha and Clayson, 2000). The final adopted value in this study complies with the above literature.
The e k in Eq. (11) is found to be 0.015. Van der Zanden et al. (2017c) suggested that e k = 0.3 based on the turbulence level at the reference level for suspended sediment concentration. However, it is also indicated that the value of e k varies greatly with distance from the bed. In the present study, the k nbp at the first model grid above the bed is used,  which is much higher than the turbulence energy at the sediment concentration reference level. The smaller value of k e is therefore expected. Theses spatial variation demonstrated the importance of representing wave breaking effects in the model as key hydrodynamic and sediment resuspension parameters. For the specific mobile bed experiment examined in this study, the wave breaking effect is included along the shoreward slope of the breaker bar with γ k linearly increasing from 0 to 10. Reniers et al. (2004) suggested value of 0.5 for γ k based on the estimated breaking turbulence given by surface roller energy dissipation towards bed surface (Roelvink and Stive, 1989). In the present study, however, the turbulence kinetic energy at the first grid point is used for the k b , which is expected to be much lower than the estimated by the surface roller energy, similar to above e k calibration. Table 1 lists the final adopted values for the model parameters as discussed above. Further details of the calibration procedures can be found in Zheng (2017).

Results
Based on above model parameters, the model results on various processes are compared with the detailed measurements of Van der Zanden et al. (2016) and Van der Zanden et al. (2017a). The errors from these tests on phase-averaged current, sediment concentration, transport rates and bed level evolution are analysed in this section, to reveal the model's effectiveness on representing wave breaking and nonlinearity to the cross-shore transport process. Fig. 6 presents the computed wave height and mean water level, which both show a good agreement with the measurements. In particular, the model accurately reproduces the break point location in this test. Further onshore in the surf zone (i.e. x > 52 m), the wave height also agrees very well with the measurements. The under-predicted wave height over the offshore slope of the bar is partly due to the strong wave reflection from the bar itself, which is not included in the present model. However, discrepancies between the computed and measured values are limited, indicating that the cross-shore behaviour of wave dynamics and set-down/set-up are well captured by the model. Fig. 7 presents the comparison of computed and measured timeaveraged cross-shore Eulerian velocity at twelve profiles where detailed measurements are available. Overall fairly good agreement is achieved between the model. simulations and the measurements. Higher above the bed at almost all profiles, the model accurately reproduces the shape of the undertow profile over the water column, which varies strongly with height above the bed. Close to the bed surface, the predicted values also follow the ACVP measurements (tiny blue circles) reasonably well. The locations of the first model grid point above the bed at all twelve profiles are indicated with the horizontal red dashed lines in Fig. 7. Table 2 lists the normalised root mean square (NRMS) errors of the computed results of various elevations above the bed at these twelve locations to show the overall prediction of velocity profiles across the depth. The absolute error (AE) and relative error (RE) for velocity at the first grid are also listed in the table to show the model's performance on the near bed flows. The good agreements in NRMS indicates that the model captured the wave driven mean flow across the bar very well. As shown in AE and RE, the near bed flow is also predicted with good accuracy for most profiles, except at profile 5 and 6 where large (about 7 cm) and rapid bed level changes occurred in the experiment and hence relatively large velocity discrepancies may thus be the result of the fixed bed assumption in the simulation. Fig. 8 compares the computed and measured (phase-averaged) TKE at the same twelve profiles. Generally good agreement is obtained over the middle and lower part of water column; while in the upper part of water column the TKE is over-predicted, especially near the surface at profiles 3-7 where intense wave breaking occurs. Similar overprediction has been reported in previous studies involving a variety of turbulence closure models, e.g. Christensen (2006), , Brown et al. (2016), Zhou et al. (2016). Different explanations for this issue also have been proposed. Larsen and Fuhrman (2018) showed that the asymptotic exponential growth rates for the TKE and eddy viscosity under surface waves leads to unrealistic turbulence energy production. A stress-limiting approach is proposed to stabilise the turbulence generation and dissipation in the outer surf zone (Larsen and Fuhrman, 2018;Larsen et al., 2020). The present model is based on the same strain tensor and dissipation concept as in other similar turbulence models, hence it is expected that the same stability issues exist. However, as the discrepancies are mostly found near the surface at the break point where their effects on the suspended sediment concentration are small (see section 4.4.4), and the overall agreement is reasonable, the results are thus considered satisfactory for the aim of the present study.

Near bed wave orbital velocity, wave skewness and asymmetry
As stated in previous section, the magnitude of the near-bed wave orbital velocity, wave skewness and asymmetry are the key hydrodynamic parameters for SANTOSS formula. The measured and predicted peak orbital velocities at the top of the boundary layer are presented in Fig. 9(a). It is clear in the figure that the velocity magnitudes are well captured at the shoreward slope (x > 55 m) of the sandbar inside the surf zone. Along the offshore slope (x < 55 m), the predicted trend agrees well with the measurements but the magnitudes are somewhat underestimated, which is partly due to the under-predicted wave height in these positions as previously shown in Fig. 6. The phase-averaged flow (undertow) velocity at the top of wave boundary layer u δ , is interpolated from model results to the height of δ wbl based on a logarithmic vertical velocity profile assumption, and compared with the measured data in Fig. 9(b). The overall agreement is good, except within the outer surf zone from x = 57 − 59 m where the model under-estimates the mean flow. In this region, as highlighted by Van der Zanden et al. (2016), the particularly strong undertow is larger in magnitude compared to the wave-induced oscillatory flow which leads to a mean velocity profile that deviates from the logarithmic profile used in the present model within the near-bed boundary layer.
The computed near-bed velocity skewness Sk(u w ) is compared to the measurements in Fig. 9(c). On the offshore slope before the sandbar crest, the model predictions (approximately 0.6) follow the experiment data reasonably well. After the wave plunge point (x = 55 − 56 m), the measurements vary strongly from the minimum value of 0.03 to the maximum of nearly 1.0 within the trough area, which is not seen in the model result. This is likely due to the complex behaviour of higher-order wave harmonics as the wave breaks and de-shoals along the shoreward slope of the sandbar (Beji and Battjes, 1993). Particularly, waves split into newly reformed surf bores and secondary crests after wave breaking. The parameterisation of Ruessink et al. (2012), which is based on wide range of field conditions, does not capture these very localised wave deformation processes. Fig. 9(d) presents the predicted and measured orbital velocity asymmetry Asy(u w ). The magnitude of the predicted Asy(u w ) is clearly underestimated by approximately 40-45% of the mearused value along the entire sandbar, which is most likely because the present laboratory experiment involved monochromatic waves, while the parameterisation of Ruessink et al. (2012) is based on irregular short waves from field measurements; the wave irregularity may lead to a lower mean wave asymmetry as observed here. Furthermore, the measured data used to calibrate the parameterisation of Ruessink et al. (2012) were time-averaged (over 15nullmin sampling intervals) and bin-averaged (over Ursell number classes), which can result in smoothing of local maximum values. It is interesting to note that Rafati et al. (2021) used the same parameterisation by Ruessink et al. (2012) in XBeach-Surfbeat for simulating cross-shore transport and sandbar migration at the DUCK field site. Their results show better performance of the parametrisation in comparison with the field measurements in energetic conditions, but for low-energetic conditions they obtained similar under-predictions of 30%-45% as in the present study at the bar crest.

Suspended sediment concentration
The computed and observed phase-averaged suspended sediment concentration profiles at the twelve locations are presented in Fig. 10, in which the red and blue circles represent the data obtained by the Transverse Suction System nozzles (TSS) and ACVP, respectively. The errors at each position are listed in Table 3. It is noted that in these results, the near bed sediment pickup is computed based on Shields parameter due to skin friction only. The enhanced pickup due to breaking turbulence is not included. Overall, the model accuracy for the suspension layer higher in the water column is considered to be good across the twelve locations, comparing with the TSS data, especially Fig. 7. Comparison of simulation results (solid black lines) with observed vertical profiles (circles) for the cross-shore wave-averaged velocity; for clarity, the ADV measured data is color coded at different profiles, while the near-bed ACVP measurements are represented by tiny blue circles. The vertical black dashed lines indicate the profile measurement locations and zero values for each profile; the horizontal red dashed lines near the bed indicate the location of the first grid above the bed used in the model.

Table 2
Normalised root mean square (NRMS) error for the cross-shore wave-averaged velocity as shown in Fig. 7, absolute error (AE) and relative error (RE) of the simulated wave-averaged velocities at the first model grid point above the bed.  from the break point (profile 3) shore-wards into the outer surf zone across the bar trough region (profile 8) where the mean relative error (MRE) is less than 30%. The good agreement can be attributed to the well captured phase-averaged current (undertow) and TKE in this region, as shown in Figs. 7 and 8, respectively, indicating that in the suspension layer, mixing due to breaking turbulence is well represented within the model. The over-predicted TKE near the water surface as shown in Fig. 8 has no apparent effect on the accuracy due to the particularly low concentration levels. In the shoaling zone at profile 1 and 2, the concentration is under-estimated, which can be explained partly by the under-prediction in wave orbital velocity and velocity skewness, as shown in Fig. 9, leading to the particularly low near bed reference concentration. Similar under-prediction in the suspension layer (− 0.9 < z/h < − 0.7) can also be seen at profile 10-12 in comparison with the TSS data. This discrepancy can be explained partly as the under-prediction in phase-averaged flow speed as shown in Fig. 9(b). Inside the boundary layer close to the bed surface, the model tends to under-predict the concentration in comparison with the ACVP data (z/h < − 0.9) from profile 3 to 8, indicating the near-bed pick up rate is underestimated in the outer surf zone. In contrast, in the inner surf zone, i.e. profile 10 to 12, the computed near bed concentration follows the ACVP data best, even though the error in the predicted phase-averaged velocity and TKE are similar to those in the outer surf zone. These differences suggest that the wave breaking turbulence can influence the near bed pick up in the outer surf zone, i.e. from the breaking point to the splash point in the present study, more apparently than that in the inner surf zone. In these area, the reference concentration calculated based on Shield parameter due to skin friction, as that in Fig. 10, is clearly not sufficient to provide required resuspension level at the bed surface in comparing with the measured data. The breaking turbulence effects as represented in Eq. (11) should be included in the model prediction.
To confirm these findings, the breaking turbulence is added to the pickup as in Eq. (11) with a calibrated e k in the surf zone (profile 3-12). Fig. 11 presents the corresponding computed suspended concentration with Eq. (11) against the laboratory measurements. The computed sediment concentrations without (broken line, C p2 ) and with (solid line, C k p2 ) inclusion of breaking turbulence are compared with the measured data at the top of the boundary layer in Fig. 11(a). A similar comparison at the bed surface (rst model grid cell above the bed) is shown in Fig. 11 (b). In Fig. 11(c), the computed concentrations are also compared with the measurements across the whole water column, to illustrate the changes in the overall profile. As expected, breaking turbulence represented by Eq. (11) increases the concentration near the bed surface and leads to a better agreement with the ACVP data in the vicinity of the plunge point (x = 56 m), with concentrations rising to approximately 80% of the measured maximum value as shown in Fig. 11(b). However, the measured data clearly indicate that the increase in sediment concentration largely focuses around profile 5-7, and the concentration for x < 55 m and x > 56.5 m largely remain low as shown in Fig. 11(b). These results suggest that the strong plunging jet is able to penetrate through water column and bring breaking turbulence to the boundary layer, but such effects focus only around plunging point, i.e., profiles 5-7 in the present case. At the other positions, the influences from the breaking turbulence is less significant and the skin friction based Shield parameter still works better. With the rise in the near bed concentration at profile 5-7, it is noted that, the inclusion of breaking turbulence also leads to a slight over estimation comparing with the measured concentration on the top of boundary layer at these positions as shown in Fig. 11(a), although the model also overestimates concentration at x < 55m and x > 57m due to the application of Eq. (11) to the whole outer surf zone. These discrepancies suggest a lower level of sediment diffusivity within the near bed boundary layer in the wave flume than the model predictions. Van der Zanden et al. (2017b) also noted this phenomena based on the measured sediment concentration in the flume, and argued that the sediment mixing level drops around the plunge point. A detailed description of breaking effects on both fluid advection and sediment diffusion process is therefore required to fully capture the near bed resuspension and mixing within the boundary layer. Higher above bed at these positions (profile 5-7) as shown in Fig. 11(c), the inclusion of breaking turbulence improves the model prediction on the suspension level, especially from z/h = − 0.8 to − 0.4, indicating the noticeable impacts of breaking turbulence on the overall sediment resuspension and mixing in the present case.

Phase-averaged sediment transport rates
The computed phase-averaged suspended load transport rates along the profile are compared with the measure-ments of Van der Zanden et al. (2017b) in Fig. 12(a).
It is noted that in Van der Zanden et al. (2017b), the measured

Table 3
Mean absolute error (MAE), mean relative error (MRE) and root-mean-square error (RMSE) of the simulated suspended sediment concentration as shown in Fig. 10 suspended load transport rate was decomposed into a 'wave bottom boundary layer' (WBL) component, covering sediment flux from a reference level of 5 mm to 2 cm above the bed, and an 'outer layer' transport component based on the flux above 2 cm. To be consistent with the present model definition, only the 'outer' suspended transport rate is adopted for the suspended load transport rate comparison, since the transport within the WBL is accounted for by the SANTOSS near-bed transport formula. Hence, the model results shown in Fig. 12(a) are obtained by integrating the sediment flux from 2 cm above the bed, i.e. q s,p = ∫ ζ c − h+2cm u l Cdz. The computed suspended load shows a clear offshore directed transport due to the strong wave driven undertow across the beach, with maxima around the plunging point (x = 56m). Overall, the calculated q s agrees well with the measurements in terms of both cross-shore distribution pattern and magnitude. Table 4 shows that the mean relative error (RE) at the 12 profiles are below 41%. Prior breaking, there is certain underestimation at at profile 2 (x = 53 m) on the offshore slope, which corresponds to the location of the under-predicted sediment concentration in suspension shown in Fig. 10. Similar reason can also explain the lower transport rate in the computed values comparing with the measurements from x = 58m shore-wards. However, the model clearly produced the overall correct offshore transport, especially the magnitude and position of the peak immediately after the plunging point due to the good representation of wave breaking induced mean flows, turbulence mixing and sediment suspension as shown in previous figures. Fig. 12(b) compares the predicted near-bed transport rate (solid line) by the SANTOSS formula (q nb,p ) with the measured q nb,m (squares). As discussed previously, the near-bed sediment transport calculated using the SANTOSS formula includes the total sediment transport flux within the wave boundary layer. To be consistent with the model, the measured near-bed sediment transport, q nb,m (Fig. 12), is taken as the sum of the measured suspended load transport from 5 mm above the bed to the top of the wave boundary layer, and the measured bedload transport in Van der Zanden et al. (2017a). Unlike the suspended load, the near-bed transport is generally positive, i.e. onshore-directed at most locations except the bar trough (x ≈ 58 m) area. This indicates that the overall transport within the boundary layer is dominated by wave nonlinearity despite the near-bed mean flows (undertow) are offshore-directed along the entire profile as shown in Fig. 7. The offshore-directed near-bed transport in the bar trough is likely due to the combination of the positive bed slope and the strong offshore-directed undertow velocities in this region (Van der Zanden et al., 2017a).
The computed near-bed transport q nb,p follows the measurements very well, increases in magnitude as the depth decreases along the offshore slope of the sandbar (x = 35.0 − 52.0 m), and reaches first peak near the breaking point at x = 52 m, then decreases towards the plunge point (x = 55 m) at the top of the sandbar. This pattern is attributed to the reduction in magnitude of near-bed orbital velocity and change in wave asymmetry after the initiation of wave breaking (see Fig. 9). However, the substantial under-prediction in orbital velocity asymmetry shown in Fig. 9(d) leads to an underestimation of the onshore transport rate peak at x = 53 m. The measured second peak in q nb,m on the shoreward slope of the bar at x = 56 m, and the negative transport in the bar trough at x = 58 m have all been well captured by the model but offset by about 0.7 m in the onshore direction compared to q nb,p . Van der Zanden et al. (2017a) indicated that the large onshore transport at x = 56 m is likely due to the slope of the bar onshore face, which approaches the natural angle of repose and leads to downward (onshore) transport by gravity. Therefore, the mismatches in the peak position in this region may attribute to the combination of under-estimation in wave asymmetry and gravity driven onshore transport which is absent in the present model (see also Section 6). The smaller peak in the measured transport rate at x = 60 m in the inner surf zone can be attributed to the increasing skewness, which is not captured by the model as shown in Fig. 9. Again, this is due to the discrepancy in the wave skewness as shown in previous section. Table 5 lists the error matrix for the comparison. Overall, a good agreement can be found in terms of the transport direction and its magnitude over the measured section. Apart from those mismatched positions, most of the REs are between 20% and 90%. Without considering the last profile (12) due to its rather small values, the average RE for the total near-bed transport rate is approximately 80%, indicating the present approach within SANTOSS formula for representing wave asymmetry and skewness, breaking and near bed streaming are suitable for the modelling of cross-shore transport process within the surf zone. Further discussions on the contributions from various processes to the near bed transport are presented in the following Section 6.
The sum of the computed suspended load and near-bed transport rates gives the predicted total load transport rate (q tot,p ), which is compared with the measured data (q tot,m ) in Fig. 12(c). Table 6 lists the errors involved at the 12 measurement locations. Overall, the predicted q tot,p follows the measured data well along the sandbar, and the errors at most locations are less than 55%. On the offshore slope, i.e. 48 m < x < 53 m, the predicted total transport is lower than measured, which is attributed to the under-prediction of onshore transport in the near-bed region as shown in Fig. 12(b). The largest discrepancy is found at profile 6 (x = 56m), on the shore-ward slope of the sandbar and the bar trough, i.e. around the plunging point, due to the over-prediction in offshore suspended load transport and under-prediction in onshore nearbed transport.

Fig. 12.
Comparison of computed suspended load transport rate, near bed transport rate and total transport rate with the measured values.

Table 4
Absolute error (AE) and relative error (RE) of the simulated suspended sediment concentration at wave boundary layer level and the suspended load transport rate as shown in Fig. 12 Table 5 Absolute error (AE) and relative error (RE) of the simulated near-bed sediment transport rate as shown in Fig. 12 Table 6 Absolute error (AE) and relative error (RE) of the simulated total sediment transport rate as shown in Fig. 12

Morphological evolution
Using the calibrated parameters shown in Table 1, the model system was subsequently applied to the experimental period from 30nullmin to 90nullmin with a mobile bed condition enabled, to test the model's performance on the overall bed evolution and its feedback to the hydrodynamics across the sandbar. Fig. 13(a) shows a comparison of the measured and predicted bed profile at t = 60nullmin, based on With the bed profile measured at t = 30nullmin. The discrepancy between the model and laboratory data are denoted as the magenta lines on the top of the figure. The measurements show a region of slight erosion in the middle part of the offshore slope (x = 46.0 − 51.0 m), significant accumulation at the upper end of the offshore and shoreward slope of the sandbar (x = 51.5 − 56.8 m) and significant erosion around the trough (x = 57.0 − 60.0 m), leading to a steepening of both the offshore and shore-ward slopes of the sandbar. The predicted bed profile agrees well with these measurements, and all the aforementioned erosion and deposition characteristics have been correctly captured. The simulation continued for further 30nullmin and the resultant bed profile is compared with the measured profile at t = 90nullmin in Fig. 13(b). From t = 60nullmin-90nullmin, both the measured and modelled profiles show a similar morphological evolution as shown in Fig. 13(a), i.e. the sandbar crest continues to grow while the bar trough continues to deepen, leading to a steeper shore-ward slope. Over both 30-60nullmin and 60-90nullmin periods, the model captured correctly the erosion on the offshore slope, deepening of the trough area, and the deposition over the crest of the bar. Comparing with the measurements, the magnitude of the error is within ±5 cm, with a root-mean-square error (RMSE) of 1.8 cm (<12% of the total bed level change), and correlation coefficient (CCF) and model skill of 0.93 and 0.92, respectively. These results demonstrate the model's effectiveness in simulating both onshore and offshore transport processes with good accuracy, and hence its capability in reproducing the trend of bar migration.

Model validation
Based on the calibrated parameters, the model system was applied to LIP11D experiments of Roelvink and Reniers (1995) as independent validations. In particular, the LIP1B and LIP1C cases were selected to test the model's performance for both erosive (LIP1B) and accretive (LIP1C) conditions, respectively.

Model setup
The numerical model covered 201m long beach, and the water depth was 4.1m. Cross-shore profile of the model domain can be seen in Fig. 14 for the main test section in both cases. Along the offshore regions, a significant wave height of 1.4m with a 5s wave peak period in LIP1B and 0.6m wave height and 8s peak period in LIP1C was imposed. The corresponding random wave spectrum were also applied as in the experiment with the total 10 frequency bins. The Dean's number is 11.2 for LIP1B and 3 for LIP1C, respectively, clearly distinguishing the erosive and accretive conditions in the two very different cases. The median grain diameter D 50 = 0.22mm, sand density ρ s = 2650nullkg/m 3 , water    c Linearly increases from 0 to 10 along the shoreward slope.
More detailed information regarding model setup are given in Table 7. The computational domain was discretised into uniform triangular meshes with horizontal resolutions of 0.5 m. In the vertical direction, total of 31 σ-layers were used to resolve the dynamic processes in the water column. The cyclic condition was applied to both lateral boundaries to eliminate the variations in the alongshore direction. The model was run for 1 h simulation period. Two morphodynamic acceleration factors, 18 and 13, were used to simulate the morphological evolution in the laboratory experiment of LIP1B and LIP1C tests, respectively. The values for the calibrated parameters were kept the same in these two validation cases, apart from the C 1 reduced slightly to 0.85 in the deep water region and H w reduced to the range between 0.25 and 0.4 to reflect the very different spilling breaking waves in these validation tests, comparing with the plunging breaking waves in the calibration.

Results
The computed results in wave propagation, fluid hydrodynamics, sediment transport and beach evolution were compared with the available measurements for both cases. Results from LIP1C are shown here for detailed discussion. Similar results were also obtained for LIP1B test. Fig. 15(a) shows the computed wave height and water level from shoaling zone (x < 130m) to the surf zone (130m < x) for the LIP1C test. Inside the surf zone, wave height decreases from around 0.6m to about 0.4m as the waves break and dissipate energy. The vertical profiles of the corresponding undertow currents at the 10 measuring positions are also presented in Fig. 15 compared with the measured data. Unlike that in the calibration test, the speed of the undertow current across the beach is fairly small, with a maximum velocity of about 0.15nullm/s in the inner surf zone (x = 152m), as shown in Fig. 15 which is typically found in the acrretive condition. Similar results were also achieved in LIP1B test, apart from strong undertow currents with a maximum value of approximately 0.6nullm/s over the bar, as under the erosive waves. The corresponding errors in the computed velocity at these positions are listed in Table 8 for both cases. The overall agreements are found very good in both tests, with the maximum NRMS less than 0.25 for LIP1B and LIP1C in most points, apart from the points close to the shoreline due to the relatively small velocity magnitude.

Sediment transport and morphological changes
The computed vertical profiles of sediment concentration in LIP1C are also shown in Fig. 16 at the 9 positions across the beach, compared with the laboratory data. In the shoaling region and outer surf zone (x < 130m), sediment concentration near the bed tends to be low, 0.5nullg/l for LIP1B and 0.2nullg/l in LIP1C. Higher values were found in the surf zone (x > 130m)with the maximum value of 0.5nullg/l and 1.5nullg/l for LIP1B and LIP1C, respectively. Error analysis for all nine measurement locations is presented in Table 9. Model results showed fair agreement with the experimental data. The mean relative error is less than 30% in LIP1B and 25% in LIP1C at most points, similar to that in the previous calibration tests. The higher divergence in LIP1B at the shoaling region is largely due to the very low concentration values at these positions. This clearly indicates that the inclusion of breaking turbulence in the model's mixing and sediment pick-up approach can improve the model's prediction accuracy and can represent the breaking wave-induced hydrodynamics and sediment resuspension reasonably well.
The computed near-bed, suspended, and total load transport rates are presented in Fig. 17 for LIP1B and Fig. 18 for LIP1C respectively. The resultant morphological evolution is shown in panel (b) in both figures, and the difference between the initial and final beach profiles from the model and measurements are also compared in panel (c). The measured total load for both cases were also shown as a comparison. In the LIP1B case, the suspended load transport shows a clearly offshore directed transport with the peak inside the surf zone at around 145 m. Its magnitude also is comparable with the peak of the near-bed transport.

Table 8
Normalised root mean square (NRMS) error for the vertical profile of cross-shore velocity as shown in Fig. 15, absolute error (AE) and relative error (RE) of the crossshore velocity at the first node above bed level for LIP1B and LIP1C.

Table 9
Mean absolute error (MAE) and mean relative error (MRE) and normalised root-mean-square (NRMS) error of the simulated suspended sediment concentration for LIP1B and LIP1C as shown in Fig. 16 The near-bed transport, as predicted by SANTOSS, is dominantly onshore directed, and the maximum is found before the breaking, x = 120m. The total transport rate obtained from the model shows a very similar change across the beach compared with the experimental data, e. g. an onshore transport prior to breaking due to the strong near-bed transport and then a change to offshore direction due to the increasing offshore directed suspended load transport. The resultant sandbar in the model result and measurement depicts a clear offshore movement due to the combination of these cross-shore transport processes. The predicted migration is approximately 10 m offshore after 18 hours of simulation time (see Fig. 17(b)), which corresponds very well with the measured values. The bed change during this period of time also agreed very well, as shown in Fig. 17(C), especially the drop of bed level between 130m < x < 140m and increase to the maximum at x = 145m and subsequent drop 150m < x < 160m. The error of the bed level change at the bar positions x = 135m and x = 145m, e.g. where the maximum erosion and accretion, are less than 10% of the measured values. The CCF and Skill values for the model prediction as shown in Table 10 are all similar to that in the validation case. As a comparison, in LIP1C case in Fig. 18(a), it is noted that the computed onshore near bed transport dominates the overall sediment motion across the beach. The offshore-directed suspended load is fairly small in its magnitude in comparison with the near bed transport. The measured total transport rate agrees with the model results in its direction, largely remains onshore directed, apart from that in the surf zone (x = 140m) where the transport reversed towards offshore. However, the measured total transport rate is stronger in its magnitude in the shoaling region up to 140m comparing with the computed values. These differences suggest that the near bed transport rate was under-predicted by the model in the shoaling region, while inside the surf zone (x > 140m), it is over-predicted. However, the good overall pattern of the computed total load transport means the bed morphology changes are still well predicted as shown in Fig. 18(b). The corresponding bed level changes over the 13 hours also are reproduced with very good accuracy as seen in Fig. 18(c), such as the slight accretion prior x = 130m and then switches to erosion up to x = 160m. Similar in LIP1B, the error of the bed level change at the bar positions x = 137m and x = 140m are within than 25% of the measured values, and the CCF and Skill as shown in Table 10 clearly indicates the good performance of the model. The differences between Fig. 17 for LIP1B and Fig. 18 for LIP1C clearly demonstrate the mechanism of bar migration under erosive and accretive conditions respectively. The offshore suspended transport process erodes the bar inside the surf zone and move it offshore. The onshore near bed transport, on the other hand, tends to accrete sediment in the shoaling region and push the bar onshore. The relative strength between these two transport processes determines the overall bar migration and beach morphology. From model development point of view, it is therefore critical to capture both onshore near bed transport and the offshore suspended transport due to various mechanisms involved. In the present study, the inclusion of SANTOSS formula for the wave nonlinearity effects on near bed transport, breaking induced sediment resuspension clearly improved the model's performance in both suspended and near bed transport predictions and hence the good agreements with the measurements.
One of the major difference in the predicted total transport rate in LIP1C case as shown in Fig. 18(a), is that the peak in the onshore direction occurs at around 135m in the measurements, while the corresponding predicted peak is shifted around 5m onshore, at around 140m. This is not the case in LIP1B in which both maxima of onshore and offshore transport were well captured by the model as shown in Fig. 17. Part of the reason for such difference in the two cases, is the model predicted wave asymmetry as presented in Fig. 19, in which the computed wave skewness and asymmetry for LIP1B and LIP1C were compared with the experiments. The model skewness and asymmetry are derived from wave orbital velocity calculated using a combination of Abreu et al. (2010) and Ruessink et al. (2012) formulas. The computed results for LIP1B showed fair agreement in the shoaling region at < x 130m. Inside the surf zone, the model rightly predicted the drop of wave asymmetry after the breaking region (130m) although overall the value is over-estimated before gradually falling close to the beach. In LIP1C simulation, however, the model predicted a nearly constant rate of drop in wave asymmetry, and failed to reproduce the large decrease from − 0.6 to − 1.0 after x = 135 − 140m as shown in the measurements. This is similar to the previous validation case where the model predicted much lower value of wave asymmetry. As discussed in the following section, the under-estimation in wave asymmetry magnitude noticeably contributes to the error in the near bed transport inside the surf zone.

Model applications
After it has been calibrated and validated extensively, the model is used to investigate the effects of wave breaking and nonlinearity on the near bed transport. In particular, the effects of four main mechanisms on the near-bed transport rate are explored: (1) bottom boundary layer streaming; (2) wave asymmetry effects; (3) bed slope effects; and (4) wave breaking effects. The experiment setup of Van der Zanden et al. (2017a) is used for the model application scenario because of the comprehensive measurements available in the data-set. A series of systematic numerical tests were conducted by individually switching off the above mechanisms within the SANTOSS formula assuming the bed is fixed, with the aim to identify the effects of the corresponding physical process on the overall near bed sediment transport rate. The predicted near-bed transport rates with one of the four mechanisms excluded are shown in Fig. 20(a)-(d), respectively. For comparison, the computed near-bed transport rate (see Fig. 12(b)) from the model simulation with the final calibrated parameters as shown in Table 1 and with all effects switched on is also included in these figures (q nb,p , the red line).
Without the near-bed streaming, the predicted near-bed q nb,ps1 in Fig. 20(a) becomes offshore-directed along the entire profile with the exception of a small area on the shore-ward slope of the bar. The magnitude of the transport rate is also much smaller than q nb,p , especially along the offshore side of the sandbar. These changes can be explained by the fact that the near-bed mean flow is offshore directed and of relatively high magnitude as shown in Fig. 9(b). The absence of wave boundary layer streaming will therefore lead to offshore-directed transport of the medium sands by the mean flow. The small section of onshore transport largely attributes to the gravitational effects along the  shore-ward slope of the bar. Rafati et al. (2021) found that the wave boundary layer streaming in the SANTOSS formula constitutes 30%-50% of the total wave driven onshore transport when tested with the XB-SB model for different erosive conditions. These results are consistent with the findings of Henderson et al. (2004) and Yu et al. (2010), who suggest that the boundary layer streaming is an important contributor to the near-bed and total transport rates, especially in the shoaling region prior to wave breaking. Fig. 20(b) shows the computed near-bed transport rate by the SAN-TOSS formula without the effect of wave asymmetry (q nb,ps2 ). Compared to the full processes driven results, q nb,p , the magnitude of q nb,ps2 drops slightly between x = 52.0 and x = 57.0 m, corresponding to the moderate increase of computed Asy(u w ) values in this region (see Fig. 9(d)). Further onshore, there is no notable difference between q nb,ps2 and q nb,p , which is explained by the very low Asy(u w ) in the computed hydrodynamic input to the SANTOSS formula. However, as discussed previously, Asy(u w ) is significantly under-predicted in comparison with the measurements. The wave asymmetry effect is therefore expected to be underestimated considerably in this comparison. In order to evaluate the wave asymmetry effects reliably, an additional simulation was conducted by adopting the measured Asy(u w ) as hydrodynamic input to the SANTOSS formula. The corresponding result q * nb,q is shown as the black dashed line in Fig. 20(b), together with the measured near-bed transport rate as symbols. Clearly q * nb,q depicts a much better agreement with the measurement than q nb,p along the upper end of the offshore-facing slope (x = 50.0 − 55.0 m), the onshore transport rate is enhanced to a magnitude that is closer to the measured value, e.g. almost 100% increase in the computed values; from the sandbar trough further onshore (x = 57.0 − 63.0 m), the transport direction is now correctly predicted as onshore, even though its magnitude is smaller than measured. These comparisons indicate that the wave asymmetry parameter Asy(u w ) in the SANTOSS formula strongly influences the near-bed transport rate by increasing the onshore transport in the surf zones where waves deform significantly, which is consistent with many previous studies. The results in previous validation case LIP1C also clearly show the under-prediction in wave asymmetry, which correlates with the under-estimation of near bed onshore transport and the shift in the position of the transport maxima around breaking region. The wave asymmetry effects are considered in SANTOSS through changes in the near-bed shear force and through modifications to the on-/offshore half cycle sediment settling times to account for the wave shape. The present results suggest that with better parameterisation on the wave asymmetry process, SANTOSS formula is therefore expected to represent wave nonlinearity effects on the near bed transport inside the surf zone with good confidence.
The computed near-bed transport rate without the bed slope effect is presented in Fig. 20(c). As expected, the magnitude of up-slope transport increases without the opposing effect of gravity while down-slope transport is reduced, although its overall effect on the net transport rate is not as significant as the previous two mechanisms. In addition, it is noted that the local peak at x = 55.3 m disappears after the bed slope effect is eliminated as discussed previously. Fig. 20(d) depicts the result without the effect of wave-breaking induced turbulence in the SANTOSS formula, i.e. γ k = 0 in Eq. (14). The magnitude of the predicted near-bed transport rate (q nb,p ) reduces along the upper part of the shore-ward slope of the bar (x = 55.0 − 56.1 m). In addition, the location where q nb,p changes direction (from onshore to offshore) shifts by about 0.5 m in the offshore direction and the peak offshore transport rate increases by almost 100%. These differences are caused by the absence of extra stirring effects due to wave-breaking induced turbulence in the near-bed region in the model, Eq. (14), which then leads to a reduction in bed shear stress and near-bed onshore transport. The results are consistent with the observations in Van der Zanden et al. (2017a) and Van der Zanden et al. (2017b) in which the near-bed onshore transport is correlated positively to the enhanced turbulence under the plunge point.

Discussion and conclusions
The present study focused on calibration and validation of a phase averaged morphological model based on migrating sandbar under breaking waves in controlled laboratory experiments. The comprehensive laboratory measurements of Van der Zanden et al. (2017a) provide opportunity to calibrate key parameters within all components of the morphodynamic model system simultaneously, from wave-induced mean flows, wave-averaged turbulence characteristics, sediment concentration to the suspended load and near-bed transport rates, and overall beach morphology. Results from the present model provide the direct confirmation on the parameterisations for the cross-shore sediment transport processes, as well as their contributions to the onshore sandbar migrations.
The results on transport rates highlight the importance of including the wave skewness, wave asymmetry, wave boundary layer streaming, gravitation driven transport and breaking turbulence, which vary in their contribution to the total transport at different regions along the sandbar. In particular, the wave driven near-bed streaming is found to contribute considerably to the onshore sediment flux in the shoaling zone and surf zone prior the plunging point in the present study. By including the measured wave asymmetry, the present model is able to reproduce the measured near-bed transport rate with excellent agreement (Fig. 20), suggests that from the breaking point shore-wards wave asymmetry dominates the onshore transport within the bed boundary layer. Around the plunging point, the wave asymmetry reduces to a much lower level, which minimises the near-bed onshore transport rate to near zero. At the onshore slope, the strong gravitational flux causes another peak in onshore transport. In the inner surf zone, the under-estimated wave skewness seems to cause an under-prediction in the onshore transport. The observed significance of wave skewness, asymmetry and wave boundary layer streaming is broadly consistent with previous studies by Hoefel and Elgar (2003), Hsu et al. (2006), Ruessink et al. (2007), Dubarbier et al. (2015) and Ferna'ndez- Mora et al. (2015). In the present study, the inclusion of wave skewness, asymmetry and near bed streaming are based on a parameterised shear force together with alterations in the period of T cu or T tu in the SANTOSS formula. The good agreements on near-bed transport rates indicates that such an approach works reasonably well in the present simulations. The local "half cycle" method in the SANTOSS formula is essential to incorporate the required intra-wave processes at the appropriate wave phase. However, it is also clear that the accuracy of the parameterisation of each process is often fundamental to the model accuracy, as demonstrated most clearly by the analysis of wave skewness and asymmetry effects shown in Fig. 20.
Wave breaking effects are included in the present model through an enhanced turbulence mixing due to the plunging jet at the water surface, by including breaking-induced turbulence on the pickup in Eq. (11) and near-bed onshore transport in Eq. (14) used by the SANTOSS formula. As shown in Fig. 20, its contribution to the near-bed transport is mostly apparent in the region between plunging point and splash point, which confirms many previous studies, e.g., Hsu et al. (2006); Lim et al. (2020). The breaking turbulence impacts on the sediment pickup and subsequent suspended load transport are revealed by using Eq. (11). Although the results are broadly similar to previous research (Zhou et al., 2016;Fernandez-Mora et al., 2017), the detailed comparison in Fig. 11, however, indicates the region at the bed surface affected is restricted to the immediate neighbouring of the plunging point. Meanwhile, better Fig. 20. The cross-shore distribution of the near-bed transport rate calculated from various sensitivity numerical tests: q nb,p (red line) is the same predicted result as shown in Fig. 12(b); q nb,ps1 , q nb,ps2 , q nb,ps3 and q nb,ps4 (black line) are the SANTOSS formula predicted results without (a) the effect of near-bed streaming, (b) the effects of wave velocity asymmetry, (c) bed slope effects and (d) wave-breaking induced turbulence. The model setup used for q * nb,q is same with that used for q nb,p , except that the measured Asy(u w ) is used as hydrodynamic input to the SANTOSS formula. Dot-dashed line in each figure represents the measured bed profile at t = 30nullmin. approach is required for sediment diffusivity tends to avoid the over-estimation in concentration at the top of the wave boundary layer. However, as most phase-averaged morphological models, the present model does not resolve the detailed near-bed turbulence structure. Consequently, the estimation on the transformation of breaking turbulence from surface to the bed surface becomes important to the good prediction of sediment re-suspension under the plunging breaking waves. From a modelling point of view, these results indicate the necessity to include wave breaking induced turbulence in the sediment pickup function to achieve a proper level of near bed resuspension, but at the same time, emphasise the correct modelling of turbulence characteristics and vertical diffusion.
Similar to what Hsu et al. (2006) argued, the present study also suggests that due to the presence of various competing transport mechanisms, a good overall model performance on morphology does not necessarily indicate a correct prediction of all underlying transport processes. For example, discrepancies in predicted wave skewness and asymmetry (see Fig. 9) clearly result in underestimated near-bed onshore transport rates (see Fig. 20), but because at the same time the offshore directed suspended load transport is under-predicted in this region (see Fig. 12), the total transport rate and the overall bed level evolution are still well predicted. This result demonstrates the complexity of cross-shore sediment transport which depends on all the processes involved. Accurate representation of each process is therefore required to ensure robustness of the model performance for a wide range of conditions. This study has highlighted the importance of inclusion and accurately representing all relevant physical mechanisms for the cross-shore sediment transport in a phase-averaged morphodynamic model. The main conclusions are summarised as follows.
1. The improved SANTOSS formula is capable of capturing the onshore transport rates in the near-bed region. The good agreement with the measured bed evolution also indicates that the present model describes the feedback between the bed level changes and hydrodynamics well, which further assures the accuracy and effectiveness of the presently adopted modelling approaches and their applicability in practical morphodynamic simulations. 2. Model results suggest the wave asymmetry has a strong influence on the near-bed onshore sediment transport, especially from breaking point shore-wards, due to the strong deformation of the waves shape.
The SANTOSS model provides good predictions of the near bed transport rate when the measured wave asymmetry is used, which indicates that accurate parameterisations of the key intra-wave processes will improve the model's performance inside the surf zone considerably. 3. Model results suggest that the near bed turbulence energy can be used to estimate the sediment pickup rate in the adjacent of the plunging point, which is substantially enhanced due to the breaking effects. However, such an approach requires detailed turbulence modelling within the wave bottom boundary layer in order to provide appropriate levels of sediment resuspension. 4. Model tests on various processes show that in the shoaling region up to the break point, the predicted near-bed transport rate is highly sensitive to the wave induced near-bed streaming. Around the bar crest, where waves have a highly nonlinear shape, results are shown to be sensitive to the wave asymmetry of the near-bed flow. Bed slope effects play an important role at both the offshore and onshore slope of the bar, but have a minor effect in the inner surf zone. 5. Sediment transport in the cross-shore direction is often a combined process involving the key transport mechanisms listed above. These mechanisms can be compensatory under certain conditions, hence it is important to validate the model skill on all key physical quantities in order to minimise model uncertainties.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.
As discussed in Section 3.2.4, the representative half cycle periods is defined based only on the wave orbital velocities in the present study. Consequently, in the presence of a mean offshore-directed undertow (u δ in Fig. 2), the onshore transport period T c is increased and the offshore transport period T t is reduced in the present method (i.e. T c and T t ) compared to the original definition of  (i.e. T ori c and T ori t ), although the magnitude and direction of transport flux are still dictated by the combined flows (u δ + u r,c , u δ − u r,t ). When mean offshore flow equals (or exceeds) the magnitude of positive orbital velocity as depicted in Figure Appendix A.1, the present calculation is still based on both the representative boundary layer flow velocity u δ + u r,c , crest period T c and the representative velocity u δ − u r,t , trough period T t (i.e., the blue lines), in contrast to the original approach that uses only the offshore peak representative boundary layer flow velocity u δ − u r,t and the wave period T (i.e., the red lines).

Appendix B. Wave skewness and asymmetry
Appendix B.1. Analytical expression of wave orbital velocity The analytical formulation of Abreu et al. (2010) that reproduces the near-bed wave orbital velocity is as follows: 1 − r cos(ωt + φ) (B.1) u w orb is the orbital velocity amplitude, w( = 2π/T p ) is the angular angular frequency, φ is the waveform parameter (− π ≤ φ ≤ 0) and r is the parameter of skewness (− 1 < r ≤ 1). The corresponding acceleration time series is given by:

Appendix B.2. Parameterizations for the parameter of waveform (φ) and skewness (r)
In order to utilize the analytical expression of Abreu et al. (2010) to estimate the near-bed wave orbital velocities, Ruessink et al. (2012) proposed a parameterisation method for calculating the parameter of waveform (φ) and skewness (r). This parameterisation method could be separated into three steps: Firstly, the Ursell number is calculated as a function of wave height (H), wave period (T p ) and water depth (h) with the following equation: in which p 1 to p 6 are six calibration parameters with advised values of p 1 = 0.0, p 2 = 0.857, p 3 = − 0.471, p 4 = 0.297, p 5 = 0.815 and p 6 = 0.672. In addition, it is noted that the two commonly used parameters to measure the skewness and asymmetry of the wave motion, the wave orbital velocity skewness (Sk(u w )) and velocity asymmetry (Asy(u w )), which are defined respectively as Sk(u w ) = (uw) 3 where σ uw is the standard deviation of orbital velocity u w , and H(u w ) is the Hilbert transform of u w .
In the third step, the parameter of skewness (r) is determined via its relation to the parameter B as and the parameter of waveform (φ) is determined via its relation to the parameter φ as φ = − φ − π 2 (B.8) Via the Eqs. (B.3) -(B.8), the parameter of waveform (φ) and skewness (r) used in the analytical expression of Abreu et al. (2010) can be determined from the wave parameters H, T p and the water depth (h). However, it is also noted that it is difficult to obtain the parameter r from a known parameter B from the Eq. (B.7). Veen (2014) thus proposed a simple fitting curve as r = 0.0517B 3 − 0.4095B 2 + 1.0853B − 0.0099 (B.9)