# Adaptation of the SANTOSS transport formula for 3D nearshore models: Application to cross-shore sandbar migration

Hassan Shafiei, Julien Chauchat, Cyrille Bonamy, Patrick Marchesiello

► To cite this version:

Hassan Shafiei, Julien Chauchat, Cyrille Bonamy, Patrick Marchesiello. Adaptation of the SANTOSS transport formula for 3D nearshore models: Application to cross-shore sandbar migration. Ocean Modelling, 2023, 181, pp.102138. 10.1016/j.ocemod.2022.102138. hal-04257166

HAL Id: hal-04257166

https://hal.science/hal-04257166v1

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Adaptation of the SANTOSS transport formula for 3D nearshore models: application to cross-shore sandbar migration

Hassan Shafieiª, Julien Chauchatª,*, Cyrille Bonamyª, Patrick Marchesielloᵇ

## Abstract

The ability to model bar migration in the nearshore zone is critically important for predicting coastal morphology. We present a cost-effective wave-averaged 3D model (CROCO) with parameterization of wave-related bedload transport (adapted SANTOSS), and evaluate/calibrate its performance in comparison with data measured in a two-dimensional wave channel. In this model, suspended-load transport (both in and above the wave bottom boundary layer) is resolved by the flow model. We assume that if the tuning parameters are not robust to varying forcing conditions, the parameterization can be considered a failure. The present work is therefore a first step towards a numerical model capable of predicting the onshore and offshore migration of a sandbar under storm and post-storm conditions using the same set of parameters. One way to achieve this goal is to treat the effects of waves and currents in separate formulations to avoid conflicts and redundancy. We show that the implementation and adaptation of SANTOSS in CROCO can achieve this goal, provided that the hydrodynamics are accurately reproduced and the separation between the wave-related bedload and current-related suspended load formulations is respected.

**Keywords:** Beach morphology, Surf zone, Bed load, Sediment load, sandbar migration, beach erosion/accretion

<footnote>ªUniv. Grenoble Alpes, CNRS, Grenoble INP*, LEGI, 38000 Grenoble, France</footnote>

<footnote>* Institute of Engineering Univ. Grenoble Alpes</footnote>

<footnote>ᵇIRD/LEGOS, Toulouse, France</footnote>

<footnote>*Correspondence to: Univ. Grenoble Alpes, LEGI, 1209 Rue de la Piscine, 38610 Gires. E-mail addresses: julien.chauchat@univ-grenoble-alpes.fr (J. Chauchat).</footnote>
# 1. Introduction

Sediment transport due to waves and currents plays a crucial role in the morphology and ecosystem of coastal areas. Waves can change coastal morphology in any direction depending on hydrodynamic and morphological conditions. Sandy beaches are among the areas most rapidly affected by various environmental conditions. They are continuously subjected to waves that result in the creation of sandbars near the shoreline. Therefore, sandbars are the primary indicator of cross-shore sediment transport in most cases and can migrate onshore or offshore. Their contribution to wave breaking is considerable, and they can therefore serve as a beach protection mechanism.

Many studies have explored the drivers of the motion of coastal sediment particles (e.g. van Rijn et al., 2013). The processes responsible for shoreward transport are wave-related (wave streaming, wave skewness and asymmetry, and Stokes drift), while the undertow is the main driver of offshore motion (Henderson et al., 2004). Wave streaming is the near-bed onshore mean Eulerian flow caused by waves propagating towards the beach (Longuet-Higgins, 1953); the Stokes drift is the difference between the Lagrangian and Eulerian flow velocities (Stokes, 1880; van Den Bremer and Breivik, 2018); and a nonlinear wave form is described by skewness (sharp high peaks and broad shallow troughs) and asymmetry (steep front faces and mild-sloped back faces) (Elgar, 1987; Henderson et al., 2004). The undertow is the wave-average seaward current that moves under the trough of the waves. It is the result of the balance between the radiation stress gradient due to breaking waves, the turbulent shear stress and the cross-shore pressure gradient due to wave set-up (Johnson, 1919; Dyhr-Nielsen and Sørensen, 1970; Hansen and Svendsen, 1985). The final balance of these wave and current forcing determines the net direction of coastal particle transport, leading to sandbar migration in many of the cases (e.g. Elgar et al., 2001; Hoefel and Elgar, 2003; Hsu et al., 2006).

One of the key phenomena influencing the role of each process mentioned above is the position of wave breaking, in particular whether the waves break onshore or offshore of the sandbar (Symonds and Bowen, 1984; Roelvink and Stive, 1989; Henderson et al., 2004). In particular, this affects the cross-shore distribution of breaker-induced undertow, which is considered in many studies as a key player on bar formation (e.g. Dyhr-Nielsen and Sørensen, 1970; Dally and Dean, 1984; Stive and Battjes, 1985; Stive, 1987). Breakers can also increase the onshore bedload transport through their effect on wave asymmetry and skewness and near-bed turbulence. However, the magnitude of this onshore forcing relative to the offshore forcing by the undertow depends on other parameters such as the breaker type (e.g. Ting and Kirby, 1994).

The shear stress caused by near-bottom flows triggers sediment mobilization, while currents and nonlinear waves move the near-bed sediments. However, there is no exact boundary between bed- and suspended-load and the bedload boundary layer is rather an arbitrary mathematical definition. The resultant balance of these transport modes determines the direction (i.e. onshore or offshore) of the bar migration. With the increase of the bed shear stress, the particles start to roll and slide (Fredsoe and
Deigaard, 1992). If the shear velocity surpasses a threshold, the particles will detach from the bed and, depending on their horizontal velocities (e.g., due to turbulent eddies) and vertical velocities (e.g., due to settling velocity), they may settle on the bed or stay suspended (Fredsoe and Deigaard, 1992).

In sediment transport modeling, the bed-elevation is obtained using the Exner equation which includes bedload transport parametrization and pickup and settling rate of the suspended sediment. The sediment concentration is obtained using a suspended tracer advection-diffusion equation coupled with the hydrodynamic model through the velocity field (Warren and Bach, 1992; Lesser et al., 2004b; Warner et al., 2008, e.g.). In suspended load models, the sediments enter/exit the flow domain through a source/sink term in the bottom layer. The bedload transport is coupled with the flow model through either near-bottom velocity or shear stress (Amoudry, 2008). The most commonly-used bedload transport models are those of Meyer-Peter and Müller (1948); van Rijn et al. (1993); Soulsby and Damgaard (2005); van der A et al. (2013), and their capabilities are studied by various researchers. The averaging approaches have been a major restriction to the applicability range of these models. For instance, Rafati et al. (2021) successfully used XBeach (Roelvink et al., 2009), a depth-averaged model, to predict the offshore bar migration during storms. However, onshore migration could not be simulated using the same set of parameters. One approach to solving this problem is to use 3D models that solve the vertical flow structure, such as Delft3D (Lesser et al., 2004b). Success may again be limited in predicting beach accretion, due to wave-averaging procedure and inability of the bedload model to represent the effects of wave asymmetry (van Rijn et al., 2011). To finally predict the correct direction of sandbar migration, the SANTOSS transport formula (van der A et al., 2013), which follows the wave half-cycle concept to account for wave asymmetry and skewness, was implemented in Delft3D (Veen et al., 2014; van der Werf et al., 2015). Their model reproduced both onshore- and offshore-directed sandbar migration, but erosion and accretion conditions had to be calibrated separately (possibly due to an overspecification of undertow transport in the bedload formulation).

The present study builds on the mixed success of Delft3D using SANTOSS, but applying it to the Coastal and Regional Ocean Community (CROCO) model (Debreu et al., 2012; Marchesiello et al., 2021a), which is based on the Regional Oceanic Modeling System ROMS (Shchepetkin and McWilliams, 2005). The first attempt to add suspended-sediment-transport equations to ROMS was done by Blaas et al. (2007) for the Southern California shelf. Warner et al. (2008) supplemented the ROMS shelf sediment transport model with a bedload formulation and then implemented wave-averaged equations (using at first a three-dimensional radiation stress formulation) to address nearshore problems (Warner et al., 2010). Next, the successful implementation of vortex-force formalism by Uchiyama et al. (2010) was adopted by many others (Bennis et al., 2011; Kumar et al., 2012) and implemented in all versions of ROMS, and most notably in CROCO (Marchesiello et al., 2015). For wave forcing, CROCO can be fully coupled to either a ray-theory spectrum-peak wave propagation model (Uchiyama et al., 2010) or a more comprehensive spectral model, i.e., the third generation wave model WAVEWATCH III (Tolman et al.,
2009).

In all ROMS-related models mentioned above, bedload is calculated using either Meyer-Peter Müller (Meyer-Peter and Müller, 1948) or Soulsby-Damgaard (Soulsby and Damgaard, 2005) models. One of the main challenges of these models is to estimate their average contribution over each wave period. This leads to a poor representation of the effects of waves on bedload transport, i.e., velocity and acceleration skewness. As a result, bedload transport is underestimated by the model and fitting parameters are used to compensate for this defect (e.g., Rafati et al. (2021)). This problem can be solved as in the SANTOSS formulation of van der A et al. (2013) by using the half-cycle concept (e.g. Dibajnia and Watanabe, 1992a; Camenen and Larson, 2005). The half-cycle approach enables the model to calculate the amount of sediment mobilized and transported in each half-cycle of the wave. Therefore, provided an accurate estimate of the wave shape is given, the model can reasonably calculate the net sand transport due to wave asymmetry and skewness as well as phase lag between mobilization and transport at the intra-wave cycle — the sand entrained by the bed shear stress during a half wave-cycle is transported in the other direction during the next half wave-cycle. This makes the model more physically sound, more robust and less dependent on the calibration of parameters. As already mentioned, the SANTOSS formulation has been implemented in Delft3D, but also more recently in ROMS-COAWST (Kalra et al., 2019), on which the implementation in CROCO is based.

Modeling the evolution of the sandbar cross-shore profile has been a popular but challenging topic among coastal engineers studying the hydrodynamics and morphodynamics of the surf zone. Older models are reviewed by Roelvink and Broker (1993), focusing on a "process-based" approach that explicitly accounts for actual physical processes. Although the various models use similar governing equations, their application are often specific due to differences in parametrization, numerical algorithms, and boundary conditions (Spielmann et al., 2004; van Rijn et al., 2013). Some specific target include wave energy conditions, bathymetric complexity or duration of the desired simulation. For example, focusing on long-term evolution, the model of Walstra et al. (2012) successfully captured the dominant mechanisms in bar growth and decay. For short-term bar migration, Spielmann et al. (2004) proposed various effective parameters on bar migration in an erosion case, while Hsu et al. (2006) developed a model more suitable for the prediction of onshore bar migration. However, all of these models were limited to either offshore or onshore migration and were heavily dependent on the calibration of empirical parameters. In contrast, Hoefel and Elgar (2003) specifically considered wave-induced forcing and successfully simulated onshore bar migration, but their model failed under offshore bar migration conditions. With more success, the models of Zheng et al. (2014) and Rafati et al. (2021) predicted both onshore and offshore bar migrations, but again they used different values of the fitting parameters to achieve this result. This brief review of the literature reveals that a comprehensive hydro-morphodynamic model capable of predicting onshore/offshore bar migration is still lacking.

The present work aims at developing and implementing a comprehensive 3D nearshore model to
predict the onshore and offshore sandbar migration under storm and post-storm conditions, without the need to modify the setting parameters. To achieve this goal, the sediment transport model must be sufficiently extensive to represent the dominant transport processes under each forcing condition. In this paper, we take a first step by evaluating the ability of the model to predict sandbar migration in a two-dimensional wave channel (2DV) environment. We will demonstrate that the implementation of the SANTOSS transport formula within CROCO is a reliable and robust solution, provided that hydrodynamics are accurately reproduced and that a proper separation between wave and current forcing in the bedload and suspended load formulations is respected.

This manuscript is structured as follows. In Section 2, the hydrodynamic and sediment transport model equations are presented. The hydrodynamics are modeled by a set of 3D wave-averaged equations with wave-current interactions, coupled with a ray-theory spectrum-peak wave propagation model. The model description for sediment transport and morphodynamic evolution follows, with suspended load transport equations and bedload wave-related transport parameterization using the SANTOSS formula- tion. In Section 3, the numerical configuration including bathymetry and boundary forcing is explained. Then, Section 4 presents the hydrodynamic calibration procedure. The main results of the model are presented in Section 5, where the main drivers of sandbar migration are discussed. The discussion is then completed in Sections 6 and 7 providing the concluding remarks of this study.

## 2. Model formulation

In this paper, the Coastal and Regional Ocean Community (CROCO) model is implemented to study the morphological modeling of sandy beaches, including hydrodynamics, waves, sediment transport, and morphology. A non-hydrostatic version of the model is available that can explicitly solve for surface waves (Marchesiello et al., 2021a,b), but we rely here on a more traditional and less computationally expensive wave-averaged approach. Eulerian wave-averaged current equations for mass, momentum, and tracers are included based on a vortex-force formalism from the asymptotic theory of McWilliams et al. (2004), plus non-conservative wave effects due to wave breaking, roller waves, bottom streaming, wave-enhanced vertical mixing and bottom drag (Uchiyama et al., 2010; Marchesiello et al., 2015). The wave information is provided by a spectrum-peak WKB wave-refraction model that includes the effect of currents on waves (rather than WW3, the spectrum-resolving wave model). The velocity field is used to carry the suspended load (Blaas et al., 2007), while the recent SANTOSS bedload formulation is implemented in this study to include nonlinear wave effects.

### 2.1. Hydrodynamic model

#### 2.1.1. Flow model

CROCO solves the three-dimensional free-surface RANS (Reynolds-Averaged Navier-Stokes) equa- tions transformed into a terrain-following coordinate system, here with hydrostatic and Boussinesq ap-
proximations as in ROMS (Shchepetkin and McWilliams, 2005; Debreu et al., 2012). The Wave Effect on Currents (WEC) is introduced in the mean flow equations using a Vortex Force formalism (VF, (McWilliams et al., 2004)) rather than the radiation-stress representation. The effects of wave averaging are divided into Bernoulli head gradients (i.e., minus the gradient of the horizontal kinetic energy density of the primary and long waves) and a vortex force (i.e., the coupling of the flow vorticity with Stokes drift). This formalism was fully implemented, including non-conservative terms of acceleration and dissipation, into ROMS by Uchiyama et al. (2010), then in CROCO by Marchesiello et al. (2015). Here, we use a simplified set of two-dimensional vertical equations, neglecting the Coriolis force — but the implementation, including the new bedload transport formula, is fully 3D and is tested in a 3D configuration in parallel to this work. After reorganising advection and vortex form terms, these equations in 2DV Eulerian and Cartesian form are:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + w \frac{\partial u}{\partial z} + \frac{\partial \phi}{\partial x} = - \frac{\partial \mathcal{K}}{\partial x} + J + F^w + \frac{\partial}{\partial z} \left[ \nu^t \frac{\partial u}{\partial z} \right], \quad (1) $$

