Modeling Streambank and Artificial Gravel Deposit Erosion for Sediment Replenishment

Sediment replenishment by artificial gravel deposits is a measure to increase sediment supply in gravel-bed rivers. Thereby, streambank erosion is the dominant process for gravel entrainment. In this contribution, we quantitatively validate a numerical morphodynamic 2D model and the relevant model approaches to reproduce non-cohesive streambank erosion. Therefore, a calibration and a sensitivity analysis of the relevant model approaches and parameters are carried out based on a reference laboratory experiment on streambank erosion in a straight channel from the literature. The relevant model approaches identified to successfully reproduce lateral streambank erosion are the gravitational bank collapse, the lateral bed slope effect on the bed load transport direction and the local bed slope effect on the critical Shields stress. Based on these findings, the numerical model was compared against data from laboratory experiments on gravel deposit erosion. Thereby, the focus was on the influence of the hydraulic discharge, the grain size distribution of the sediment and the geometrical quantities of the gravel deposits, such as the width, height and length of the deposit. It is shown that the dynamics of the erosion process were well reproduced by the numerical model using non-uniform sediment. Furthermore, the erosion rates were in good agreement with the laboratory experiments, except for the initial phase of the experiments, where the erosion rates were highest and settling of the gravel deposit was observed in the laboratory experiments. Overall, the numerical model proved to be a suitable tool to predict the erosion process of artificial gravel deposits, and hence, can be recommended for the design of sediment replenishment measures.


Introduction
Many trained gravel-bed rivers suffer from a distinct sediment deficit. In most cases, the sediment continuity is interrupted by dams, run-of-river (ROR) hydropower plants (HHPs) or sediment traps. Consequently, degradation and armoring of the river bed has been widely observed. Sediment replenishment by artificial gravel deposits is a possible measure to increase sediment supply in gravel-bed rivers. Thereby, the main objectives are to reestablish bed load transport, to enhance aquatic habitat for macroinvertebrates, plants and fish and to protect infrastructure by preventing further degradation [1][2][3][4][5]. Meanwhile, this technique has been applied worldwide, for example for the rivers Reuss [6], Limmat, and Aare in Switzerland [7], for the rivers Moosach, Lech, and Isar in Germany [8], for the Trinity River in the U.S. State of California [9] and for the Nunome River and several other rivers in the Kizu River basin in Japan [10]. Usually, additional sediment is obtained from sediment traps [7], stone quarries and gravel pits [7] or is excavated upstream of an ROR-HHP [6]. The transport to the replenishment site is mostly done by heavy trucks or a conveyor belt [1]. Finally, the gravel is placed in the river by a front loader, excavator or from a belt. In some cases, gravel is directly entrained by the flow. In most cases, however, artificial gravel deposits are placed in the form of gravel bars along the river embankment (see Figure 1). Eventually, the gravel will be entrained during flood events, where streambank erosion is the relevant fluvial process. The erosion process at the embankment of the gravel deposit is similar to the erosion process at river embankments. There is still a need for research to better understand the erosion process of non-cohesive embankments and to validate the numerical model approaches to reproduce this dynamic process. The correct reproduction of the non-cohesive streambank erosion process is essential for modeling artificial gravel deposit erosion. The material used for sediment replenishment is usually non-cohesive gravel with low organic content to minimize the negative ecological impact of fine sediments and turbidity, which can damage fish eggs and cause clogging of the river bed.
Various numerical models based on a cross-sectional approach have been proposed in the literature to compare to Ikeda's experimental data (e.g., [11][12][13][14][15]). Thereby, the main uncertainties were related to the appropriate choice of the the angle of repose, the Ikeda parameter, the bed load formula [11] and the critical Shields stress [13]. However, only a few numerical 2D simulations can be found to quantitatively compare with Ikeda's experimental data (e.g., [16]). In our opinion, this experiment presents a meaningful validation case for numerical modeling of non-cohesive streambank erosion. Furthermore, this allows for a detailed sensitivity analysis of the relevant model parameters assuming a heterogeneous embankment. In contrast, for cohesive streambanks, the natural variability of streambank and erodibility properties is usually high [17]. Hence, the input parameters for modeling river bank stability (factor of safety) are subject to a certain degree of uncertainty [18].
For numerical modeling, the gravitational bank collapse of non-cohesive sediment is often taken into account using a geometrical approach, where the collapsing bank material is transported in the direction of the steepest slope to the neighboring computational cells, until the angle of repose is reached again [16,19]. Furthermore, the angle of repose can be distinguished for submerged and dry cells. Hence, a larger angle of repose (steeper slope) can be assigned to dry cells to mimic the apparent cohesion effect [20]. In contrast, model approaches for streambank erosion of cohesive soils are generally based on the stability analysis with regard to a particular failure mechanism [21][22][23][24][25].
In the first step, our goal is to quantitatively validate the relevant model approaches to reproduce non-cohesive streambank erosion using the data from Ikeda's laboratory experiments on non-cohesive streambank erosion in straight channels [26]. The relevant model approaches for non-cohesive streambank erosion are identified, in particular the (i) local bed slope effect (correction of the critical Shields stress), (ii) lateral bed slope effect (bed load transport direction) and (iii) gravitationally-induced bank collapse. We highlight the importance of the interaction of all three model concepts to capture the dynamics of the lateral erosion process. Furthermore, the analysis further includes the quantification of the influence of relevant model parameters, such as the bed load pre-factor, the pre-factor of the lateral bed slope effect (also referred to as the Ikeda parameter) and the angle of repose for the bank collapse approach.
In a second step, the numerical model is compared against Friedl's laboratory data on the erosion process of non-cohesive artificial gravel deposits [27], where the focus was on the effect of the geometrical properties of the gravel deposit, such as width, height, length and geometry, hydraulic load (discharge) and grain size distribution. In addition to the bed load transport approach used for Ikeda's experiments, we apply a non-uniform bed load transport approach and highlight the effect of grain sorting for artificial gravel deposit and streambank erosion. A validated numerical model can support the planning of sediment replenishment measures to better anticipate the effect of hydrology and the optimal placement of artificial gravel deposits.
This leads to the research question of validating the key model approaches and governing parameters required for modeling the process of non-cohesive streambank erosion in straight channels and the erosion of artificial gravel deposits. This article is structured as follows: First, we describe the numerical 2D model focusing on the relevant model approaches for sediment transport, the two laboratory experiments and the setup used for numerical modeling. Then, we compare the simulation results with the reference data from both laboratory experiments. Finally, in the discussion, we compare our results with numerical simulations on streambank erosion from the literature along with a conclusion on the presented findings.

