**2. Case Study**

A problematic reach of the upper Hungarian Danube reach (river km (rkm) 1796–rkm 1794, Figure 1) has undergone major morphological changes during the last decades. Many studies presented in the literature [13–19] that because of installations of many river regulation measures (e.g., groin fields, ripraps and hydropower plant at rkm 1819) in the last decades, intensive gravel formations [19], important bed level incision [20] and bed armoring processes [16,19] could be detectable, mainly in the main channel [21]. In contrast, the bed content is much finer in the groin fields, causing siltation and erosion of the finer sediments during flood waves [19,22].

**Figure 1.** The sketch of the investigated Danube study reach (**left**) and grain size distributions taken from the investigated reach (**right**). The characteristic water discharges are mean flow *Qm* = 2000 m3/s, range of bankfull discharge *Qbf* = 4300–4500 m3/s, 2, 10 and 100-year flood event *Q2* = 5950 m3/s, *Q*<sup>10</sup> = 7950 m3/s and *Q*<sup>100</sup> = 10,400 m3/s [23].

The river reach can be characterized by the following parameters: the channel bed width at mean water-stage ranges between 150 m and 350 m [19] with the average water surface gradient of 0.0002–0.00025. As Figure 1 shows, the river section is regulated by conventional structures, such as groins and the banks are protected against erosion by ripraps. Also, sidearms, islands, gravel bars, confluence zone can be observed, which refer to the complex topography of the river section. Bed material samples were taken from the main stream, groin fields and gravel bars. Some of the grain size distributions can be seen in Figure 1 right. The figures refer to very diverse spatial bed contents (0.32 mm < *d*<sup>50</sup> < 70.5 mm, where *d*<sup>50</sup> is the median grain size). Such a wide dispersion of the bed content is a unique feature of the Danube River (~rkm 1600–rkm 1800); at the Lower Austrian Danube, (~rkm 1885, 90 km upstream), the Danube flows through a gravel bed, where *d*<sup>50</sup> is 21.1 mm without any finer fractions [24]. In turn, the middle Hungarian Danube (~200 km downstream) has a typical sandy bed with *d*<sup>50</sup> < 0.05 mm [25]. The complexity of the topography and bed content suggest spatially and temporally varied sediment transport nature [22]. That is in some places the gravel, elsewhere the sand transport dominates [26]. This kind of individual complexity was presented e.g., by the field measurements of Török and Baranya [19], or in [27,28].

Based on field measurements, the researchers could make approximate estimates regarding the ongoing local- and reach-scale morphodynamic processes. For instance, the essential bed changes in the last decades caused important water management related problems and also difficulties in navigation, which could be detected with the field measurements [19,20,27–29]. However, the studies emphasize that additional studies are required for even more reliable morphodynamic investigation. Computational modeling offers an alternative for analyzing and predicting the recent and expected processes [30,31]. That is, the numerical modeling based investigation could be an important complementary manner with meaningful added value.

For this reason, the reliable calculation of the morphological changes is a major interest to researchers and the application of a 3D computational fluid dynamics (CFD) sediment transport model became justified. However, the choice of the applied sediment transport model was not obvious. Many formulas can be found in the literature (e.g., [1–4,32]). Most of them are developed focusing on a given morphodynamic process (e.g., Wilcock and Crowe model: bed armoring [8]; van Rijn model: sand erosion and deposition [33] etc.). However, in the case of the examined river section, spatially and temporally varied sediment transport nature occurs. That is, none of the existing sediment transport formula are expected to operate reliably for both the sand and coarse bed material, within a given river reach.

#### **3. Materials and Methods**

#### *3.1. Introducing the Combining Sediment Transport Calculation Method*

A novel combined approach [12] of the van Rijn and the Wilcock and Crowe bed load sediment transport formulas were applied because of the complex and spatially varied bed material and dominant sediment transport nature. From now, the Wilcock and Crowe formula will be indicated with *W&C*, while the van Rijn will be with *vR*. A brief summary of these models can be read in the Appendix A. The combined manner was already presented and validated against laboratory measurements [12]. In that study, laboratory experiments were used, which were carried out in an 11 m long, 1 m wide straight flume with non-uniform bed material. A groin was installed 5 m far from the outlet, which caused scouring, bed armoring and sand aggradation processes in the channel [34].

Török and Baranya [22,26] pointed out a novel decision criterion, which is a suitable method for indicating whether the sand or rather the sand transport dominates locally. The herein presented combining method was based on this statement. Namely, if the shear Reynolds number (*Re\**) is below 300, the sand transport is prevalent. Otherwise, if the *Re\** occurs above 400 the gravel transport dominates. Based on these, the combined calculation method said that the local bed load sediment transport rate was calculated as the following (*qbi*,*W*&*<sup>C</sup>* is the sediment transport rate calculated by *W&C* and *qbi*,*vR* is the rate by *vR*):

