The Effect of Sediment Transport Models on Simulating River Dune Dynamics

River dunes, dynamic bedforms in the river bed, limit navigable depths during low flows and increase bed roughness. To predict the navigable depth or where maintenance dredging is needed, a numerical dune development model can be a powerful tool. To study the effect of sediment transport on dune shape and propagation, four different sediment transport models were applied in an existing dune development model. Each sediment transport model was able to simulate dune propagation, while only sediment transport models based on the shear stress reshaped the river dunes. The tested sediment transport models can simulate dune celerity similar to observations and realistic, though different, dune shapes for low and median discharges. Implementation of a gravitational bed slope effect combined with a critical shear stress results in low angle dunes, which are representative for river dunes during low river flows. Sediment transport models with spatial relaxation, also result in low angle dunes. However, the relaxation parameters need to be redefined for low flow situation to prevent transition to upper stage plane bed at too low flow velocities. Further analysis of the resulting dune shapes shows that the sediment transport model determines the dune shape in terms of slope angles, while the dune height is related to the total transport capacity.


Introduction
River dunes are dynamic periodic bedforms in the river bed (Figure 1) that propagate and change shape with changing flow conditions.These river dunes can cause bottlenecks for inland navigation during low flows and increase bed roughness (Van Rijn, 1984b).Predictions of dune development in rivers can help to efficiently plan fairway maintenance and predict flood levels.Numerical dune development models, which are able to quickly calculate changes in bed morphology, can support such predictions (Paarlberg et al., 2010).Computational power is ever increasing (Moore, 1965), which allows for increasing the complexity of numerical morphodynamic models.However, to make valuable predictions, a relatively simple and therefore fast dune development model is still needed to make the output relevant for operational decisions.recent studies have shown that the dune shapes found in large rivers show different behavior than in flume experiments: dune height relative to the water depth is smaller in deep flows (Bradley & Venditti, 2017), lee slopes in deep rivers tend to have small slope angles (Cisneros et al., 2020), and dune length becomes longer during low flows (Lokin et al., 2022).Until now, dune models have mainly been applied on flume conditions (Paarlberg et al., 2009;Shimizu et al., 2009) and under normal to high flows (Paarlberg et al., 2010;Van Duin et al., 2021), while application and validation for low flows are not yet present in literature.
Several numerical dune development models exist, with varying complexity in the flow and transport solvers.Ranging from numerical models solving the two-dimensional shallow water equations in the vertical plane [2DV] (Giri & Shimizu, 2006;Nelson et al., 2011;Paarlberg et al., 2009), to direct numerical simulation with particle-flow interaction for sediment transport (Nabi et al., 2013).For dune development forecasts, a model should be able to simulate dune shape and propagation accurately and much faster than real time.A more complex turbulence closure model, such as large eddy simulations, or a direct numerical simulation increases the level of detail in the resolved flow field.This also helps to better resolve changes in bed morphology, but takes multiple days to resolve laboratory scale dunes (Nabi et al., 2013), even with the currently available computer power.On the contrary, by using a constant eddy viscosity, simulations of the development of a river dune over a couple of days, take mere minutes to hours and also results in realistic dunes (Paarlberg et al., 2009;Van den Berg et al., 2012).Therefore, 2DV models using a constant eddy viscosity are promising, for real-time forecasting.In this study, we adopted the dune development model developed by Paarlberg et al. (2009) and extended by Van Duin et al. (2017).This model has been validated for high flow situations (Paarlberg et al., 2010) and upper stage plane bed [USPB] (Van Duin et al., 2021), but this model was not yet extended or tested for low flow conditions.
During low flows, when bed load transport is dominant, dune shape depends on transport capacity through the ratio between shear stress and critical shear stress, with a change in shape at the offset of a mixed sediment transport state (Bradley & Venditti, 2019).Dunes wash out toward USPB, when most of the bed load is brought into suspension, during extreme high flows (Naqshband et al., 2014).This transition has been simulated by parameterizing the effect of suspended sediment with a step length model (Naqshband et al., 2016;Van Duin et al., 2017, 2021).This raises the question on the influence of this parameterization on the dune shape and propagation during low flow conditions.The objective of this study is to assess the influence of sediment transport mechanisms through applying different sediment transport models on dune shape and development in a 2DV numerical morphodynamic model for low flow conditions.The dune shape and celerity are compared to observed dunes in the Waal River (Lokin et al., 2022) focusing on low and median flow conditions.