Numerical Model
The software BASEMENT was developed at the Laboratory of Hydraulics, Hydrology and Glaciology (VAW) of ETHZurich. For this study, we use the numerical 2D model of BASEMENT [28]. Water flow is modeled by the depth-averaged 2D shallow-water equations (SWE) consisting of the continuity and momentum equations. Morphological changes are based on the 2D Exner equation with closures to determine bed load transport (e.g., Equation (1)) and the Hirano approach [29] using the active layer concept in the case of non-uniform sediment. Figure 2 gives an overview of the numerical model and the relevant submodules. The topography is discretized with an unstructured triangular grid. Based on this, the SWE are solved with a finite volume method (cell-centered), while the Exner equation is spatially discretized on median dual cells (node-centered, [20]). Time integration is performed with an explicit Euler scheme. Thereby, first the hydrodynamics are calculated, and based on this, the bed load transport is determined (de-coupled approach).
In the following, we focus on the relevant bed load transport approaches, which are used in this study for modeling non-cohesive streambank erosion. The absolute value of the bed load transport q b = | q b | is calculated for uniform and unimodal sediment with the general form of the Meyer-Peter and Mueller bed load formula [30,31] as: where α = bed load pre-factor, ε = bed load exponent, g = gravitational acceleration, s = ρ s /ρ w = 2.65 (with ρ s = sediment density and ρ w = water density), d m = mean grain size diameter, θ = Shields stress, also denoted as dimensionless bottom shear stress (θ = | u| 2 /(c 2 f g(s − 1)d m ), with u = (u, v) T = depth-averaged flow velocity vector and c f = the dimensionless Chézy friction coefficient according to Equation (2) [32]), and θ c,δ = the critical Shields stress of the mean grain size diameter corrected for the local bed slope effect (Equation (6)). Note that the original pre-factor and exponent have values of α = 8 and ε = 1.5, respectively [30], but have been adjusted later to α = 4.93 and ε = 1.6, respectively [31].
The bed load transport formula for uniform sediment (Equation (1)) can be extended to non-uniform sediment as: where, in addition to Equation (1), f i = volume fraction of the grain class i of a total of n grain classes, d i = grain size diameter for grain class i, θ i = the Shields stress of grain class i (θ i = | u| 2 /(c 2 f g(s − 1)d i )) and ζ i = the hiding function of Ashida and Michiue [33]. The latter is a slight modification of Egiazaroff's hiding function [34] for smaller relative grain sizes (d i /d m < 0.4). A similar approach was proposed by Ribberink extending Equation (1) to non-uniform sediment using Egiazaroff's hiding function [35]. Note that in the case of non-uniform sediment, hydraulic friction (Equation (2)) is affected by grain sorting effects, where the equivalent sand roughness is calculated as k s = (2.1-2.6)d 90 for 2.0 ≤ σ g ≤ 2.5 [36], where d 90 is dynamically determined based on local grain size composition (d 90 = characteristic diameter of the grain size distribution where 90 percent is finer, σ g = geometric standard deviation of the grain size distribution).
The critical Shields stress θ c for a flat bed is determined based on Yalin and da Silva's approximation of the Shields diagram [37] as: where D * denotes the dimensionless grain diameter defined as: where ν denotes the kinematic viscosity (ν = 10 −6 m 2 s −1 at temperature T = 20 • C). The effect of local bed slope (subscript δ) on the critical Shields stress is taken into account using Van Rijn's approach as: where γ denotes the angle of repose and δ l and δ t are the bed slope in streamwise and transverse directions, respectively [38]. The effect of lateral bed slope on the direction of bed load transport (see Figure 3) is considered as follows: where ϕ b = bed load transport direction with respect to the main flow direction ϕ q , n q = normal vector of main flow q = (uh, vh) T (with u and v denoting depth-averaged flow velocity in the x-and y-direction, respectively), where n q points in the downhill direction ( s · n q < 0), s = ( ∂z b ∂x , ∂z b ∂y ) = bed slope (with z b = bed elevation), N l = the coefficient of the lateral bed slope approach, also referred to as the Ikeda parameter, and M l = the exponent of the lateral bed slope approach. Typically, values for N l range between 1.43 and 2.63, while M l has a value of 0.5 [39][40][41][42][43]. After the bed load transport is calculated, the bank collapse is applied based on a geometrical approach, where bank failure occurs if the local bed slope exceeds a predefined critical angle for dry and wet elements, γ dry and γ wet , respectively. The material is passed to the next lower cells until the local bed slope becomes smaller than the critical angle [19,20]. The effect of apparent cohesion can be reproduced to some degree by considering a larger angle for the material than its angle of repose. The angle of repose of a sediment mixture depends on grain size, grain shape, gradation, compaction, saturation and material and may range from 30-42 degrees for non-cohesive sediments [19].