$$q\_{hi} = \begin{cases} \begin{array}{c} q\_{hi, \mathcal{W} \& \mathcal{C}} \text{ if } \mathcal{R} \mathbf{c}^\* > 400 \\ q\_{hi, \mathcal{rR}} \text{ if } \mathcal{R} \mathbf{c}^\* \le 300 \\ \end{array} \\ \begin{array}{c} \text{else} \end{array} \end{cases} \tag{1}$$

where

$$f = \frac{1}{100}(R\varepsilon^\* - 300). \tag{2}$$

These models do not contain any variable to calibrate, the constants (e.g., the constant von Karman, or the lift coefficient) and base equations (e.g., the original Shields curve [35] was used in the van Rijn model) were defined as they are published and recommended in the original papers [8,33].

In addition to the bed load transport estimation, the suspended sediment transport was calculated in each computational grid according to the suspended *vR* formula [36]. Thus, the bed load was calculated according to Equations (1) and (2), while the suspended load was estimated by the van Rijn equation.

### *3.2. Applied 3D Flow Model*

The numerical model used in this study [37,38] solved the 3D Reynolds-averaged Navier–Stokes (RANS) equations with the *k-*ε turbulence closure (see e.g., [39]) by using a finite-volume method and the (semi-implicit method for pressure linked equations) SIMPLE algorithm [40,41] on a 3D non-orthogonal grid. At the boundaries, where the fluid flow cannot be considered as a free turbulence zone, the wall law was applied for the velocity profile calculation [42]. The momentum equations were in the complete form, describing the hydrodynamic effects in all directions. The roughening impact of the vegetation in the flood plain area was described as an energy loss term in the Navier–Stokes equations [43], and could be specified for each cell. Using this option, the effect of the vegetation was taken into account as a drag-effect.

In order to eliminate the boundary effect, the computational grid was longer in both upstream (rkm 1801) and downstream (rkm 1793.5) direction than the investigated ~ rkm 1795 and rkm 1799 river reach. The applied grid can be seen in Figure 2.

**Figure 2.** The computational grid of the investigated Danube study reach.

The study side was discretized with 355 cells in the streamwise direction and 150 cells in the lateral direction, respectively, resulting in the streamwise direction in an average resolution of 18 m and transversely 5 m in the main channel, while 13 m in the floodplain area. Vertically 11 layers were defined, at every ten percent of the depth, resulting in a maximum of 585,750 computational cells. The resolution of the grid was set based on previous 3D model studies [37,44–48] and Section 3.2.2 demonstrates its applicability.

The bed material was discretized by five fractions, which were: *d*<sup>1</sup> = 0.3125 mm, *d*<sup>2</sup> = 1.25 mm, *d*<sup>3</sup> = 5.7 mm, *d*<sup>4</sup> = 16.2 mm and *d*<sup>5</sup> = 56.57 mm. The ripraps and groins were characterized by *d* = 300 mm. According to field measurement based considerations, the active layer thickness was set to 0.5 m.

#### 3.2.1. Parameterization

For instance, Parker mentioned [11] that the bed material of most river reaches is less complex and the grain sizes happen in a narrow range. In turn, in rare cases—e.g., the herein studied river section—the occurring grain sizes cover a significantly wider range (silt–gravel), resulting in a very complex spatial distribution. The bed material cannot be supposed as spatially uniform because of this (as many studies do at most river sections [44,46]), which makes the allocation of the bed material a less obvious method. According to Baranya [49], a relation can be stated between the calculated local bed shear stress value by 3D flow model at the mean-water stage and the local *d50*. Thus, applying the fitted function, a transitional and continuous *d50* map can be estimated based on the calculated bed shear stress distributions. This method was used based on 33 bed material samples [19,26]. The standard deviation of the calculated *d50* to the measured *d50* was 3.2 mm.

As the boundary conditions for the RANS equations, at the inflow boundaries the water discharge, at the outflow boundary, the water level was set. The discharge time series (Figure 3) and the water levels of the Danube were defined based on the measured time series. Additionally, the flow discharge time series of the Mosoni-Danube was set based on a 1D numerical Danube model [50].

**Figure 3.** Discharge time series at rkm 1801 for the period October 2012–October 2014. The different colored periods marked with letters a–h display the flood waves (including the 100-year historical flood wave from 2013, e) which exceed the bed-forming flow discharge (*Q* > 2100 m3/s).

An essential part of the model setup is the correct set of the inflow sediment rate. For this purpose, the flow discharge dependence of the suspended load [51] and the bed load [19] functions were used.

