Nonhydrostatic Numerical Modeling of Fixed and Mobile Barred Beaches: Limitations of Depth-Averaged Wave Resolving Models around Sandbars

: Along sandy coastlines, submerged, shore-parallel sandbars play an essential role in shoreline morphology by dissipating wave energy through depth-induced wave breaking. While wave breaking and sediment transport around sandbars are complex three-dimensional (3D) processes, shoreline morphology is typically simulated with depth-averaged models that feature lower computational demand than do 3D models. In this context, this study examines the implications of depth-averaging the ﬂ ow ﬁ eld and approximating the breaking process in nonhydrostatic models (e.g., XBeach nonhydrostatic) for the hydro-and morphodynamic processes around sandbars. The implications are drawn based on reproducing large-scale experiments of a barred beach pro ﬁ le using the single-layer (XBNH) and the reduced two-layer (XBNH + ) modes of XBeach. While hydrodynamic processes were predicted with high accuracy on the sandbar ’ s seaward side, wave heights were overpredicted on the bar ’ s landward side. The overestimation was due to the simpli ﬁ ed reproduction of the complex breaking process near the sandbar ’ s peak, particularly in terms of the generated turbulence in the water column. Moreover, the velocity pro ﬁ le with a strong undertow could only be represented in a simpli ﬁ ed way even using the two-layer mode XBNH + , thus resulting in inaccurate predictions of sediment loads around the sandbar. A parametric study is performed, and it revealed which model parameters control the simulation of the wave-breaking process. Thus, wave height predictions could be improved by tuning the energy-dissipation parameters. However, ﬂ ow velocities and morphodynamic predictions could not be improved accordingly. Thus, this study identi ﬁ es possible hydrodynamic model improvements, such as incorporating a roller dissipation model. Moreover, it improves understanding of key drivers and processes that should be included in nonhydrostatic depth-averaged models to simulate morphological changes around sandbars more ef ﬁ ciently. DOI: 10.1061/(ASCE)WW.1943-5460.0000685 . This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.