$$ \frac{\partial \phi}{\partial z} + \frac{g\rho}{\rho_0} = - \frac{\partial \mathcal{K}}{\partial z} + K, \quad (2) $$

$$ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \quad (3) $$

where t, x and z are time and cross-shore and vertical directions respectively; (u, w) are the Eulerian mean wave-averaged velocities in (x,z), and are equal to Lagrangian velocities minus stokes drift (u^st, w^st); g and φ are, respectively, the acceleration due to gravity and the dynamic pressure normalized by density ρ₀; ν^t is the turbulent viscosity that includes wave effects; K is the Bernoulli head and (J, K) are the horizontal and vertical components of vortex force terms (see McWilliams et al. (2004); Uchiyama et al. (2010) for details). F^w represents wave-induced non-conservative forces.

The right-hand side of these equations are the WEC terms: Bernoulli head gradients and vortex force, plus non-conservative forces (wave-induced acceleration and vertical mixing). However, in the actual implementation, the terms are rearranged, particularly the advection and vortex force terms that produce flux-form advection terms with Lagrangian advection velocities (Uchiyama et al., 2010). F^w is the sum of wave-induced acceleration processes (Equation (4)): wave-breaking and roller acceleration (B^b, Equation (5)) and wave-induced bottom streaming (B^wd, Equation (6)):

$$ F^w = B^b + B^{wd}, \quad (4) $$

$$ B^b = \frac{(1 - \alpha_r) \epsilon^b + \epsilon^r}{\rho_0 \sigma_p} f^b(z), \quad (5) $$

$$ B^{wd} = \frac{\epsilon^{wd}}{\rho_0 \sigma_p} f^{wd}(z), \quad (6) $$
where $\epsilon^b$, $\epsilon^r$ and $\epsilon^{wd}$ are wave dissipation rate due to wave breaking, surface rollers and wave bottom drag respectively and are calculated using Equations (7), (8), and (9). $f^{wd}(z)$ and $f^b(z)$ are vertical distribution functions. $\alpha_r$ (roller-fraction in CROCO, Uchiyama et al. (2010)) is the fraction of wave dissipation energy converted into rollers (i.e. surface bores of broken primary waves travelling onshore). $\rho_0$ and $\sigma$ are the mean water density and the angular wave peak frequency respectively (we will note $f_p = \sigma_p/2\pi$ the wave peak frequency). $\epsilon^b$, $\epsilon^r$ and $\epsilon^{wd}$ are given by:

$$ \epsilon^b = \frac{3\sqrt{\pi}}{16} \rho g f_p B_b^3 \frac{H_{rms}^3}{D} \left[ 1 + \tanh \left\{ 8 \left( \frac{H_{rms}}{\gamma_b D} - 1 \right) \right\} \right] \left[ 1 - \left\{ 1 + \left( \frac{H_{rms}}{\gamma_b D} \right)^2 \right\}^{-2.5} \right], \quad (7) $$

$$ \epsilon^r = \frac{0.1gE^r}{c}, \quad (8) $$

$$ \epsilon^{wd} = \frac{1}{2\sqrt{\pi}} \rho f_w |u_{orb}|^3, \quad (9) $$

The parameterization of $\epsilon^b$ follows Church and Thornton (1993), where $H_{rms}$ is the root mean square wave height (computed by the wave model in Sec. 2.1.2), and $D$ is wave-averaged thickness of the water column (i.e., the sum of resting depth of the ocean and local wave-averaged sea level, see Uchiyama et al. (2010)). $\gamma_b$ and $B_b$ (gamma_tg and B_tg in CROCO) are the empirical wave breaking parameters which need to be tuned. $\gamma_b$ is the ratio of wave height to water depth, and $B_b$ accounts for the type of breaker. In Equation (8), $c$ is the amplitude of phase velocity of the primary wave and $E^r$ is the roller energy related to the roller action density $\mathcal{A}^r = E^r/\sigma_p$ (Svendsen, 1984) calculated by solving a prognostic equation in the wave model (Eq. 17 in Sec. 2.1.2). In Equation (9), $f_w$ and $u_{orb}$ are wave friction factor and bottom wave orbital velocity respectively that are calculated in Equations (13) and (14) respectively.

In Uchiyama et al. (2010), the turbulent viscosity $\nu^t$ in the vertical mixing term is computed using a modified K-profile parameterization (KPP, Large et al. (1994)) to account for breaker-induced turbulence. CROCO has the additional option of a Generic Length Scale (GLS) approach (Warner et al., 2005). In the present study, we chose the $k-\omega$ GLS model, solving the closure equations for turbulent kinetic energy (TKE) and dissipation rate $\omega$. The introduction of wave breaking into this model follows Kumar et al. (2012), where the injection of surface TKE is given by surface boundary conditions in the closure equations. In the surf zone, only a fraction of the breaker dissipation produces turbulence (here we assume 5%) while the rest goes to the wave rollers. The surface boundary condition for $\omega$ is a function of the surface mixing length, which is taken as $0.5H_{rms}$.

The bottom drag (divided by $\rho$) is traditionally introduced as a bottom boundary condition of the vertical mixing term. Here, it is a combined wave-current shear stress formulated as in Soulsby (1995):

$$ \tau_{bot}^{cd} = \tau_c \left( 1 + 1.2 \left( \frac{\tau_w}{\tau_w + \tau_c} \right)^{3.2} \right), \quad (10) $$
$\tau_c$ is the shear stress due to currents:

$$ \tau_c = \rho \frac{\kappa^2}{\log^2 \frac{z_m}{z_0}} u_{z=z_m}^2, \qquad (11) $$

where $\kappa$ is the von Karman constant, $z_m$ the height of the first grid point above the bed, and $z_0$ (Zob in CROCO) is an empirically determined bedload roughness height (Blaas et al., 2007). $\tau_w$ is the shear stress due to waves:

$$ \tau_w = \frac{1}{2} \rho f_w u_{orb}^2, \qquad (12) $$

$$ f_w = 1.39 \left( \frac{\sigma_p z_0}{u_{orb}} \right)^{0.52}, \qquad (13) $$

$$ u_{orb} = \frac{\sigma_p H_{rms}}{2 \sinh kD}. \qquad (14) $$

### 2.1.2. Wave model

CROCO couples a spectrum-peak wave model with its flow model. It conserves the monochromatic wave action $\mathcal{A} = E/\sigma_p$ ($E = \frac{1}{8}\rho g H_{rms}^2$ is the depth-integrated wave energy) and cross-shore wave number $k_x$ (in 1D):

$$ \frac{\partial \mathcal{A}}{\partial t} + \frac{\partial (\mathcal{A} c_{g,x})}{\partial x} = - \frac{(\epsilon^b + \epsilon^{wd})}{\sigma_p}, \qquad (15) $$

$$ \frac{\partial k_x}{\partial x} + c_{g,x} \frac{\partial k_x}{\partial x} = -k_x \frac{\partial u}{\partial x} - \frac{\kappa \sigma_p}{\sinh 2\kappa D}, \qquad (16) $$

where $c_{g,x}$ is the cross-shore wave group velocity (considering the Doppler shift due to currents on the group velocity of the primary carrier waves). The prognostic equation for roller action is:

$$ \frac{\partial \mathcal{A}^r}{\partial t} + \frac{\partial (\mathcal{A}^r c_x)}{\partial x} = \frac{\alpha_r \epsilon^b - \epsilon^r}{\sigma_p}, \qquad (17) $$

where $c_x$ is the phase speed of the primary wave in x direction (again considering the Doppler shift due to currents). $\alpha_r$, the roller fraction, is generally taken in the range [0-1], but we noticed that with values greater than 1, an additional breaking delay is generated similar in essence to the mechanism parametrized in Delft3D by Roelvink et al. (1995). $\epsilon^b$, $\epsilon^r$ and $\epsilon^{wd}$ are parametrized as shown in Equation (7), (8) and (9). More details are provided in Uchiyama et al. (2010) and Weir et al. (2011).

### 2.2. Sediment-transport model

#### 2.2.1. Suspended-load transport

In the sediment transport model, the suspended sediment concentration is solved with an advection-diffusion equation (Blaas et al., 2007; Warner et al., 2008). Equation (18) solves the wave-averaged
sediment concentration C for each sediment class (multiple sediment classes can be implemented):

$$ \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + (w - W_s) \frac{\partial C}{\partial z} - E_s = -u^{st} \frac{\partial C}{\partial x} - w^{st} \frac{\partial C}{\partial z} - \frac{\partial}{\partial z} \left[ \nu_c^t \frac{\partial C}{\partial z} \right]. \quad (18) $$

Here again, this equation is implemented in CROCO in flux form by using the non-divergence of both Eulerian and Stokes drift velocities. No explicit diffusion is required since a third-order upstream biased advection scheme is used to reduce both numerical dispersion and diffusion and improve the effective resolution of the model (Soufflet et al., 2016). $W_s$ is the settling velocity that can be computed from the empirical formula of Soulsby (1997a) (Blaas et al. (2007)). $\nu_c^t$ is the turbulent tracer diffusivity including wave effects. $E_s$ represents the upward erosion flux (kg/m²/s) modeled as:

$$ E_s = (1 - p)E_0 \max \left( \frac{\tau_{sf} - \tau_{CE}}{\tau_{CE}}, 0 \right), \quad (19) $$

where the surface erosion rate $E_0$ and the critical bed shear stress for erosion $\tau_{CE}$ are parameters to be tuned (respectively, ERATE and TAU\_CE in CROCO)(Smith and McLean, 1977; Soulsby and Whitehouse, 1997), p is the bed porosity and $\tau_{sf}$ represents the bottom skin stress calculated as $\tau_{bot}^{cd}$ (Equation (10)), but with a skin roughness increased by ripples in the bar trough area.

### 2.2.2. Morphodynamics

The cross-shore bed evolution (variation in time of $z_b$, the height of the bed, Equation (20)) is calculated from the cross-shore divergence of sediment fluxes (Exner equation, Exner (1925, 1931)), which results from the difference between erosion and sedimentation of suspended sediments (calculated as explained in Section 2.2.1), as well as from the divergence of bedload fluxes $q_{b,slope}$ (which needs to be evaluated separately; see Section 2.2.3).

$$ \frac{\partial z_b}{\partial t} = -\frac{f_{mor}}{1-p} \left( \frac{\partial q_{b,slope}}{\partial x} - W_s \frac{\partial C}{\partial z} + E_s \right). \quad (20) $$

Equation (20) accounts for: a morphological acceleration factor $f_{mor}$ introduced by Lesser et al. (2004a); and the effect of bed slope on the bedload fluxes, introduced as a flux correction (Lesser et al., 2004a; Warner et al., 2008). The concept of morphological acceleration is based on the fact that morphodynamic changes are slower than hydrodynamic ones (van Rijn (1993)). In this case, the bed evolution can be accelerated without affecting the hydro-morphological solution. The avalanche module used to reproduce the gravitational flow of sand when the bottom slope exceeds the critical slope angle is presented in Equation (21)(Lesser et al., 2004a; Warner et al., 2008):

$$ q_{b,slope} = q_b \left( \frac{0.65}{(0.65 - \tan \beta) \cos \beta} \right), \quad (21) $$

where $q_{b,slope}$ is the modification of $q_b$ by considering the effect of the bed slope $\beta = \tan^{-1}(dz_b/dx)$. The value 0.65 is derived from the consideration of an angle of repose of 33°.
### 2.2.3. Bedload transport: an adapted SANTOSS formulation

The sediment flux $q_b$ can be calculated using different bedload models already implemented in CROCO, which are described in Warner et al. (2008), e.g., those of Meyer-Peter Müller (Meyer-Peter and Müller, 1948) or Soulsby and Damgaard (Soulsby and Damgaard, 2005). These models are more suited to continental shelf problems. In this paper, the near-bed total transport formulation SANTOSS proposed by van der A et al. (2013), suitable for nearshore applications, is adapted to CROCO, a depth-resolving wave-averaged 3D flow model. This adaptation is necessary because at least part of the current-related transport mechanisms in SANTOSS are already implemented in CROCO, and are thus duplicated in the bedload transport. In the problem addressed in this paper, we will see that double counting the effect of currents in SANTOSS leads to unrealistic results, while only marginal errors are imparted by neglecting these currents in SANTOSS and relying on the suspended load model for the current-related transport. Therefore, here, the current-related transport is solved by the flow model and the remaining physics (i.e., the bedload transport) is parameterized in SANTOSS. This formulation is implemented as in Kalra et al. (2019) and validated against a flume experiment. It is based on the half-wave cycle concept proposed by Dibajnia and Watanabe (1992b) that captures wave velocity and acceleration skewness and the effect of phase lag between mobilization and transport (Dohmen-Janssen et al., 2002).

In this section, we present our adapted version of the SANTOSS model. The modifications mainly concern the removal of the wave-averaged current effects to keep only the nonlinear wave effects in the formulation. In doing so, we follow up on the observation by Veen et al. (2014) that the full SANTOSS model would tend to double count the transport by mean currents in the wave boundary layer if a suspended load transport model is used in combination. In the case of onshore bar migration for example (i.e., LIP1c), this double counting prevents the replication of the onshore bar migration observed in the experiments and even reverses the migration (see Figure 13). In the brief presentation of the model below, we retain the terms associated with wave-averaged currents for completeness (they are implemented in practice with the option to disable them or not), but focus on wave effects.

The method to obtain sediment transport under asymmetric and skewed waves can be divided into three major steps (van der A et al., 2013; Kalra et al., 2019). In the first step, the asymmetric waveform based on the Ursell number is evaluated using wave statistics. The Shields parameter for each half cycle of the wave form is computed in the second step. Finally, a phase lag is estimated from the velocity field that determines the amount of bedload transported in the half cycle following the mobilization half cycle. The non-dimensional bedload transport rate $\Phi$ is thus given by:

$$ \Phi = \frac{q_b}{\sqrt{(s-1)gd_{50}^3}} = \frac{1}{T} \left[ \frac{\theta_c}{|\theta_c|^{1/2}} T_c \left( \Omega_{cc} + \frac{T_c}{2T_{cu}} \Omega_{tc} \right) + \frac{\theta_t}{|\theta_t|^{1/2}} T_t \left( \Omega_{tt} + \frac{T_t}{2T_{tu}} \Omega_{ct} \right) \right], \quad (22) $$