The riverbed topography of the main river channel was available from 2012 and 2014. The initial bed geometry was set according to the map from 2012. The calculated bed change map could be prepared for this two-year-long period, which included the historical flood wave from 2013 (Figure 3). This bed change map was used as a benchmark for the validation purpose (Figure 4).

**Figure 4.** Measured bed changes for the period October 2012–October 2014.

Four regions in the river reach were highlighted by green rectangles and an ellipse, marked by A, B, C and D. In these places, the following bedforms and morphodynamic processes were detected by field measurements [16,19]. At region A and D, the blue spots refer to a pronounced scouring downstream of the groins. In region B a groin field could be found. Here, local bed changes took place, both scouring processes (blue spots) and sediment depositions (brown spots). The ellipse (C) and the brown spot in region D show the places where gravel bars were located. As these phenomena basically determine the reach-scale morphodynamic processes, a key question is whether the novel sediment transport calculation manner introduced a more reliable estimation of them.

The numerical simulation of the 2 years, 722 days long period demands very large computational capacity. According to the preliminary estimation calculations, the simulation of one model variant for such a long time period would take around half a year. Therefore, to reduce the duration of the simulation, only the periods that exceeded the bed-forming flow discharge (*Q* > 2100 m3/s) [19,24] were simulated. It meant that the inflow boundary condition was set to a continuous flow discharge series, including the highlighted periods in Figure 3a–h, respectively and omitting the intermediate periods. In this range, 66% of the annual bed load amount passed [23]. In turn, it is emphasized that the ignored 34% of the annual bed load yield was still significant. That is, the numerical model neglects the simulation of bed changes that take place during the lower water regime. However, many studies (including papers regarding the investigated river reach) display that the major bed changes such as scouring, bar formations and flushing of the groin fields are expected rather during floods [16,19,24,52]. That is, the mean and lower water stages are less important in the view of bed changes. Based on the above, the bed changes caused by the eight flood waves (Figure 3a–h, a total of 211 days) approach well the real two years changes. Hence, the calculated bed changes were comparable to the measured.

#### 3.2.2. 3D Flow Model Validation

The herein applied 3D CFD model was already adapted and validated for the investigated Hungarian reach of River Danube, which was published in previous research works, e.g., [45,46,53,54]. Those studies have already demonstrated the reliable application of the 3D flow model. Regardless of these, the flow model validation was elaborated for the peak of the historical flood wave in 2013. The cross-sectional acoustic Doppler current profiler (ADCP) flow measurements by the North-Transdanubian Water Directorate regarding the peak stage of the flood wave were used as benchmark flow values. Figure 5 shows the measured (left) and the calculated (right) cross-sectional velocity distributions, from exactly the same points. However, as the ADCP instrument is not able to measure at the direct water - and bed surfaces, the measured cross-sections are smaller slices, respectively. Furthermore, regarding the comparison it is important to mention that the moving boat ADCP measurements record the momentary flow conditions, which include the effect of turbulence, resulting in velocity pulsation in 0.1 m/s order of magnitude [25,53,55]. In turn, the RANS model calculates the time-averaged velocity values.

Comparing the measured and calculated cross-sections, the following remarks can be stated. The velocity values were in the same ranges (0–3 m/s). The highest velocities (yellow and red spots) were calculated at the same place as the cross-sections than in the real case, which underlined the reliable estimation of the main stream. The locations of the lower velocities (blue and light green spots) were calculated as trustworthy also. That is the calculated flow pattern could be realized reliable at the groin fields (e.g., at the right-bank sides of cross-section III, IV and V, Figure 5) and gravel bars also (e.g., at the right-bank sides of cross-section VII). Finally, the root-mean-square deviation (RMSD) between the measured and calculated velocities could be calculated for all the seven cross-sections. These values can be seen in Table 1.

**Figure 5.** Measured and calculated cross-sectional horizontal velocity distributions.

**Table 1.** The root-mean-square deviation (RMSD) between the measured and calculated velocities, the maximum cross-sectional horizontal velocity and the percentage deviation of them for the seven cross-sections.


Considering the maximum cross-sectional horizontal velocities (Table 1, third row), the RMSD values displayed that the deviation of the flow calculation was below ~12% of the maximum values. The only two sections, where the inaccuracy happens ~15% were cross-sections V and VI. In these cross-sections, the higher deviation suggested that the numerical model estimated a narrower main stream as the measurements display. Around these cross-sections, the varying geometry increasec the transverse flow structure. However, a known limitation of the time-averaged RANS modeling was just this phenomenon, which is the underestimation of the secondary flow strength [56–58]. Next to this section, the flow model could be considered as validated for higher flood waves too, taking into account the available measurement data. For a more comprehensive validation, fixed boat ADCP measurements would have been expedient, which could be used to obtain time-averaged velocity values and also to estimate local bed shear stresses. The model validation for the 100-year flood wave should be evaluated in this light.