Introduction
Coastal zones are attractive locations for tourism, recreation, economic activities, and global trade exchange (Sardain et al. 2019).
However, global warming and sea-level rise stress many coastal ecosystems by intensifying coastal erosion and flooding events (Elsayed and Oumeraci 2016;Nguyen et al. 2020).Currently, 24% of the world's sandy beaches are eroding at rates exceeding 0.5 m/yr and 28% are accreting, while 48% are stable (Luijendijk et al. 2018).Submerged sandbars (Fig. 1) represent an important morphological feature that has implications for human activities in coastal zones (e.g., tourism and surfing) as well as for coastal protection against coastal erosion (Kuang et al. 2019;Poate 2011).Sandbars dissipate wave energy through depth-induced wave breaking and reduce wave period through nonlinear dispersive wave transformation in the surf zone (Cohn et al. 2014;Kuznetsova and Saprykina 2019;Mulligan et al. 2019).Therefore, they influence cross-shore sediment transport processes, affect the shoreline position, and mitigate coastal erosion (Eichentopf et al. morphological processes (e.g., sediment stirring, suspension, and transport) as well as subsequent morphological evolution (e.g., bed changes and bar migration).
Around sandbars, the time-averaged velocity varies over the depth due to the strong depth variations in momentum (i.e., wave shear stress).Thereby, flow velocity profiles can be simplified as shown in Fig. 1 by an onshore component near the surface and an offshore component (i.e., return flow or undertow) closer to the bed (van der A et al. 2017).These two velocity components represent two counteracting sediment stirring drivers that control bar formation and its net migration.The latter components are often generated as follows: 1.The net onshore mass flux (i.e., Stokes drift) is induced by wave orbital velocity skewness and asymmetry as well as by surface rollers (bores) induced by intense wave breaking over the bar.2. The offshore undertow then balances the net onshore mass flux (van der A et al. 2017;Dubarbier et al. 2015;Rafati et al. 2021;Svendsen 1984;Wenneker et al. 2011;Xie et al. 2017;van der Zanden 2016;van der Zanden et al. 2016van der Zanden et al. , 2017c, b;, b;Zheng et al. 2014).
The onshore sediment stirring driver dominates on the offshore side of the bar and favors shoreward bar migration (Kim et al. 2017;van der Zanden et al. 2017b), while the offshore driver dominates on the onshore side of the bar and leads to offshore bar migration (van der Zanden et al. 2016).The net migration of a sandbar in a specific time interval is often determined by the imbalance between the effects of these two drivers during the same interval (Hewageegana and Canestrelli 2021;Melito et al. 2020; Vidal-Ruiz and de Alegría-Arzaburu 2020; Yuhi et al. 2016).
Several (open-source) numerical models have been developed to assess nearshore hydrodynamics and associated sediment transport and bed morphology.Among these modeling tools is the processbased model XBeach (Roelvink et al. 2009).Due to its depth-averaged nature, two-dimensional (2D) XBeach is often more computationally efficient than are 3D models, particularly for large-domain and region-scale wave models (Gharagozlou et al. 2020).XBeach can be operated in two main hydrodynamic modes (Elsayed 2017;Roelvink et al. 2018): (1) a hydrostatic wave averaging mode that can either be run under stationary or transient (surfbeat) conditions (Roelvink et al. 2009), and (2) a nonhydrostatic (hereafter XBNH) mode that resolves wave-by-wave flows and surface elevation variations due to short waves (Smit et al. 2010).Both modes utilize the same morphodynamic module, while the hydrodynamics are represented differently.
In the surf-beat mode of XBeach (hereafter XBSB), infragravity (long) waves of a frequency range between 0.005 and 0.05 Hz are solved together with currents using the nonlinear shallow water equations (NLSWEs).Short waves of frequency >0.05 Hz are solved separately in XBSB using wave and roller-action balance equations, and then shortwaves' forces exerted on the water surface are incorporated into the former NLSWEs as a source term (Roelvink et al. 2009).The separation in solving long and short waves in the latter mode significantly reduces the computational time, but phases of shortwaves are not resolved.Concurrently imposed intrawave processes, expressing themselves through wave skewness, asymmetry, and breaking, are therefore not physically resolved but implicitly integrated in XBSB, using empirical formulations that may require parameterization for each modeling case (Cohn et al. 2014;Elsayed and Oumeraci 2017;Kalligeris et al. 2020;Rafati et al. 2021).Therefore, resolving shortwaves is essential for predicting hydrodynamics and net migration of sandbars as intrawave processes affect the aforementioned net onshore mass flux (Fernández-Mora et al. 2015;Hewageegana and Canestrelli 2021;Hoefel and Elgar 2003).
XBNH resolves shortwave dynamics, including wave skewness, asymmetry, and breaking.In addition, it predicts the breakpoint accurately, which is an advantage over XBSB and Boussinesq-type models (e.g., Cienfuegos et al. 2010;David et al. 2017; Kirby 2016; Fig. 1.Schematic representation of nearshore hydromorphodynamic processes, including submerged sandbars' role in dissipating wave energy through depth-induced wave breaking.Waves propagate to the shore through four subsequent zones: shoaling, breaking, inner surf, and swash zones.In the shoaling zone, wave shapes become highly skewed and asymmetric until they break over the bar producing bores (rollers) in the inner surf zone.Then in the swash zone, waves produce swash motions in the form of runup and rundown.Waves breaking over the bar are often accompanied by wave overturning, high turbulence, strong return flow (undertow), and high sediments suspension.The point in between onshore and offshore mass fluxes represents the inflection point.(Adapted from van der Zanden 2016.)Madsen et al. 1991;Nwogu 1993;Wenneker et al. 2011), which generally need a separate breaking model to initiate the breaking process.However, the simulation of intrawave processes in XBNH requires a significantly higher spatial resolution and associated smaller time step, thereby increasing the computational demand.Nevertheless, XBNH's computational demand is still much lower than that of full 3D numerical models, which, besides their unfeasibility for large areas and region-scale simulations, do not necessarily provide better predictions (Lashley et al. 2020).
To account for the interactions between hydrodynamic and morphologic processes in the surf zone in a phase and depth resolving fashion, Deltares has recently released a new version of XBeach (version 1.23.5446,known as XBeachX), which has the option of subdividing the water depth into two layers by means of a fraction α (=h down /h); see Fig. 2.This option (hereafter XBNH+) provides quasi-3D computations to represent, for instance, the undertow in the surf zone utilizing the reduced two-layer approach of Cui et al. (2014).Because of the vertical discretization in two layers, frequency dispersion is computed more accurately than using XBNH of a single layer (De Ridder 2018;De Ridder et al. 2020).
Although the reduced two-layer mode may better represent the onshore and offshore drivers, the approximations and assumptions in XBNH+ (Fig. 2) still impose unavoidable limitations and inaccuracies on the simulation of complex processes around the sandbar.For instance, to account for the wave-breaking process, XBNH+ (and also XBNH) utilizes an edge-based compact difference scheme to approximate the vertical gradient of the nonhydrostatic pressure.Thus, wave breaking is regarded as a subgrid process where propagating waves are allowed to steepen until the front face is almost vertical, as depicted in Fig. 2.This affects the computation of the detailed wave-breaking processes (e.g., spilling/overturning) and excludes relevant wave-roller dynamics.Hence, this approximation may have significant implications for the simulation of hydrodynamic and morphodynamic processes around submerged sandbars.Moreover, being based on the reduced two-layer concept, XBNH+ provides implicit computations of layer-averaged flow velocities (u up and u down , as depicted in Fig. 2) and does not calculate the full velocity variation over the water depth.Layer-averaged flow velocities may be calculated as a function of the depth-averaged velocity (u) and the velocity difference between the layers (dU) (Fig. 2).The sediment transport calculations are thereby still based on the depth-averaged velocity (u) as is the case in XBNH, which may lead to unrealistic predictions of the transport rates and evolving morphology, leading to inaccurate predictions of sandbar migration.It is currently unknown what the implications of these model approximations are for the simulation of wave propagation, energy dissipation, and sediment transport in the sandbar area, as well as their effects on the morphological evolution of the sandbar.Identifying these limitations is crucial for further model development.Moreover, it is essential for the broad community using XBNH (or any similar nonhydrostatic depth-averaged model) to consider such model limitations in model applications to real-world case studies so that uncertainty related to the model structure can be evaluated.From this perspective, this study aims to identify and quantify how the approximations in XBNH+ and XBNH affect the simulation of hydro-and morphodynamic processes around sandbars.
To date, a few studies have focused on the hydrodynamic processes simulated by XBNH+ (e.g., Klaver et al. 2019;Lashley et al. 2020;De Ridder 2018;De Ridder et al. 2020) and found a good prediction capability.However, a comprehensive study on the simulated hydrodynamic and morphodynamic processes around sandbars is still lacking.Moreover, the morphological processes in XBNH have not yet undergone extensive validation (Roelvink et al. 2018;Ruffini et al. 2020;Rutten et al. 2021).Only a few studies have applied XBNH to investigate nearshore morphodynamic processes (e.g., Daly et al. 2017;Jongedijk 2017;Kuznetsova and Saprykina 2019;Ruffini et al. 2020;Sanuy and Jiménez 2019;van Weeghel 2020).Among these studies, Kuznetsova and Saprykina (2019) was the only study to address barred beach modeling by applying XBNH to artificial beach profiles.The latter study Fig. 2. Schematic representation of approximated nearshore hydrodynamic processes in the reduced two-layer nonhydrostatic mode of XBeach.The water depth h is divided into two water layers of depths h up and h down for the upper and lower water layers, respectively.Intrawave processes (e.g., waves' skewness and asymmetry) are resolved.Wave breaking is resolved, but wave spilling/overturning and roller accompanying breaking are missing.The reduced two-layer option enables implicit calculations of the layers' velocities u up and u down as a function of the depth-averaged velocity (u) and the velocity difference between the layers (dU).
investigated the effect of a sandbar location on the dissipation of wave energy and the cross-shore transport of sediment but did not provide a validation of the model.The remaining studies mainly focused on the swash zone, showing conflicting findings.Daly et al. (2017) and Jongedijk (2017) showed that XBNH reasonably simulates beach erosion and accretion in the swash zone after fine-tuning to relevant sediment transport parameters (e.g., bed slope and dilatancy effects of De Vet 2014).However, Sanuy and Jiménez (2019), Ruffini et al. (2020), andvan Weeghel (2020) showed high accuracy of XBNH to simulate hydrodynamic swash-swash interactions and worse performance in capturing intrawaves sediment transport processes.Thus, to the authors' best knowledge, there is no comprehensive study that has validated the ability of XBNH or XBNH+ to simulate hydrodynamic and associated morphodynamic processes in the surf zone or at least examined the effect of the model approximations, mentioned previously, on hydro-and morphodynamic processes around sandbars.
In this context, the SINBAD data set (van der A et al. 2017;van der Zanden 2016;van der Zanden et al. 2016van der Zanden et al. , 2017avan der Zanden et al. , b, c, d, 2019) ) provides an opportunity to investigate roles of these outlined shortcomings thoroughly, based on detailed measurements of hydrodynamics and subsequent morphodynamics around a sandbar.This data set describes large-scale flume experiments in which hydrodynamic processes near a fixed sandbar (van der A et al. 2017) and related morphodynamic processes near a mobile sandbar (van der Zanden 2016;van der Zanden et al. 2016van der Zanden et al. , 2017b, c, d) , c, d) were investigated.Therefore, the detailed objectives of the present paper are: 1.To validate XBNH and XBNH+ with the SINBAD data set and examine their capabilities and limitations in simulating hydrodynamic and morphodynamics processes near sandbars; 2. To examine the capability of the quasi-3D computations in XBNH+ in reproducing the undertow and the associated sediment transport processes while finding the optimal layer fraction for such purposes; 3. To investigate the effect of the wave-breaking approximation and the omitted wave overturning, as well as the depth averaging of flow velocity, on hydrodynamic and morphological processes near sandbars; 4. To examine the parameters that control the simulation of wave breaking and relevant energy dissipation in XBNH and XBNH+; and 5. To identify missing processes in XBNH and XBNH+ necessary for proper modeling of hydromorphodynamic processes around sandbars.The remainder of this paper is organized as follows: first, a brief summary of the SINBAD data set, XBNH and XBNH+, and the numerical model setup for the data set are given in the section "Methods and Data Set."The outcomes of reproducing the data set using XBNH and XBNH+ are then presented, followed by a discussion of the findings and suggestions for further model developments.Finally, the conclusions of the study are outlined.