Streambank Erosion in a Straight Trapezoidal Channel
Ikeda investigated lateral erosion in a straight channel using non-cohesive sandy material [26]. The laboratory flume was 15 m long, 0.5 m wide and had a bed slope of 0.215%. A constant discharge of 4.13 l/s was applied at the upstream boundary (Run No. 17 in [26]). Data of the exact sediment feed are not available, but sediment was added in the laboratory experiment at the upstream boundary in order to minimize erosion at the inlet. However, erosion at the inlet could not be completely prevented because there was no embankment upstream, which could have contributed to sediment supply. The mobile bed corresponded to a half symmetrical model of a trapezoidal channel, where the left side was limited by a vertical glass wall (see Figure 4). The duration of the laboratory experiment was twelve hours. The topography was measured in the two cross-sections CS9 and CS11, 9 m and 11 m from the inlet, respectively. The bed and the bank material consisted of a well-graded sand with a median diameter of d 50 = 1.3 mm and a geometric standard deviation of σ g = 1.3. The tendency towards meandering was prevented by the fixed glass wall [26]. The numerical setup was based on the laboratory scale of Ikeda's experiment. The topography was first discretized with a computational grid using approximately 35,000 elements. In addition to this medium discretized grid, a finer and coarser grid were evaluated in the scope of the sensitivity analysis ( Table 1). The initial topography of the numerical model was created based on the given geometry (see Figure 4b) and deviated only slightly from the original topography of the laboratory experiments (root mean square error of RMSE = 0.0031 m and RMSE = 0.0025 m for CS9 and CS11, respectively).
At the upstream hydraulic boundary, a constant discharge of 4.13 l/s was defined assuming uniform flow with the same bed slope as in the channel. Similarly, uniform flow was defined for the downstream hydraulic boundary condition. The initial hydraulic conditions consisted of a steady flow field. The roughness of the channel bed and the embankments were defined with an equivalent sand roughness of k s = 3 mm using Einstein's logarithmic friction law according to Equation (2) valid for hydraulic rough boundaries. Note that this results in a similar dimensionless friction coefficient c f as with the logarithmic friction law of Yalin and da Silva [37] (relative difference < 1.5 % with k s = 3 mm and water depth h = 60 mm). The initial water depth of h 0 = 60.2 mm matched fairly well with the measured water depth of h 0,re f = 61.0 mm. Note that during the simulation, the water table slightly increased by 2.4 mm, which was found to be in good agreement with the increase of 2.5 mm in Ikeda's experiment [26].
Sediment was added at the upstream boundary according to the local transport capacity, whereas a transparent boundary condition was applied downstream. The uniform bed load transport was modeled according to Equation (1) using a mean grain size of d m = 1.3 mm. The latter corresponds to the median grain size according to Ikeda's laboratory experiment [26] (d m = d 50 ).