Sediment Transport Models
Many sediment transport models have been developed over the past decades, all of which are based on different assumptions and fitted to sediment transport under various conditions (Ackers & White, 1973;Bagnold, 1966;Einstein, 1950;Engelund & Hansen, 1967;Meyer-Peter & Müller, 1948;Nakagawa & Tsujimoto, 1980;Van Rijn, 1984a).The equations are developed based on the assumption that flow parameters, such as flow velocity and shear stress, are the driving forces in enabling sediment movement and are therefore a proxy of the total transport.This section shortly reviews four transport models, to better understand how they influence dune development through their assumptions and limitations: (a) Engelund and Hansen (1967), (b) Meyer-Peter and Müller (1948), (c) Meyer-Peter and Müller with linear relaxation to incorporate spatial lag in sediment transport (Tsujimoto et al., 1990), and (d) the pick-up and deposition model by Nakagawa and Tsujimoto (1980).
These four transport models are chosen because they cover bed load and mixed load transport, which are the predominant transport modes during low flows.Also, the variety of the parameterization of the driving force, through both shear stress and flow velocity is covered.Lastly, because these transport models are commonly used in dune development models (Giri & Shimizu, 2006;Shimizu et al., 2009;Van Duin et al., 2017).Here q s is used for any sediment transport, regardless whether total, bed load or suspended transport is calculated by the model.Engelund and Hansen (1967) [EH] have formulated a total transport model considering both bed load and suspended load transport.This transport model resolves the total sediment transport using the depth averaged flow velocity, U(x), and the median grain size, D 50 , as main parameters.It does not include incipient motion.The model has been derived based on flume experiments and has a validity range for sandy bed material, with a D 50 ranging from 63 μm to 2 mm.

Engelund and Hansen (1967)
The implementation of the EH model follows Equation 1: g is the gravitational constant.Δ is the dimensionless relative sediment density (ρ s − ρ w )/ρ w with ρ s the density of the sediment and ρ w the density of the water.u * is the shear velocity equal to  √  , with h the reach averaged water depth and i the reach averaged bed slope.The spatial variability of the total transport is only proportional to the local depth averaged flow velocity squared, U 2 (x), this makes the total transport sensitive to small variations in the flow velocity.In time the transport indirectly depends on the average water depth through the shear velocity.All other parameters are uniform and steady over the computational domain.

Meyer-Peter and Müller (1948)
The bed load model of Meyer-Peter and Müller (1948) [MPM] is one of the oldest and most widely used sediment transport models.The model only resolves bed load transport, which is determined based on the bed shear stress and a criterion for incipient motion trough the critical shear stress.In this study, the applied formula includes the gravitational bed slope effects (Sekine & Parker, 1992): with β is a calibration factor calculated by   =  Δ with m = 4, τ(x) the local shear stress, τ c the critical shear stress, η the bed slope parameter, which ensures the down slope preference of particles and is inversely related to the angle of repose, ϕ, (Sekine & Parker, 1992) and the local bed slope with z b the bed elevation and x the horizontal distance.
The relation between the resulting bed shear stress  (() − ()) and the bed load transport is non-linear, firstly trough the factor n = 3, and through the Heaviside step function  (() − ()) , that ensures zero sediment transport when τ(x) < τ c (x).The critical shear stress, including the gravitational bed slope effects is defined as follows: The initial critical shear stress, τ c,0 , is derived from the Shields parameter for a flat bed,   *  , by: (4)

Meyer-Peter and Müller With Linear Relaxation
Assuming that sediment transport needs some distance to adapt to non-uniform flow, the equilibrium transport as calculated with the MPM model is adjusted with a so called relaxation [MPM-LRE] (Tsujimoto et al., 1990).
The relaxation is calculated with the relation proposed by Tsujimoto et al. (1990): where q s,e is the equilibrium transport as calculated by the MPM model (Equation 2), q b the transport with relaxation and Λ the mean step length.This mean step length is defined as the average distance traveled by a sediment particle when it is entrained from the bed (Einstein, 1950), and can be calculated as follows: where α is the non-dimensional step length.Many models have been proposed to determine the value for α (Sekine & Kikkawa, 1992;Shimizu et al., 2009).In this paper the model as proposed by Van Duin et al. (2021) is used, because it has been adjusted for deep rivers.The value of α is defined as: where α min is the base value of the non dimensional step length, equal to 50, and h ref is a calibration factor set at 0.1166.Value τ* = 0.5 is related to the transition from dunes to USPB (Van Duin et al., 2021).

Nakagawa and Tsujimoto (1980)
The pick-up and deposition model of Nakagawa and Tsujimoto (1980) [NT] determines a transport gradient based on the probability of a particle being entrained from the bed, p s (x), and the probability of settling of that particle at distance Δx downstream of the pick-up point, p d (x).This model therefore accounts for the hopping of sediment particles over the bed. 10.1029/2023WR034607 5 of 17 The pick-up formula (Equation 9), F 0 is a calibration parameter which is determined at a value of 0.001 in this study, to match the dune celerity of the simulated dunes to the celerity in the Waal River.All stresses are defined in their non-dimensional form, see Equation 4. Like in the MPM model, the critical shields stress is corrected for the bed slope.
The deposition formula (Equation 10) states that the probability of deposited sediment at location x is the integral of the probability of the sediment being deposited at a distance s downstream of the point where it has been picked up.Therefore, first the pick-up at all locations, p x (x − s), has to be known.f(s) is the exponential distribution function of the deposition given by: Both the MPM-LRE model and NT model have been used to study the transition to USPB during high flows (Van Duin et al., 2017;Yamaguchi et al., 2019).