Methods and Data Set
Experimental SINBAD Data Set The SINBAD data set includes measurements of a fixed and a mobile barred bed experimental campaign performed in the Maritime Research and Experimentation Wave Flume (CIEM) at Universitat Politècnica de Catalunya (UPC), Barcelona.This large research facility is 100 m long, 3 m wide, and 4.5 m deep, with a wedge-type wave paddle of a maximum stroke of 2 m.The campaign was part of the SINBAD project (Sand Transport under Irregular and Breaking Waves Conditions) that aimed at providing new insights into the wave-breaking hydrodynamics (van der A et al. 2017) and morphodynamic evolution of a sandbar under energetic breaking waves (van der Zanden et al. 2016(van der Zanden et al. , 2017c)).Thus, this experimental campaign consisted of two similar experiments (Fig. 3) in which (1) a mobile medium-sand bed was installed (van der Zanden et al. 2016), and (2) the bed was thereafter fixed with a concrete layer to provide insights into the intrawave velocities and turbulence.The latter is difficult to achieve under mobile bed conditions because of practical complications associated with sediment suspension (van der A et al. 2017).Hereafter, these experiments are referred to as the "mobile-bed" and "fixed-bed" experiments.
As shown in Fig. 3, the sandbar had a 1:12 offshore slope, a 0.6 m elevation above the landward runnel, and a 1:7 landward slope.The bathymetry was followed by a 10 m-long 1:125 slope leading toward a fixed 1:4 sloped beach profile.The mobile-and fixed-bed experiments were conducted with a still water level (SWL) of 2.55 m and 2.65, respectively.In both experiments, regular plunging waves were generated with a wave height of H = 0.85 m and a period of T = 4 s.Thus, the bar slope resulted in a surf similarity parameter (ξ) of ξ = 0.44.Table 1 summarizes the main test parameters for both fixedand mobile-bed experiments, as respectively outlined in van der A et al. (2017) and van der Zanden et al. (2016van der Zanden et al. ( , 2017d)).
The fixed-bed experiment tests contained 123 runs of 38 min while 12 mobile bed runs were performed.The idea behind the significant number of model runs was to get a data set of high repetition rates and to take measurements at many locations along the flume using a movable frame, where some measuring devices were installed.Each run of the mobile-bed experiment consisted of six periods of 15 min of waves (90 min in total).The bed profile was measured after every two periods (i.e., each 30 min).The flow velocities and wave heights were measured in each run using the devices shown in Fig. 3. Wave heights were measured using resistive wave gauges (RWG) and pressure transducers (PT) installed on the flume wall PT wall and a movable frame PT frame .Flow velocities were measured over the water depth using Nortek Vectrino Acoustic Doppler Velocimeters (ADV), Laser Doppler Anemometer (LDA) located as shown in Fig. 3(b).Sediment concentrations were measured using a transverse suction system (TSS) and acoustic concentration and velocity profiler (ACVP).Fig. 3 shows the number and locations of the installed RWG, PT, ADV, LDA, TSS, and ACVP.
Since a hydrodynamic equilibrium was first reached after 300 s (=75 waves) in the mobile-bed experiments, the first 300 s of the measurements were discarded (following van der Zanden et al. 2016).For the fixed-bed experiments, the first 400 s were also disregarded (=100 waves) for the same reason (following van der A et al. 2017).Hence, 10 min (150 waves) and 31.3 min (470 waves) were used to make a quantitative comparison with the processed measurements of the SIN-BAD data set.More details on the treatment of mobile-and fixed-bed experimental results are also described in van der Zanden et al. (2016) and van der A et al. (2017), respectively.

Nonhydrostatic Mode of XBeach
Accounting for intrawave processes (e.g., wave skewness and asymmetry) in modeling the hydromorphodynamic processes around sandbars is essential for an accurate simulation of the net onshore flux and subsequent intrawave sediment transport (Hewageegana and Canestrelli 2021).Therefore, XBNH and XBNH+ are used in this study.The detailed formulations of XBNH and XBNH+ are, respectively, available in Smit et al. (2010) and De Ridder et al. (2020).This section provides only a brief overview of XBNH and XBNH+, including relevant equations to this study and descriptions of the model parameters used later in the discussion section.XBNH and XBNH+ solve the relevant hydrodynamic processes driving sediment transport using the NLSWEs, including a formulation for the depth-averaged normalized dynamic (nonhydrostatic) pressure similar to a prototype version of the SWASH model of Zijlema et al. (2011).The NLSWE for a single cross-shore profile reads as where x and t = horizontal spatial and temporal coordinates; η = free surface elevation above an arbitrary horizontal plane; u = depth-averaged cross-shore horizontal velocity; h = total water depth; ρ = density of water; g = gravitational acceleration; n = Manning bed friction factor; ν h = horizontal viscosity; and q = depth-averaged dynamic pressure normalized by the density.The normalized dynamic pressure q is derived similarly to the one-layer version of the SWASH model (Smit et al. 2013(Smit et al. , 2014;;Zijlema et al. 2011), in which the depth-averaged dynamic pressure is computed from the mean of the dynamic pressure at the surface and at the bed.The dynamic pressure at the surface is assumed to be zero and changes linearly along the depth.In analogy to SWASH, XBNH applies a conservation of momentum approach to a breaking wave, thereby treating it as a hydraulic jump.This allows the model to capture the macro-scale effects of wave breaking on hydrodynamics in the surf zone without modeling turbulence (Mulligan et al. 2019).The horizontal viscosity ν h allows for lateral mixing in the horizontal circulations due to sub-grid eddies.It is computed in XBeach (XBNH and XBSB) using the Smagorinsky (1963) model.
Because of the 2DH nature of XBeach, vertical mixing is not modeled,  which means that a part of the turbulent viscosity (i.e., the vertical viscosity ν v ) is omitted.The horizontal eddy-viscosity ν h is given assuming alongshore uniformity and omission of alongshore current, as where c s = Smagorinsky constant ranging between 0 and 1. XBeach defaults to c s = 0.1 and can be adjusted with the XBeach's keyword (KW) nuh.The quantity Δx = computational grid size.To improve the computed location and magnitude of wave breaking, XBNH applies the hydrostatic front approximation (HFA) of Smit et al. (2013), in which the pressure distribution under breaking bores is assumed to be hydrostatic.The HFA includes an additional horizontal viscosity to prevent the generation of high-frequency noise in the wave profile due to the discrete activation of the HFA.These frequencies are introduced because the model must adapt to the enforced hydrostatic pressure distribution.This additional viscosity ν HFA h is expressed as where L mix = typical horizontal length scale over which mixing occurs and is computed as a function of the friction coefficient μ and the local water depth h (i.e., L mix = μh).The value of μ ranges between 0.75 and 3 (Smit et al. 2013) and has a default value of 1 in XBeach, which can be modified internally using the KW breakvisclen.Eq. ( 4) can be combined with Eq. ( 3) to read in which β = tuning factor to increase the viscosity during wave breaking.The value of β ranges between 1 and 3 and has a default value of 1.5 that can also be modified internally using the KW breakviscfac (Deltares 2015).The HFA is initiated once the rising rate of water surface exceeds a predetermined value of the maximum wave steepness criterion δ = ∂η/∂t, which is a model parameter ranging between 0.3 and 0.8.By default, XBNH considers hydrostatic bores when δ = 0.4 and the value of δ can be changed using the KW maxbrsteep.A variation in δ either causes premature or delayed initiation of energy dissipation, resulting in over-or underprediction of wave heights in the surf zone.According to Smit et al. (2013), predicted wave heights by nonhydrostatic computations are sensitive to the selected value of the additional viscosity ν HFA h (controlled by μ or β).They are also sensitive to the value of δ as an indicator of how well the model can represent the wave shape for highly nonlinear waves.
With the new release (XBeachX), a two-layer mode (XBNH+ of De Ridder et al. 2020) can be run to perform quasi-3D computations with the advantage of simulating sediment transport and bed morphology.As a result, the computation of the layer-averaged velocities u up and u down (Fig. 2) become possible to account, for instance, for onshore mass transfer in the upper layer (of a thickness h up ) and the resulting undertow in the lower layer (of a thickness h down ).Because of being based on the reduced two-layer approach of Cui et al. (2014) with minimized Poisson equations formulation to reduce the computational demand, XBNH+ does not explicitly compute the layer-averaged velocities u up and u down .Instead, it computes the velocity differences between the two layers dU and the depth-averaged velocity u, which can be converted into layers' velocities using where α = layer fraction (=h down /h, as shown in Fig. 2) and known in XBeach using the KW nhlay.Eq. ( 6) considers the mass conservation along the cross-shore profile by considering that the onshore flow in the upper layer is equal to the return flow in the lower layer (i.e., h up • u up = h down • u down ).In combination with the depth-averaged velocity u, u up , and u down can be computed based on the layer fraction α following Eq.( 6), so that the capability of XBNH+ to reproduce the undertow can be examined.
As illustrated in Fig. 1, propagating waves are shoaling toward the shoreline, evolving toward steeper front faces, and eventually become highly skewed, forming a sawtooth shape due to the transferring energy from the primary wave components to their higher harmonics (Elsayed and Oumeraci 2017;Hoefel and Elgar 2003;Martins et al. 2020).Increasing asymmetry leads to wave breaking in a plunging, spilling, surging, or collapsing way, depending on wave and bed profile characteristics.Because the water surface may overturn, the breaking process is a complex phenomenon and, in addition, air-water interaction occurs (Rainey 2007).Wave breaking increases turbulent kinetic energy throughout the water column (van der Zanden et al. 2016).As a result, additional impact and turbulence vortices are transferred into the water column, particularly in the case of plunging breakers.More specifically, the curling wave front transforms into a jet that impinges the water surface and penetrates the water column (van der Zanden et al. 2019) and may increase sediment stirring and suspension (Lim et al. 2020;van der Zanden et al. 2017b).
Though breaker-induced turbulence is not physically modeled in either XBNH or XBNH+, its effect on sediment stirring is implicitly integrated based on a parametrized model suggested by van der Zanden et al. (2017c).This model accounts for the turbulent kinetic energy in the sediment stirring calculations and has proven to enhance the sediment stirring, suspension, and accretion in the sandbar zone.A depth-averaged form of the latter turbulence-effect model was recently integrated into the calculations of the equilibrium sediment concentration in XBNH and XBNH+ by Jongedijk (2017) so that the sediment concentration in the water column may be increased based on the value of the depth-averaged turbulence kinetic energy k turb .The latter energy is calculated as a function of the depth-averaged velocity u and a calibration parameter β d that can be adjusted using the KW betad.The following section discusses the numerical setup of the fixed and mobile bed experiments using XBNH and XBNH+.

Numerical Model Setup
Hydrodynamic and morphodynamic simulations are used to numerically reproduce the fixed and mobile bed profiles shown in Fig. 3.These results were then compared with the available experimental data (van der A et al. 2017;van der Zanden et al. 2016van der Zanden et al. , 2017c)).Table 1 presents the main physical test conditions and boundaries used to set up the fixedand mobile-bed numerical models.Since the experimental fixed-bed tests (the 123 tests) and the mobile-bed tests (the 12 tests) were conducted with identical boundary conditions, numerical setups using single sets of initial and boundary conditions were sufficient.However, to examine the optimal value of the layer fraction (α) along the sandbar zone, four numerical runs were performed for the fixed bed conditions (Table 2) with (1) one single layer, (2) two layers of a 0.15% fraction (i.e., α = 0.15), (3) two layers with default layer fraction of 0.33 (i.e., α = 0.33), and (4) two layers with an equal fraction (α = 0.50).The layer fraction that provided the best results in terms of water surface, wave height, and undertow profile was used for the morphodynamic simulation, for which two runs (one with a single layer and one with two layers) were performed.
To numerically reproduce the SINBAD experiments, a 1D profile was generated using the bed levels shown in Fig. 3.A total of 514 and 856 horizontal grid cells were used for the hydrodynamic and morphodynamic models.The hydrodynamic and morphodynamic model bed levels were z = −2.65 m and z = −2.55m, since z = 0 m was defined at the SWL and positive upward.
For the purpose of future reproduction, Table 3 summarizes the boundary conditions.A stationary wave boundary of constant wave energy is used to represent the regular waves at the offshore wave boundary.The lateral wave boundary conditions are considered by the default of XBeach (i.e., Neumann boundary condition), which means that the alongshore gradient of the wave energy is zero, assuming homogeneity in the alongshore direction.To solve the NLSWEs numerically, an initial offshore flow boundary condition was required.Therefore, the still-water level in the flume (0.00 m) was assigned with the KW front.From the lateral sides, the lateral flow boundary condition was used as a wall to prevent water discharge from or to the computational domain, similar to the flume sidewalls in the experiments.
Table 4 presents further model parameters.The simulation time was based on the experiments (i.e., 38 min for fixed bed and 90 min for mobile bed).The simulation starts at t 0 = 0 s with an output time step equals 0.1 s (tintg = 0.1).To only consider model outputs after the hydraulic equilibrium was reached (i.e., after 400 and 300 s for the fixed and mobile bed, respectively), the KW tstart was defined by these values.Furthermore, a uniform bed roughness was defined by a Manning coefficient n of 0.025 m −(1/3) •s, based on the reported value for the bed roughness in the SINBAD data.For the mobile bed simulation, D 50 , D 90 , and the morphological acceleration factor were set as presented in Table 4.