Erosion of Artificial Gravel Deposits
Friedl performed laboratory experiments on the erosion of artificial gravel deposits in a trapezoidal flume with fixed bed and fixed embankments [27,44]. The experimental setup is shown in Figure 5 for both the trapezoidal and rectangular geometry of the gravel deposit. The experiments were inspired by the Reuss River below the ROR-HPP Bremgarten-Zufikon in Switzerland at a scale of 1:25 based on Froude similarity. The flume had a length of 35 m, a channel width of W = 2 m and a bed slope of S = 0.172%. The parameters examined were the water discharge Q, deposit width W d , deposit height H d , deposit length L d , mean grain size diameter d m , deposit geometry and geometric standard deviation of the grain size distribution σ g ( Table 2).  The experiments lasted between 0.5 and 3 h depending on the hydraulic load and the streambank erosion. This corresponds to a duration between 2.5 and 15 h in prototype scale. The topography was measured with a Laser Distance Measurement Sensor (LDMS) with an accuracy of <1 mm [44] in cross-sections at a distance of half a meter at irregular time intervals. Based on this, the downstream sediment supply Q b,sup was evaluated for 4-5 time steps depending on the experimental run. The downstream sediment supply Q b,sup was defined as the sediment supply to the downstream river reach just below the gravel deposit (Cross-section B-B in Figure 5). However, the experimental values contain a certain degree of uncertainty, because the bulk density ρ b = 1370 kg/m 3 from the initial gravel deposit differed from the bulk density of the material deposited in the near-bank zone (ρ b = 1800 kg/m 3 ) [27]. Furthermore, settling of the gravel deposit was observed at the beginning of most of the laboratory experiments. The work in [27] corrected the values of Q b,sup for the effect of settling during the first time step. Nevertheless, the experimental data of Q b,sup provide an idea of the overall erosion process and are used in the following for a comparison with the numerical simulations.
The numerical simulations were performed in prototype scale (1:1). The topography was parameterized as shown in Table 2. The unstructured and triangular computational mesh consisted of approximately 30,000 cells. In order to capture the lateral erosion dynamics, the spatial resolution of the gravel deposit was 15-20 cells in the lateral direction. This results in computational times of approximately real-time speed.
A constant water discharge was defined at the upstream hydraulic boundary assuming uniform flow with the same bed slope as in the channel. Furthermore, uniform flow conditions were defined at the downstream boundary. The hydraulic initial condition was defined based on a steady flow field. Hydraulic friction was accounted for using the logarithmic friction law according to Equation (2).
No sediment was added at the upstream boundary, while a transparent sediment boundary condition was applied at the downstream boundary. The initial topography corresponded to the initial geometry of the laboratory experiment. The gravel deposits were erodible, whereas the bed and the embankments were non-erodible. The non-uniform bed load transport formula was used according to Equation (3). The grain size distribution used in the laboratory experiment was discretized for the numerical simulations using seven grain classes according to Parker's procedure [45] (Figure 6). Similarly, grain size distributions were discretized for the coarse and fine sediment, Experiment (Exp.) Nos. 25 and 28, respectively ( Table 2). The thickness of the active layer was set to a constant value of L a = 2d 90,0 , where d 90,0 denotes the initial characteristic diameter of the grain size distribution where 90 percent is finer (d 90,0 = 68 mm, Figure 6).  [27], which are used for validation of the numerical model ( Figure 7). The simulated depth-averaged flow velocities are compared to the flow velocities from the ADV measurements at a relative height of 0.368 of the local water depth. At this height, the flow velocity corresponds to the depth-averaged flow velocity [37]. The lateral distribution of the depth-averaged flow velocities was reproduced fairly well at all cross-sections. The ADV measurements indicate that the flow is attached to the embankment apart from the lower end of the deposit, where recirculation occurs and 3D effects influence the flow separation zone. Based on the hydraulic calibration, the calibration of the morphodynamic model was performed. Good agreement with the laboratory experiments was found for a bed load pre-factor of α = 6.0, a bed load exponent of = 1.5 and an Ikeda parameter of N l = 2.0 (Equation (7)). Furthermore, the roughness of the gravel deposit and the area with deposited material was taken into account due to local grain-size composition using k s = 2.5d 90 , where d 90 denotes the local grain size diameter in the active layer where 90 percent is finer. This approach is in good agreement with the findings of [36] suggesting k s = σ 1.05 g d 90 = 2.22d 90 (with σ g = 2.14 according to Figure 6) and the more general range of 2.1d 90 ≤ k s ≤ 2.6d 90 (for 2.0 ≤ σ g ≤ 2.5). Based on the observations from the laboratory experiments, the angle of repose for dry and wet elements γ dry and γ wet was 40 • and 30 • , respectively. The initial bank slope of the gravel deposits was 35 • , and bank collapse was immediately observed in the experiments [44], resulting in smaller bed slopes for the submerged material and steeper slopes for the dry material due to apparent cohesion. Note that all values of the calibrated parameters are within the application range of the respective approaches according to the literature (see Section 2.1).

Evaluation of Morphodynamic Simulations
In the present paper, the numerical simulations were evaluated using the Brier Skill Score (BSS, Equation (8)) and the Root Mean Square Error (RMSE).
In Equation (8), MSE denotes the mean square error, which is defined as: where -brackets denote the mean value, Y a set of N predictions (model) y 1 , y 2 , .., y N , X a set of N observations x 1 , x 2 , .., x N , where each element of Y can be compared to X at the same place in space and time, and B a set of N observations of the initial topography b 1 , b 2 , .., b N . Note that linear interpolation between observations might be necessary for evaluating numerical results. The following classification for the evaluation using BSS (range given in square brackets) is adopted [46]