Dune Model Description
The The model consists of three modules, first a flow module, second a sediment transport module and third a bed update module (Figure 2a).These modules operate in a decoupled way.The flow field and shear stresses are calculated based on the input bed and the periodic boundary conditions.The resulting water depth, flow velocity and shear stress, serve as input for the sediment transport module, from which the output is used for the bed update.
The model domain consists of one single dune (Figure 2b), with periodic vertical boundary conditions.The flow conditions, transport rate and bed level at the upstream boundary are equal to those at the downstream boundary.This results in a virtually infinite row of identical dunes.This setup implies that dune length of the simulated dune is equal to the imposed model domain length.We chose a constant domain length during the calculations to reduce calculation times, even though the dune length can be assessed separately with a stability analysis (Camporeale & Ridolfi, 2011;Hulscher, 1996).

Flow Module
The flow field over the dune is calculated using the two-dimensional shallow water equations in the vertical plane [2DV], assuming hydrostatic pressure and a constant eddy viscosity.This restricts the applicability of the model to small Froude numbers.The horizontal axis, x, is aligned to the average bed slope, i.The vertical axis, z, is defined perpendicular to the bed average bed slope (Figure 2b).The free surface is defined as the average water depth, h, combined with the local variations, ζ.At the free surface the boundary conditions are: no flow through the free surface and zero stress.At the bed, z b , the boundary conditions are: no flow trough the bed and a partial slip condition.This partial slip condition is needed for turbulence closure in combination with the constant eddy viscosity.The derivation of the equations and the numerical implementation of the flow model are described in Van den Berg et al. (2012).
The assumption of a constant eddy viscosity results in a relatively small computational load.It only resolves a parabolic flow profile and neglects any flow separation at the dune lee and energy losses through flow expansion.
Low-angle dunes are dominant during low to median flow in rivers (Cisneros et al., 2020;Lokin et al., 2022) and flow separation is related to steep lee slope angles, with angles larger than 10° (Kwoll et al., 2017;Lefebvre et al., 2016).Therefore, flow separation is not the most relevant process to reproduce dune dynamics.
The flow velocity field is resolved using the water depth, while the specific discharge is used as an input value.
An iterative process determines the water depth and flow field combination that matches the specific discharge (Paarlberg et al., 2009).Because the model is not specifically calibrated to reproduce the water depths found in the field but to study the dune shapes, this routine results in a water depth that may not be equal to the water depth found in the field.Nevertheless, the resulting water depths are fairly close to the observed water depths (Table 1).

Sediment Transport and Bed Updating
The sediment transport rates are determined with the transport models described in Section 2. The bed evolution resulting from gradients in the sediment transport is calculated with the Exner equation (Equation 12), including the porosity, ɛ.

Modeling Approach
Dunes were simulated for a range of representative discharges covering the extreme low flows, up to mean flows in the Waal River (Table 1).The specific discharges (q Waal ) presented in Table 1 are representative for the mean discharge of the given discharge range (Q Waal ), and are extracted from the specific discharge at the river axis in a 2DH model (Rijkswaterstaat, 2017) which is used by the Dutch Ministry of Infrastructure and Water Management for water level calculations.The water depths (h Waal ) are derived from the gauging station Dodewaard (Figure 1b) and are defined as the average water depth for the given discharge range.These water depths are compared to the depths determined by the model (h mod ).
The simulations started from an undulation with a single sinusoidal shape with an amplitude of 0.1 m from which the dunes could grow to an equilibrium under steady forcing.The equilibrium is reached when the dune height increase is smaller than 1 mm per day.The domain length was set at 100 m for all discharge settings, as the observed dune length in the Waal river varies within 20% of 100 m (Table 1).The numerical grid size in x-direction is always 0.5 m.In the vertical direction the number of cells is always 25, equally distributed over the water depth.This grid ensures a minimal aspect ratio, cell height over length, of 0.3 for the simulation with the smallest water depth.
All simulations covered 100 days using a timestep of 10 s, this took 6-10 hr on a laptop (Intel(R) Core(TM) i7-10750H), depending on the chosen transport model.Simulations were terminated earlier when the simulated bed became numerically unstable, or when the dune height reached a value smaller than 10% of the initial bed which occurs when the simulation settings are representative for USPB.

Dune Dimensions and Celerity From the Simulations
The simulated dune celerity is defined as the horizontal displacement of the dune crest per day, averaged over the last 10 days of a simulation.The dune crest location is defined as the location of the maximum bed elevation,   max .The dune height is the vertical distance between the dune crest and trough, which is equivalent to the difference between the maximum and minimum value of z b (t).For the simulated dunes, the relative amplitude of the higher order wave numbers in the amplitude spectrum, k = 1/λ, indicates the deformation.
The simulated sediment fluxes are analyzed to assess the part of the calculated sediment transport capacity that contributes to the propagation of the dunes.By combining the simulated dune celerity and height, the sediment flux related to dune propagation is estimated by Equation 13 (McElroy & Mohrig, 2009;Simons et al., 1966).These estimated fluxes are compared with the fluxes determined by the model.