Assessment of Model Performance
In addition to the qualitative representation of the comparison between the predicted and measured hydrodynamics and morphodynamics, the root mean square error (RMSE), the bias, and the Brier skill score (BSS) were used to quantitatively assess the model performance.These measures compared observed (V m ) and simulated (V p ) outcomes as The angled brackets in Eq. ( 7) indicate the average of N readings or cross-shore locations.RMSE and Bias = commonly used statistical indicators to assess the performance of numerical models.In practice, bias is a measure of the overestimation or underestimation of the model.BSS = predictive skill relative to a baseline prediction V i , which is a frequently used indicator for the morphodynamic prediction skill (Elsayed and Oumeraci 2017).In this study, the RMSE and bias are used to assess the model's hydrodynamic performance, and the RMSE and the BSS are used to assess the morphodynamical performance.According to van Rijn et al. (2003), model performance with the BSS can be classified as follows: BSS < 0 equals bad fit, BSS = 0-0.3equals poor fit, BSS = 0.3-0.6 equals reasonable/fair fit, BSS = 0.6-0.8equals good fit, and BSS = 0.8-1.0equals excellent/best fit.

Results and Initial Discussion
This section presents the comparison of the simulated and observed hydrodynamic and morphodynamic processes.The observations of the fixed-bed and the mobile-bed experiments are separately used to validate the simulated hydrodynamic and morphological processes.

Hydrodynamic Processes
The following subsections investigate the ability of XBNH and XBNH+ to reproduce the fixed-bed experiments.First, the focus is on the layer fraction's effect on simulating the instantaneous wave shape, wave propagation process, and water-surface elevation.
The wave height, wave set-down and setup, and the depth-averaged flow velocities are then considered.

