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

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, 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 fixed- and mobile-bed experiments, as respectively outlined in van der A et al. (2017) and van der Zanden et al. (2016, 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 SINBAD 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.

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 aswhere 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 $\overline{q}$ = depth-averaged dynamic pressure normalized by the density. The normalized dynamic pressure $\overline{q}$ is derived similarly to the one-layer version of the SWASH model (Smit et al. 2013, 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, aswhere 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 ${\nu }_h^{\mathrm{HFA}}$ is expressed aswhere 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 readin 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 ${\nu }_h^{\mathrm{HFA}}$ (controlled by μ or $\beta$). 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 usingwhere α = 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+.

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. 2016, 2017c). Table 1 presents the main physical test conditions and boundaries used to set up the fixed- and 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.55 m, 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 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₅₀, D₉₀, and the morphological acceleration factor were set as presented in Table 4.

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.3 equals poor fit, BSS = 0.3–0.6 equals reasonable/fair fit, BSS = 0.6–0.8 equals good fit, and BSS = 0.8–1.0 equals excellent/best fit.

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.

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.

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.

As a first attempt to improve XBeach capabilities in terms of breaker-induced 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 $(\beta )$ [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 $\beta$ (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 $\beta$, δ, 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).

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 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. 2020, 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).