Validation Data
The simulated dune dimensions and celerity are compared with observed dunes from the Waal River.Lokin et al. (2022) have determined these dune parameters along the river axis for a wide range of flow conditions.The dune dimensions used here are the dune height (δ), dune length (λ) and celerity (c), as derived by Lokin et al. (2022).All dune parameters were determined in the middle of the fairway which is similar to the main channel axis.To compare the observed dune height, length and celerity with the modeled dunes, the average of the observed values within a discharge range was determined (Table 1).
For a quantitative dune shape comparison, a Fourier analysis was applied to the observed dunes in a smaller stretch the Waal River (Netherlands) for three periods with distinctive discharge regimes: (a) extreme low flows during the summer and autumn of 2018, (b) median flows during the spring of 2013 and, (c) decreasing flows from high to median during the winter and spring of 2018 (Figure 3).Within the observations analyzed in Lokin et al. (2022), a stretch of 1 km has been selected and studied over five consecutive, biweekly, bed elevation measurements.This 1 km stretch was followed using a Lagrangian window, such that the stretch always starts at the same dune.The amplitude spectra resulting from the Fourier transformation indicate the amplitudes belonging to the dominant wave length in the dune field and show whether the dune field has one or multiple dominant wavelengths.The amplitude spectra from the simulated dunes are compared with the spectra of the observed dunes to assess the dune dimensions.
Analysis of the amplitude spectra shows that the dunes change little in shape during the extreme low flows with a discharge varying between 500 and 750 m 3 /s (Figure 4a).The peak of the amplitude spectra decreases a little from 0.18 m for the 6 August 2018 profile to 0.14 m on 5 October 2018 (Figure 4b).The wave numbers for this peak correspond to wavelengths of 143-167 m.This indicates that dune length increases and the average dune height decreases over time while the extreme low flow persists.A small peak that stands out in the amplitude spectrum is the hump at k = 0.2 m −1 (Figure 4b), this corresponds to secondary bedforms superimposed on the primary dunes.These secondary bedforms may be a means of sediment transport during low flow (Zomer et al., 2021), but are not simulated in the model.
Several pronounced peaks are visible in the amplitude spectra of the dunes during median flow (Figures 4c  and 4d).The most pronounced peak is located at a wavelength of 125-143 m (k = 7.0e−3 to 8.0e−3) with an amplitude varying between 0.27 and 0.33 m.Next to the main peak, both in the higher and lower wave numbers pronounced amplitude peaks are found.The peak at the lower wave number corresponds to a wave length that is twice the wave length of the main peak.While the higher order peak corresponds to a wavelength of 90 m.This indicates that both the main wave length is present but that also a second wave length, although decreasing in height over time, is present in the bed.
A clear shift in the peak location in the amplitude spectrum is visible when the flow decreases from high to median flow (Figures 4e and 4f).The wave length shifts from 90 m just before the peak discharge to 66 m just after the peak discharge, and to the 143 m in the last profile in this period (22 March 2018).This shows that during the peak of the flood wave dune lengths still decrease, and increase again after the flood wave has passed.
While the individual dunes in all profiles in Figures 4a, 4c, and 4e have different wavelengths, the length corresponding to the maximum amplitude varies around 100 m.All spectra show a clear peak at the main dune length and some smaller peaks at higher wave numbers, corresponding to higher order waves present in the bed causing deformation of the primary dune.
Next to the amplitude spectra, a mean dune shape was determined for the low and median flow conditions.First, individual dunes were identified in the 5 measurements of each discharge regime, by selecting all crests and clipping the dunes on the surrounding troughs.Afterward, the individual dunes were shifted such that the crests were located in at the same vertical and horizontal position.Finally, the average, 25th and 75th percentiles of these dune levels at each position in the horizontal were determined to derive the dune shape.These mean dunes were compared with the simulated dunes.

Dune Height
The dune height for each combination transport model-discharge setting is shown in Figure 5a.For the simulations with the EH transport model, the dune height remains constant.Contrary to the dune height behavior in the field observations, the simulated dune height (δ at t = 100 days) decreases with increasing discharge in the MPM and MPM-LRE simulations.This decrease is most pronounced with the MPM-LRE model.Only for the NT model, the dune height increases slightly with the discharge.
Regarding the dune height over water depth relation, the trend in the simulation is similar to the trend found in the observations (Figure 5b).As the water depth and flow velocity increase the dune height becomes smaller relative to the water depth.This is mainly due to the differences in the simulated water depth, as the trend in dune height development in the simulations is opposite to the observations.In both, simulations and observations, the water depth over dune height is similar to values derived by previous research (Bradley & Venditti, 2017;Yalin, 1964).
A period of 75-90 days with a continuously low discharge is possible (Van Brenk et al., 2022), but the same amount of days with a high discharge in a perennial, rain dominated river is not very likely.Dune height development during the rising limb of a flood wave shows hysteresis with the changes in discharge (Bagnold, 1966;Warmink, 2014).This hysteresis can be simulated with the MPM and NT models as the dune height increases more slowly during low flows than during high flows (Figures 5d and 5e).Because the dune heights in the simulations with both the MPM and NT models grow to a fairly constant value, adjustments are needed to simulate the dune height decrease during the falling limb of a flood wave.