Water Surface Elevation
The simulated water-surface elevation at the time at which the hydraulic equilibrium was reached (i.e., at t = 400 s), as well as half a wave period afterward (i.have the same boundary and initial conditions, the two-layer runs yield a wavelength that differs from that obtained with a single layer, as shown in Figs.4(a and b).As a result, the single-layer run produces different wave propagation, particularly on top and onshore of the bar.According to Cui et al. (2014), the reduced twolayer mode (i.e., XBNH+) results in a more accurate wave period (and, hence, wave length) than does the single-layer option.This difference in wavelength/wave period explains why the breaking point locations differ for different numbers of computational layers.
The number of layers influences the instantaneous water surface, yet the layer fraction has almost no effect on wave propagation.Moreover, the single-layer model run produced instabilities, as shown in Fig. 4, where sudden bed level gradients are situated (i.e., at x = 34 m).The two-layer mode improves the modeling of frequency dispersion where steep bathymetric gradients and sudden bed changes are present, thus reducing instabilities in two-layer runs.This means that improving the frequency dispersion by increasing water layers enhances model behavior at locations of sudden bed changes (see also, e.g., Klaver et al. 2019).
To quantify the model performance, the water surface elevation (η) is averaged over the wave period at the locations shown in Fig. 4(c) (i.e., Sections d-k).As a result, the outcomes of the four model runs are compared with the experimental results of van der A et al. ( 2017) at these locations, as shown in Fig. 4(d).
The RMSE and bias of the comparison are outlined in Table 5.
Table 5 presents that XBNH and XBNH+ accurately predicted the phase-resolving water-surface elevation offshore of the bar, but their predictive capability decreased above the bar crest.The single-layer run by XBNH provides the lowest RMSE of the watersurface elevation offshore of the bar [which is also visible in Sections d-f] and the largest RMSE on top of the bar [Sections g-j].The two-layer run of 33%-layer fraction provides the lowest RMSE and bias on top of the bar.Thus, it provides the best predictions of the water surface at that location.It is concluded that the predictive capabilities of XBeach in terms of the water surface predictions are enhanced by the two-layer option (i.e., XBNH+).

Wave Heights and Setup
The predicted wave heights by the four model runs are compared with observed data in Fig. 5(a).In general, the model predicted the wave height accurately between the wave maker and the bar.However, on top and landward of the bar, XBNH and XBNH+ underestimate the energy dissipation.As a result, the wave energy predicted by them, and hence the wave heights, are generally higher than the observations.This may be explained by the simplified representation of the wave-overturning process in XBNH and XBNH + and the omitted roller-energy dissipation behind the breaker bar.The underestimation of the energy dissipation decreases for XBNH +; this may be due to the improved modeling of the frequency dispersion.Table 6 lists the RMSE and bias of the wave heights along the entire x-axis.The RMSE values of the wave height are approximately identical for the two-layer runs and were lower than those of the single-layer run.The latter reflects that the layer fraction does not affect the predicted wave heights.However, the positive bias values confirm that the model overestimates the wave heights.The layer fraction of 33% provides the lowest RMSE and bias values, thereby providing the best predictions of the wave heights.
The maximum and minimum water levels computed by XBeach are also compared with the observations [Fig.5(b)].The model accurately predicts the minimum water level.However, the underestimation of the energy dissipation, as mentioned earlier, leads to a higher maximum water level on top and landward of the bar.
The predicted time-averaged mean water level (MWL) is compared with the measured values, as shown in Fig. 5(c).Offshore and onshore of the bar, wave set down and setup developed due to the wave breaking.The model generally captures the timeaveraged MWL, as stated in Table 6, using the RMSE and bias values.Furthermore, Fig. 5(c) shows that the prediction of the MWL is fairly insensitive to the layer fraction (Table 6).However, on the bar's offshore side, the predicted MWL showed a better agreement with observations.Onshore of the bar, the single-layer run provided a better prediction of the MWL.The clear lag between modeled and observed setup, particularly for the two-layer runs, can be explained by the omitted surface roller that often captures a strong onshore momentum stored at the surface after breaking.Moreover, the higher predicted setup onshore of the bar may be explained by the less-predicted value of undertow (e.g., Gallop et al. 2020), as discussed in the following section.

Depth-Averaged and Depth-Resolved Flow Velocities
As shown in Fig. 3, the SINBAD data set's velocity measurements were carried out at several depth locations using ADVs and LDAs at 12 cross-shore locations at the bar area, between x = 49 and x = 65 m.The measurements are time and depth averaged to compare them with the computed time-and depth-averaged velocities (u), as shown in Fig. 6(a).The negative values of u, particularly over the bar crest, indicate the dominating offshore directed undertow, which drives continuous erosion of the beach with offshore transport.Fig. 6(a) also shows that computed u-values by XBNH and XBNH+ are well predicted offshore of the bar.On top of the bar, where the breaker-induced turbulence governs, XBNH and XBNH+ underestimate the depth-averaged velocities.Again, this underestimation may be a result of the limitation of XBeach to fully incorporate the complex wave-breaking process over the bar and the accompanying strong undertow.Landward of the bar, model predictions were reasonable.Fig. 6(b) shows that although the time-and depth-averaged velocities computed by the four model runs were identical, the velocity difference dU between the layers was different due to the difference in layer fractions.The velocity difference dU was considered zero in the case of the Table 5. RMSE and bias of the water surface elevations η predicted at the locations of measurements [i.e., Sections d-k in Fig. 4

Metric
Layer fraction single-layer run.In the case of different layer fractions, the values of dU and hence the layer-averaged velocities u up and u down vary along the profile because of differences in layer thicknesses.Using Eq. ( 6), the predicted layer-averaged velocities were computed and compared with the observed depth-resolved velocities at the 12 cross sections, as shown in Fig. 6(d).To perform a fair comparison between layers' averaged velocities and the measurements, observed depth-resolved velocities are also averaged over the thickness of the layers of 33% and 50% fractions.Outcomes of the comparison show high differences between computed and measured layer velocities, particularly at the sections located directly at the bar's onshore slope (e.g., h-j), where the undertow is significant.The latter proves that two layers are generally insufficient to fully capture the velocity variations over the depth and to predict the strong undertow as well as the location of the inflection point (Fig. 1) between the onshore-and offshore-directed momentum fluxes.Moreover, the inflection point's position varies vertically from one section to another.Therefore, a fixed-layer fraction  along the complete profile, as present in XBNH+, may limit the model performance; it eventually violates the physical processes with complex depth-varying velocity distributions often observed in ridge-runnel situations.

Morphodynamic Processes
This subsection investigates XBNH and XBNH+ performances to reproduce the mobile-bed experiments and morphological evolution of the sandbar.A comparison between the simulated profile evolution, sediment transport rates, and suspended sediment concentrations is made with those observed during the mobile-bed experiment.The outcomes of the single-layer run (Run 1), as well as the run of the 33% fraction run (Run 3), are presented, as the 33% fraction run had proven to be closest to the water level and velocity observations.

Sandbar Evolution and Time-Averaged Sediment Transport Rates
The difference in simulated and observed sandbar evolution is shown in Fig. 7(a).In general, the outcomes show similar patterns but largely different magnitudes.When compared with the observed sandbar evolution, the simulated morphological evolution of the model with a single layer and two layers showed a difference in the evolution of the sandbar after 30 and 60 min as well as after 90 min.
During the mobile-bed experiment, the sandbar amplitude (vertical distance from the sandbar crest to trough) increased from 0.39 to 0.74 m as a result of an increasing bar crest elevation and deepening of the bar trough.The bar crest increase was fed from one side by the dominant onshore-directed bedload driven by velocity skewness and asymmetry on the outer bar slope.The onshore-directed bedload indicated in Fig. 7(b) by the positive time-averaged bedload, which eroded the outer bar slope and accreted the sandbar crest (for more details, see van der Zanden et al. 2017b).Due to the wave asymmetry, onshore orbital velocities (during the wave crest phase) were larger than offshore orbital velocities (during the wave trough phase), thereby steering net bedload onshore.After wave breaking, the bedload increased again in the onshore direction on the sandbar's landward slope [Fig. 7(b)].This was explained by a gravitational contribution due to the landward sloping bar, that is, the bed-slope effect (van der Zanden et al. 2017b).On the landward side of the bar, the negative time-averaged suspended load [Fig.7(b)] dominated and was directed offshore due to the strong undertow, as indicated in Fig. 1 and shown through the observed velocities in Fig. 6(a).Hence, the suspended load contributed to both sandbar trough deepening and sandbar crest growth (van der Zanden et al. 2017b).This accumulation of sediment on the bar's crest due to onshore-directed bedload and offshore-directed suspended load was accompanied by an onshore bar crest migration of 34 cm.
While the onshore migration of the bar peak was reasonably predicted (50 and 30 cm with the single-layer and two-layer models), both XBNH and XBNH+ simulated a slowly decreasing sandbar crest elevation and sandbar trough infilling.As presented in Table 7, which compares outcomes of the single and two layers with observations, the RMSE values between cross-shore locations from x = 44 to 64 m increased to more than 10 cm during the 90 minutes of simulation.The latter indicates the incapability of XBNH and XBNH+ to reproduce the observed bar evolution properly.The run of two layers shows lower RMSE values, and its prediction is, therefore, theoretically better.However, the negative BSS scores indicate the "poor" performance of both XBNH and XBNH+, which was not overcome by the extra layer in XBNH+.
The incapability to simulate the sandbar's growth accurately is a typical limitation of current sediment transport formulations under breaking wave conditions, even after extensive calibration (van Rijn et al. 2013;Walstra et al. 2012;van der Zanden et al. 2017c).The sediment transport equations used in XBNH and XBNH+ are not intended for intrawave models, particularly under skewed and symmetric waves, and can cause deviations in the results (Mancini et al. 2020(Mancini et al. , 2021;;Ruffini et al. 2020;van Weeghel 2020).XBNH+ significantly underpredicts the offshore-directed undertow velocities, which would have resulted in additional offshore sand transport (Fig. 6).This returns (1) to the incapability of the model to simulate the wave overturning and the induced turbulence as well as the surface roller and (2) to the sediment transport formulation.Sandbar trough formation occurs when energetic waves break and stir sediment up into the water column.The latter indicates that present suspended sediment transport formulations are unable to correctly compute sediment suspension under these conditions (van der Zanden et al. 2017c).Hence, the overall amount of sediment being transported was underestimated.The simulated total load was solely onshore directed [see the positive values in Fig. 7(b)] and was 98% based on the suspended load's contribution.This contrasts with the experimental observations where no significant onshore-directed bedload was observed.The latter means that none of those processes (i.e., the wave asymmetry effect, the wave breaking effect, or the bed-slope effect) influenced the bed load in the simulation.b) the time-averaged total sediment load (q tot ), suspended load (q sus ), and bedload (q bed ).At the sandbar trough, the measured offshore directed suspended load during the experiment was more than the measured onshore directed bedload, and therefore the total load was offshore directed.The suspended load was offshore directed at the sandbar trough due to the large undertow velocity combined with the large suspended-sediment concentration (due to wave-breaking turbulence).In XBNH and XBNH+, however, the suspended load was onshore directed.This can be explained by the underestimation of both the undertow and the time-averaged suspended-sediment concentration.While depth/layer averaging leads to an underestimated undertow, the lack of wave-breaking turbulence to stir up sediment causes the underestimation of the suspended-sediment concentration.This process is described in the next subsection, focusing on comparing the intrawave suspended load.