where $s = \rho_s/\rho$ ($\rho_s$ is the density of sand) and $d_{50}$ is the median grain diameter. $T$, $T_c$, $T_t$, $T_{cu}$ and $T_{tu}$ are the wave period, duration of wave crest half cycle, duration of wave trough half cycle, duration of
255

accelerating flow within the crest half cycle and duration of accelerating flow within the trough half cycle respectively (see Figure 1 for definitions), $\theta_c$ and $\theta_t$ represent the Shields numbers associated with the wave crest and trough half cycles. The sand load transported during the crest period is the combination of $\Omega_{cc}$ (mobilized during the crest period) and $\Omega_{tc}$ (mobilized during the trough period). Similarly, $\Omega_{tt}$ and $\Omega_{ct}$ are the sand load transported during the trough period (mobilized during the trough and crest periods respectively).

260

<figure><img src="image_2.png" /><figcaption>Figure 1: Definition sketch for the time series of near-bed velocity in the wave direction. The parameters $T_c$ and $T_t$ are the positive (peak) and negative (trough) flow durations. Similarly, $T_{cu}$ and $T_{tu}$ are the flow acceleration durations in the positive and negative x-directions respectively (van der A et al., 2013). The vertical axis is the free-stream velocity normalized by its maximum value.</figcaption></figure>

The sand load transported during each half-cycle is conventionally modeled according to a power law of Shields number:

$$ \Omega_i = \max \left( 11 \left( |\theta_i| - \theta_{cr} \right)^{1.2}, 0 \right), \quad (23) $$

where $\theta_{cr}$ is the critical Shields number (Soulsby (1997a)) and, hereafter, the subscript "i" is either "c" for crest or "t" for trough half cycles. To determine $\Omega_{ct}$ and $\Omega_{tc}$, i.e., the portion of the bedload remaining in suspension to be transported in the next half cycle, a phase lag parameter is evaluated (see van der A et al. (2013) for details). The Shields number for the peak or trough ($\theta_i = \theta_t$ or $\theta_c$) is calculated according to:

$$ \theta_i = \frac{\frac{1}{2} f_{w\delta i} |u_{i,r}| u_{i,r}}{(s-1)g d_{50}}. \quad (24) $$

$u_{i,r}$ is the representative cross-shore combined wave-current velocity at trough or crest half cycles calculated as:

$$ u_{i,r} = \frac{\hat{u}_i}{\sqrt{2}} + |u_\delta|, \quad (25) $$

where $\hat{u}_i$ is the peak crest or trough orbital velocities, $u_\delta$ is the steady current velocity at the top of the wave boundary layer. $f_{w\delta i}$ is the linear wave-current friction factor at crest or trough calculated by
(Ribbonink, 1998):

$$f_{w\delta i} = \frac{\hat{u}}{u_{\delta} + \hat{u}} f_{wi} + \frac{u_{\delta}}{u_{\delta} + \hat{u}} f_{\delta} \quad (26)$$

where $\hat{u}$ is the representative orbital velocity amplitude for the whole flow cycle (given by $\hat{u} = \sqrt{2}u_{orb}$). In Equation (26), $f_{\delta}$ is the current-related friction factor dependent on a current-related roughness $k_{s\delta}$ and $f_{wi}$ is the wave friction factor, calculated separately for the crest and trough half-cycles and depends on a wave-related roughness $k_{sw}$ (see Appendix A of van der A et al. (2013) for more detail). If the representative orbital excursion amplitude $\hat{a} = \hat{u}T/2\pi$ is large enough (i.e., greater than 1.587$k_{sw}$):

$$f_{wi} = 0.00251 e^{5.21 \left[ \left( \frac{2T_{iu}}{T_i} \right)^{2.6} \frac{\hat{a}}{k_{sw}} \right]^{-0.19}}, \quad (27)$$

otherwise, $f_{wi} = 0.3$.

If $u_{\delta}=0$ ($u_{i,r} = \hat{u}_i/\sqrt{2}$ and $f_{w\delta i} = f_{wi}$), the effect of currents is completely removed from the bed-load transport calculation. This choice represents our reference configuration to avoid double counting the transport by wave-averaged currents. The sensitivity to this modeling choice will be shown in the discussion Section 6.

Finally, following Kalra et al. (2019), after calculating $\Phi$, we apply to $q_b$ a bedload factor $f_{bld}$ (bed-load\_coeff in CROCO). This factor allows us to adjust the relative contribution to sediment transport of wave-induced bedload compared to the suspended load transported by mean currents. In this way, we have a better control of the antagonistic mechanisms that govern onshore and offshore transports respectively.

### 3. Numerical configuration

In this section, the initial and boundary conditions as well as the numerical settings of the simulations are presented. The LIP1b (erosion) and LIP1c (accretion) cases of the LIP11D experiments (Roelvink and Reniers, 1995) are used to evaluate and validate the model. LIP1b and LIP1c experiments are conducted sequentially, meaning that the final bathymetry of LIP1b is the initial bathymetry of LIP1c. The numerical domain is 201 meters long and water depth is 4.1 meters. The median grain diameter $d_{50}$ =0.22 mm, sand density $\rho_s$=2650 kg/m³, water density $\rho_0$=1000 kg/m³, and settling velocity $w_s$ mm/s (calculated using Soulsby formula). The wave boundary conditions including significant wave height and period, the Dean number and duration of wave generation are provided in Table 1.

Two sets of instruments were used for the measurements: fixed and movable by a carriage. In this paper, the measurements obtained from the velocity meter and suction tubes deployed on the carriage are used to evaluate the modelled cross-shore and vertical profiles of velocity and concentration respectively. During the LIP11D experiment, the carriage has been moved cross-shore in each wave hour; therefore, each vertical profile is obtained at different wave hour during the simulation. However, the simulation results presented in this paper are extracted at half of the total duration of the simulation (i.e. at 9 and
<table><caption>Table 1: Wave boundary conditions ($H_s$ and $T_p$), Dean number $\Omega$ and duration of LIP1b and LIP1c experiments</caption><thead><tr><th>Parameter</th><th>LIP1b</th><th>LIP1c</th></tr></thead><tbody><tr><td>$H_s$ [m]</td><td>1.4</td><td>0.6</td></tr><tr><td>$T_p$ [s]</td><td>5</td><td>8</td></tr><tr><td>$\Omega$</td><td>11.2</td><td>3</td></tr><tr><td>Duration [hour]</td><td>18</td><td>13</td></tr></tbody></table>

6.5 hours for LIP1b and LIP1c respectively). It should be noted that this assumption is a source of errors between the model simulation and the measurements.

As shown in Figure 2, LIP1b experiment starts with a developed sandbar. After 8 hours, the bar moves offshore by about 4 meters with a slight increase of its height. 10 hours later, it moves 5 meters further offshore, the crest becomes sharper and about 5 cm higher. During this erosion phase, the sand seems to be transported mostly from the bar trough and upper shoreface to the offshore side of the bar. 18 hours after the start of the experiment, wave conditions transition to the accretion phase (Table 1), resulting in a 4 m shoreward migration of the bar after 7 hours and a 7 cm increase in bar height. However, over the next 6 hours, the onshore migration slows down, while the bar crest continues to rise by about 5 cm. During the accretion phase, the sand seems to be transported mostly from the offshore side of the bar to its lee side. Note that, in both cases, well-developed ripples are apparent on the bar through, increasing the roughness of the bottom at this point.

The numerical domain is discretized using a uniform grid with horizontal resolution of 1.5 m — an optimal choice according to a sensitivity analysis (not shown). Longshore gradients are neglected, so there is only one grid cell in this direction. With terrain following coordinates, the number of vertical layers is 20 throughout the domain; therefore, the heights of the cells in the deep region and around the bar are about 21 cm and 5 cm respectively.

In the wave model, the offshore wave height is forced at the model boundary with values provided in Table 1. The significant wave height $H_s$ is 1.4 m and 0.6 m for LIP1b and LIP1c cases respectively (the model uses $H_{rms} = H_s/\sqrt{2}$) and the peak wave period is 5 and 8 s respectively. The resulting Dean number $\Omega = H_s/T_pW_s$ (Dean, 1973, 2018) clearly differentiates the erosion and accretion conditions from the experiments. Apart from the forcing conditions, all other wave model parameters are the same for both cases (see next section). For the sediment transport model, the morphological factors are 18 and 13 for the LIP1b and LIP1c cases respectively. Therefore, the simulation time is only 1 hour in both cases, while the real time is accelerated differently without any consequence on the results (as confirmed by experiments without morphological acceleration, not shown here). There are two main calibration parameters to be tuned in the suspended load model: the critical bed shear stress $\tau_{CE}$ and erosion rate $E_0$. $\tau_{CE}$ is set to
<figure><img src="image_3.png" /><figcaption>Figure 2: Initial, intermediate and final bathymetry of LIP1b and LIP1c</figcaption></figure>

0.18 N/m², following the empirical relation of Soulsby and Whitehouse (1997) (depending only on grain size), which controls the initiation of sand erosion in the suspended load model. The range of erosion rate values tested in this paper is 0 < E₀ < 0.025 kg/m²/s. Note that the grain roughness z₀s = d₅₀/12 (znots in CROCO) is used for the sediment bed roughness. The bedload model described in Section 2.2.3 has its own empirical formulations for sand mobilization based on grain size and wave statistics. We only used one calibration parameter, the bedload factor, which was tested in the range 0 < f_bld < 1.

## 4. Hydrodynamic calibration

In this section, before the full hydro-morphodynamic simulation, the hydrodynamic model is independently calibrated. The sediment model is inactive and intermediate bathymetries (i.e., the bathymetry after 8h for LIP1b and 7h for LIP1c cases) are used throughout the simulations. First, the best-fit solutions are presented. Then, the most sensitive hydrodynamic parameters are tested and discussed. To avoid overload, some of the figures are presented in the "Appendix: hydro-sensitivity". The calibrated hydrodynamic parameters will be used in the morphodynamic calibration (Section 5).

### 4.1. Reference solutions

The best parameter values are chosen as those that provide the best model-data comparison both using the normalised Root-Mean-Square Error (nRMSE, Equation 28) value (Roelvink et al., 2009) and qualitatively through figures presented in this paper, with major emphasis on the position of the bar
crest in both the erosion and accretion phases. In Equation 28, *m* stands for measured value and *c* for computed value, ⟨●⟩ represents an averaging operator and *σ* the associated standard deviation. Note that the presented nRMSE values of velocity are the linear average of the nRMSEs of each vertical velocity profile.

$$nRMSE = \frac{\sigma_{m-c}}{\max(\sigma_m, \langle m \rangle)} \qquad (28)$$