Celerity and Sediment Flux
For all transport models, the simulated dune celerity has the same order of magnitude as the observed dune celerity (Figure 6a), even though the transport models were not specifically calibrated for the use on dunes.However, the simulated sediment flux varies with multiple orders of magnitude between the transport models (Figure 6b).Important to note is the sediment flux simulated with the NT model, which is the pick-up rate only.This value is two orders of magnitude larger than the sediment flux with the other three models.The largest part of the pick-up is deposited within the same grid cell, so it does not contribute to bed development.
To test whether the calculated sediment flux contributes to the propagation of the dunes, the simulated sediment flux and the sediment flux were also estimated with Equation 13.For the transport models MPM and MPM-LRE, these fluxes are quite similar (Figure 6b).This indicates that the full sediment flux contributes to the dune propagation.For the MPM-LRE simulations with a discharge above 6.5 m 3 /s, part of the sediment stays does not contribute to dune propagation and stays in suspension.At these discharges, the relaxation becomes more important compared to the calculated bed load transport.The relaxation mimics the effect of suspended transport, which is not contributing to dune propagation and initiates the transition to USPB.
The EH transport model simulates dunes that do not grow and have the smallest celerity.Spatial variability in sediment flux in the EH model only depends on the spatial variability in the mean flow velocity through q s (x) ∝ u(x) 2 .As a result, the maximum transport coincides with the dune crest throughout the simulation causing the displacement of initial bed without deformation (Figure 7).The difference between the sediment flux determined by the model and calculated with Equation 13, results from the small variability in flow velocity compared to the dune averaged flow velocity.The mean flow velocity is of order 1 m/s, while the spatial variation over the dune is two orders of magnitude smaller.Only the variation around mean flow velocity contributes to the propagation of the dunes.The mean flow velocity contribution to the sediment flux can be seen as suspended load or an active layer of bed load, which does not contribute to bed changes.Therefore, a sediment transport model based on the mean flow velocity is able to simulate dune propagation.The amplitude spectra corresponding to the equilibrium dune shape with in gray the amplitude spectra from the observations.

Dune Shape
Where the dune height and celerity are important parameters for the prediction of navigable depths during low flows, the equilibrium dune shape is important to assess whether the model is capable of simulating dunes that resemble observed dunes.Each transport model is able to deform the dunes.Even for the EH transport model a small deformation of the initial sinusoidal undulation is visible.The lee slope steepens in both the median and low flow case.The MPM, MPM-LRE and NT transport models all show large deformations of the initial bed, resulting in different but realistic dune shapes (Figure 7).
Low lee slope angles are one aspect of these realistic dune shapes.The steepest slopes are simulated with the MPM model.The equilibrium dunes simulated with the MPM model have a maximum lee slope angle of 15.1 and 15.4° for q = 3.1 and q = 7.4 m 2 /s.In both cases the steepest slope is located at the trough.Here the transport capacity reaches a minimum due to a balance between bed slope effect for the critical shear stress and the imposed shear stress by the flow.By adding linear relaxation the maximum lee slope angles decrease to 5.7 and 3.1° for the low and median flow simulations respectively.These maximum lee slope angles are located on the lower half of the lee slope, and the location of this maximum moves toward the dune crest with higher discharges.
The NT model results in maximum lee slope angles of 12.2 and 8.1° halfway the lee slope.Although the lee slope angles are larger than found in the observations (Figure 7), in none of the cases the slope angles reach the angle of repose and an avalanching scheme is therefore not needed.Also the location, on the lower part of the dune lee is in line with findings in other rivers (Cisneros et al., 2020).
When considering the dune shapes obtained with the MPM model, the location of the maximum lee slope angle coincides with the location with the smallest shear stress.This coincidence is related to the slope correction on the critical shear stress, which has also the largest effect at the location with the steepest slope.The steep slope results in a smaller critical shear stress needed to maintain sediment transport, while the shear stress imposed by the flow also decreases over the lee slope.These two effects are in balance in the equilibrium dune, resulting in a smaller lee slope angle than the angle of repose.The steepest slope is then found at the point with the smallest shear stress.These results suggest that the local slope effect on sediment transport contributes to the formation of a low angle dune in bed load dominated sediment fluxes.
The comparison of amplitude spectra of the simulations with those of the observations (Figures 7b and 7d) shows that the amplitude of the main component of the simulated dune (k = 10 −2 m −1 ) after 100 days is two times higher than the observed dune during low flow.Except for the EH transport model, where the amplitude spectrum for the low flow situation is similar to the observed spectrum.However, for this case the initial amplitude is similar to the dune height in the field data.The relative values of the amplitudes of the first and second order components of the spectra for the MPM, MPM-LRE and NT models during median flow follow the observed amplitude spectrum.
The amplitude spectra for both simulated discharges show similarities with the observed spectra.The major difference is that the simulation only allows wave numbers complying to k = n/L, where n is any positive integer, due to the cyclic boundary conditions.Therefore the two dominant peaks in the measurements, one at k = 7.0e − 3 m −1 and the other at k = 1.1e − 2 m −1 cannot be simulated by the model, regardless of the chosen transport model.