Intrawaves Suspended-Sediment Transport
Since the suspended load is the dominant transport mode at the sandbar trough, this section evaluates the model's capability to predict the sediment suspension in more detail.To this end, the first 15 min of the simulated and observed sediment concentrations in the sandbar area were compared.Fig. 8 shows the measured and simulated normalized suspended sediment concentration c sus [-] and normalized horizontal flow velocity u [-].The normalization was performed over all the 12 cross-shore locations in Fig. 8(a), including the variation within the wave phase.Therefore, the cross-shore maximum and minimum c sus and u, as well as their relative variation over the wave phase, were compared.The normalization is defined as where c sus = suspended-sediment concentration; and u max and u min = maximum and minimum phase-averaged velocities measured or predicted at a specific cross section.The values of maximum and minimum concentrations are defined in a similar way.The phase variation of the horizontal flow velocity and the suspended concentration were averaged over 100 waves, starting after 300 s in the lab experiments (i.e., when the hydrodynamic equilibrium was reached).The experimental velocities were depth-averaged from three ADVs that have multiple positions over the depth.Moreover, the suspended-sediment concentrations were measured by the ACVP and averaged over a thin layer between 0.005 and 0.1 m from the bottom.From the simulations, the depth-averaged velocity and the depth-averaged total suspended-sediment concentrations were extracted.
One of the main differences between the simulation and the experiment is the location of the concentration peak.For instance, for the simulation, c sus = 1 can be seen in Section d, while for the experiment, this is achieved in Section g.This means that, during the experiment, the relatively large suspended-sediment concentrations were located at the sandbar trough (around x = 56 m).During the simulations, these were present around the bar crest (x = 53-55 m).This difference in location of maximum suspended-sediment concentration can be related to the effect of the wave-breaking turbulence on the suspension (i.e., entrainment) of sediment, which is not captured in XBNH and XBNH+.The models instead determine the sediment entrainment rate by the depth-averaged horizontal flow velocity u, which reaches its maximum on the bar crest [see Sections c-f in Fig. 8(b)].As a result of the advection and diffusion of the suspended sediment [implemented in XBNH and XBNH+ by the depth-averaged advection-diffusion equation of Galappatti and Vreugdenhil (1985)], erosion and deposition take place at different cross-shore locations.remained approximately constant during the wave phase.This means that the onshore and offshore transport rates were relatively balanced during the wave crest and wave trough phases.At the bar's landward slope, observed sediment concentrations were high due to the strong undertow, which stirred and suspended more sediment [see Sections f-h in Fig. 8(b)].Moreover, the wavebreaking turbulence in this zone enhanced sediment stirring and suspension.However, sediment concentrations varied considerably in the simulation during the wave phase [see Sections b-m in Fig. 8(b)].During the wave-crest phase, simulated-sediment concentrations were higher (following the increase in horizontal flow velocity), while during the wave-trough phase, the sediment concentrations decreased.As a result, the net simulated suspended sediment flux was directed onshore due to the difference in relative flow velocity and suspended sediment concentration between the wave-crest and wave-trough phases.

Turbulence Effect on Sediment Stirring and Bar Migration
The mismatching observed and modeled suspended sediment concentrations can partly be attributed to the omission in the models of wave-breaking turbulence on sediment stirring and suspension.To test this hypothesis numerically, Jongedijk's (2017) model extension (mentioned in the model description section) is considered in this section by activating the KW lwt to include the effect of the depth-averaged turbulence kinetic energy (k turb ) on sediment stirring in XBNH+.
Run 3 (two layers, α = 0.33) was replicated considering five values of the calibration parameter β d , and the sensitivity of the simulated sandbar evolution to β d was drawn in Fig. 9(a).Although a wide range of β d -values was tested (i.e., 0, 1, 10, 100, 150), no enhancement in predicting the bar migration was achieved.These findings agree with those of van Weeghel (2020), who concluded that accounting for the turbulence effect in XBNH does not provide much better predictions of sediment suspension and transport.The depth-averaged velocity at the sandbar trough zone (55 < x < 60) remains inaccurately predicted and tending toward zero [Fig. 6(a)].As a result, the predicted k turb at the bar trough is also low.Therefore, sediment stirring on the bar's onshore side is not supported by enhanced turbulence.This is further shown in Fig. 9(b), which compares predicted and measured k turb , when betad equal to 0, 100, and 150, respectively.
The measured depth-averaged k turb is the average of measured turbulent kinetic energy at three different relative depths ζ measured from the bed using the three ADVs of multiple positions over the depth (i.e., at ζ = 0.11, 0.38 and 0.85) along the 90 min simulation time (see van der Zanden et al. 2017b).A betad value of 150 provided comparable values with the measured k turb .However, predicted k turb remained near to zero-value at the zone of the bar trough (55 < x < 60) because of the significantly underestimated depth averaged velocity at this zone [Fig. 6(a)].This further indicates that the low hydrodynamic predictivity of XBNH+ is still the reason for the poor predictivity of the sandbar migration, even when the turbulence effect on sediment transport is considered.Hence, in addition to the limitation due to depth-averaging of the flow velocities and the inability to simulate overturning of the breaking waves, the sediment stirring by wave-breaking turbulence is underestimated.Effects attributed to these processes are, hence, insufficiently captured in XBNH+ and are seen as the main limitation causing the observed difference in suspended sediment transport and net bar migration.