Streambank Erosion in a Straight Channel
The numerical model was calibrated using the medium-sized computational grid ( Table 1). The following parameters were found for the calibrated model: bed load pre-factor α = 1.6, lateral transport factor N l = 1.4, and angle of repose for dry and wet material γ dry =34 • and γ wet =30 • , respectively. The calibrated numerical model agreed fairly well with Ikeda's reference data over the entire time span of twelve hours for both cross-sections CS9 ( Figure 8) and CS11 (Figure 9). This was further quantified by means of goodness-of-fit measures, such as the Brier Skill Score (BSS, Equation (8)) and the RMSE. As a result, the BSS indicated an increasing agreement over time for both cross-sections. For instance, in cross-section CS9, the BSS improved from a good agreement for the first two time steps to an excellent agreement for the last three time steps of the experiment. The RMSE increased during the first three time steps and decreased during the last two time steps. However, the evaluation for cross-section CS11 was even slightly better. In this case, the BSS indicated a good agreement for the first time step and an excellent agreement for the last four time steps. Furthermore, the continuous decrease of the RMSE indicates an increasing agreement over the course of the simulation.
A sensitivity analysis was performed for the governing model parameters. Thereby, the focus was put on the deposition height at the bank toe (d) and the amount, or width, of the lateral erosion (e) of the first node above the water surface including a margin of ten percent of the freeboard (see Figure 4c). In order to refer to Ikeda's experiment, deposition d and lateral erosion e were compared to the reference data d re f and e re f , respectively. Thus, a perfect model fit would result in values of d/d re f = e/e re f = 1.0 (e.g., the straight black line in Figure 10). The following sensitivity analysis refers to the simulated topography of the cross-section CS11.
Firstly, the three main model approaches were evaluated: (i) local bed slope effect, (ii) lateral bed slope effect and (iii) bank collapse. Based on the calibrated model, these approaches were deactivated separately to analyze their influence ( Figure 10). For instance, without the local bed slope effect, the progress of lateral erosion was too slow. Furthermore, without considering the lateral bed slope effect, there was insufficient lateral transport from the streambank towards the middle of the channel. Finally, the bank collapse was relevant not only for the dislocation and transport of material, but also essential for lowering the dry nodes at the top of the streambank. The lowered nodes eventually became wet and exposed to the flow and hence were further eroded by fluvial entrainment. Without this model approach, lateral erosion was not possible (Figure 10). Based on these findings, it can be concluded that all three model approaches (i-iii) are necessary to capture the dynamics of streambank erosion.   Secondly, the sensitivity was evaluated for the (i) bed load pre-factor α, (ii) lateral bed slope factor N l (Ikeda parameter), (iii) angle of repose γ for bank collapse and (iv) grid resolution. The bed load pre-factor α turned out to be a sensitive parameter with regard to the lateral erosion e and deposition in the channel d ( Figure 11). Compared to the calibrated model, larger values of α increased the amount of lateral erosion e and vice versa (Figure 11, left). Furthermore, larger α increased deposition d during the first two hours due to increased lateral sediment entrainment. After twelve hours, however, the deposition d was smaller due to generally increased transport capacity (Figure 11, right). The lateral bed slope factor N l was a sensitive parameter, as well. This refers in particular to the amount of lateral erosion and to a lower degree to the deposition ( Figure 12). Compared to the calibrated model, larger values of N l increased the amount of lateral erosion e and vice versa (Figure 12, left). Hence, deposition d increased, as well, due to increased lateral sediment entrainment (Figure 12, right).
The angle of repose for dry and wet material (γ dry and γ wet ), which were relevant to trigger bank collapse, were also sensitive parameters concerning lateral erosion e and deposition d (Figure 13). In order to simplify the sensitivity analysis, identical values for critical angles were chosen for dry and wet material (γ dry = γ wet ). Compared to the reference run (γ dry = γ wet = 30 • ), larger critical angles led to a decreased amount of lateral erosion e and vice versa (Figure 13, left). A similar, although less pronounced trend was observed for the deposition d (Figure 13, right).  Finally, the grid resolution affected the lateral erosion process in particular during the initial phase of the experiment (Figure 14). Compared to the medium and fine grid, the lateral erosion was apparently underestimated when using the coarse grid, e.g., for time t = 0.5 h, t = 2.0 h and t = 4.0 h (Figure 14, left). Furthermore, the simulation with the coarse grid tended to underestimate the deposition in the channel, in particular during the first four hours. However, discrepancies almost vanished at t = 12.0 h (Figure 14, right). The simulation with the fine mesh did not significantly improve the results compared to the simulation with the medium grid (Figure 14, left).
The lateral erosion slowed down over the course of the simulation. However, at the end of the experiment (after twelve hours), bed load transport was not yet in equilibrium over the entire channel length. The sediment outflow was still four-times larger than the sediment inflow (transport capacity at the upper boundary). However, the sediment outflow and inflow were steadily converging over time, which indicates that the development towards an equilibrium channel was still ongoing.

Erosion of Artificial Gravel Deposits
Based on the above findings gained from the simulation of Ikeda's experiments, the numerical model was first calibrated as described in Section 2.2.2 and then tested against data from Friedl's laboratory experiments on the erosion process of artificial gravel deposits [27,44]. In the following, we first compare the simulation results with the reference data of a single experiment based on two cross-sections. However, for a more compact comparison of different simulations, we focus on the downstream sediment supply Q b,sup as an integrated quantity.