The Influence of Transport Models on Dune Shape
To assess the influence of sediment transport mechanisms on dune development, dunes were simulated using four different sediment transport models.The results show that all transport models based on bed shear stresses are able to simulate dune propagation.
Through the dependency of the total sediment transport on the flow velocity in the EH-model, only a small part of the total sediment flux contributes to the propagation and deformation of the dunes.The EH model was able to simulate dune propagation but not dune growth or decay, because of the coincidence of the point of maximum sediment flux with the dune crest, in line with the hypothesis on dune propagation by Naqshband et al. (2017).
The MPM model performs best in terms of propagation speed and adaptation times of the dunes.Due to a balance between the shear stress and the bed slope effects, the lee slope angle remains smaller than the angle of repose, but larger than values found in rivers which are often smaller than 10° (Cisneros et al., 2020;Lokin et al., 2022).The location of the maximum lee slope angle is found at the dune toe, this is the location with the smallest bed shear stress.This balance limits development of the lee slope toward the angle of repose as long as the time-averaged shear stress is positive.This positive shear stress can only be achieved if no persistent flow separation zone is present at the lee slope.
Linear relaxation in the MPM-LRE model decreases the maximum lee slope angle and shifts the location of the maximum lee slope angle toward the crest.However, it also introduces parameterized suspended bed load, which contributes to the transition to USPB even at low flow velocities.The relaxation can enhance small lee slope angles, but needs to be implemented with care to prevent undesirable behavior, such as too much suspension, due to large step lengths.The same holds for the step length in the NT model.While the flattening occurs at a higher discharge, that discharge is still smaller than the expected discharge where USPB is expected in the Waal River (Van Duin et al., 2021).
In summary, processes included in the different transport models can help the understanding of dune shapes found in rivers.Firstly, a small lee slope angle may persist due to a balance between the time averaged bed shear stress over the lee slope and gravitational slope effects.This balance holds as long as the time averaged bed shear stress is positive.Secondly, the spatial lag between variations in flow and sediment transport decreases the lee slope angle, resulting in realistic low angle dunes as observed in the data during low flows.We note that the way this process is parameterized also mimics bed material brought into suspension to enable USPB.Therefore, the parameters in the step length model applied in MPM-LRE and NT need to be reconsidered for the application in low flow situations.Especially since the step length introduces diffusion which is seen during low flow in the Waal River (Lokin et al., 2022), further calibration of the dimensionless step length parameter is needed such that the dune development model can mimic this diffusive behavior.

Controls on the Dune Height
Two aspects stand out for the MPM, MPM-LRE and the NT transport models regarding the simulated dune height.First, the simulated equilibrium dunes are higher than the observed dunes.Second, the simulated equilibrium dune height decreases with discharge, and therefore also with water depth, for both the MPM and MPM-LRE models for all discharges and for the NT model discharges larger than 6 m 2 /s.Overall, the dune height varies little between discharges.These findings suggest that parameters, other than the chosen sediment transport model, control the dune height in the dune model.Nevertheless, the transport model determines the simulated dune shape in terms of slope angles, location of the maximum lee slope angle and smoothness of the dune crest and trough.
In river dune and marine sand wave models with periodic boundary conditions, dunes tend to grow to an equilibrium length that fills the complete domain (Campmans et al., 2018;Warmink et al., 2014), regardless of the chosen initial condition.The longest possible dune in the model domain will eventually become dominant.The domain length in this paper was kept constant at 100 m, which is the average dune length in the Waal River (Table 1).Changing the length to 80 and 120 m in the in the MPM and NT models results in the same trend in dune height (Figures 8a and 8c), but with ≈10% dune height difference, compared to a change in length of 20%.
Only the MPM-LRE model shows a different trend for dune height (Figure 8b).This is related to the step length in the transport model, which is set constant and independent of the domain length.As the domain length, and therefore the dune length, decreases and the step length remains the same, more sediment picked-up from the dune lee surpasses the dune trough leading to USPB (Naqshband et al., 2014(Naqshband et al., , 2017)).Therefore, the influence of the domain length on dune shape is fairly small.
The MPM and NT models result in similar dune heights.Also, these dune heights are larger than the dune heights in the observations.The sediment transport in both models depends on the relation between the flow induced shear stress and the critical shear stress, resulting in a transport capacity.The shear stress is calculated in the flow solver, which is based on the assumptions of a constant eddy viscosity and hydrostatic pressure.These assumptions are valid for steady uniform flow, which is not realistic in flow over a dune, with strong flow expansion on the lee slope.The simulated lee slope angles can induce intermittent flow separation (Best & Kostaschuk, 2002;Kwoll et al., 2016Kwoll et al., , 2017;;Lefebvre et al., 2016), increasing flow resistance, leading to smaller flow velocities and shear stresses and therefore decrease the transport capacity.In a bed load dominated system, the dune height is related to the sediment flux (Equation 13).Including those processes in the flow module may result in a more realistic dune height at the cost of significant increase of computational load.
This hypothesis was tested by increasing the critical shear stress in the MPM model.This increase is equivalent to a smaller shear stress due to a decreased flow velocity.The shear stress and critical shear stress balance each other in the transport capacity trough   ∝ ( − )  .The relation is non-linear, because n = 3 in the MPM transport model, the total transport decreases quite rapidly when  ( − ) goes to zero.When the dune becomes higher, sediment transport on the crest increases.Also the trough becomes deeper such that no transport is possible and it will fill up again limiting the growth at small τ − τ c .At higher  ( − ) the dune height increases up to a point that the water depth becomes growth limiting.Increasing the critical shear stress by a factor 2 results in smaller dune heights, that increase with increasing discharge (Figure 8a).
The effect of the increased critical shear stress on dune shape and sediment flux for both low and median flows is seen in Figures 8d and 8e.Within the low flow simulation with increased critical shear stress, the threshold of incipient motion is just surpassed and the maximum sediment transport flux is a factor 6 smaller than in the lower critical shear stress case.For the median flow simulation, the maximum sediment transport flux only decreases by a factor 2 for the same increase of critical shear stress.The result is a smaller dune height during low flow than during median flow, which corresponds to field observations.Also, the simulated dune shapes correspond better to the observed dune shape (Figures 8d and 8e) Furthermore, the dune celerity decreases, from 2.8 to 1.5 m/day for the low flow case, which is just below the 25th percentile of the observed celerity (1.7 m/day) and from 5.6 to 4.3 m/day in the median flow case, which is just above the 25th percentile of the observed celerity (4.1 m/day) (Lokin et al., 2022).Therefore, accounting for energy losses due to dune shape in a parameterized critical shear stress, can result in a more realistic dune height, without increasing the computational load of the model.