Discussion around Tuning of Breaking Parameters and Possible Model Improvements
This analysis has revealed new insight pertaining to inherent limitations of XBNH and XBNH+ to simulate sandbar morphology, in particular: (1) lacking detail of simulating wave-breaking processes such as wave overturning, subsequent air-entrainment, and induced surface roller, and (2) shortcomings in reproducing the breakerinduced vertical turbulence and subsequent vortices and circulation in combination with the strong undertow.The first model limitation results in an underestimation of the rate of energy dissipation during the breaking process, leading to an overestimation of the wave height between the sandbar and the shoreline.Subsequently, this affects the predictions of the water-surface elevation and wave runup.The second model limitation results in an underestimation of the depth-averaged velocities over the sandbar and significant differences between velocities in XBNH+'s layers and measured depth-varying velocities.As a result, both XBNH and XBNH+ poorly predict the bed and suspended loads as well as the morphological evolution of the sandbar.Even when considering the effect of turbulence on sediment stirring in a parametrized way [as discussed in Jongedijk (2017)], the underlying predicted flow velocities were not yet accurate enough for a fair estimation of the turbulence effect.As a first attempt to overcome the presented limitations and enhance energy-dissipation rate during the breaking process, the following subsection examines a set of hydrodynamic parameters that affect wave-breaking simulation.

Role of Breaking Parameters
As a first attempt to improve XBeach capabilities in terms of breakerinduced turbulence and energy-dissipation rate, a parametric study is performed aimed at those parameters controlling the breaking process.A careful parameterization could help until more elaborate model assumptions and structure arise.To that end, four model parameters are considered: (1) the extra breaking-viscosity parameter (β) [Eq.( 5)], (2) Manning coefficient (n) [Eq.( 2)], (3) the maximum wave steepness criterion (δ), and (4) the horizontal viscosity ν h controlled by the Smagorinsky constant C s [Eq.( 3)].These parameters might theoretically contribute to underestimating the rate of energy dissipation in XBNH and XBNH+.Reasons for the parameter selection are outlined more closely in Table 8.To highlight the potential of these parameters, without performing a full sensitivity analysis for them as empirical tuning coefficients, five more model runs were performed by XBNH+, considering two layers of 33%-layer fraction.Thus, the new model runs are identical to Run 3 (Table 2) but with changes in these parameters' values.These runs are labeled Run 5 to Run 9, with details presented in Table 8.
To visualize the enhancement of the prediction capability compared with the baseline model run (i.e., Run 3), the instantaneous (at t = 400 s) and time-averaged wave heights are plotted in Fig. 10.Increasing the local breaking viscosity by doubling the value of β (Run 5) leads to more energy dissipation and thus dramatically reduces the overestimated wave heights compared with Run 3.This enhancement in model prediction can be noticed by comparing the RMSE and bias of the wave height (RMSE_H and Bias_H) for Run 3 and Run 5, as shown in Fig. 11.Doubling the bed friction, as the case in Run 6, also leads to slight reductions of the wave heights due to the friction-induced dissipation, thus reducing the RMSE and bias of the predicted wave heights compared with Run 3 (Fig. 11).Decreasing the maximum wave steepness criterion δ from 0.4 to 0.3 (Run 7) enhances the wave heights' prediction capability as it controls wave shapes.Because the magnitude of ν h depends on the gradients in the velocity field as well as the Smagorinsky constant C s according to Eq. ( 3), increasing the value of C s introduces extra dissipation and thus more accurate prediction of wave heights.Thus, it might be a substitute for the omitted vertical mixing in XBeach, which requires, for instance, a standard k−ɛ turbulence model (Launder and Spalding 1974) to be accounted for, similar to the SWASH model.Setting the Smagorinsky constant C s = 1 (Run 8) enhanced the energy dissipation dramatically, leading to significantly improved wave-height predictions (Fig. 10).The latter is reflected by much lower RMSE and bias values of the wave height (Fig. 11).Enhanced energy dissipation with Run 8 proves that the underestimation of energy dissipation in XBeach relates mostly to the model's inability to capture the details of the breaking process, including the vertical mixing.In Run 9, the combined effect of parameters β, δ, and C s is considered.The outcomes of Run 9 showed the lowest RMSE and smallest absolute bias value (Fig. 11).The parameters of the most significant influence on the energy dissipation are C s and β, which are measures of the local and turbulence viscosities, as shown in Fig. 11.
Although tuning the previous parameters adjusts the wave height by increasing the rate of energy dissipation, the RMSE and bias values for the MWL increase, as shown in Fig. 11, compared with the reference model run (i.e., Run 3) before the tuning.The latter is further confirmed by the clear overestimation of the MWL accompanying the breaking parameters tuning, as shown in Fig. 12(a).Moreover, these parameters' tuning does not enhance the predictive capability of the depth-averaged velocity over the bar, as shown in Fig. 12(b).In contrast, it decreases its prediction quality.Decreased prediction quality of MWL and flow velocity with tuning the breaking parameters means that parameters tuning negatively affects MWL and the flow velocities, which are essential contributors to sediment transport.Therefore, such tuning does not provide an improved prediction of the sandbar evolution that is basically driven by the velocity distribution throughout the water column.For this reason, the results of sandbar migration after the parameters tuning are not presented here.
Based on the outcomes of breaking parameters' analysis, it can be concluded that limitations of XBNH and XBNH+ to reproduce hydro and morphodynamic processes around sandbars are not related to tuning the breaking parameters rather than the underprediction or omission of some physical processes, in which tuning the breaking parameters failed to fully substitute their effects.Turbulence-induced mixing and surface roller represent examples of omitted processes, which are responsible for energy dissipation behind the bar and underestimated undertow (Rafati et al. 2021;Smit et al. 2010).

Toward Improved Sandbar Modeling
The previous analysis revealed that the key step toward improved predictive capabilities with respect to morphodynamics in XBNH and XBNH+ should start with improving the accuracy of the   8).
hydrodynamic modeling.This study aimed only to (1) draw the implications of model simplifications (e.g., simplified breaking without wave overturning) and assumptions (e.g., depth or layer averaging) for hydro-and morphodynamic processes around sandbars, and (2) report the current model limitations accordingly.Therefore, further model extensions to account for underpredicted/omitted processes are beyond this paper's scope.Nevertheless, suggestions for model improvements are presented in this section as an outlook for this study.
The first suggestion, which would not necessarily cover the case of plunging waves, is to increase the number of water layers so that the vertical structure of wave shear stress and wave-induced undertow is better captured, and the vertical wave turbulent mixing can be included.Eventually, the higher vertical resolution would result in a 3D model (similar to the SWASH model), and any gain in accuracy must be carefully discussed considering additional computational demand.A second more feasible way forward might be to integrate further processes to the governing equations.These are, for example, the surface roller, which would (1) increase the energy dissipation after the breaking, leading to better prediction of wave heights behind the sandbar, (2) capture a strong onshore momentum that reduces the clear lag between modeled and observed wave setup, as shown in Figs.5(c) and 12(a), and (3) result in more-intense undertow (see, e.g., Rafati et al. 2021;Wenneker et al. 2011;Xie et al. 2017).The surface-roller effect can be incorporated into the momentum equation of NLSWEs [Eq.( 2)] as a source term force in a way similar to XBSB, Wenneker et al. (2011), or Rafati et al. (2021).
From the perspective of morphodynamics, it is first crucial to develop/consider sediment transport equations that incorporate the intrawave effects of skewed and asymmetric waves on bedload transport (e.g., Mancini et al. 2020Mancini et al. , 2021)).In this way, the prediction of the dominant transport on the offshore bar slope can be improved.Second, improvements relevant to sediment entrainment predictions by energetic turbulent vortices and their effect on keeping sediment in suspension are necessary.Sediments must be kept in suspension to simulate offshore directed transport, as it is normally observed under energetic breaking wave conditions (Christensen et al. 2019;Lim et al. 2020;Mieras et al. 2019).