In Figure 3, the top panel compares the vertical profiles of velocity across the beach, the middle panel shows the cross-shore profile of $H_{rms}$ and wave set-up and the bottom panel shows the near-bottom velocity (at 10 cm above the bed). This Figure is obtained by using $\gamma_b = 0.4$, $B_b = 0.6$ and $\alpha_r = 1.5$. The bottom roughness in the offshore and onshore side of the bar are set to 0.002 and 0.01 m respectively. The list of all the parameters tested is presented in Table 2; their corresponding normalised RMSE values of undertow and wave height are also provided. It is important to note that using the same settings for the Soulsby transport model (Soulsby and Damgaard (2005), the bar migrates offshore in both cases (see the respective results in D.1).

On the offshore side of the bar in LIP1b ($60 < x < 100$), the incident waves are already influenced by the bottom and begin to break. Therefore, the wave height decreases as the waves propagate toward the shoreline and the velocity near the bottom (used as a proxy for undertow) is significant. On the other hand, in the offshore region of LIP1c, the incident waves are not large enough to start the breaking process. Consequently, the wave height, wave set-up, and undertow remain unaffected in this region.

Over the bar, a strong undertow is generated in both cases (stronger in LIP1b). This is due to strong breaking on the bar that transfers momentum, and also causes a drop in wave height and wave set-up. Since LIP1b has more incident momentum, more of it is transferred on the bar than in LIP1c, leading to a larger cross-shore pressure gradient and stronger undertow. The model predicts the wave height and wave set-up well in both cases. The undertow is also fairly well predicted, but the results are more sensitive to the wave-breaking parameters as will be shown.

On the onshore side of the bar, in both cases, hydrodynamic changes are generally less important. The undertow and wave height tend to decrease while the wave set-up increases, although the presence of a shallow terrace, beginning at about $x=160$ m, generates a secondary breakpoint. The model follows the observed trend, but it seems to underestimate the undertow, particularly in LIP1b.

In LIP1c, the model shows a secondary undertow maximum associated with the secondary breakpoint, which is absent in the flume experiment. This may reflect a lack of synopticity in the measurements — the currents have been measured at different times corresponding to different bathymetries — or a modeling issue in the shallow nearshore region (e.g., overestimation of wave-setup).

As can be observed, the undertow is underestimated in 1b and overestimated in 1c. Although the results of each case can be improved by calibrating the wave breaking parameters independently, this is the best result we could achieve upon numerous tests using the same parameters for both cases. It
365

370

should be reminded that the lack of synopticity in the measurements is another source of error (due
to the movement of the bar during the measurement). In the LIP experiments, each velocity profile is
measured at different wave hours while we have fixed the bathymetry using the bed profile measured in
the middle of the experiment. For example, activating the sediment transport model, strengthens the
near-bed velocity to -0.29 *m/s* after 5 hours from the initial time (the same wave hour at which the
measurement is conducted). This will change the best choice of wave breaking parameters. Although
some studies present simulation results according to the corresponding wave hour of each vertical profile
(e.g., Marchesiello et al. (2021b)), we present our results obtained in the middle of the simulation time
period because they are satisfactory enough.
<figure><img src="image_4.png" /><figcaption>Figure 3: Model-data comparison of cross-shore velocity profiles, wave height $H_{rms}$, wave set-up (eta), and near-bottom current (undertow) using the same tuning parameters for both LIP1b and LIP1c (reference simulations). Wave set-up (eta) is shown using dashed line.</figcaption></figure>
4.2. Sensitivity to wave breaking parameters $\gamma_b$ and $B_b$

Wave breaking parameter $\gamma_b$ controls the cross-shore position of the breakpoint. Figure 4 illustrates the sensitivity of model results to this parameter. As expected from Equation (7), when $\gamma_b$ is increased, the waves tend to break closer to shore and release their energy in shallower water. This generally increases the wave height and decreases the wave set-up in the surf zone.

On the offshore side of the bar, LIP1b is particularly sensitive to $\gamma_b$ because waves start to break earlier in this case due to higher offshore wave height. The top panel of Figure 4(a) shows that at $x = 65$ m, an increase of $\gamma_b$ results in a decrease of breaker-induced current shear. What remains is the Eulerian return flow, which verifies the Lagrangian continuity equation (non-divergent Lagrangian flow). Due to the delayed breaking with a higher $\gamma_b$ in LIP1b, there is a large increase in wave breaking on the bar with an increased undertow. As a result, the breaking pattern in LIP1b and LIP1c becomes more similar. On the contrary, for lower $\gamma_b$, the breaking occurs earlier and the offshore points have higher current shear, while the onshore points have reduced wave height and undertow.

On the bar, the effects of increasing $\gamma_b$ is opposite for the erosion and accretion cases: in LIP1b, higher $\gamma_b$ moves the breaking point from offshore to the bar crest. The higher the $\gamma_b$ is the stronger the undertow becomes. In contrast, in LIP1c with a higher value of $\gamma_b$, wave breaking is stronger inshore of the bar and on the terrace slope, which is the sensitive area for LIP1c. We choose $\gamma_b = 0.4$ as the best compromise in both cases.

The second tuning parameter in the wave dissipation equation is $B_b$ (Equation (7)). It controls the amount of wave energy dissipated during breaking. The effects of increasing $B_b$ is similar to decreasing $\gamma_b$ throughout the domain and in both cases, although the model is less sensitive to $B_b$ than $\gamma_b$. The sensitivity of the model to this parameter is shown in Appendix B.1. $B_b = 0.6$ is chosen as the best value to be used in the following.
<figure><img src="image_5.png" /><figcaption>Figure 4: Model sensitivity to the breaking parameter $\gamma_b$ for experiments LIP1b and LIP1c, showing the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow.</figcaption></figure>
4.3. Sensitivity to the bottom roughness length $Z_{0b}$

Another important source of dissipation is the bottom friction (Equation (10)). The logarithmic law of the wall requires for the bottom stress a value for the bottom roughness length $Z_{0b}$ (Equations (11) and (13)). This coefficient affects both current- and wave-induced friction. However, the friction force due to waves is about one order of magnitude larger than that due to mean currents in the present case. In this model, due to the existence of well-developed ripples in the trough area between the bar crest and the terrace (Figure 2), higher roughness is considered for the onshore side of the bar than that of the offshore side. The offshore-onshore transition point is set to be at $x = 140m$ in both LIP1b and LIP1c cases. To avoid sharp changes of roughness, we use a hyperbolic tangent function.

The sensitivity of the model to this coefficient is relatively moderate. It is illustrated for both erosion and accretion cases in Figure 5 (for offshore roughness) and Appendix B.2 (for onshore roughness). As expected, $Z_{0b}$ mostly affects the near-bottom velocity: higher values lead to lower velocities. It should be noted in Figure 5 that the sensitivity of $H_{rms}$ in LIP1c is due to the fact that, since the wave energy is moderate, the dissipation by bottom friction is comparable to that of the wave breaking in this case. The best value for offshore $Z_{0b}$ (2 mm) is chosen by giving more weight to the offshore near-bottom velocity. Inshore of the bar, the near-bottom velocity is over- and under-estimated in LIP1b and LIP1c respectively. Therefore, the best compromise is found to be with $Z_{0b}$=10 mm for the onshore side of the bar.
<figure><img src="image_6.png" /><figcaption>Figure 5: Model sensitivity to the offshore bottom roughness $Z_{0b}$ for experiments LIP1b and LIP1c, showing the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow.</figcaption></figure>
4.4. Sensitivity to bottom streaming and roller fraction

Two other parameters in Equation (4) for which sensitivity analysis has been performed are bottom
streaming (Equation (6)) and roller fraction (Equations (5) and (17)). As the sensitivity of the results
to these parameters are weaker than those to γb and offshore Z0b, the results are shown in the appendix
section rather than in the main text (Appendices B.3 and B.4). Bottom streaming is a wave-averaged
Eulerian shoreward flow caused by phase differences in horizontal and vertical orbital velocities occurring
in the wave boundary layer (Longuet-Higgins (1953) and Xu and Bowen (1994)). Therefore, this term
tends to counteract the undertow near the bottom. Note that the original SANTOSS formulation also
accounts for wave streaming, but we have neglected it here since it is already accounted for in the flow
model (Uchiyama et al., 2010).

The roller fraction α_r is the fraction of energy transferred from the waves to the rollers during wave
breaking, thus expanding the surf zone shoreward compared to a local transfer at the breakpoint. With-
out roller (α_r=0), the undertow is localized on the bar crest. Using a value of α_r higher than 1 is
unconventional as it implies getting some of the roller energy from the mean flow. However, it acts as an
additional breaker delay (Roelvink et al. (1995)) that seems particularly relevant to the breaking process
on a sandbar. As we observe in Figure B.4, in the measurement, wave breaking effect on the near-bed
velocity is spread throughout the inner surf zone while in the model it is more spatially concentrated over
the primary bar; this may be due to under-prediction of turbulence in the flow model (not shown here).
Using a roller fraction greater than unity inserts additional energy into the bores, which represents this
dissipation delay in the model. For this reason, we decided to challenge the roller model and choose a
value of 1.5 that is greater than unity.
<table><caption>Table 2: Parameters studied in the hydrodynamic calibration (the chosen values in each sensitivity test is shown in bold)</caption><thead><tr><th rowspan="2">Run</th><th colspan="5">hydrodynamic parameters</th><th colspan="2">LIP1b</th><th colspan="2">LIP1c</th></tr><tr><th>γ<sub>b</sub></th><th>B<sub>b</sub></th><th>Z<sub>0</sub> [m]<br>(offshore)</th><th>Z<sub>0</sub> [m]<br>(onshore)</th><th>α<sub>r</sub></th><th>nRMSE_H<br>[-]</th><th>nRMSE_U<br>[-]</th><th>nRMSE_H<br>[-]</th><th>nRMSE_U<br>[-]</th></tr></thead><tbody><tr><td>1</td><td>0.3</td><td>0.6</td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.063</td><td>0.224</td><td>0.101</td><td>0.91</td></tr><tr><td>2</td><td><b>0.4</b></td><td>0.6</td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.048</td><td>0.258</td><td>0.075</td><td>1.13</td></tr><tr><td>3</td><td>0.5</td><td>0.6</td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.066</td><td>0.32</td><td>0.073</td><td>1.31</td></tr><tr><td>4</td><td>0.4</td><td>0.5</td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.057</td><td>0.29</td><td>0.071</td><td>1.25</td></tr><tr><td>5</td><td>0.4</td><td><b>0.6</b></td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.048</td><td>0.258</td><td>0.075</td><td>1.13</td></tr><tr><td>6</td><td>0.4</td><td>0.7</td><td>0.002</td><td>0.01</td><td>1.5</td><td>0.049</td><td>0.237</td><td>0.084</td><td>1.03</td></tr><tr><td>7</td><td>0.4</td><td>0.6</td><td>0.0005</td><td>0.01</td><td>1.5</td><td>0.052</td><td>0.218</td><td>0.091</td><td>0.71</td></tr><tr><td>8</td><td>0.4</td><td>0.6</td><td><b>0.002</b></td><td>0.01</td><td>1.5</td><td>0.048</td><td>0.258</td><td>0.75</td><td>1.13</td></tr><tr><td>9</td><td>0.4</td><td>0.6</td><td>0.005</td><td>0.01</td><td>1.5</td><td>0.045</td><td>0.36</td><td>0.06</td><td>1.44</td></tr><tr><td>10</td><td>0.4</td><td>0.6</td><td>0.002</td><td>0.008</td><td>1.5</td><td>0.048</td><td>0.245</td><td>0.074</td><td>1.12</td></tr><tr><td>11</td><td>0.4</td><td>0.6</td><td>0.002</td><td><b>0.01</b></td><td>1.5</td><td>0.048</td><td>0.258</td><td>0.075</td><td>1.13</td></tr><tr><td>12</td><td>0.4</td><td>0.6</td><td>0.002</td><td>0.02</td><td>1.5</td><td>0.052</td><td>0.295</td><td>0.081</td><td>1.15</td></tr><tr><td>13</td><td>0.4</td><td>0.6</td><td>0.002</td><td>0.01</td><td>1</td><td>0.05</td><td>0.265</td><td>0.074</td><td>1.11</td></tr><tr><td>14</td><td>0.4</td><td>0.6</td><td>0.002</td><td>0.01</td><td><b>1.5</b></td><td>0.048</td><td>0.258</td><td>0.075</td><td>1.13</td></tr><tr><td>15</td><td>0.4</td><td>0.6</td><td>0.002</td><td>0.01</td><td>2</td><td>0.04</td><td>0.251</td><td>0.076</td><td>1.16</td></tr></tbody></table>

## 5. Morphodynamic simulations: sandbar migration

Here we pursue the objective of evaluating the model to predict accretion and erosion conditions with the same set of parameters. This capability of the model is tested by doing sensitivity analysis to two tuning parameters following an iterative approach. Due to the disregard of current effects in the bedload formulation (i.e., our adaptation of the SANTOSS formulation to a 3D flow model), this iterative calibration process has become simpler: while the current-driven offshore bar migration is controlled by the erosion rate in the suspended load formulation, the control of wave-driven onshore bar migration is achieved by the bedload factor. Considering the onshore-directed migration of the bar, we expect the sandbar migration to be controlled by bedload transport in the accretion case. Consequently, in the first step of our iteration procedure, we focus on the model sensitivity to the bedload factor in LIP1c, with no suspended load ($E_0 = 0 \text{ kg/m}^2/\text{s}$). Second, using the bedload factor selected in the first step, we find the most realistic erosion rate for LIP1b where the suspended-load transport is dominant. Then, the same procedure as the first step is repeated using the updated value of the erosion rate. This iteration is repeated until the differences in two successive iteration cycles become small enough. Finally, the values
by which the model predicts both onshore and offshore sandbar migration are chosen as the best-fit case.
More details are presented in the next section. Note that the critical bed shear stress for erosion is not
considered as a tuning parameter but is calculated using the Soulsby formulation as in Blaas et al. (2007).

5.1. Sensitivity to the bedload factor

As mentioned above, we start the morphodynamic calibration procedure by considering zero erosion
rate $E_0 = 0 \text{ kg/m}^2\text{/s}$ (equation (19)). Therefore in this case, the bed height in Equation (20) depends
only on near-bed transport. Figure 6 illustrates the sensitivity of the model to the bedload factor with
zero suspended load. In this figure, each part consists of two panels: the top panel shows the modeled final
bed elevation and compares it with the measured initial and final bed levels and the bottom panel shows
the differences between the initial and final bed levels ($z_b^i$ and $z_b^f$ respectively) for both measurement and
simulation.

In LIP1c, the bar moves onshore (figure 6(b)) which is consistent with the flume experiment. With
higher bedload factor $f_{bld}$, a larger mass of sand is moved from the offshore to the lee side of the bar,
which leads to both onshore migration and amplitude increase of the bar. The same phenomenon occurs
for the terrace around $x = 160 \text{ m}$ which tends to form a secondary inner bar in the model. The bar also
moves toward the coast in LIP1b, which is contrary to the experiment (Figure 6(a)). In this case, as
explained by Marchesiello et al. (2021b), wave asymmetry and skewness are stronger in LIP1b than in
LIP1c, resulting in greater wave-driven onshore sand transport. This shows the considerable importance
of the suspended load for the offshore migration of the sandbar (see Section 5.2).

Note that a sensitivity numerical experiment (not shown) without the phase lag between mobilization
and transport in the SANTOSS formulation (setting $\Omega_{ct}$ and $\Omega_{tc}$ to zero in Eq. 22) indicates that the
main driver of wave-related transport in these experiments is not the phase lag process (because the sand
is not fine and settling velocity is high), but the net result of wave asymmetry or skewness. In what
follows, we will consider wave nonlinearity as the dominant process of wave-related sediment transport
along with wave streaming effects.

Consequently, we start the iteration process by considering $f_{bld}=0.6$ in our sensitivity test to erosion
rate (explained in the next section). In the second cycle of our iteration procedure (after choosing the
best erosion rate for LIP1b $E_0=0.02 \text{ kg/m}^2\text{/s}$), we reconnect the suspended load transport in LIP1c. The
results are not drastically changed but the undertow, even under these moderate conditions, requires a
compensating increase in $f_{bld}$ from 0.6 to 0.8 for better comparison to the data ($f_{bld}=0.8$ is chosen as the
final value due to its slightly better results, compare Figures 7 and C.2). Yet, in LIP1c, the $f_{bld} = 0$ case
shown in Figure C.1(b) confirms that suspended load transport plays a minor role in sandbar migration,
although it contributes more significantly to terrace shape.

In LIP1b, the addition of suspended load transport is so overwhelming that it changes the direction
of the bar migration, while decreasing the sensitivity to the bedload factor (compare Figure C.1(a) with
485

Figure 6(a)). The wave-driven onshore bar migration becomes a minor process in comparison with the current-driven offshore migration associated with suspended load. However, neglecting bedload transport results in an overestimation of the offshore bar migration.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Model sensitivity to the bedload factor $f_{bld}$ with erosion rate $E_0=0$ kg/m$^2$/s in LIP1b (a) and LIP1c (b). In each case, the upper panel shows the bed level for the model and data in initial and final states; the lower panel shows the difference between initial and final bed levels ($z_b^i$ and $z_b^f$ respectively).</figcaption></figure>
## 5.2. Sensitivity to erosion rate

The second step in each cycle of the calibration procedure is to adjust the erosion rate to obtain the most accurate offshore bar migration possible, retaining the best bedload factor from the first step ($f_{bld}=0.6$). Figure 7(a) illustrates how the prediction of the offshore bar migration in LIP1b is improved by considering the suspended load. The vertical sand concentration profiles of LIP1b is shown in Figure 8. These profiles are obtained in the middle of the simulation period while the corresponding measurements had been acquired with a wider spread in the sampling period. Again, the lack of synopticity is a source of discrepancy between the model and the experiments since the change in bathymetry during the experiment affects the hydro-sedimentary fields. However, in general, the sand concentration and its vertical gradients are well simulated, confirming the ability of the 3D flow model to simulate mixing processes. In LIP1c, as already discussed in Section 5.1, the low sand concentration and low undertow make $E_0 \ kg/m^2/s$ have even less impact on the bar, and its most notable effect is on the terrace. Note that, as expected, with $E_0 = 0 \ kg/m^2/s$, the concentration is zero throughout the domain during the simulations of LIP1b and LIP1c (not shown).

Figure 8 shows that, on the offshore side of the bar, stronger erosion rate $E_0$ results in an increase in suspended sand concentration on the bar, as expected from the sand source term in Equation (18). Therefore, the undertow transports a larger amount of sand offshore, leading to a greater offshore bar migration without notable change in the bar amplitude. On the other hand in LIP1c, the erosion rate does not affect the position of the bar crest, although it affects the bar shape and associated distribution of sand around the crest. One reason is that the undertow in LIP1c is weaker and restricted to the inshore side of the bar as also observed by Marchesiello et al. (2021b) using a wave-resolving model. Vertical turbulent mixing is also lower, which does not promote suspended sediment transport and limits it to a local area (the inshore side of the bar). Finally, as mentioned, the suspended sediment transport is dominant in LIP1b, making it less sensitive to the bedload factor.
<figure><img src="image_9.png" /><figcaption>Figure 7: Model sensitivity to erosion rate $E_0$ with bedload factor $f_{bld}=0.6$ in LIP1b (a) and LIP1c (b) (same panels as Figure 6)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Model sensitivity of vertical sand concentration profiles to the erosion rate in the erosion case LIP1b, with bedload factor $f_{bld}$=0.6 - the top panel illustrates the cross-shore position of the profiles in comparison to the intermediate bathymetry after 8 hours</figcaption></figure>

### 5.3. Optimization of the morphological calibration

This section proves the ability of the model to predict onshore and offshore sandbar migration using the same model tuning parameters. Choosing the most suitable values for erosion rate $E_0$ and bedload factor $f_{bld}$ is tricky because deviation of the simulation from measurement varies spatially across the beach. A classical metric to evaluate the quality of a morphological prediction is the BSS Sutherland et al. (2004) defined as follows (Brier et al., 1950):

$$BSS = 1 - \frac{Var(c - m_0)}{Var(m_f - m_0)}, \qquad (29)$$
where *Var* is the variance of the quantity between brackets, *m*₀ stands for initial measured bed profile, *m*_f for final measured bed profile and *c* for computed bed profile. The BSS relates the variance of the difference between data and model to the variance of the data: BSS=1 means perfect skill, BSS=0 means no skill, BSS<0 means model is worse than no change scenario Roelvink et al. (2009).

Considering error metrics, nRMSE and BSS, (see Table 3) as well as the position of bar crest and maximum bed evolution (see Figure 9) altogether, we chose the best values to be 0.02 kg/m²/s and 0.6 for erosion rate and bedload coefficient respectively. Although using $f_{bld} = 0.8$ (i.e., run#23) predicts the cross-shore position of the bar slightly better, it leads to negative value of BSS in LIP1c. The reason is that the measured bed evolution is small on the offshore toe of the primary bar and in the shallow region while our model has the largest errors in those regions. Consequently, despite the fact that both cases predict the bar position well, run#22 is preferred considering that it has the higher BSS for both LIP1c (BSS=0.142) and LIPB1b (BSS=0.142). Respectively, these values are considered to be reasonable/fair and good-excellent prediction according to Sutherland et al. (2004). Note that in Figure 9 the bottom panel illustrates the vertical concentration profiles at their corresponding position in the domain.

Both the offshore and onshore bar migration (LIP1b and LIP1c respectively) are generally well predicted by the model. In LIP1b, the offshore side of the bar follows the same bed profile as in the experiment, but the model overestimates the accumulation of sand in the lee side. This excessive sand deposition is reflected in an overestimation of sand concentration near the bottom. This overestimation also appears in LIP1c both on the bar and the terrace. As mentioned earlier, one reason may be the lack of synopticity in the data. For example, in LIP1b, the near-bottom concentration (5 cm above the bed) at x = 138 m decreases from 6.9 g/l (Figure 9(a)) to 3.6 g/l when viewed at 5 h from the initial state (the same wave hour at which the concentration is measured). However, another reason in low energy conditions, especially in the shallow waters of the terrace, may be the overestimation of the undertow (as a result of the overestimation of wave height; see Section 4.1) that increases the bed shear stress and erosion flux $E_s$.
<figure><img src="image_11.png" /><figcaption>(a) 1b</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) 1c</figcaption></figure>

Figure 9: Reference morphodynamic simulations in LIP1b (a) and LIP1c (b) (same panels as Figure 6 with the addition of a lower panel presenting the vertical sand concentration profiles in various cross-shore positions).
# 6. Discussion

In this section, we will further explain the improvements this study provides to our understanding of sandbar transport modeling. First, we will use the maximum bed-evolution peak as a proxy to sandbar migration to better assess the sensitivity of the results to bedload and suspended load coefficients. Next, we will use a sand transport budget to assess the contributions of bedload and suspended load (Section 6.2) and evaluate the role of wave nonlinearity parameters in the onshore transport process (Section 6.3). Then, the improvements associated with the modification of the SANTOSS formulation will be discussed. Finally, we will test the performance of the model with the same parameters on another laboratory experiment of cross-shore sand transport.

## 6.1. Model robustness

To better assess the robustness of the onshore/offshore bar migration prediction, the result of 60 numerical simulations are presented in Table 3, corresponding to different values of bedload factor $f_{bld} \in \{0, 0.25, 0.4, 0.6, 0.8, 1\}$ and erosion rate $E_0 \in \{0, 1, 1.5, 2, 2.5\} \times 10^{-2} \text{ kg/m}^2/\text{s}$ (the same values are used for 1b and 1c cases). Moreover, the nRMSE values (see Equation 28) of concentration profiles as well as the normalised nRMSE and BSS of bed evolution (at $100 < x < 175$) are presented. $nRMSE_C$ is the linear average of the nRMSE values of each vertical concentration profile and, in general, it increases with the erosion rate parameter. However, this did not affect the procedure of choosing the best values since the best cases are chosen considering the final bed profile rather than concentration profiles. On the other hand, the lowest $nRMSE_{bed}$ values are obtained when $f_{bld}=0.4, 0.6$ or $0.8$ which is consistent with our best choice considering the relatively small variation of nRMSE values. Moreover, BSS is closer to unity with higher erosion rates.

Figure 10 shows the distribution of bar height evolution (maximum value of $z_b^f - z_b^i$ on the y-axis) versus cross-shore bar migration (horizontal displacement $x_{bar}^f - x_{bar}^i$ on the x-axis). In Figure 10(a) for LIP1b, the symbol is in the upper left quadrant corresponding to an offshore bar migration with an increase in bar height evolution, while in Figure 10(b) for LIP1c, the symbol is in the upper right quadrant corresponding to an onshore bar migration also with an increase in bar height evolution. In these figures, each symbol type represents a simulation performed with a given value of the bedload factor $f_{bld}$ while each color corresponds to a given value of the erosion rate $E_0$. For the LIP1b erosion case, Figure 10(a) shows that, with the exception of $E_0 = 0 \text{ kg/m}^2/\text{s}$, all model results are in the correct quadrant and, as expected, both the offshore bar migration and bar height evolution increase with erosion rate. As the erosion rate increases from $10^{-2}$ to $2.5 \times 10^{-2} \text{ kg/m}^2/\text{s}$, the sensitivity of the model prediction to bedload coefficient decreases significantly. For the LIP1c accretion case, Figure 10(b) shows a weaker sensitivity than in 1b case. Consistently, in 1c also higher erosion rate leads to less onshore migration, but with the decrease of bar height evolution.
It is interesting to note that in both the erosion and accretion cases, the increase in the bedload factor
moves the bed shoreward and upward. In the absence of suspended load in the erosion case, in 1b, the
increase in the bedload factor tends to increase only the height of the bar with little shoreward movement
while, in 1c, lower erosion rate leads to higher sensitivity to the bedload coefficient.

Although the two selected parameters move the sandbar in opposite directions, increasing one does not have the same effect as decreasing the other. The main difference is in their effects on sand accumulation. For example, considering the case of ($f_{bld}$, $E_0$)=(0.8, 0.015) (with the symbol △ in blue color), in LIP1b, increasing $E_0$ from 0.015 to 0.02 kg/m²/s (with the symbol △ in cyan color) increases the offshore migration of the bar while accumulating sand on the bar crest (bar sharpening). On the other hand, decreasing $f_{bld}$ from 0.8 to 0.6 (with the symbol ∇ in blue) also leads to an increase in the offshore bar migration but lowers the bar height evolution during the migration process (bar smoothing). This is because wave-related onshore transport, even when dominated by current-related offshore transport, nevertheless contributes to the convergence of sand on the bar. However, when bedload transport dominates, decreasing the bedload coefficient and increasing the erosion rate have the same effects with more sensitivity on bedload coefficient.
<table><caption>Table 3: Values of erosion rate and bedload factor presented in this study along with their respective marker used in Figure 10; the error metrics are also provided. The case with bold-faced values (Run#22) is chosen as the best case</caption><thead><tr><th rowspan="2">Run<br>Id</th><th rowspan="2">f<sub>bld</sub><br>[kg.m<sup>-2</sup>.s<sup>-1</sup>]</th><th rowspan="2">E<sub>0</sub><br>marker</th><th colspan="3">LIP1b</th><th colspan="3">LIP1c</th></tr><tr><th>nRMSE<sub>C</sub></th><th>nRMSE<sub>bed</sub></th><th>BSS<sub>bed</sub></th><th>nRMSE<sub>C</sub></th><th>nRMSE<sub>bed</sub></th><th>BSS<sub>bed</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>0</td><td>o</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2</td><td>0.25</td><td>0</td><td>+</td><td>0.465</td><td>1.27</td><td>-0.63</td><td>0.543</td><td>0.756</td><td>0.279</td></tr><tr><td>3</td><td>0.4</td><td>0</td><td>x</td><td>0.465</td><td>1.51</td><td>-1.27</td><td>0.543</td><td>0.718</td><td>0.238</td></tr><tr><td>4</td><td>0.6</td><td>0</td><td>∇</td><td>0.463</td><td>1.82</td><td>-2.3</td><td>0.543</td><td>0.879</td><td>-0.150</td></tr><tr><td>5</td><td>0.8</td><td>0</td><td>Δ</td><td>0.463</td><td>2.12</td><td>-3.43</td><td>0.543</td><td>1.15</td><td>-0.848</td></tr><tr><td>6</td><td>1</td><td>0</td><td>*</td><td>0.463</td><td>2.36</td><td>-4.47</td><td>0.543</td><td>1.39</td><td>-1.54</td></tr><tr><td>7</td><td>0</td><td>0.01</td><td>o</td><td>0.423</td><td>0.634</td><td>0.600</td><td>0.908</td><td>1.07</td><td>-0.117</td></tr><tr><td>8</td><td>0.25</td><td>0.01</td><td>+</td><td>0.413</td><td>0.926</td><td>0.158</td><td>0.717</td><td>0.801</td><td>0.235</td></tr><tr><td>9</td><td>0.4</td><td>0.01</td><td>x</td><td>0.418</td><td>1.15</td><td>-0.304</td><td>0.733</td><td>0.717</td><td>0.288</td></tr><tr><td>10</td><td>0.6</td><td>0.01</td><td>∇</td><td>0.435</td><td>1.39</td><td>-0.911</td><td>0.752</td><td>0.752</td><td>0.126</td></tr><tr><td>11</td><td>0.8</td><td>0.01</td><td>Δ</td><td>0.458</td><td>1.57</td><td>-1.42</td><td>0.76</td><td>0.917</td><td>-0.274</td></tr><tr><td>12</td><td>1</td><td>0.01</td><td>*</td><td>0.466</td><td>1.78</td><td>-2.10</td><td>0.77</td><td>1.14</td><td>-0.882</td></tr><tr><td>13</td><td>0</td><td>0.015</td><td>o</td><td>0.677</td><td>0.512</td><td>0.74</td><td>1.13</td><td>1.12</td><td>-0.221</td></tr><tr><td>14</td><td>0.25</td><td>0.015</td><td>+</td><td>0.644</td><td>0.677</td><td>0.553</td><td>0.873</td><td>0.859</td><td>0.154</td></tr><tr><td>15</td><td>0.4</td><td>0.015</td><td>x</td><td>0.661</td><td>0.834</td><td>0.322</td><td>0.895</td><td>0.772</td><td>0.222</td></tr><tr><td>16</td><td>0.6</td><td>0.015</td><td>∇</td><td>0.674</td><td>1.02</td><td>-0.008</td><td>0.917</td><td>0.748</td><td>0.166</td></tr><tr><td>17</td><td>0.8</td><td>0.015</td><td>Δ</td><td>0.707</td><td>1.21</td><td>-0.42</td><td>0.93</td><td>0.858</td><td>-0.135</td></tr><tr><td>18</td><td>1</td><td>0.015</td><td>*</td><td>0.741</td><td>1.39</td><td>-0.881</td><td>0.947</td><td>1.021</td><td>-0.571</td></tr><tr><td>19</td><td>0</td><td>0.02</td><td>o</td><td>0.905</td><td>0.669</td><td>0.555</td><td>1.35</td><td>1.191</td><td>-0.369</td></tr><tr><td>20</td><td>0.25</td><td>0.02</td><td>+</td><td>0.957</td><td>0.585</td><td>0.660</td><td>1.44</td><td>0.931</td><td>0.0457</td></tr><tr><td>21</td><td>0.4</td><td>0.02</td><td>x</td><td>0.986</td><td>0.61</td><td>0.632</td><td>1.05</td><td>0.841</td><td>0.116</td></tr><tr><td><b>22</b></td><td><b>0.6</b></td><td><b>0.02</b></td><td>∇</td><td>0.947</td><td>0.711</td><td>0.504</td><td>1.07</td><td>0.779</td><td>0.142</td></tr><tr><td>23</td><td>0.8</td><td>0.02</td><td>Δ</td><td>0.974</td><td>0.835</td><td>0.318</td><td>1.09</td><td>0.849</td><td>-0.097</td></tr><tr><td>24</td><td>1</td><td>0.02</td><td>*</td><td>1.04</td><td>1.01</td><td>0.030</td><td>1.12</td><td>0.958</td><td>-0.403</td></tr><tr><td>25</td><td>0</td><td>0.025</td><td>o</td><td>1.17</td><td>1.03</td><td>-0.038</td><td>1.52</td><td>1.27</td><td>-0.548</td></tr><tr><td>26</td><td>0.25</td><td>0.025</td><td>+</td><td>1.19</td><td>0.874</td><td>0.238</td><td>1.63</td><td>1.01</td><td>-0.099</td></tr><tr><td>27</td><td>0.4</td><td>0.025</td><td>x</td><td>1.21</td><td>0.804</td><td>0.355</td><td>1.72</td><td>0.910</td><td>0.008</td></tr><tr><td>28</td><td>0.6</td><td>0.025</td><td>∇</td><td>1.25</td><td>0.758</td><td>0.427</td><td>1.22</td><td>0.813</td><td>0.102</td></tr><tr><td>29</td><td>0.8</td><td>0.025</td><td>Δ</td><td>1.29</td><td>0.756</td><td>0.437</td><td>1.24</td><td>0.821</td><td>-0.008</td></tr><tr><td>30</td><td>1</td><td>0.025</td><td>*</td><td>1.35</td><td>0.819</td><td>0.351</td><td>1.27</td><td>0.922</td><td>-0.312</td></tr></tbody></table>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Distribution of the bar height evolution versus bar migration (between final and initial bed profiles) for LIP1b and LIP1c, and a range of erosion rate $E_0 \in [0; 0.025] \text{ kg.m}^{-1}.\text{s}^{-2}$ and bedload factor $f_{bld} \in [0; 1]$.</figcaption></figure>
## 6.2. Contribution of bedload and suspended load to total sand transport

In this section, a sediment transport budget allows us to separate the transport analysis of the bedload (by the adapted SANTOSS formulation) and the suspended load transport (by the flow model) to improve our understanding of the simulated sandbar migration. Figure 11 illustrates the contribution of each mechanism to total sand transport ("transport" here means "flux divergence") over the sandbar for both erosion and accretion cases. In this figure, sediment is removed/accumulated in regions with positive/negative signs. In both cases, bedload transport is negative on the offshore side of the bar and positive in the lee side. Conversely, suspended-load divergence is positive on the offshore side of the bar and negative on lee side (opposite signs to bedload transport around the sandbar). As expected, the total transport follows the same pattern as the suspended load in 1b and the bedload in 1c, as in each case there is a dominant mechanism. Both transport mechanisms are stronger in 1b than in 1c. Note that these patterns compare well with a similar analysis performed with a wave-resolving model (Marchesiello et al., 2021b), although for 1c, the present wave-averaged model gives lower values for suspended load and higher for bedload. The reason may be due to differences in the bottom shear stress in the two models. Concerning the erosion case, 1b, the morphological evolution results from the balance between bedload and suspended load which have nearly the same order of magnitude. On the bar crest suspended load contribution is positive while bedload is negative with slightly smaller value leading to a net accumulation of sand. In the inner bar, suspended load contribution is negative while bedload is positive with slightly smaller value leading to a net erosion of sand. Overall, the morphodynamic evolution during erosion case is controlled by suspended load with bedload playing a secondary role by opposing the effect of suspended load. In the accretion case, suspended load has almost no effect. The morphodynamic evolution is almost entirely controlled by bedload.
<figure><img src="image_15.png" /><figcaption>Figure 11: Contribution of bedload and suspended-load divergence to total sediment transport over the outer bar.</figcaption></figure>
### 6.3. Near-bed wave skewness and asymmetry

In Section 6.2, we discussed the contribution of bedload in the sediment budget. In this section, the cross-shore profile of near-bed wave nonlinearity, which is the main driver of bedload transport is explained (Figure 12). In this figure, both normalised and absolute values of wave skewness and asymmetry are illustrated (for both cases of 1b and 1c). In our model, normalised plots are calculated using the formulations proposed by Ruessink et al. (2012) and the absolute values are obtained by multiplying the normalised values by local root-mean-squared velocity $U_{rms}$. The absolute values therefore show larger differences between 1b and 1c than the normalized values.

Ruessink et al. (2012) used the Ursell number (Doering and Bowen, 1995) to describe the variability of wave skewness and asymmetry as a function of wave length, wave height, and water depth ($Ur \propto H\lambda^2h^{-3}$). According to Ruessink et al. (2012), asymmetry and skewness reach their maximum at about $Ur=2$ and 10, respectively (by definition, the negative values of asymmetry indicate that the waves pitch forward). Therefore, when the Ursell number exceeds the value of 2, the skewness decreases while the asymmetry continues to increase. This is consistent with our results, as shown in Figure 12 (the right axes of the bottom panels show the cross-shore profile of Ursell number).

The waves propagating from deep water begin to become skewed (starting from $x = 20~m$ for both cases) before they become asymmetric (starting from $x = 40~m$ for both cases). However, in the shoaling zone, cross-shore changes of absolute skewness is smaller than those of asymmetry. Over the bar, skewness slightly decreases while asymmetry reaches its maximum with a relatively sharp change. This emphasizes the effective role of asymmetry on bedload as the cross-shore positions of maximum bed evolution and asymmetry peak coincide. In addition, it suggests that the magnitude of absolute skewness and asymmetry in bedload transport plays a stronger role than their cross-shore distribution (the same behavior is observed in the wave-resolving model Marchesiello et al. (2021b)).
<figure><img src="image_16.png" /><figcaption>Figure 12: Cross-shore wave skewness and asymmetry profiles of erosion (top panel) and accretion (bottom panel) cases.</figcaption></figure>
6.4. *Comparison of full and adapted SANTOSS formula*

In their publication, van der A et al. (2013) state that the SANTOSS formula is built to model wave-
and current-driven sand transport in the wave boundary layer, including both bedload and suspended load,
i.e. the total transport. For non-breaking waves and relatively weak mean currents, sand transport can be
assumed to occur entirely within the wave boundary layer and the transport formula therefore describes
the total transport rate for these conditions. However, for stronger conditions, if there is a significant
amount of sediments in suspension above the sheet flow layer or the bedload layer, a separate suspended
load model based on an advection-diffusion equation for the mean concentration evolution C(x, z, t) needs
to be solved (e.g. Soulsby, 1997b). To avoid double counting the near-bed sand transport due to mean
currents, van der A et al. (2013) recommend canceling the sediment transport by the suspended load
model in the wave boundary layer. In the present work, we take a slightly different approach in which
the suspended load model is active even in the near-bed region (including a representation of bottom
streaming). With the typical vertical resolution we use (Δz<sub>b</sub> ≈ 3cm above the bar crest), current-related
transport is assumed to be essentially resolved by the suspension model. Under typical flow conditions
above the bar crest (LIP1b), the undertow generates a typical Shields number θ ≈ 0.8 and suspension
number W<sub>s</sub>/u<sub>∗</sub> ≈ 0.4 (u<sub>∗</sub> is the bed friction velocity). For such flow conditions, current-driven bedload
only represents about 5% of the total current-driven transport Roberts et al. (2003), which justifies our
modeling approach. In LIP1c, wave-related transport becomes dominant and even a higher percentage
of bedload transport by currents remains a secondary process. To verify our model and justify our
hypothesis, a sensitivity analysis to the sediment transport model parameters with the full SANTOSS
formula is presented below.

Figure 13(b) illustrates the behavior of the full SANTOSS formulation for predicting sandbar migration in the accretion case (1c) using the same calibration parameters as with the adapted formulation proposed here. In this case, the full formulation predicts a seaward migration of the bar of about 5 m instead of a the observed shoreward migration of 3 m. It also stretches the terrace offshore contrary to observation. Even without suspended load (i.e., E₀ = 0), the full SANTOSS formula is not able to reproduce the onshore bar migration. In the absence of suspended load, the bedload transport by currents in SANTOSS is the only mechanism responsible for offshore bar migration. Contrary to the flume experiment, the onshore transport due to wave skewness and asymmetry is thus much weaker in the model than the current-related offshore transport (compare the green and magenta lines in the Figure 13(b)). Note that it may be possible to achieve onshore migration in 1c by reducing the contribution of currents in the full SANTOSS (using a scaling factor "alphac" implemented in CROCO, as in COAWST, and not discussed in this paper) thereby minimizing the undertow effect. However, perhaps as long as the near-bed resolution is sufficient, it may be more realistic – and less reliant on another unknown parameter – to separate the undertow and wave contributions in a wave-averaged depth-resolving model, as suggested in this paper.

In the erosion case 1b (Figure 13(a)), maintaining the same tuning parameters as those obtained
for the adapted SANTOSS formula, the full formula leads to an overestimation of the offshore sandbar migration. A more accurate prediction of offshore transport can only be made when the suspension model is disabled (i.e., $E_0 = 0$). This result confirms our hypothesis that using current-related sediment transport in SANTOSS duplicates the contribution of the undertow since in our simulation the suspended load is resolved. Incidentally, similar results can be obtained in 1b by adjusting the erosion rate in the suspended load model and turning off bedload transport, reflecting the dominance of suspended load transport in this case.

The results presented in this section illustrates one of the difficulties in deriving a sand-transport formula for depth-resolving wave-averaged models: accounting for waves and currents in two separate models. This difficulty is less noticeable in depth-averaged models where the suspended load is strongly parameterized, but at the cost of reduced genericity, i.e., a higher need for calibration. In depth-resolved models, most of the suspended load is resolved and, as explained above, at least under the flow conditions studied here, this is the dominant contribution from current-driven transport. However, in the general case, our current model may underrepresent some of the near-bed transport by the currents. Further work is needed to address this issue and this is left as a prospect.
<figure><img src="image_17.png" /><figcaption>Figure 13: Failure of the full SANTOSS formulations to predict onshore sandbar migration. The values in the legend correspond to (bedload coefficient, ERATE).</figcaption></figure>
6.5. Further model verification

Finally, the model is tested against the SANDS experiments (Eichentopf et al., 2018) without further calibration. The results for both erosion and accretion cases are presented in Figure 14. In the top panel, the model correctly predicts the dominant offshore direction of sandbar migration and the cross-shore position of the sand pickup and amalgamation. In this case the wave height and period are (H_s, T_p)=(0.53 m, 4.14 s). In the accretion case, (H_s, T_p)=(0.32 m, 5.44 s), which results in onshore sand transport as illustrated in the lower panel. Here again, the model reproduces the observed accretion in a realistic way.

<figure><img src="image_18.png" /><figcaption>Figure 14: validation of the model using the SANDS experiments</figcaption></figure>

In summary, the present simulations describe a sandbar that is driven onshore by wave-related bedload transport and offshore by current-related suspended-load transport. Calibrating these two processes separately makes it possible to show that in the erosive case, the model is robust to changes in the bedload factor (provided that the erosion rate coefficient is within a realistic range of values). Similarly, in the accretion phase, the cross-shore position of the sandbar is well predicted and not overly sensitive to the bedload factor in a given range (e.g., from 0.6 to 1). This suggests that the model captures enough mechanisms not to strongly rely on the calibration parameters. Such robust model behavior in predicting bar migration is essential for predicting beach morphodynamics on time scales of days to
weeks, including upstate and downstate sequences. This is of course beyond the scope of this study, but
the results presented here are a first step on which to build future research.

7. Conclusion

In this paper we applied a depth-resolving wave-averaged morphodynamic model to the sandbar migration in the nearshore. The novelty of the approach lies in the usage of the half-wave cycle concept sand transport formula, namely SANTOSS (van der A et al., 2013), to predict on/offshore sandbar migration depending on the wave conditions, *i.e* accretion/erosion cases.

We studied separately the effects of the main transport drivers: undertow, wave nonlinearity and phase lag between mobilization and transport, which allowed us to better understand the respective contributions of each mechanism. In the present model, we have separated the current-driven suspended-load and the wave-driven near-bed load contributions to the total sand transport. This has improved the understanding of their functioning leading to either onshore or offshore sandbar migration without tuning the empirical parameters differently for the two cases. In the erosive case (1b), both mechanisms exist, but suspended load transport by the current is dominant. In contrast, in the accretion case (1c), near-bed transport by waves is dominant. Isolating the effects of velocity skewness/asymmetry on onshore transport allows us to observe the coincident position of maximum bed evolution and a distinct peak in wave asymmetry. This suggests that the magnitude of the wave nonlinearity, rather than its distribution, is the dominant factor. On the other hand, the model shows negligible sensitivity to phase lag between mobilization and transport at different wave half-cycles. This was expected since, in this study, the beach is made of medium rather than fine sand for which phase-lag effects at the laboratory scale are not expected. Additionally, wave streaming slightly increases the onshore migration of the bar in both cases (1b and 1c) and generates slightly taller bar crests.

All of these results were obtained by disabling the contribution of currents to bedload in the wave
boundary layer, a limitation that needs to be better evaluated and understood. Nevertheless, our present
adaptation of the SANTOSS formulation for 3D wave-averaged models allows, for the first time, to predict
the onshore/offshore sandbar migration under storm and post-storm conditions using the same set of
parameters. As a major perspective of this work, the model should now be applied to the morphological
evolution of natural 3D beaches (i.e. adding longshore variations) to validate the good behavior observed
for the cross-shore migration of sandbars.

Acknowledgements

This work was performed under the auspices of the DGA by Shom and Univ. Grenoble Alpes, CNRS,
Grenoble INP<fnref n="1" />, LEGI, 38000 Grenoble, France under the contract MEPELS and Protevs II. It also

<footnote n="1">¹ Institute of Engineering Univ. Grenoble Alpes</footnote>
730

received support from a consortium of French research agencies, as part of CROCOs development project (Insu GdR-2014 named CROCO). We would like to thank Dano Roelvink for sharing the LIP data. Apart from these, all data were acquired by the authors and the CROCO source code is freely available at www.croco-ocean.org. Both laboratory and modeling data are available upon request.

735

Over the last two decades, an important research effort has been supported by the Shom, financed by DGA, through the ECORS (2007-2012) and MEPELS (2017-2024) projects in order to better understand and model sandy beach morphodynamics.
Appendix A. Nomenclature

<table><thead><tr><th>Symbol</th><th>Description</th><th>Dimensions</th></tr></thead><tbody><tr><td>â</td><td>representative orbital excursion amplitude</td><td>m</td></tr><tr><td>A</td><td>wave action density</td><td>kg m<sup>2</sup>/s</td></tr><tr><td>A<sup>r</sup></td><td>roller action density</td><td>kg m<sup>2</sup>/s</td></tr><tr><td>B<sub>b</sub></td><td>empirical wave breaking parameter for the type of breaker</td><td>—</td></tr><tr><td>B<sup>b</sup></td><td>body force due to breaking and roller accelerations</td><td>m/s<sup>2</sup></td></tr><tr><td>BSS</td><td>Brier Skill Score</td><td>—</td></tr><tr><td>B<sup>wd</sup></td><td>body force due to wave bottom drag</td><td>m/s<sup>2</sup></td></tr><tr><td>c</td><td>computed value</td><td>variable</td></tr><tr><td>c<sub>g,x</sub></td><td>cross-shore wave group velocity</td><td>m/s</td></tr><tr><td>c<sub>x</sub></td><td>cross-shore phase velocity of the primary wave</td><td>m/s</td></tr><tr><td>C</td><td>sediment Reynolds-averaged, wave-averaged concentration</td><td>kg/m<sup>3</sup></td></tr><tr><td>d<sub>50</sub></td><td>median grain diameter</td><td>m</td></tr><tr><td>D</td><td>wave-averaged thickness of the water column</td><td>m</td></tr><tr><td>E</td><td>depth-integrated wave energy</td><td>kg m<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>E<sub>s</sub></td><td>erosion mass flow rate</td><td>kg/m<sup>2</sup>/s</td></tr><tr><td>E<sub>0</sub></td><td>empirical erosion rate coefficient</td><td>kg/m<sup>2</sup>/s</td></tr><tr><td>E<sup>r</sup></td><td>roller energy density</td><td>kg m<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>f<sub>bld</sub></td><td>bedload coefficient</td><td>—</td></tr><tr><td>f<sub>mor</sub></td><td>morphological acceleration factor</td><td>—</td></tr><tr><td>f<sub>p</sub></td><td>spectral wave peak frequency</td><td>1/s</td></tr><tr><td>f<sub>w</sub></td><td>wave-friction factor</td><td>—</td></tr><tr><td>f<sub>wi</sub></td><td>wave-friction factor for either the crest or trough half-cycles</td><td>—</td></tr><tr><td>f<sub>δ</sub></td><td>current-related friction factor</td><td>—</td></tr><tr><td>f<sub>wδi</sub></td><td>linear wave-current friction factor at crest or trough</td><td>—</td></tr><tr><td>F<sup>w</sup></td><td>sum of wave-induced acceleration processes</td><td>m/s<sup>2</sup></td></tr><tr><td>F<sup>b</sup>(z)</td><td>vertical distribution function of breaking waves and rollers</td><td>—</td></tr><tr><td>F<sup>wd</sup>(z)</td><td>vertical distribution function of turbulent wave bottom boundary layer</td><td>—</td></tr><tr><td>g</td><td>gravitational acceleration</td><td>m/s<sup>2</sup></td></tr><tr><td>h</td><td>water depth</td><td>m</td></tr><tr><td>H</td><td>wave height</td><td>m</td></tr></tbody><tfoot><tr><td colspan="3">Continued on next page</td></tr></tfoot></table>
continued from previous page

<table><thead><tr><th>Symbol</th><th>Description</th><th>Dimensions</th></tr></thead><tbody><tr><td>H<sub>s</sub></td><td>significant wave height</td><td>m</td></tr><tr><td>H<sub>rms</sub></td><td>root mean square wave height</td><td>m</td></tr><tr><td>H<sub>*</sub></td><td>local wave height</td><td>m</td></tr><tr><td>i</td><td>subscript representing either 'c for crest' or 't' for trough</td><td>—</td></tr><tr><td>IV</td><td>initial value</td><td>variable</td></tr><tr><td>J</td><td>horizontal component of vortex force</td><td>m/s<sup>2</sup></td></tr><tr><td>k/TKE</td><td>turbulent kinetic energy</td><td>m<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>k<sub>sw</sub></td><td>bed roughness height for waves</td><td>m</td></tr><tr><td>k<sub>sδ</sub></td><td>bed roughness height for current</td><td>m</td></tr><tr><td>k<sub>x</sub></td><td>cross-shore local wave number</td><td>1/m</td></tr><tr><td>K</td><td>vertical component of vortex force</td><td>m/s<sup>2</sup></td></tr><tr><td>κ</td><td>Bernoulli head</td><td>m<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>m</td><td>measured value</td><td>variable</td></tr><tr><td>m<sub>0</sub></td><td>initial measured bed profile</td><td>m</td></tr><tr><td>m<sub>f</sub></td><td>final measured bed profile</td><td>m</td></tr><tr><td>nRMSE</td><td>normalised root-mean-square error</td><td>—</td></tr><tr><td>p</td><td>bed porosity</td><td>—</td></tr><tr><td>q<sub>b</sub></td><td>cross-shore sediment flux without avalanche module</td><td>kg/(m<sup>2</sup>s)</td></tr><tr><td>q<sub>b,slope</sub></td><td>cross-shore sediment flux with avalanche module</td><td>kg/(m<sup>2</sup>s)</td></tr><tr><td>RMSE</td><td>root-mean-square error</td><td>—</td></tr><tr><td>s</td><td>sand/water density ratio</td><td>—</td></tr><tr><td>t</td><td>time</td><td>s</td></tr><tr><td>T</td><td>wave period</td><td>s</td></tr><tr><td>T<sub>c</sub></td><td>duration of the crest half cycle</td><td>s</td></tr><tr><td>T<sub>cu</sub></td><td>duration of accelerating flow within the crest half cycle</td><td>s</td></tr><tr><td>T<sub>iu</sub></td><td>duration of accelerating flow within either the crest or trough half cycle</td><td>s</td></tr><tr><td>T<sub>tu</sub></td><td>duration of accelerating flow within the trough half cycle</td><td>s</td></tr><tr><td>T<sub>i</sub></td><td>duration of either the crest or trough half cycle</td><td>s</td></tr><tr><td>T<sub>p</sub></td><td>peak period</td><td>s</td></tr><tr><td>T<sub>t</sub></td><td>duration of the trough half cycle</td><td>s</td></tr><tr><td>T<sub>tu</sub></td><td>duration of accelerating flow within the trough half cycle</td><td>s</td></tr></tbody></table>

Continued on next page
continued from previous page

<table><thead><tr><th>Symbol</th><th>Description</th><th>Dimensions</th></tr></thead><tbody><tr><td>u</td><td>Eulerian cross-shore velocity</td><td>m/s</td></tr><tr><td>u<sub>i,r</sub></td><td>cross-shore combined wave-current velocity at trough or crest half cycles</td><td>m/s</td></tr><tr><td>u<sub>orb</sub></td><td>bottom orbital velocity</td><td>m/s</td></tr><tr><td>u<sub>rms</sub></td><td>root-mean-squared velocity</td><td>m/s</td></tr><tr><td>u<sub>w</sub>(t)</td><td>time-varying free-stream orbital velocity</td><td>m/s</td></tr><tr><td>u<sub>δ</sub></td><td>steady current velocity at the top of wave boundary layer</td><td>m/s</td></tr><tr><td>û<sub>i</sub></td><td>peak crest or trough orbital velocities</td><td>m/s</td></tr><tr><td>û</td><td>representative orbital velocity amplitude for the whole flow cycle</td><td>m/s</td></tr><tr><td>u<sub>∞</sub><sup>f</sup></td><td>free-stream velocity</td><td>m/s</td></tr><tr><td>u<sup>st</sup></td><td>cross-shore stokes velocity</td><td>m/s</td></tr><tr><td>U<sub>m</sub><sup>f</sup></td><td>maximum free-stream velocity</td><td>m/s</td></tr><tr><td>Ur</td><td>Ursell number</td><td>—</td></tr><tr><td>u<sub>*</sub></td><td>bed friction velocity</td><td>m/s</td></tr><tr><td>Var</td><td>Variance</td><td>—</td></tr><tr><td>w</td><td>Eulerian vertical velocity</td><td>m/s</td></tr><tr><td>w<sup>st</sup></td><td>vertical stokes velocity</td><td>m/s</td></tr><tr><td>W<sub>s</sub></td><td>sand settling velocity</td><td>m/s</td></tr><tr><td>x</td><td>cross-shore direction</td><td>m</td></tr><tr><td>x<sub>maxΔz</sub><sup>f</sup></td><td>cross-shore position of the maximum bed evolution at the end of the simulation</td><td>m</td></tr><tr><td>x<sub>maxΔz</sub><sup>i</sup></td><td>ross-shore position of the maximum bed evolution at the start of the simulation</td><td>m</td></tr><tr><td>z</td><td>vertical direction</td><td>m</td></tr><tr><td>z<sub>0</sub></td><td>empirically determined bedload roughness height</td><td>m</td></tr><tr><td>z<sub>0b</sub></td><td>bottom roughness length</td><td>m</td></tr><tr><td>z<sub>0s</sub></td><td>grain roughness length</td><td>m</td></tr><tr><td>z<sub>b</sub></td><td>bed height</td><td>m</td></tr><tr><td>z<sub>m</sub></td><td>height of the first grid point above the bed</td><td>m</td></tr><tr><td>z<sub>b</sub><sup>f</sup></td><td>bed level at the end of the simulation</td><td>m</td></tr><tr><td>z<sub>b</sub><sup>i</sup></td><td>bed level at the start of the simulation</td><td>m</td></tr><tr><td>α<sub>r</sub></td><td>fraction of wave energy converted into roller energy</td><td>—</td></tr></tbody></table>

Continued on next page
continued from previous page

<table><thead><tr><th>Symbol</th><th>Description</th><th>Dimensions</th></tr></thead><tbody><tr><td>β</td><td>local bottom slope</td><td>degrees</td></tr><tr><td>γ<sub>b</sub></td><td>ratio of wave height to water depth (empirical wave breaking parameter)</td><td>—</td></tr><tr><td>Δx</td><td>cross-shore grid size</td><td>m</td></tr><tr><td>δ</td><td>wave boundary layer thickness</td><td>m</td></tr><tr><td>ε<sup>b</sup></td><td>wave dissipation due to wave breaking</td><td>kg s<sup>3</sup>/m<sup>2</sup></td></tr><tr><td>ε<sup>r</sup></td><td>wave dissipation due to wave conversion into roller</td><td>kg s<sup>3</sup>/m<sup>2</sup></td></tr><tr><td>ε<sup>wd</sup></td><td>wave dissipation due to bottom viscous drag</td><td>kg s<sup>3</sup>/m<sup>2</sup></td></tr><tr><td>η and<br>eta</td><td>wave set-up</td><td>m</td></tr><tr><td>θ</td><td>Shields number</td><td>—</td></tr><tr><td>θ<sub>c</sub></td><td>Shields number associated with the wave crest half cycles</td><td>—</td></tr><tr><td>θ<sub>cr</sub></td><td>critical Shields number</td><td>—</td></tr><tr><td>θ<sub>i</sub></td><td>Shields number associated with either the wave crest or trough half cycles</td><td>—</td></tr><tr><td>θ<sub>t</sub></td><td>cross-shore Shields number associated with the wave trough half cycle</td><td>—</td></tr><tr><td>κ</td><td>Von Karman constant</td><td>—</td></tr><tr><td>λ</td><td>wavelength</td><td>m</td></tr><tr><td>ν<sup>t</sup></td><td>turbulent viscosity that includes wave effects</td><td>m<sup>2</sup>/s</td></tr><tr><td>ν<sub>c</sub><sup>t</sup></td><td>turbulent tracer diffusivity including wave effects</td><td>m<sup>2</sup>/s</td></tr><tr><td>ρ</td><td>total density of water</td><td>kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>mean density of water</td><td>kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>s</sub></td><td>sand density</td><td>kg/m<sup>3</sup></td></tr><tr><td>σ</td><td>standard deviation</td><td>—</td></tr><tr><td>σ<sub>p</sub></td><td>spectral wave peak angular frequency</td><td>rad/s</td></tr><tr><td>τ<sub>c</sub></td><td>current-only skin stress</td><td>N/m<sup>2</sup></td></tr><tr><td>τ<sub>CE</sub></td><td>critical bed shear stress for erosion</td><td>N/m<sup>2</sup></td></tr><tr><td>τ<sub>sf</sub></td><td>combined wave-current bottom friction stress including the ripples</td><td>N/m<sup>2</sup></td></tr><tr><td>τ<sub>w</sub></td><td>wave-only skin stress</td><td>N/m<sup>2</sup></td></tr><tr><td>τ<sub>bot</sub><sup>cd</sup></td><td>combined wave-current bottom friction stress</td><td>N/m<sup>2</sup></td></tr><tr><td>φ</td><td>normalized dynamic pressure</td><td>m<sup>2</sup>/s<sup>2</sup></td></tr></tbody><tfoot><tr><td colspan="3">Continued on next page</td></tr></tfoot></table>
continued from previous page

<table><thead><tr><td>Symbol</td><td>Description</td><td>Dimensions</td></tr></thead><tbody><tr><td>Φ</td><td>non-dimensional net transport rate</td><td>—</td></tr><tr><td>ω</td><td>turbulent dissipation rate</td><td>kg s<sup>3</sup>/m<sup>2</sup></td></tr><tr><td>Ω<sub>cc</sub></td><td>sand load entrained during the wave crest period and transported during the wave crest period</td><td>m<sup>2</sup>/s</td></tr><tr><td>Ω</td><td>Dean number</td><td>m<sup>2</sup>/s</td></tr><tr><td>Ω<sub>ct</sub></td><td>cross-shore sand load entrained during the wave crest period and transported during the wave trough period</td><td>m<sup>2</sup>/s</td></tr><tr><td>Ω<sub>i</sub></td><td>cross-shore sand load entrained during either the wave trough or crest period</td><td>m<sup>2</sup>/s</td></tr><tr><td>Ω<sub>tc</sub></td><td>cross-shore sand load entrained during the wave trough period and transported during the wave crest period</td><td>m<sup>2</sup>/s</td></tr><tr><td>Ω<sub>tt</sub></td><td>cross-shore sand load entrained during the wave trough period and transported during the wave trough period</td><td>m<sup>2</sup>/s</td></tr></tbody></table>
Appendix B. Detailed results of the sensitivity analysis to hydrodynamic parameters

<figure><img src="image_19.png" /><figcaption>Figure B.1: Model sensitivity to the breaking parameter $B_b^{52}$ for experiments LIP1b and LIP1c, showing the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure B.2: Model sensitivity to the onshore bottom roughness $Z_{0b}$ for experiments LIP1b and LIP1c, showing the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure B.3: Model effects of wave streaming on the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow for experiments LIP1b and LIP1c.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure B.4: Model sensitivity to the roller fraction parameter for experiments LIP1b and LIP1c, showing the cross-shore current profiles, wave height $H_{rms}$, wave set-up (eta), and undertow.</figcaption></figure>
Appendix C. Detailed results of the sensitivity analysis to morphodynamic parameters

<figure><img src="image_23.png" /><figcaption>Figure C.1: Model sensitivity to the bedload factor $f_{bld}$ with erosion rate $E_0=0.02 \text{ kg/m}^2/\text{s}$ in LIP1b (a) and LIP1c (b) (same panels as Figure 6)</figcaption></figure>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

Figure C.2: Model sensitivity to erosion rate $E_0$ with bedload factor $f_{bld}=0.8$ in LIP1b (a) and LIP1c (b) (same panels as Figure 6)
740 Appendix D. sandbar migration prediction using Soulsby and Damgaard (2005) formulations

<figure><img src="image_26.png" /><figcaption>Figure D.1: Illustration of the overestimation of offshore sand transport using Soulsby and Damgaard (2005) formulations</figcaption></figure>
Amoudry, L., 2008. A review on coastal sediment transport modelling.

Bennis, A.-C., Ardhuin, F., Dumas, F., 2011. On the coupling of wave and three-dimensional circulation models: Choice of theoretical framework, practical implementation and adiabatic tests. Ocean Modelling 40 (3-4), 260-272.

745

Blaas, M., Dong, C., Marchesiello, P., McWilliams, J. C., Stolzenbach, K. D., 2007. Sediment-transport modeling on southern californian shelves: A roms case study. Continental shelf research 27 (6), 832-853.

Brier, G. W., et al., 1950. Verification of forecasts expressed in terms of probability. Monthly weather review 78 (1), 1-3.

Camenen, B., Larson, M., 2005. A general formula for non-cohesive bed load sediment transport. Estuarine, Coastal and Shelf Science 63 (1), 249-260.

750

Church, J. C., Thornton, E. B., 1993. Effects of breaking wave induced turbulence within a longshore current model. Coastal Engineering 20 (1), 1-28.
URL http://www.sciencedirect.com/science/article/pii/037838399390053B

Dally, W. R., Dean, R. G., 1984. Suspended sediment transport and beach profile evolution. Journal of waterway, port, coastal, and ocean engineering 110 (1), 15-33.

755

Dean, R., 2018. Physical modelling of littoral processes. In: Physical modelling in coastal engineering. Routledge, pp. 119-140.

Dean, R. G., 1973. Heuristic models of sand transport in the surf zone. In: Proc. 1st Australian Conf. Coastal Eng. pp. 209-214.

760

Debreu, L., Marchesiello, P., Penven, P., Cambon, G., 2012. Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean Modelling 49, 1-21.

Dibajnia, M., Watanabe, A., 1992a. Sheet flow under nonlinear waves and currents. In: Coastal Engineering 1992. pp. 2015-2028.

Dibajnia, M., Watanabe, A., Jan. 1992b. Sheet flow under nonlinear waves and currents. Coastal Engineering Proceedings 1 (23).

765

Doering, J., Bowen, A., 1995. Parametrization of orbital velocity asymmetries of shoaling and breaking waves using bispectral analysis. Coastal engineering 26 (1-2), 15-33.

Dohmen-Janssen, C., Kroekenstoel, D. F., Hassan, W. N., Ribberink, J. S., 2002. Phase lags in oscillatory sheet flow: experiments and bed load modelling. Coastal Engineering 46 (1), 61-87.

770

Dyhr-Nielsen, M., Sørensen, T., 1970. Some sand transport phenomena on coasts with bars. In: Coastal Engineering 1970. pp. 855-865.
Eichentopf, S., Cáceres, I., Alsina, J. M., 2018. Breaker bar morphodynamics under erosive and accretive wave conditions in large-scale experiments. Coastal Engineering 138, 36–48.

775

Elgar, S., 1987. Relationships involving third moments and bispectra of a harmonic process. IEEE transactions on acoustics, speech, and signal processing 35 (12), 1725–1726.

Elgar, S., Gallagher, E. L., Guza, R., 2001. Nearshore sandbar migration. Journal of Geophysical Research: Oceans 106 (C6), 11623–11627.

780

Exner, F., 1931. Zur dynamik der bewegungsformen auf der erdoberfläche. Ergebnisse der kosmischen Physik 1, 373–445.

Exner, F. M., 1925. Über die wechselwirkung zwischen wasser und geschiebe in flussen. Akad. Wiss. Wien Math. Naturwiss. Klasse 134 (2a), 165–204.

Fredsoe, J., Deigaard, R., 1992. Mechanics of costal sediment transport. word Scientific.

785

Hansen, J. B., Svendsen, I., 1985. A theoretical and experimental study of undertow. In: Coastal Engineering 1984. pp. 2246–2262.

Henderson, S. M., Allen, J. S., Newberger, P., 2004. Nearshore sandbar migration predicted by an eddy-diffusive boundary layer model. Journal of Geophysical Research: Oceans 109 (C6).

Hoefel, F., Elgar, S., 2003. Wave-induced sediment transport and sandbar migration. Science 299 (5614), 1885–1887.

790

Hsu, T.-J., Elgar, S., Guza, R., 2006. Wave-induced sediment transport and onshore sandbar migration. Coastal Engineering 53 (10), 817–824.

Johnson, D. W., 1919. Shore processes and shoreline development. John Wiley & Sons, Incorporated.

Kalra, T., Sherwood, C., Warner, J., Rafati, Y., Hsu, T.-J., 2019. Investigating bedload transport under asymmetrical waves using a coupled ocean-wave model. pp. 591–604.

795

Kumar, N., Voulgaris, G., Warner, J. C., Olabarrieta, M., 2012. Implementation of the vortex force formalism in the coupled ocean-atmosphere-wave-sediment transport (coawst) modeling system for inner shelf and surf zone applications. Ocean Modelling 47, 65–95.

Large, W. G., McWilliams, J. C., Doney, S. C., 1994. Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Reviews of geophysics 32 (4), 363–403.

800

Lesser, G., Roelvink, J., van Kester, J., Stelling, G., Oct. 2004a. Development and validation of a three-dimensional morphological model. Coastal Engineering 51 (8–9), 883–915.
Lesser, G. R., Roelvink, J. v., van Kester, J. T. M., Stelling, G., 2004b. Development and validation of a three-dimensional morphological model. Coastal engineering 51 (8-9), 883-915.

Longuet-Higgins, M. S., 1953. Mass transport in water waves. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 245 (903), 535-581.

Marchesiello, P., Auclair, F., Debreu, L., McWilliams, J., Almar, R., Benshila, R., Dumas, F., 2021a. Tridimensional nonhydrostatic transient rip currents in a wave-resolving model. Ocean Modelling 163, 101816.

Marchesiello, P., Benshila, R., Almar, R., Uchiyama, Y., McWilliams, J. C., Shchepetkin, A., 2015. On tridimensional rip current modeling. Ocean Modelling 96, 36-48, waves and coastal, regional and global processes.

Marchesiello, P., Chauchat, J., Shafiei, H., Almar, R., Benshila, R., Dumas, F., 2021b. 3d wave-resolving simulation of sandbar migration.

McWilliams, J. C., Restrepo, J. M., Lane, E. M., 2004. An asymptotic theory for the interaction of waves and currents in coastal waters. Journal of Fluid Mechanics 511, 135.

Meyer-Peter, E., Müller, R., 1948. Formulas for bed-load transport. In: IAHSR 2nd meeting, Stockholm, appendix 2. IAHR.

Rafati, Y., Hsu, T.-J., Elgar, S., Raubenheimer, B., Quataert, E., van Dongeren, A., 2021. Modeling the hydrodynamics and morphodynamics of sandbar migration events. Coastal Engineering 166, 103885.

Ribberink, J. S., Jul. 1998. Bed-load transport for steady flows and unsteady oscillatory flows. Coastal Engineering 34 (1-2), 59-82.

Roberts, J. D., Jepsen, R. A., James, S. C., 2003. Measurements of sediment erosion and transport with the adjustable shear stress erosion and transport flume. Journal of Hydraulic Engineering 129 (11), 862-871.

URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9429%282003%29129%3A11%28862%29

Roelvink, D., Reniers, A., Van Dongeren, A., De Vries, J. V. T., McCall, R., Lescinski, J., 2009. Modelling storm impacts on beaches, dunes and barrier islands. Coastal engineering 56 (11-12), 1133-1152.

Roelvink, J., Reniers, A., 1995. LIP 11D delta flume experiments: a dataset for profile model validation. WL/Delft Hydraulics.

Roelvink, J., Stive, M., 1989. Bar-generating cross-shore flow mechanisms on a beach. Journal of Geophysical Research: Oceans 94 (C4), 4785-4800.
Roelvink, J. A., Broker, I., 1993. Cross-shore profile models. Coastal Engineering 21 (1-3), 163-191.

Roelvink, J. A., Meijer, T. J. G. P., Houwman, K., Bakker, R., Spanhoff, R., 1995. field validation and application of a coastal profile model. In: ASCE (Ed.), Proceedings of Coastal Dynamics 95. pp. 818-828.

835

Ruessink, B., Ramaekers, G., Van Rijn, L., 2012. On the parameterization of the free-stream non-linear wave orbital motion in nearshore morphodynamic models. Coastal engineering 65, 56-63.

Shchepetkin, A. F., McWilliams, J. C., 2005. The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean modelling 9 (4), 347-404.

840

Smith, J. D., McLean, S. R., 1977. Spatially averaged flow over a wavy surface. Journal of Geophysical Research (196-1977) 82 (12), 1735-1746.

Soufflet, Y., Marchesiello, P., Lemarié, F., Jouanno, J., Capet, X., Debreu, L., Benshila, R., 2016. On effective resolution in ocean models. Ocean Modelling 98, 36-50.

Soulsby, R., 1995. Bed shear-stresses due to combined waves and currents. Advances in coastal morphodynamics.

845

Soulsby, R., 1997a. Dynamics of marine sands.

Soulsby, R., 1997b. Dynamics of Marine Sands: A Manual for Practical Applications. Thomas Telford, London., London.

Soulsby, R., Whitehouse, R., 1997. Threshold of sediment motion in coastal environments.

850

Soulsby, R. L., Damgaard, J. S., Aug. 2005. Bedload sediment transport in coastal waters. Coastal Engineering 52 (8), 673-689.

Spielmann, K., Astruc, D., Thual, O., 2004. Analysis of some key parametrizations in a beach profile morphodynamical model. Coastal Engineering 51 (10), 1021-1049.

Stive, M. J., 1987. A model for cross-shore sediment transport. In: Coastal Engineering 1986. pp. 1550-1564.

855

Stive, M. J., Battjes, J. A., 1985. A model for offshore sediment transport. In: Coastal Engineering 1984. pp. 1420-1436.

Stokes, G. G., 1880. On the theory of oscillatory waves. Transactions of the Cambridge philosophical society.

860

Sutherland, J., Peet, A., Soulsby, R., 2004. Evaluating the performance of morphological models. Coastal engineering 51 (8-9), 917-939.
Svendsen, I. A., 1984. Mass flux and undertow in a surf zone. Coastal Engineering 8 (4), 347–365.

Symonds, G., Bowen, A., 1984. Interactions of nearshore bars with incoming wave groups. Journal of Geophysical Research: Oceans 89 (C2), 1953–1959.

Ting, F. C., Kirby, J. T., 1994. Observation of undertow and turbulence in a laboratory surf zone. Coastal Engineering 24 (1-2), 51–80.

Tolman, H. L., et al., 2009. User manual and system documentation of wavewatch iii tm version 3.14. Technical note, MMAB Contribution 276, 220.

Uchiyama, Y., McWilliams, J. C., Shchepetkin, A. F., 2010. Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone. Ocean Modelling 34 (1-2), 16–35.

van Den Bremer, T., Breivik, 2018. Stokes drift. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376 (2111), 20170104.

van der A, D. A., Ribberink, J. S., van der Werf, J. J., O'Donoghue, T., Buijsrogge, R. H., Kranenburg, W. M., 2013. Practical sand transport formula for non-breaking waves and currents. Coastal Engineering 76, 26 – 42.

van der Werf, J., Veen, R., Ribberink, J., van der Zanden, J., 2015. Testing of the new santoss transport formula in the delft3d morphological modeling system. In: The Proceedings of the Coastal Sediments 2015. World Scientific.

van Rijn, L., 1993. Principles of sediment transport in rivers, estuaries and coastal seas. Aqua Publications, Amsterdam, The Netherlands.

van Rijn, L., Tonnon, P., Walstra, D., 2011. Numerical modelling of erosion and accretion of plane sloping beaches at different scales. Coastal Engineering 58 (7), 637–655.

van Rijn, L. C., Ribberink, J. S., Van Der Werf, J., Walstra, D. J., 2013. Coastal sediment dynamics: recent advances and future research needs. Journal of Hydraulic Research 51 (5), 475–493.

van Rijn, L. C., et al., 1993. Principles of sediment transport in rivers, estuaries and coastal seas. Vol. 1006. Aqua publications Amsterdam.

Veen, R., van der Werf, J., Ribberink, J., van der Zanden, J., 2014. Improving wave-driven cross-shore sand transport modelling in delft3d using the new santoss transport formula. In: NCK-Days 2014: Preparedness.

Walstra, D., Reniers, A., Ranasinghe, R., Roelvink, J., Ruessink, B., 2012. On bar growth and decay during interannual net offshore migration. Coastal Engineering 60, 190–200.
Warner, J., Sherwood, C., Arango, H., Signell, R., 2005. Performance of four turbulence closure methods implemented using a generic length scale method. Ocean Modelling 8, 81–113.

895

Warner, J. C., Armstrong, B., He, R., Zambon, J. B., 2010. Development of a coupled ocean–atmosphere–wave–sediment transport (coawst) modeling system. Ocean modelling 35 (3), 230–244.

Warner, J. C., Sherwood, C. R., Signell, R. P., Harris, C. K., Arango, H. G., 2008. Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Computers and geosciences 34 (10), 1284–1306.

900

Warren, I., Bach, H., 1992. Mike 21: a modelling system for estuaries, coastal waters and seas. Environmental Software 7 (4), 229–240.

Weir, B., Uchiyama, Y., Lane, E., Restrepo, J., McWilliams, J., 2011. A vortex force analysis of the interaction of rip currents and surface gravity waves. Journal of Geophysical Research: Oceans 116 (C5).

Xu, Z., Bowen, A., 1994. Wave-and wind-driven flow in water of finite depth. Journal of physical oceanography 24 (9), 1850–1866.

905

Zheng, J., Zhang, C., Demirbilek, Z., Lin, L., 2014. Numerical study of sandbar migration under wave-undertow interaction. Journal of Waterway, Port, Coastal, and Ocean Engineering 140 (2), 146–159.