Erosion Process
The erosion of the gravel deposit and the successive deposition downstream was analyzed based on the laboratory Experiment No. 8 ( Table 2). The erosion process of the gravel deposit is shown for the cross-section in the middle of the gravel deposit ( Figure 15). Compared with the laboratory experiment, the erosion process was captured very well by the numerical model over the whole time span of ten hours (0.85 ≤ BSS ≤ 0.996, 0.028 ≤ RMSE ≤ 0.065). Furthermore, the bed surface texture at the toe of the embankment continuously coarsened over the course of the simulation, expressed by the mean grain diameter in the active layer d m compared to the initial mean grain diameter d m,0 in Figure 15. Obviously, the cross-section became wider due to lateral erosion. Consequently, bed shear stresses decrease near the embankment. This leads to a more selective bed load transport at the embankment, i.e., the finer grain sizes are more likely transported than the coarser grain sizes. During the stages of coarsening, a mobile armor develops, while a static armor is formed in the limit of vanishing bed load transport [47]. This grain sorting effect slightly stabilized the embankment and reduced streambank erosion compared to the simulation using uniform sediment (Figure 16).
In the laboratory experiments, the sediment eroded from the gravel deposit was eventually deposited downstream. A flow separation was observed at the end of the gravel deposit. Hence, a recirculation zone established downstream of the gravel deposit [44]. Despite some limitations in reproducing recirculation and 3D flow structures, the depth-averaged numerical model was able to reproduce the evolution of the sediment deposition dynamics well ( Figure 17).  Figure 5) over a period of 600 min using non-uniform sediment (axes normalized, Y = y/W d , Z = z/H d ).

Sediment Supply to the Downstream Reach
For the numerical simulations, the downstream sediment supply Q b,sup was determined based on the change in volume in the flume upstream from Cross-section B-B (see Figure 5) averaging over time intervals of 1500 or 3000 seconds and over the same irregular time intervals as used in the laboratory experiments [44]. As already mentioned in the description of the laboratory experiments, it is important to note that the reference data of the downstream sediment supply Q b,sup contain a certain amount of uncertainty and are used for comparison to give an idea of the overall erosion process. The downstream sediment supply Q b,sup from numerical simulations is further presented using a higher temporal resolution to show more clearly the effect of the varying parameters. Before discussing the results, the effects of a number of important parameters on the downstream sediment supply are investigated by sensitivity analysis based on the laboratory experiments.

Influence of Gravel Deposit Geometry and Grain Size Distribution
In the following, we investigate the influence of the gravel deposit width, height, length and geometry (trapezoidal connected to shore, rectangular off shore) based on the experimental setup. Firstly, the influence of the deposit width W d was studied in terms of a variation of the relative deposit width W d /W = 0.2, 0.25 and 0.3 ( Table 2). The supply rates obtained from the numerical model corresponded fairly well to those found in the laboratory experiments (Figure 18a). However, differences were observed mainly for the initial phase of the experiments. The influence of the deposit height H d was analyzed in terms of a variation of the relative deposit height H d /h 0 = 0.5, 1.0 and 1.5 (see Table 2). The gravel deposit was submerged during the laboratory experiment of runs No. 11 and 14 with H d /h 0 = 1.0 and H d /h 0 = 0.5, respectively. The erosion processes of the submerged and the dry gravel deposit observed in the laboratory experiment were slightly different. The non-submerged gravel deposits were eroded predominantly by lateral erosion, while the submerged gravel deposits experienced fluvial entrainment also on the top of the gravel deposit. The erosion rates from the numerical simulations were in good agreement with the laboratory experiments (Figure 18b). Furthermore, the influence of the deposit length L d was analyzed in terms of a variation of the relative deposit length L d /W = 1.0, 2.0, 4.0 ( Table 2). The erosion rates obtained from the numerical simulations agreed well with the erosion rates from the laboratory experiments (Figure 18c). Finally, the influence of the grain size distribution of the sediment mixture was investigated using (i) a fine mixture with a mean grain size d m = 25.0 mm, (ii) a medium mixture with d m = 37.5 mm and (iii) a coarse mixture with d m = 50.0 mm ( Table 2). In the beginning, the downstream sediment supply Q b,sup was in a similar range for all mixtures (Figure 18d). At this point, the Shields stress was sufficiently high to transport all grain classes. However, after three hours, differences became more apparent. Over time, the Shields stress on the embankment reduced, and grain sorting effects became relevant, especially towards the end of the simulation. The numerical model was able to reproduce the general trend of a slower erosion process and hence smaller Q b,sup for coarser sediment mixtures. This effect can be attributed to both grain sorting effects and larger grain size diameters.