#### **4. Results**

#### *4.1. Comparison of the Calculation Methods*

In order to confirm the better operation of the novel *Re\** dependent combined method, the bed change calculation by the van Rijn, Wilcock and Crowe and the combined method were compared.

The simulations were performed only for the d and e flood waves because of the significant computational time (see Figure 3d,e). The initial model setup (the flow field, water levels and the bed material) was given for each run from the results of the model runs for the first three flood waves by the *Re\** dependent combined method. The flood wave d was a relatively low (peak is around 3460 m3/s), but durable (~2 months long) flood wave, while the e was the historical one with a peak higher than 10,000 m3/s. Thus, the comparative analysis presented the operational characteristics of the sediment transport models for both the durable lower and also for the extreme water regimes. As a benchmark, the measured bed change map indicated the extent of the possible bed changes. However, the measured and calculated maps could not be compared directly, because the measured belonged to the whole two-year-long period (Figure 3).

The calculated bed change maps for the flood wave d (Figure 3) are presented in Figures 6–8.

**Figure 6.** Calculated bed changes by the van Rijn (*vR*) formula for a 2.5 month-long period (Figure 3).

**Figure 7.** Calculated bed changes by the Wilcock and Crowe (*W&C*) formula for a 2.5 month-long period (Figure 3).

**Figure 8.** Calculated bed changes by the combined method for a 2.5 month-long period (Figure 3).

As the results show, the *vR* model estimated unrealistic changes both spatially and in magnitude. The unrealistically huge bed changes suggest that the *vR* model did not seem to be an appropriate model for the given Danube reach, particularly not for bed change calculation in the main channel.

The *W&C* sediment transport model estimated a more stable bed surface than the *vR* (Figure 7). In this particular case, the bed surface seemed so resistant that the mean flow field was too weak to cause any significant bed changes. The motion of the very fine, basically suspended inlet load was calculated by the *W&*C model as bed load. Therefore, that part of the inlet sediments settled progressively along the channel. However, because of the quite low suspended load [51], the bed level rise caused by sedimentation was negligible (<0.005 m).

Figure 8 shows the bed changes calculated by the combined method. The red lines illustrate the borderline which separated the areas where the *vR* or the *W&C* model was activated in the initial moment of the model run. Accordingly, it can be seen that the *vR* formula was invoked at the near-bank areas, at the groin fields and also at a smaller part of the Vének lower gravel bar. At these regions, more significant (~0.05 m) sedimentation was estimated. That is, at these less hydraulically rough parts of the river bed, the deposition of both the finer bed load and suspended load were expected, which areas can be detected by the *Re\** [22]. In turn, in the navigation channel, no considerable bed change happened. According to the suspended form of the *vR* formula, the finer suspended load passes over the calculation domain, while the bed surface remains still, calculated by the *W&C* formula. This assumption is consistent with the conclusions of the field measurements [19]; the main channel seems to be armored enough to be resistant at the mean water regime.

The other flood wave, for which the comparative analysis of the three calculation methods (*vR*: Figure 9, *W&C*: Figure 10 and combined: Figure 11) was established was the historical flood wave from 2013 (Figure 3). Regarding this hydrological case, the *vR* model estimated also an unrealistic bed change map (Figure 9.). This was mainly true for the main channel. There, the motion of the coarser grains was probably overestimated, resulting in huge erosions and depositions. In turn, at the near-bank regions, at gravel bars and at the groin fields, the changes seemed to be partly in the expectable order of magnitude. However, it is clearly visible that the *vR* formula is not an appropriate choice for the morphological change calculation of such a complex river reach.

**Figure 9.** Calculated bed changes by the *vR* formula for the historical flood wave (Figure 3).

**Figure 10.** Calculated bed changes by the *W&*C formula for the historical flood wave (Figure 3).

**Figure 11.** Calculated bed changes by the combined method for the historical flood wave (Figure 3).

The *W&C* model calculated more realistic bed changes (Figure 10), especially in the main channel. However, the measurements at the lower gravel bar and also at the whole main channel showed significant (Δ*z* > ±0.2 m) changes. Based on these, the *W&C* formula likely overestimated the stability