Summary and Concluding Remarks
This study examined the effect of 2D nonhydrostatic model assumptions (e.g., depth averaging) and simplifications (e.g., omitting wave overturning and surface roller) on the hydro-and morphodynamic processes around sandbars.The XBeach modes XBNH of a single water layer and XBNH+ of reduced two-water layers were used as tools for the latter purpose.The SINBAD data set (van der A et al. 2017;van der Zanden et al. 2016van der Zanden et al. , 2017avan der Zanden et al. , b, c, 2019) ) was used to assess the predictive skills of XBNH and XBNH+.The simulated hydrodynamic and morphodynamic processes in XBNH and XBNH+ were compared with the outcomes of these experiments.
The outcomes of the numerical simulations showed that XBNH and XBNH+ accurately predicted the phase-resolving watersurface elevation on the offshore side of the bar.The predictive skills decreased on top of the bar because of the simplified representation of the complex wave-breaking process.Fundamental processes, such as overturning, air-entrainment, and wave-generated turbulence, are not considered because of the model's 2D nature.XBNH and XBNH+, therefore, underestimated the rate of energy dissipation by wave breaking induced by the sandbar.As a result, wave heights after wave breaking were overestimated.XBNH+ showed slightly better results due to the more accurate computation of the frequency dispersion.The simulated instantaneous watersurface elevation was sensitive to the number of layers but insensitive to the layer fraction (i.e., the relative depth of the two layers).The two-layer run of 33%-layer fraction provided the best predictions of the water surface over and onshore of the bar.The results also showed that XBNH+ does not provide a better prediction of the undertow than XBNH, though discretizing the water depth into two water layers, which is an advantage over XBNH.Two layers are found insufficient to fully capture the velocity variation over the water depth.Using three different water layer fractions (α = 0.15, 0.33, and 0.50) showed that the position of the inflection point between the onshore-and offshore-directed momentum fluxes is variable along the cross-shore.Thus, selecting a fixed layer fraction in the complete model may limit the model performance.
In terms of morphodynamics, XBNH and XBNH+ simulated a slowly decreasing sandbar crest and no trough deepening.These observations contrasted with those of the experimental campaign.
The "bad" model skill was a result of (1) the implemented sediment transport equations that do not sufficiently incorporate the intrawave effect of skewed and asymmetric waves, (2) the underestimation of the undertow, and (3) the lacking relation between wave-breaking turbulence and its effect on keeping sediments in suspension.As a result, both XBNH and XBNH+ provided modest prediction skills of bed and suspended loads as well as of the bed morphology and bar migration.
The analysis showed that the limited representation of the wavebreaking process in single-and two-layer models has several implications for the simulation of surf zone processes.These include the overestimation of the wave height and setup landward of the sandbar and the inability to simulate bar trough deepening.To improve such a model's predictive skill in the surf zone, the option to incorporate a surface-roller model may (1) reduce the overpredicted wave setup by increasing surface momentum, (2) increase the energy dissipation behind the sandbar, and (3) enhance the offshoredirected undertow.Another improvement could be the increase in the number of depth layers and an improved representation of the wave-breaking-induced turbulence.Increasing the number of layers would optimize the inflection point's location and result in a better prediction of the undertow.The representation of breakerinduced turbulence, including its effect on sediment entrainment and trapping, could significantly affect morphodynamic predictions in the surf zone and require further research.
The reported model limitations of XBNH and XBNH+ may be generalized to any similar nonhydrostatic 2D model (e.g., REEF3D-SFLOW of Wang et al. 2020).The model user should be aware of such model limitations before applying it to real-world problems.The study assessed the model limitations in the case of a plunging breaker.With spilling waves, for instance, the discussed limitations might have much less effect because there is much less turbulence generated by the breaking wave.

Fig. 3 .
Fig. 3. Overview of the SINBAD experiments in the CIEM flume: (a) setup of the fixed and mobile bed experiments showing the location of the Resistive Wave Gauges (RWG), fixed on the flume's wall Pressure Transducers (PTwall) and mobile on a movable frame Pressure Transducers PTframe; (b) zoom in sandbar zone showing positions of Acoustic Doppler Velocimeters (ADVs), Laser Doppler Anemometer (LDA), Acoustic Concentration and Velocity Profiler (ACVP), and Transverse Suction System (TSS); and (c) close-up of instrumentation positions at one of the mobile frames located around x = 58 m. (Adapted from van der Aet al. 2017;van der Zanden et al. 2016van der Zanden et al. , 2017a, c.) , c.)

Fig. 4 .
Fig. 4. Comparison of measured and predicted water levels (η): (a) instantaneous predicted water levels for the considered four model runs at the time of reaching the hydrodynamic equilibrium (t = 400 s); (b) the same at t = 402 s; and (c) bed profile and locations of Sections d-k, where head measurements are available.Plots in (d) present averaged-over-wave-period T measured versus predicted water levels at locations of head measurements for the considered four model runs.

Fig. 5 .
Fig. 5. Predicted versus measured data for Runs 1-4 (Table 2) along the bed profile: (a) time-averaged mean wave heights; (b) maximum and minimum time-averaged water levels; (c) time-averaged mean water levels (MWL); and (d) bed profile and still water level.

Fig. 6 .
Fig. 6.Comparison between measured and computed flow velocities: (a) time-averaged depth-averaged velocity (u); (b) time-averaged differences between layers velocities (dU); (c) bed profile and locations of Sections d-o, where velocity measurements are available; and (d) time-averaged depth-resolved velocities at Sections d-o showed in (c).Sections are noted in brackets; [d] = Section d, etc.

Fig. 8
Fig. 8.Comparison of the suspended sediment concentrations: (a) profile evolution in the first 15 minutes and locations of Sections b-m, where sediment concentrations are compared; and (b) normalized sediment concentrations c sus and velocities u at Sections b-m.Sections are noted in brackets; [b] = Section b, etc.

Fig. 9 .
Fig. 9. Turbulence effect on sediment stirring and bar migration.Comparisons between: (a) observed and predicted bed-levels in the bar zone considering different values of the calibration parameter betad; and (b) predicted and measured time and depth-averaged turbulence kinetic energy (k turb ).

Fig. 11 .
Fig. 11.RMSE and Bias values (in cm) for the wave height and the MWL for the runs with adjusted model parameters (Table8).

Table 1 .
Summary of physical test conditions and boundaries for the fixed and mobile bed experiments H = 0.85 m and T = 4 s (regular waves) H = 0.85 m and T = 4 s (regular waves)

Table 3 .
Summary of boundary conditions assigned to the SINBAD numerical test Neumann Offshore flow Nonhydrostatic weakly reflective flow boundary Front: nonreflective boundary, defined by KW front = nonh_1d; back: nonreflective boundary, defined by KW back = abs_1d Lateral flow Wall boundary considers no flow through the flume sides Defined using the KWs: left = wall and right = wall

Table 2 .
Model runs for the numerical reproductions of the SINBAD data set

Table 4 .
Main parameters used to reproduce the SINBAD data set numerically

Table 6 .
RMSE and bias (in cm) of the wave height and mean water level predictions by XBNH and XBNH+ along the entire x-axis J. Waterway, Port, Coastal, Ocean Eng.

Table 7 .
Comparison of sandbar evolution as well as statistical indicators for the performance of the simulations J. Waterway, Port, Coastal, Ocean Eng.

Table 8 .
Characteristics of extra model runs to investigate the parameters that control the breaking process and energy dissipation in XBNH and XBNH+ s Respectively 1.5, 0.4 and 0.1 Respectively 3, 0.3, and 1 Combination of Run 5, 7, and 8