Influence of Discharge
The influence of the discharge Q was investigated with experiment Runs No. 7, 8 and 9 ( Table 2). Note that the gravel deposit was submerged in the case of the laboratory Experiment No. 9 (Q = 703 m 3 /s). The downstream sediment supply rate was most different from the reference sediment supply rate during the first hour (see Figure 19a). On the one hand, the erosion in the laboratory experiments was influenced by settling of the gravel deposit, especially during the initial phase. Therefore, the evaluation of the experiments is subject to a certain amount of uncertainty. On the other hand, the numerical model is not necessarily able to reproduce the initially rapid erosion process. However, for the rest of the experiment, the numerical model is in good agreement with the laboratory experiment predicting significantly smaller erosion rates for smaller discharges.
For most simulations, the gravel deposit geometry was trapezoidal (connected to shore) and only rectangular (off shore) for three simulations. The sediment supply rates from numerical simulations are compared to the laboratory experiments in Figure 19b. Obviously, the rectangular gravel deposits eroded much faster compared to the trapezoidal gravel deposits. This is mainly due to the fact that the erosion takes place on both sides of the rectangular gravel deposit. Hence, the sediment supply was found to amount to around twice the sediment supply of the trapezoidal gravel deposit. Note that for the larger discharge (Q = 703 m 3 /s), no data are available because the erosion process in the laboratory experiments was very fast.

Discussion
In the following, we first discuss the effect of grid resolution and compare our simulation results from Ikeda's laboratory experiment with other numerical models from the literature. Finally, we highlight the grain sorting effects on the erosion process of artificial gravel deposits comparing with observations from laboratory experiments and field studies from the literature.

Effect of Grid Resolution
The quality of numerical 2D simulations depends on the spatial resolution of the computational grid, among others. In any case, the local bed slope is controlled by the local grid resolution. For example, morphological changes on a different grid size ∆y result in different bed slope angles γ ( Figure 20). Obviously, all relevant model approaches, such as the (i) bank collapse, (ii) lateral bed slope effect and (iii) local bed slope effect, depend on the local bed slope. Consequently, the calibration of the governing parameters corresponds to a certain grid resolution. Generally, a fine grid resolution is required to capture the dynamics of the lateral erosion process. The finite volume method implies that the results converge to a solution as grid resolution gets finer. However, a comparison between simulations with different grids remains difficult due to the dependency of the calibrated model approaches on the local grid resolution, such as the local bed slope effect, the lateral bed slope effect and the bank collapse. Therefore, a certain set of model parameters implicitly refers to a certain spatial resolution of the computational grid. Further research is necessary to find new consistent approaches showing convergence, e.g., using adaptive grids [48].

Comparison with Numerical Models from the Literature
Ikeda's experiment was used by various authors to validate mathematical and numerical models on equilibrium channels based on a cross-sectional approach [11][12][13][14][15]. These approaches are based on the same equations as adopted in the presented numerical model. However, they neglect longitudinal variations to reduce the complexity and focus on a single cross-section [15]. 2D models should be able to reproduce lateral erosion processes taking longitudinal variation into account, such as for natural rivers. For instance, Abderrezzak et al. used Ikeda's experiments to validate a numerical 2D model using van Rijn's bed load formula [16]. Their results were evaluated with a Brier Skill Score (BSS) ranging between 0.78 and 0.94 over the whole simulation, with BSS = 0.91 after twelve hours. These results are comparable to ours, but the present results show a slightly better agreement with the laboratory experiments and also reflect the temporal dynamics very well. For instance, the evaluation after twelve hours shows an excellent agreement with BSS = 0.99 and 1.0 for CS9 ( Figure 8) and CS11 (Figure 9), respectively. Abderrezzak et al. [16] performed a sensitivity analysis by testing the influence of another bed load formula (Meyer-Peter and Mueller), a finer and coarser mesh, a decreased critical angle for bank failure, an increased effect of bed slope in the flow direction on the bed load magnitude and of an approach to account for the local bed slope effect on the critical Shields stress. These results were compared to a reference simulation run only considering the first four hours of the experiment. The effect using another bed load formula was significant, which is in agreement with the present study showing a strong influence of the bed load pre-factor (α in Equation 1). A comparison of the grid dependency of the two numerical models is rather difficult due to the dependency on the local grid resolution of all relevant model approaches, as discussed above. However, the results of [16] seem to support the findings that a set of calibrated parameters is related to a certain grid resolution. The decrease of the critical angle for bank failure (from 40 • -38 • and 35 • ) was reported to have only a minor effect on the pattern of bed changes. However, this is not consistent with our findings, which indicate an increased lateral erosion when decreasing the critical angles (γ dry and γ wet ) from 38 • -30 • (see Figure 13). The main differences between the two numerical models are based on the approaches to account for bed slope effects. The work in [16] corrected both the magnitude of the bed load transport due to bed slope in the flow direction and the bed load transport direction due to lateral bed slope effect, both approaches proposed by Koch and Flokstra [16,49]. Because the bed slope in the flow direction is rather small in the laboratory experiments (0.215%), [16] found no significant influence when increasing the effect of bed slope in the flow direction on the magnitude of the bed load transport. In general, the local bed slope effect on the critical Shields stress not only affects incipient motion, but also the magnitude of the bed load transport. Although [16] used a different bed load formula and another approach to account for the local bed slope effect, they mentioned that this did not significantly improve the results. In contrast, we found that the local bed slope effect on the critical Shields stress significantly improved the temporal dynamics of the lateral erosion process. Finally, we investigated the effect of the lateral bed slope effect on bed load direction according to Equation (7) using an approach based on contributions of Ikeda and Parker [39,40,50], later adopted by Talmon [43]. The numerical simulations presented in [16] were performed with different parameter values (M l = 1.0, N l = 3.14; with θ c = 0.045) leading to a stronger lateral bed slope effect compared to the simulations presented herein. However, it is difficult to compare lateral bed slope effects for different numerical models because this approach adds diffusion to the numerical model. The total numerical diffusion in a model is composed of both the numerical diffusion due to the discretization scheme and the diffusion induced by the lateral bed slope effect [51][52][53]. As shown in the present sensitivity analysis, an increased Ikeda parameter N l (increased lateral diffusion) generally leads to increased lateral erosion and hence to larger deposition in the middle of the channel. A general approach for the lateral bed slope effect was proposed by Schmautz [54] distinguishing different types of grain movement, such as rolling and sliding (N l = 1.45, M l = 0.5) and saltation (N l = 0.75, M l = 0.25). Recent laboratory experiments from Baar et al. [55] suggest that the lateral bed slope effect depends on grain size, sediment mobility (θ/θ c ) and secondary flow intensity in the case of river bends.