Toward Dune Forecasting
The model setup presented here can be used to obtain an estimation of the mean dune shape and celerity.As the trend in dune height resembles the observations well and dune celerity has the correct order of magnitude, even with little calibration on the critical shear stress.For application on different rivers, the choice of sediment transport model depends on the dominant transport mechanism based on the D 50 and flow regime.For the flow module some calibration is needed to obtain a good relation between flow depth and dune height.
For dredging operations and shipping, the dimensions of interest are those of the dune or dunes that will be the highest and limit the water depth and not specifically the mean dune height.A first step toward such predictions is to replace the initial condition in the model from the sinusoidal undulation by a measured profile over a longer domain such as presented in Figure 7 and force the model with the measured water depth.With such settings, interaction between different length scales can be simulated for periods up to a couple of weeks ahead.Running the simulation until an equilibrium dune is obtained will not be fruitful, because eventually a longer domain will be filled with only one dune (Campmans et al., 2018), since splitting does not occur in the model (Warmink et al., 2014).The tendency of the presented model to let individual dunes merge can also give insight into which dunes are more likely to survive the merging process during the falling limb of a flood wave.
Starting with a measured bed as initial condition still simplifies the river bed with 3D bedforms (Figure 4c) into a 2D line.However, the main component of the propagation direction of river dunes is in the same direction as the governing flow (Zomer et al., 2021, Figure 6).In rivers with one main channel, such as the Waal River, the mean flow direction is parallel to the river axis and the dune profiles are also taken along this axis.Therefore, the simplification to 2D profiles is not unrealistic.To determine the development of the 3D dune field, several profiles can be taken over the width and used in the model.Due to gradients in flow depth and velocity over the cross-section, three dimensional development of the dune field can then be estimated.This is an important issue for future research.

Conclusion
Four different sediment transport models have been applied on a 2DV dune development model to study the influence of different processes on the dune height, celerity and shape for low to median flow conditions.Without calibration, all tested transport models simulate dunes that propagate with a dune celerity comparable to observations and most models show realistic dimensions.When calibrated, the presented dune development can predict dune height and celerity during low to median flows.
The sediment transport model based on the flow velocity can simulate propagation of the dune without deformation.The sediment transport models based on shear stress are all able to simulate representative dune shapes.Due to the slope correction for the critical shear stress, lee slope angles smaller than the angle of repose can be achieved.This is related to a balance between the lee slope angle and the shear stress over the lee slope.Adding a form of relaxation on the equilibrium transport even further decreases the lee slope angle.
While the dune shape can be explained based on the parameters included in the different transport models, the dune height is controlled by the transport capacity of the flow.The dune height decreases with increasing water depth and transport capacity.The simulated dune height is only weakly determined by the separately set dune wavelength.