of the channel. A remarkable bed level increase can be pointed out, which moderated towards the downstream direction. These bed level changes could be explained by the settling of the inlet finer load, which could not be taken into account as suspended load. The publication of Török et al. [59] pointed out that in the case of mixed bed content, the Shields diagram predicts lower critical bed shear stress for the bed load of the finer, sand particles, than the reference shear stress of the *W&C* model. Accordingly, the *W&C* model estimates respectively higher stability for the sand particles, than the Shields diagram and thus the *vR* model. This leads to such kind of deposition pattern in the main stream which is unrealistic compared to the measured bed changes map [19].

The combined method predicts (Figure 11) more significant bed changes, compared to the one resulted by the *W&C* model. The red borderline suggested that as at the flood wave d, the combined method calculated the sediment transport by the *vR* formula in the near-bank regions. During the flood wave e, the remarkable erosion at the groin field B meant that the groin field got flushed and the earlier deposited finer sands got eroded. Also, notable changes took place at the vicinity of the gravel bar at region D. Here, the widening of the downstream sides of the gravel bar can be seen, in accordance with the measured bed change map in Figure 4. Around the groin pair at the left bank, on the opposite side of the gravel bar, the blue spots refer to erosion. This process was probably the result of a similar, flushing process like the one which happened at the upstream groin field. In turn, in the main and navigation channel, no considerable bed change happened. According to the suspended form of the *vR* formula, the finer suspended load passed over the calculation domain, while the bed surface remained still, calculated by the *W&C* formula. The conclusions of the field measurements [19] referred also to the resistant main channel.

The results indicate that the interactions between the sediment transport at different channel sections (groin fields, gravel bars and main channel) could not be estimated by the *vR* or *W&C* formulas. However, the expedient combination of them gives an opportunity to deal with the interaction-mechanism between the local- and reach-scale processes.

Analyzing the calculated bed changes in the marked boxes, the following assumptions can be stated. In this part, the results of the *vR* model were skipped because of the unrealistic bed change calculations. At region A, the real bed level deepening (Figure 4) could not be reproduced by any method. In region B, only the combining method was able to predict significant depositions and erosions. Accordingly, the depositions probably took place during the lower water regime, while the bed level incision occurred during the flood waves. Thus, the measured bed level changes during the two years were formed most likely indeed during the whole two-year-long period. In region C, all two model results suggested that the gravel bar was in a stable state. Besides region A, the bed level increase in the main channel between region C and D (Figure 4) could not be pointed out by any sediment transport formula. Finally, in region D the model results showed that the widening of the gravel bar and also the erosion at the vicinity of the groin pair occurred rather during the higher flood wave than during the slighter flood waves, or mean water regime.

#### *4.2. Measured Data-Based Verification of the Combined Method*

A quantitative assessment of the tested sediment transport formulas was performed based on the results of the comparative analysis. First, in accordance with the conclusions of the field measurement based investigation, it was assumed that the measured erosion and deposition took place mainly during higher water regimes in region D, when the flow discharge was higher than the bed forming discharge [19,23]. Thus, the measured and calculated data could be considered as indirectly comparable in this region. Therefore, from the measured bed level change maps, the total volume of the erosion and deposition could be calculated. Furthermore, counting the number of days of the higher water levels, the average daily rate of the volume of both the erosion and deposition could be estimated. Likewise, based on the calculated bed change maps and knowing the duration of the historical flood wave, the daily average volume changes of the deposition and erosion could be estimated.

Table 2 shows data about these volumes. The table presents the ratio of the calculated volume to the measured, regarding the deposition and erosion separately, for each sediment transport formula. A value of 1 would indicate a perfect match to the measured volume change. Table 2 shows that the *vR* model was the one which overestimates most both the deposition and erosion volumes. The *W&C* model calculated the deposition rate more accurately, but still indicated more than the measured. This was partly explained by the lack of the suspended sediment calculation. In turn, the *W&C* model estimated negligible low erosion. In total, the more reliable results were provided by the combined method. With this, the erosions at the near-bank parts were calculated more accurately by the *vR* model, resulting in sediment feed for region D. Thus, because of the capturing of the coming sediments, the widening of the deposition could be better represented. As the *Re\** dependent criterion activated the *vR* formula at the groin pair, the bed level incision at its vicinity was also estimated better.

**Table 2.** The average daily volume changes for calculated (Δ*Vc*) and measured (Δ*Vm*) volume ratio values Δ*Vc*/Δ*Vm* for region *D*. A value of 1 would indicate the perfect match.


Even though the measured and calculated bed changes could not be compared directly, the nature, the magnitudes and the locations of the remarkable bed changes suggested the greater aptitude of the combined method.

As it was discussed in Section 3.2.1, the calculated bed changes were elaborated only for the higher flow discharges (>2100 m3/s, Figure 3), so the results could not be compared with the measured changes directly. Therefore, to achieve a notionally common scale, the bed change values were normalized. That is both the measured and calculated bed changes got divided by the highest bed level decrease or increase the value of the main channel:

$$
\Delta z\_{\text{norm}} = \left\{ \begin{array}{c} if \, \Delta z > 0 \to \frac{\Delta z}{\Delta z\_{\text{max}}} \\ \epsilon lse \to \frac{\Delta z}{|\Delta z\_{\text{min}}|} \end{array} \right\}. \tag{3}
$$

Thus, the occurring values develop between −1 and 1, where 1 indicates the maximum deposition height along the river reach, while −1 belongs to the biggest erosion (Δ*zmax*, *meas* = 1.5 *m*, Δ*zmin*, *meas* = 1.8 *m*; Δ*zmax*, *calc* = 0.5 *m*, Δ*zmin*, *meas* <sup>=</sup> 0.25 *<sup>m</sup>*). The measured bed changes of both the erosions and depositions were consequently higher than the calculated. That is, the numerical model underestimated the magnitudes of the bed changes. This could be partly explained by the ignoring of the third part of the annual bed load in the numerical model estimation.

By the comparison of the measured (Figure 12) and calculated normalized bed changes (Figure 13) the following remarks can be stated. The modeled bed changes did not represent the remarkable erosion at region A. It is noted that this difference also meant a sediment supply loss for the downstream in the model calculation. Since the model did not manifest any bed level decrease, the bed material was probably finer here than it was set in the model. At region B, the measured scours appeared in the calculated results (blue spots). However, not as concentric scours, but rather as lengthwise formations. The numerical model represented depositions close to the measured magnitude at region B also. However, their location iwa not accurate; the brown spots occur between the groins, instead of in the front of the groins, like in Figure 12.

**Figure 12.** Normalized bed changes of the measured values (Figure 4) regarding the whole period October 2012–October 2014.

**Figure 13.** Normalized bed changes of the calculated values. The calculation was elaborated for the eight flood waves in the period October 2012–October 2014 (Figure 3).

At the gravel bar in region C, negligible bed changes were measured. Likewise, the numerical model predicted a stable bed surface. Finally, the combined method pointed out the growing of the lower part of the Vének lower gravel bar (region D, Figure 13). The measured bed changes indicated two, separable depositions. These double depositions were also represented by the model results. Like the measured changes, the model also calculated scouring in front of the left bank groin pair. However, the extension of it was considerably lower than the scale of the measured deepening. In turn, the lengthwise deposition between the groins (between regions C and D) was also indicated in both measured and calculated maps. However, the calculated deposition occurred in a significantly lower range and formed at the right bank side and not in the main stream. Also, an important match was that neither the measured nor the calculated values suggested any essential large-scale bed changes in the main channel.

A significant difference between the measured and calculated bed changes happened in region A. Here, the effect of the potential error in the initial bed material was further examined. An investigation

was performed, which was based on the assumption that the initial bed material was set inaccurately around region A. The bed material samples around *d*<sup>50</sup> ≈ 0.01 m. Considering the grain-size distributions [19], a still realistic, but considerably finer bed material was presupposed. Therefore, the model was set up by 30% lower *d*50, which was *d*<sup>50</sup> ≈ 0.007 m. With this only one difference, the model was run for the historical flood wave. Figure 14 presents the bed changes at region A, and at the downstream of it.

**Figure 14.** Calculated bed changes by the combining method for the historical flood wave. The model in the middle figure was set up with the initial, while in the right figure with finer bed material.

The right side of Figure 14 shows the bed changes in the case of finer bed material. It can be seen that the 30% decrease of the *d50* resulted in major erosion at region A. Considering the measured changes in the left figure it is obvious that the decreasing of the *d50* led to a better match to the real bed changes. The lower row of Figure 14 represents the bed changes at the downstream. Here, important deposition formations could be measured (left Figure 14) in front of the right bank groin pair (region B) and also in the main stream, between the two gravel bars (between region B and D). These changes could not be represented by the original model setup (middle Figure 14). In turn, in the case of finer bed material (right Figure 14), the model predicted important depositions in these regions. The extension of the deposition in front of the groin pair (region B) was very similar to the measured. Also, the lengthwise extension of the deposition downstream (between regions B and D), in the main stream were also reproduced. However, the location of it was not correct. It seems that the model underestimated the crosswise sediment transport, which is a known limitation of the Reynolds averaged description of the flow field [56–58,60]. Concluding, the herein presented investigation suggested that the bed material at region A was finer than the predicted *d50* allocation for the original model setup. With this assumption, the deposition nature at the downstream has become also detectable.

#### **5. Discussion**