Grain Sorting Effects
The main driver for the erosion process of artificial gravel deposits is streambank erosion, whereby the streambank is defined as the interface of the gravel deposit and the water phase. Grain sorting effects seem to play a minor role in the overall erosion process, especially for the initial phase. However, grain sorting seems to become increasingly more relevant over the course of the simulation. Indeed, at the end of the simulation, the basal area at the streambank experienced a coarsening of the bed material in the active layer (see Figure 15). Similar coarsening effects at the near-bank zone were observed both at laboratory and field scale. For instance, Requena [56] observed coarsening of the bed surface texture at the near-bank zone in laboratory experiments on lateral erosion of gravel-bed rivers. Such effects can also be observed in natural rivers, for instance at the river widening Mittlere Aue of the Töss River in Switzerland. At this place, VAW conducted a monitoring project of the temporal development of the river widening on behalf of the Office of Waste, Water, Energy and Air (WWEA) of the Canton Zurich, Switzerland [57]. The river widening is triggered by an artificial rigid structure consisting of large boulders to divide and direct the flow towards the river banks and, hence, induce bank erosion. However, coarser fractions of the collapsed bank material were not entrained by the flow and remained at the near-bank zone. This somewhat stabilized the streambanks and damped further lateral erosion during subsequent flood events. Hence, these observations support our findings based on numerical simulations, where armoring at the near-bank zone slowed down the lateral erosion process and hence fluvial entrainment of the gravel deposit.

Conclusions
We present a numerical 2D model for the morphodynamic erosion process of both streambanks in straight channels and artificial gravel deposits consisting of non-cohesive sediment. The model was calibrated and validated based on laboratory experiments on streambank erosion in a straight channel for uniform sediment [26]. We analyzed the sensitivity of the relevant model approaches and important model parameters for streambank erosion. Above all, we were able to highlight the importance of the combined use of the relevant model approaches, such as the (i) bank collapse for gravitational transport, (ii) lateral bed slope effect on the bed load transport direction and (iii) local bed slope effect on the critical Shields stress.
Based on these findings, we calibrated and validated the numerical model to reproduce laboratory experiments on the erosion process of artificial gravel deposits consisting of non-uniform sediment [27,44]. The numerical 2D model proved to be a suitable tool to locally estimate erosion rates of artificial gravel deposits and is useful for sediment replenishment planning. To use the model as a design tool for practical applications, the following points should be considered: • Bank collapse should make use of reasonable assumptions for the critical angles of the dry and wet material γ dry and γ wet , respectively; • The lateral bed slope effect on the bed load transport direction required a reasonable value for the Ikeda parameter in the range of 1.4 ≤ N l ≤ 2.0; • The local bed slope effect on the critical Shields stress should be taken into account, e.g., using Van Rijn's approach ( [38], Equation (6)); • The effect of grid resolution should be investigated in a sensitivity analysis; • A high resolution of the computational grid for the gravel deposit is required, e.g., 15-20 cells over the width of the gravel deposit; • Consider grain sorting effects on hydraulic friction in the case of non-uniform sediment: dynamically determine the equivalent sand roughness based on the actual local grain size composition, e.g., k s ≈ 2.5 d 90 ; • Use the non-uniform bed load formula to account for sorting effects, e.g., Equation (3), with the bed load pre-factor and bed load exponent within a reasonable range (4.93 ≤ α ≤ 8.0, 1.5 ≤ ≤ 1.6); • Overall, the erosion of the gravel deposit is sensitive to the water discharge and hence to the flood hydrograph used for modeling.
The presented numerical model proved to be a powerful tool to reproduce the relevant processes related to the investigated reference experiments. Nevertheless, further research and development of the numerical model is needed to improve the understanding of the streambank erosion of natural rivers. For instance, relevant processes include the curvature effect in river bends, the effect of vegetation, the role of cohesive embankments and the effect of weakening and weathering processes on streambank stability.