Data Availability Statement
Model documents, in-and output and post-processing scripts are published in the 4TU.ResearchData repository (Lokin et al., 2023).Field observation data are published in the 4TU.ResearchData repository (Lokin, 2022).Discharge and water level data can be downloaded via https://waterinfo.rws.nl/#!/kaart/Afvoer/Debiet___20Op-pervlaktewater___20m3___2Fs and https://waterinfo.rws.nl/#!/kaart/Waterhoogten/Waterhoogte___20Opperv-laktewater___20t.o.v.___20Normaal___20Amsterdams___20Peil___20in___20cm respectively (in Dutch), data can be downloaded by selecting measuring station Tiel Waal from the list "Uit lijst" under download more data "Download meer data".This research is part of the research program Rivers2Morrow (2018-2023).Rivers2Morrow is financed by the Dutch Ministry of Infrastructure and Water Management.All measurement data were made available by Rijkswaterstaat.Our words of gratitude for collecting and sharing these data go out to technical staff of Rijkswaterstaat.We thank ir. J.W.C.A. Lokin for his help with updating the model code, and dr.ir.Anouk Bomers and Jae Evans for their feedback on the first manuscript.We also thank dr. Alice Lefebvre, dr.Christian Winter and an anonymous reviewer for their constructive reviews.

Figure 1 .
Figure 1.(a) Example of river dunes in measured bed elevation of the Waal River (Netherlands) on 24 August 2018, the vertical lines refer to kilometers within analyzed data used in this paper.The black boxes in panels (b) and (c) indicate the location of the shown bed elevation data.
setup of the Dune Development model [DuDe] as presented by Paarlberg et al. (2009) was used as a basis for this study.DuDe is derived from a process-based morphodynamic sand wave model as developed by Németh et al. (2006) and Van den Berg et al. (2012), which again is based on the numerical model ofHulscher (1996).Van Duin et al. (2017) andVan Duin et al. (2021) have extended the model to simulate USPB.

Figure 2 .
Figure 2. (a) Diagram of the morphodynamic loop as used in the model.The bed topography is the input and each timestep, Δt, this topography is updated by following the morphodynamic loop.(b) Schematic of the model setup, the box illustrates the domain between the upstream and downstream boundaries.u, v and q s are the flow velocities in the x and z direction and the sediment transport capacity respectively.The light gray arrow indicates the periodic boundary conditions.

Figure 3 .
Figure 3. Discharge of the Waal River for the three analysis periods.(a) The discharge in m 3 /s for the whole period from 2010 until 2020, with the yellow band indicating a period with persistent flow just above median flow, and the green band indicating a period with high flows that fall to median flows, and the purple band persistent extreme low flows.(b) Zoom on the median flows, (c) zoom on the high to median flows, and (d) zoom on the low flows.The gray lines indicate bed measurement days shown in Figure 4.

Figure 4 .
Figure 4. Dune profiles and their amplitude spectrum corresponding to the periods shown in Figure 3. (a) Dune profiles during extreme low flows.(b) Fourier transform of the dune profiles during extreme low flows.(c) Dune profiles during median flows.(d) Fourier transform of the dune profiles during median flows.(e) Dune profiles during high flow decreasing to median flow.(f) Fourier transform of the dune profiles during high flow decreasing to median flow.The line coloring in the right column corresponds the lines with the same shade in the left column.

Figure 5 .
Figure 5. Dune heights in the simulations.(a) Equilibrium dune heights for all simulations, the scattered dots are the observed dune heights.(b) Water depth over dune height in the simulations, the scattered dots are the observed dune heights over the measured water depth.(c)-(f) Dune height development in time in the simulations.The markers indicate the moment the equilibrium height is reached.

Figure 6 .
Figure 6.Dune celerity in the simulations.(a) Observed celerity (Lokin et al., 2022) and simulated dune celerity, crest displacement per day.(b) Sediment flux including the pores, the full lines with the filled markers indicate the simulated sediment fluxes and the dashed lines with the open markers indicate the flux estimated with Equation 13 based on simulated dune height and celerity.The gray band indicates the estimated fluxes based on the dune height and celerity values in Table 1 ±50% uncertainty.The total flux calculated with the NT model, right vertical axis, is 2 orders of magnitude larger than the fluxes calculated by the three other models, left vertical axis.

Figure 7 .
Figure 7. Dune shapes at the last time step of the simulations.(a), (b) Equilibrium dunes shapes for median and low flow (a: q = 3.1 m 2 /s, (b): q = 7.2 m 2 /s) in the colored lines, the gray lines indicate the mean dune shape with the 50% interval in the light gray area from the observed dunes in of Figure4 and (c), (d).The amplitude spectra corresponding to the equilibrium dune shape with in gray the amplitude spectra from the observations.

Figure 8 .
Figure 8. Simulated dune heights with different model domain lengths with the MPM (a), MPM-LRE (b) and NT (c) models and dune shapes for the low (d) and median (e) flow situations with different critical shear stresses, simulated with the MPM model.τ* c indicates the dimensionless critical Shields stress, which is used as input parameter for the simulations.The simulated dune height with the increased critical Shields stress for the MPM model is also shown (a).
(Lokin et al., 2022)ns indicate values related to the flow in the Waal River, Q Waal is the discharge range, q Waal is the specific discharge at the river axis based on the mean of the discharge range and h waal is the mean water depth for the given discharge range.hmod is the water depth calculated by the model for the initial condition.λdata, δ data and c d ata, are the dune length, dune height and celerity based on the data(Lokin et al., 2022)for the indicated discharge range.