There are existing proven models in the literature that work reliably for given morphodynamic cases (e.g., the van Rijn model for hydraulic smooth regimes, which mainly occurs in clear sand bed; or the Wilcock and Crowe model for hydraulic rougher regimes, which develops in coarser bed surfaces). However, there are river reaches, within which hydraulic smoother and rougher zones can form, resulting in a completely different type of sediment transport, such as bed armoring and silt aggradation parallelly. For such cases, Török et al. developed a novel sediment transport calculation method and validated it against laboratory measurements.

In this study, its validation against field measurements was introduced. The results showed that the combined application of the Wilcock and Crowe and the van Rijn models could significantly increase the precision of numerical bed change calculations. The results showed that the combined method calculated both erosion and deposition the most reliably. It estimated the volume of the scouring by order of magnitude more precisely than the other two models. Additionally, the results also evinced that the location, extension and shape of the bed formations (e.g., scours, deposition, bar evolvement) could be calculated more reliably by the novel combined approach. These results were consistent with the laboratory measurements based on validation of the calculation method.

The sediment transport interaction between the sand and gravel dominated parts of a river reach played an essential role in the bed changes. For instance, the flushing of the silt dominated groin fields can feed the downstream gravel bars, leading to the growth of the bars. The application of the combined approach really makes sense, within river reaches, where such well-separated morphodynamic situations prevail (e.g., sand aggradation in groin fields; bed armoring in the main stream). In such cases, the local-scale morphodynamic processes are calculated by the proven models, while the interaction between them leads to a more accurate reach-scale calculation. The comparison of the numerical simulations to the measured data presented that such reach-scale interactions can be taken into account by the novel combining method.

#### **6. Conclusions**

The validation against both laboratory and field measurements suggests that not necessarily the development of completely new models can lead to the evolution of the sediment transport calculation. This study highlights that the optimally combined use of sediment transport models can be a promising alternative in the sediment transport calculation. In this context, it is important to state that regardless of whether a sediment transport model is applied in combination or simply, the determination of its applicability limits by any morphodynamic variables in a well-defined way is vital. In this study, a shear Reynolds number based description was applied which allowed for the authors an easily described combination criteria. However, the essence of the combining theory does not insist on any variable for the applicability limit description. However, the herein presented example highlights the need for such a kind of description method and also the potential of the parallel use of any sediment transport models.

**Author Contributions:** Conceptualization, methodology, software, validation, writing—original draft preparation, visualization, G.T.T., resources, data curation, G.T.T. and S.B.; formal analysis, writing—review and editing, S.B.; supervision, project administration, funding acquisition, S.B. and J.J.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by MTA TKI (Office for Research Groupsof the Hungarian Academy of Sciences). Support of grant BME FIKP-VÍZ by EMMI (Ministry of Human Capacities) is also kindly acknowledged. The authors acknowledge the funding of the OTKA (Hungarian Scientific Research Fund) FK 128429 grant. The research was partly supported by the ÚNKP-19-4 New National Excellence Program of the Ministry for Innovation and Technology. The last author acknowledges the funding from the János Bolyai fellowship of the Hungarian Academy of Sciences.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

#### *Appendix A.1 Wilcock and Crowe Bed Load Transport Formula*

The well-known surface-based bed load formula of *W&C* [8] was developed for mixed gravel/sand sediments. The gravel ranged in size from 2.0 to 64 mm, the sand from 0.5 to 2.0 mm. They conducted laboratory experiments in a tilting laboratory flume, in which both the water and the sediments were recirculated. Five different sediment mixtures were investigated, varying the sand content of the sediment mixtures from 6.2% up to 34.3% [60]. The mixture grain-size distribution used in the experiments of Török et al. [34] and of Wilcock et al. [60] shows a good match. In case of the experiments of Wilcock et al. [60], *dm* grain size of the initial bed materials was in the range of 4.1–10.5 mm, whereas *dm* was 8.7 mm the measurements of Török et al. [34]. The *W&C* formula is therefore expected to be capable of describing the gravel/sand bed load transport processes in our experiments. In the case of such non-uniform bed material, the interaction between (hiding and exposure) the particles of different sizes plays an essential role in the stability of the sediments and can result in bed armoring. The *W&C* formula is reported to be capable of predicting the transient conditions of bed armoring or scouring processes. The preliminary numerical simulations of Török et al. [46] also supported this statement.

The transport formula uses the collapse similarity hypothesis, which estimates the fractional sediment transport rate (*wi\**) as the function of the ratio of the shear stress (τ) and the reference shear stress (τ*ri*). The reference shear stress τ*ri* is the value of τ at which *wi\** is equal to a small but already perceptible value *wi\** = 0.002. That is the reference shear stress has a similar physical meaning to the critical shear stress (τc), which indicates the initiation of motion of particles. The dimensionless transport rate is obtained according to the following equation:

$$\mathcal{W}\_i^\* = \begin{cases} 0.002 \theta^{7.5} & \text{if } \theta < 1.35\\ 14 \Big( 1 - \frac{0.894}{\theta^{0.5}} \Big)^{4.5} & \text{if } \theta \ge 1.35 \end{cases} \tag{A1}$$

where

$$
\theta = \pi/\tau\_{\dot{n}}.\tag{A2}
$$

The volumetric transport rate per unit width *qbi* can be expressed from the following formula:

$$\mathcal{W}\_i^\* = \frac{(s-1)\mathcal{g}q\_{hi}}{\mathcal{F}\_i \mu\_\*^3},\tag{A3}$$

where *Fs* is the proportion of sand in surface size distribution.

For the estimating of τ*ri*, the *W&C* model suggests that the stability of the sediment particles depends on the sand content. Thus τ*ri* can be calculated by the following steps:

$$
\tau\_{rm}^\* = 0.021 + 0.015e^{-20F\_S},\tag{A4}
$$

where τ\* *rm* is the dimensionless reference shear stress for the mean size of the bed surface. The reference shear stress for the mean size can be obtained from the particle Froude number equation:

$$
\tau\_{rm}^\* = \frac{\tau\_{rm}}{\mathcal{g}(\rho\_S - \rho\_W)d}. \tag{A5}
$$

Then, the last step is the estimation of the reference shear stress for the particle sizes separately:

$$\frac{\pi\_{ri}}{\pi\_{rm}} = \left(\frac{D\_i}{D\_{sm}}\right)^b,\tag{A6}$$

where

$$b = \frac{0.67}{1 + e^{(1.5 - \frac{D\_i}{D\_{\rm sun}})}}.\tag{A7}$$

#### *Appendix A.2 Van Rijn Bed Load Transport Formula*

The development of the *vR* formula is based on laboratory and field measurements. Van Rijn presented a mathematical method for the computation of bed load transport rate. The formula is basically the numerical solution of the equations of motions for a solitary particle. The mathematical model was calibrated based only on the experiments of Fernandez Luque [61,62], where the motions of uniform particles (*d* = 0.18 mm) under steady flow (*u\** = 0.04 m/s) were monitored.

The lifting term of the equation of motion for a particle is derived as [63]:

$${}\_{L}F\_{L}(shar) = \alpha\_{L} \rho \nu^{0.5} D^{2} v\_{r} \left(\frac{\partial u}{\partial z}\right)^{0.5},\tag{A8}$$

where *FL* is the lift force, α*<sup>L</sup>* is the lift coefficient, ν is the kinematic viscosity coefficient, ρ is the density of the fluid, *vr* is the relative particle velocity and ∂u/∂z is the velocity gradient. The mathematical method uses the turbulent wall law Equation (A9) for calculation of the vertical flow velocity distribution:

$$
\Delta U(z) = \frac{u\_\*}{\kappa} \ln \left( \frac{z}{z\_0} \right) \tag{A9}
$$

where *U(z)* is the horizontal velocity as the function of the *z*, *z* is the level above the bed, *u\** is the bed shear velocity, κ is constant von Karman (κ = 0.41), and *z0* is the zero-velocity level above the bed.

The volumetric bed load transport rate is defined as:

$$\frac{q\_b}{[(s-1)g]^{0.5}D\_{50}^{1.5}} = 0.053 \frac{T^{2.1}}{D\_\*^{0.3}}\,\text{\AA}\tag{A10}$$

where *qb* is the bed load transport per unit width, *D50* is the particle size, *s* is the specific density, *g* is the acceleration of gravity, *T* is the transport stage parameter, can be represented as:

$$T = \frac{\left(u\_\star'\right)^2 - \left(u\_{\star,cr}\right)^2}{\left(u\_{\star,cr}\right)^2},\tag{A11}$$

where *u\*,cr* is the critical bed shear velocity according to Shields [35], *u* <sup>∗</sup> <sup>=</sup> *g*0.5/*C u* is the bed shear velocity related to grains, where *C´* is the Chézy-coefficient related to grains and *u* is the mean flow velocity.

In the Equation (A10), *D\** is the particle parameter defined as:

$$D\_\* = D \text{50} \left[ \frac{(s-1)g}{\nu^2} \right]^{1/3} \tag{A12}$$

Finally, the *vR* formula was verified based on 56 field measurements and 524 flume data. The bed load formula is recommended to use for particles in the range of 200–2000 μm (i.e., below 2 mm). As such the formula is not expected to reliably calculate the motion of gravels and the interaction between particles of different sizes and to describe such a process, like bed armoring. However, the model can presumably describe the transport of finer grains.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
