carried out, it is far from evident what sort of empirical sediment transport formula to use. There can be found formulas for finer fractions (e.g., [20]), as well as for coarser fractions (e.g., [21]), both having limitations when considering mixed-size bed material. However, no such straightforward method can be found in the literature, which suggests the locally dominant bed load transport mode and thus the more suitable formula for nonuniform bed material. An attempt is made in this study to fill this gap, providing a methodology that works more accurately compared to earlier published methods.

## 2. Materials and Methods

### 2.1. Problem Statement

Significant variability of the dominant grain size in the bed load can take place in nonuniform bed material, e.g., in the case of the investigated Danube reach. Figure 2 shows the density function of a bimodal bed material sample, collected from the lower side bar (see Figure 1, green columns, left vertical axis). The black line (and the right vertical axis) represents the critical dimensionless bed shear stress values (or critical Shields number, Equation (1)), as the function of the grain size, calculated by the Shields curve [22]. If the dimensionless bed shear stress exceeds the critical value, the given particle is in motion, otherwise it stays calm.

$$ \tau^* = \frac{\tau}{g(\rho_s - \rho_w)d'} \qquad (1) $$

where $τ^*$ is the Shields number, $τ$ is the bed shear stress, $g$ is the acceleration due to gravity, $ρ_s$ is the sediment density, $ρ_w$ is the water density and $d$ is the grain diameter.

<figure><img src="image_112.png" /><figcaption>**Figure 2.** The density function of a bed material sample collected from a side bar (green curve). The Shields curve-related critical bed shear stress values are depicted by the black line.</figcaption></figure>

As the fraction rates indicate, sand and gravel also appear well separately, both with significant amount (>30%). The grain size distribution indicates a $d_{50}$ of 7 mm, which would refer to gravel transport. However, examining the critical Shields numbers, it turns out that the dominant bed load grain size seems to be inconstant, varying depending on the flow conditions. If the dimensionless bed shear stress is in a lower range ($τ^* < 0.05$, or equivalently $τ < 2 N/m^2$), only sand (<2 mm) movement can be expected. In turn, increasing dimensionless bed shear stress ($τ^* > 0.05$, or equivalently $τ > 2 N/m^2$) results in the motion of even the coarser gravels (>2 mm), while the motion of the finer sand fraction tends to be rather suspended, like rolling [8], resulting in even coarser bed load. Concluding the above, the $d_{50}$ grain size of the bed material itself is not able to indicate the current dominant grain size in the bed load, especially not in the case of nonuniform bed material under varying flow conditions (for further similar study cases see, e.g., [23,24]).
When implementing computational morphodynamic models in cases like the above-presented mixed-size bed material, the variable dominant fraction in the bed load can lead to inappropriate application of sediment transport models. Nevertheless, the numerical simulation of morphodynamic processes can contribute to the better performance of river engineering activities, e.g., when planning restoration measures, and it is considered a widely applied investigation tool with many benefits. However, the reliable numerical description of the interaction between river hydrodynamics and the morphological features is a challenging task in rivers with uniform bed material composition, but shows even higher difficulties in the case of mixed-size bed material conditions (e.g., [25–29]). When the grain sizes cover a smaller range in the bed material, the potential erosion results also in a quite uniform grain size range of the bed load transport. For such conditions, many bed load sediment transport models have already been developed and successfully validated [22,30–33]. On the other hand, the choice of the bed load sediment transport formula is still a challenge when the characteristic grain sizes can significantly change both in time and space [23,25,34–36]. Accordingly, the most suitable numerical formula can change in time and space. For instance, as many studies present [34,37–44], the van Rijn model [20] could provide better estimation of sediment transport in sand-dominated zones, whereas the Wilcock and Crowe model [21] suits gravel-dominated mixed-size cases better. Török et al. [34] also proved that the combined and parallel application of these two models provides better results when simulating sediment transport in mixed-size bed material. In the referred study, a physical indicator was defined as a parameter to decide which formula is to be applied. It is, however, not evident what this parameter should be. In such cases, the prediction precision of the dominant fraction in the bed load transport plays an emphatic role in the accuracy of the sediment transport calculation possibilities.

## 2.2. Existing Method for Bed Load Classification

One applicable method, the Shields–Parker river sedimentation diagram, can be used as a comprehensive bed-material and sediment transport nature classification tool for a particular river reach, with uniform bed content [7]. The diagram indicates the Shields number ($\tau^*$, or equivalently $\tau^*_{bf}$, Equation (1)) based on the $d_{50}$ specific mean grain size versus the so-called explicit particle Reynolds number ($Re_p$, Equation (2), where $d = d_{50}$ in the case of the Shields–Parker diagram). The diagram includes the well-known Shields curve (continuous black curve), which shows if the given particle is in motion (points above the curve), or remains still on the bed surface (points below the curve) [22].

$$Re_p = \frac{\sqrt{gRdd}}{v}, \qquad (2)$$

where $R = (\rho_s - \rho_W)/\rho_W$, $\rho_s$ is the water, $\rho_W$ is the sediment density, and $v$ is the kinematic viscosity of the fluid.

Parker stated [8] that alluvial rivers can usually be divided into two types based on the function of particle Reynolds number ($Re_p$, Equation (2)): sand-bed and gravel-bed rivers. According to Parker, a distinction can be defined based on $Re_p$, which actually depends only on $d_{50}$ (therefore it can, in fact, be considered as the dimensionless substitution of the grain size). That is, considering that $d = 2$ mm is the threshold value between the sand and gravel fractions, the critical value for $Re_p$ also can be calculated according to Equation (1), where $d = 2$ mm results in $Re_p = 360$ ($Re_p > 360 \rightarrow$ gravel bed, $Re_p < 360 \rightarrow$ sand-bed), which critical value is depicted by a vertical narrow in Figure 3, within the $Re_p$ range of 100–1000.
<figure><img src="image_113.png" /><figcaption>Figure 3. Shields–Parker diagram [7].</figcaption></figure>

Consequently, the method basically yields a grain size-based classification method for near-uniform bed materials. This simplified approach is not capable of classifying characteristic local morphodynamic processes for complex situations in terms of river bed morphology and inhomogeneous bed material composition [7] (e.g., in the case of bed material in Figure 2).

The basic idea here is to keep the simplicity of the classification method shown by Parker, but at the same time to involve a new indicator, which accounts for the local hydrodynamic effects that better characterize the interaction between the water flow and the river bed, i.e., the bed shear stress. Thus, based on literature review, field and laboratory measurements, a conjecture is introduced in this paper. For a better understanding of the herein-proposed methodology, the most relevant findings of the related publications are summarized in the following.

## 2.3. Introducing the Novel Method for Bed Load Classification

Shields published a well-known study regarding sediment transport investigation in 1936 [22]. The research work investigated the critical state of the bed-forming particles, applying a mathematical model of forces acting on a solitary particle. Shields pointed out that the rate of the acting forces is not constant, but the relation of the effective force of the flow to the resistance of a grain can be expressed with a universal function of the shear Reynolds number ($Re^*$, Equation (3)). The interpretation of the shear Reynolds number is expressed in Equation (3): the ratio of grain size to the thickness of the laminar boundary layer.

$$Re^* = \frac{u_*d}{\nu} = \frac{d}{v/u_*} = \frac{\text{Grain diameter}}{\text{Thickness of viscous layer}} \quad (3)$$

where $u^*$ is the shear velocity ($u_* = (\tau_p)^{0.5}$, where $\tau$ is the local bed shear stress). The following table demonstrates the difference between the shear Reynolds number and the explicit particle Reynolds number. The constant and the variable parameters in Equations (2) and (3) are depicted in Table 1. The last column indicates the ratio of $Re^*$ to $Re_p$. The main outcome of the table is that $Re^*$ implicates not only the grain size ($d$), but also a term regarding the acting forces of the flow ($u^*$). This remark confirms the $Re^*$-based investigation of Shields, instead of the $Re_p$-based method.
<table><caption>**Table 1.** The constant and the variable parameters in the equations of the shear Reynolds number (Re*) and the explicit particle Reynolds number (Re<sub>p</sub>) (Equations (2) and (3)).</caption><thead><tr><th></th><th>Constants<sup>1</sup></th><th>Variables</th><th>Re*/Re<sub>p</sub></th></tr></thead><tbody><tr><td>Re*</td><td>v</td><td>d, u*</td><td></td></tr><tr><td>Re<sub>p</sub></td><td>v, g, R</td><td>d</td><td>u*/√d</td></tr></tbody></table>

The Shields curve, actually, shows the relation between the shear Reynolds number and the acting forces [22,45]. The investigation proves that the acting forces and the stability of a grain depend not on the grain sizes, but on the hydraulic regime formed directly above the grains [7,46,47] (Considering the possible development of different kinds of turbulent flow pattern, e.g., in sand-bed rivers bed-form-induced turbulent flow, while in gravel-bed rivers, friction-induced turbulent pattern can occur [20,32,48–52]). It leads to two important consequences: (i) as the effective acting force depends not only on the grain size, but also on the acting forces of the flow, it is suggested that the sediment transport natures need to be classified based on the shear Reynolds number, too; (ii) due to this, the validation limits of the sediment transport-related mathematical models (that is, the sediment transport models) should be defined as the function of the shear Reynolds number.

The investigation of the applicability of the van Rijn bed load transport formula [20] also enhances the benefits of the shear Reynolds number-based investigation, against the grain-size-based (or characteristic grain-size-based). As one of the most widely tested and applied sediment transport formulas for river engineering applications, the van Rijn bed load transport equation [20] considers the application of the Shields curve. The transport model was developed and validated for sand transport measurements, and, therefore, the application limit was defined accordingly ($d < 2$ mm) [7]. However, the validity range of the base mathematical model can be defined by the shear Reynolds number, and the grain size is not a suitable variable for this reason. Namely, the drag coefficient ($C_D$), the lift coefficient ($\alpha_L$), or the effect of the spinning motion depends on the shear Reynolds number [53–56], and the recommended values of these variables in the van Rijn model were calibrated and validated for lower shear Reynolds numbers ($Re^* < 400$) [20,57,58]. Therefore, the determination of the applicability of the van Rijn model is reasonable as the function of the shear Reynolds number, which is the hydraulic smoother ($Re^* < 400$) regime.

At the same time, the lower $Re^*$ range implies lower bed shear velocity ($u^*$) and smaller particles ($d$) (see Equation (3)), assuming the domination of the less-resistant finer fractions in the bed load and stable coarser grains. This is proved by the about 600 measurement data [20,57,58]. Also, Parker mentions [8] that, at rougher hydraulic regime ($Re^* > 1000$), the motion of the sand fractions is rather suspended like rolling. Consequently, finer fractions disappear from the bed load and the coarser gravels remain and become dominant. Accordingly, the van Rijn bed load model is not applicable in rougher hydraulic ($Re^* > 1000$) regime.

The laboratory experiments, where coarser ($d_{50} > 5.3$ mm) and nonuniform bed material conditions were investigated by Wilcock et al. [26], are also taken into consideration here. Figure 4 shows the rate of the gravel load to the sand load as the function of the shear Reynolds number. Within any mixtures, the higher the shear Reynolds number, the coarser the composition of the bed load. It is visible that the gravel load is equal to the sand load when the $Re^*$ is between 300 and 400, around 350 (Figure 4). Above this, the gravel transport is dominant. In turn, below 350 the sand transport prevails, which underlines the earlier assumption, i.e., the fine sediment bed load transport takes place and dominates in the smoother hydraulic range. Consequently, despite that the $d_{50}$ characteristic grain size indicates only gravel transport, sand-dominated sediment transport was also found.

<footnote n="1">¹ The values of the constant terms are: v = $10^{-6} \frac{m^2}{s}$, g = $9.81 \frac{m^3}{s^2}$, R = 1.65.</footnote>
<figure><img src="image_114.png" /><figcaption>**Figure 4.** Calculated rates of the gravel load to the sand load as the function of the shear Reynolds number, based on the laboratory experiments of Wilcock et al. [26].</figcaption></figure>

Based on the above-described summary, the following findings are emphasized

* Shields has already presented [22] that the shear Reynolds number ($Re^*$) is a suitable variable (instead of the grain sizes) to define the relation of the acting forces and thus to classify the sediment transport in the near-surface zone. In other words, the effective force of the flow (together with the resistance of the grain) must be taken into account for the proper investigation of the stability and mobility of bed-forming particles.

* The applicability of both mathematical models of Shields [22] and van Rijn [20] can be described by the shear Reynolds number, but not by grain sizes. This also supports the fact that the shear Reynolds number-based classification of the sediment transport is expected to be more reliable than a grain-size-based method.

* The investigation of the measurements, for which the van Rijn bed load formula was validated and calibrated [20,57,58], suggests that the sand bed load can occur at lower $Re^*$ range, $Re^* < 400$. Complementing this by the statement of Parker [8], it is assumed that in rougher hydraulic regime ($Re^* > \sim 1000$), the sand transport appears only in suspended form, while the bed load consists mainly of gravel.

* The $Re^*$-based assessment of the laboratory experiments of Wilcock et al. [26] shows that the $Re^*$ range ~3–400 seems to be a critical zone, where the yield of the sand and gravel are roughly the same. A characteristic grain size cannot refer to such phenomena.

Accordingly, the following conjecture was stated: The $Re^*$ parameter is a more suitable parameter for the prediction and distinction of the sand- or gravel-dominated bed load transport, than either simply a characteristic grain size ($d_{90}$ or $d_{50}$), or the $Re_p$. Moreover, $Re^*$ between ~300 and 400 seems to be a critical zone: above-gravel, below-sand domination can be expected in the bed load transport.

## 2.4. *Description of the Validation Measurements*

In order to investigate accuracy of the conjecture, an assessment considering characteristic local bed grain sizes, local bed shear stress values as well as calculated Reynolds numbers was performed using datasets from recent laboratory and field experiments of the authors. First, the results from the laboratory experiments of Török et al. [59] were used. In this experiment, morphological processes (scouring, bed armoring, and aggradation) were investigated around a single groin, using mixed-size bed material with an initial $d_{50}$ of 5.16 mm. The local grain size distributions were determined by an automated image-based grain detection software tool BASEGRAIN [60]. Local bed shear stress values were estimated according to the TKE-Method [61] using the near-bed point-velocity measurements,
carried out with a 3D Acoustic Doppler Velocimeter (ADV) [62]. The measurements were carried out at quasi-equilibrium bed geometry. Then, negligible bed changes were measured in the flume ($\Delta z < \pm 1$ mm/h); however, weak bed load transport ($Q_b < 0.3$ kg/h) of mainly the sand particles ($d_{50} = 1.3$ mm) could still be observed at the outlet of the flume, where the $d_{50}$ of the bed material ranged between 4 and 6 mm [59].

Second, field data from a section of a large river was also assessed. The field experiments were carried out in the Hungarian section of the Danube River at mean flow regime, where the main characteristic morphological parameters were already introduced in the Introduction point. Forty-seven bed material samples were taken by a drag-bucket sampler, for which the grain-size distributions of the samples were determined by sieving analysis [63]. Parallel to the local bed material samplings, time-averaged flow velocity profiles were measured by an acoustic Doppler current profiler (ADCP) (WorkHorse Rio Grande 1200 kHz) at the sampling locations. Based on these long-term, fixed-boat measurements, the local bed shear stress values could be estimated using the turbulent wall law [64]. No bed load sampling could be performed during the two measurement campaigns. However, the bed load rating curve at the closest monitoring station (6 km upstream) [17] represents that low, but not negligible, bed load transport took place during both measurements ($Q_b < 0.1$ kg/s).

The data regarding the field and laboratory measurements can be found in Appendix A (Table A1).

## 2.5. Measurement Data Processing for Validation Purpose of the Stated Conjecture

Field and laboratory data from the above-introduced reach of the Danube River in Hungary and laboratory experiments were conducted to confirm the stated conjecture. In order to investigate the shear Reynolds number-based classification potential, the samples were sorted into three classes according to their origin place (Figures 5 and 6): the black squares refer to the sandy bed material samples, which were collected from the near-bank zones. Here, the bed material consists predominantly of silt and sand, gravel transport can be expected only during flood waves. The red stars originate from the side bars (e.g., bed material illustrated in Figure 2). As it was introduced in Chapter 2 and based on Figure 2, the dominant grain size in the bed load is variable (both sand and gravel can dominate) and it is strongly flow-dependent. Finally, the samples from the main channel are depicted by the blue diamonds. Based on the bed material analysis, only gravel transport can occur at this region of the channel, because sand can be vanishingly found. However, as the measurements indicate stable bed surface, pure gravel transport can be expected at mean flow conditions. The investigation examines whether the shear Reynolds number is a proper variable to separate the different classes of the points.

<figure><img src="image_115.png" /><figcaption>**Figure 5.** The bed shear stress as the function of the $d_{50}$ grain size, regarding field and laboratory measurements [17,64].</figcaption></figure>
<figure><img src="image_116.png" /><figcaption>**Figure 6.** (A) Shields number as the function of the shear Reynolds number, based on laboratory and field data [17,59,64]. Shields curve is indicated with black line. (B) Shields number as the function of the explicit particle Reynolds number, based on laboratory and field data [17,59,64]. Shields curve is indicated with black line.</figcaption></figure>

Furthermore, the comparison of the shear and particle Reynolds number-based investigations was also examined. First, the scattering of the three point classes as the function of *Re*\* and *Re*<sub>p</sub> were compared. Based on this comparison, a more accurate separation method becomes outlined visually.

Finally, in order to perform a quantitative assessment on the performance of the *Re*\*- and *Re*<sub>p</sub>-based approaches, a statistical method was applied, which provides information on the representativity of the Reynolds number ranges for the characteristic bed load transport processes. Here, the log-normal distributions of the *Re*\* and *Re*<sub>p</sub> values were calculated, separately for the three above-distinguished groups. The reason for using the log-normal distribution was the fact that asymmetric distributions were expected (based on the scattering of the points), where *Re* numbers cannot take negative values [65]. The probability density and the cumulative distribution functions are described in the Appendix A. Based on the cumulative distribution functions, the *Re*\*- and *Re*<sub>p</sub>-based separatenesses of the different point groups become qualitatively characterizable.

### 3. Results

First, the field measurement data were analyzed in Figure 5, separately. The scattering of the points reflects a clear tendency: the bed shear stress increases with higher grain sizes. Nevertheless, when indicating a separation line, which represents a *Re*\* of 350, the points from the three well-distinguishable regions indeed show different features. Points representing the shallower, sandy regions fall below the line. As most of these sampling locations can be characterized with mean grain sizes in the gravel range (*d*<sub>50</sub> > 2 mm), the Shields-Parker diagram would consequently suggest clear gravel transport. However, considering the relatively low bed shear stress values (<2 N/m²), the Shield curve predicts only sand motion (see the findings based on Figure 2, where τ\* = 0.05 → τ ≈ 2 N/m²). Indeed, the *d*<sub>50</sub> of the bed material samples at the outlet of the flume (4–6 mm) would refer to gravel domination in the bed load, but the corresponding *d*<sub>50</sub> of the bed load samples show clear sand domination (1.3 mm) [59].

Second, the points distributed along the line belong to the bars. At this part of the points, the domination of sand or gravel cannot be well distinguished, but most probably a similar amount of them is transported. Accordingly, two phenomena can occur: (i) the bed shear stress values are low and the bed materials are quite uniform and coarse, resulting in resting particles (red stars, τ < ~1 N/m²); (ii) higher flow transport capacity can cause almost the same order of magnitude in the sand and gravel discharge, if the bed material contains both fractions (red stars, where τ > ~2 N/m² and see also Figure 4).

Finally, the blue diamonds located above the line represent the stable [17,59,64] main stream with coarse grains. Because of the armored bed surface, weak sediment transport occurs during mean flow
discharge. Only higher flood waves can result in significant bed load transport [64], containing of local bed material, i.e., only gravel grains.

The constant particle Reynolds number of 360 values are indicated with a vertical dotted line. Unlike the separation line of $Re^* = 350$, the particle Reynolds number does not distinguish the different point classes.

Assessing the relationship between local Shields number and the shear Reynolds number for both the laboratory and field experiments (Figure 6A), the statements in the previous paragraphs are corroborated regarding the $Re^*$-based classification of the landforms. Again, the $Re^* \approx 300-400$ range outlines as a suitable indicator for the identification whether the gravel or the sand particles dominate.

Contrarily, no clear distinction of the dominating morphodynamic processes can be made if the particle Reynolds number is considered (Figure 6B). Most of the points, representing all three different transport modes, range between $1000 < Re_p < 10,000$, clearly below the critical condition provided by Shields. Accordingly, a stable bed surface is expected both in the laboratory and in the field study. However, this is only the case for the armored sections of the study sites (blue diamonds) [17,59,64], whereas the bed load measurements of the laboratory experiments and the bed load rating curve at the study site in the Danube River demonstrate that weak bed load transport took place in both cases. Moreover, the bed load material investigation of the laboratory experiment suggested that the bed contained mainly sand fractions. Overall, despite the fact that the Shields curve would suggest stable bed, weak and sand-dominated bed load transport was monitored during the field and laboratory measurements.

Finally, the results of the log-normal-distribution-based probability examination are presented. As a summary of the calculations, the probability data was summarized for the six log-normal distributions, considering $Re^* = 300$, $Re^* = 400$, and $Re_p = 360$ as critical values.

Table 2 shows that 5% of the near-bank points are in the $Re^* > 300$ range, 60% of the gravel bar points take $Re^*$ between 300 and 400, while only 3% of the points from the main channel are expected below 400. It can also be seen that, below the critical particle Reynolds number value ($d = 2$ mm; $Re_p = 360$), only the fine sediment transport nature can be expected. However, according to the log-normal distribution, 85% of the near-bank points are over 360. Thus, the $Re_p = 360$ indicates well that below this value mainly finer sediment transport is expected. In turn, the sediment transport nature cannot be determined obviously above this border.

<table><caption>**Table 2.** Probability values of the shear Reynolds number and the particle Reynolds number, calculated based on the fitted log-normal distribution for grouped data.</caption><thead><tr><th rowspan="2">Landform Classes</th><th colspan="2">Re*</th><th colspan="2">Re<sub>p</sub></th></tr><tr><th>&lt;300</th><th>&lt;&gt;</th><th>&gt;400</th><th>&lt;360</th><th>&gt;360</th></tr></thead><tbody><tr><td>Sand aggradation, near-bank points</td><td>95%</td><td>5%</td><td>0%</td><td>15%</td><td>85%</td></tr><tr><td>Side bar</td><td>5%</td><td>60%</td><td>35%</td><td>0%</td><td>100%</td></tr><tr><td>Main channel, bed armor</td><td>0%</td><td>3%</td><td>97%</td><td>0%</td><td>100%</td></tr></tbody></table>

The most important outcome of the probability analysis is that the $Re^*$ between 300 and 400 is indeed a critical range: below 300 the sand transport dominates, above 400 the gravel motion dominates.

The probability values indicated for the three main morphodynamic processes emphasize the suitability of $Re^*$ for the determination of the locally dominant grain size range in the bed load transport. On the other hand, the application of $Re_p$ suggests no clear classification for the same categorization, when spatially varied bed material is present, such as at the investigated river reach.
## 4. Discussion

Previous studies report that the dominant sediment transport nature in rivers with uniform bed material can be reliably determined as the function of the so-called explicit particle Reynolds number ($Re_p$), but yields less accurate estimations in cases when the bed material composition shows strong variability even along shorter reaches.

Based on the available and accessible data, namely 70 bed material samples and related local bed shear stress values gathered both from recent own laboratory and field experiments, it could be confirmed that, instead of the utilization of the explicit particle Reynolds number ($Re_p$), or the application of characteristic grain sizes, e.g., $d_{50}$ or $d_{90}$, the shear Reynolds number ($Re^*$) is a more adequate parameter to assess the locally dominant sediment transport nature. The probability calculations indicate that the sand-dominant bed load is estimable by 95% accuracy, while the gravel majority transport by 97%. In contrast, the $Re_p$-based former method shows significantly greater inaccuracy. Although the gravel-dominated locations are estimated precisely, the method calculates the prediction of fine sediment majority by 85% inaccuracy.

Besides the authors' own experimental data, it could also be shown that the $Re^*$-based approach works reliably for the laboratory and field data based on which the widely used van Rijn [20,57,58] and the Wilcock and Crowe [26] bed load formulas were validated. The method is able to indicate the fine particle domination for all of the validation and calibration data of the van Rijn model correctly. Moreover, the transitional range between the sand and gravel majority becomes detectable in the case of sand–gravel mixture bed material measurements, with which the Wilcock and Crowe model was validated.

As another example, the authors assessed the datasets, consisting of 45 coupled bed shear stress —$d_{50}$ characteristic grain size value pairs from field measurements, published by Mueller et al. [66]. In that study, the median surface grain sizes of the gravel-bed streams and rivers were reported to vary between 0.027 and 0.21 m, which indicates quite coarse bed material and armored bed surface. Accordingly, the calculated $Re^*$ are consequently high, all estimated $Re^*$ values exceed 400.

Based on the data assessment introduced above, the following classification could be set up:

* $Re^* < ~300 \rightarrow$ sand transport dominates;

* $Re^* > ~400 \rightarrow$ gravel transport dominates;

* $~300 < Re^* < ~400 \rightarrow$ gravel accumulating and side bar formation is expected.

Besides the promising validation results, the theory must be further proven by direct bed load field measurements. Also noted, that the local bed load content depends not only on the local bed material, but also on the arriving bed load material [67]. This obvious effect is not taken into account by the introduced method. Additionally, the benchmark measurements were carried out at weak bed load discharges. Thus, the validation needs to be extended to a wider range of the flow regime. However, the authors believe that the presented classification indicates meaningful progress in the investigation methods of the nonuniform bed load.

## 5. Conclusions

The herein-presented classification method of the locally dominant fraction in bed load transport can contribute to better implementation of different sediment transport investigation methods applied in large rivers. For instance, a well-known issue related to field sediment transport monitoring is the high uncertainty of direct bed load sampling methods. Having information on the local bed material and on the flow field, the offered $Re^*$-based approach can suggest where and what sort of sampling techniques would be the most suitable to collect reliable sediment information.

The results can also contribute to the development of improved computational modeling tools. Indeed, instead of applying one specific empirical sediment transport formula in a simulation to calculate local bed load transport, the $Re^*$-based approach can be utilized to distinguish between
several formulas, each of them having a certain application range. Thus, the combined use of several formulas becomes possible.

**Author Contributions:** G.T.T. conceived the study, collected the data, conducted the data analysis, preformed the statistical analysis and drafted the manuscript. G.T.T. and S.B. discussed the results. S.B. worked on subsequent drafts of the manuscript. S.B. and J.J. participated in the design and coordination of the study.

**Funding:** This research was supported by the ÚNKP-17-3 and 17-4 New National Excellence Program of the Ministry of Human Capacities, which is highly appreciated, and was partly supported by MTA TKI of the Hungarian Academy of Sciences. Support of grant BME FIKP-VÍZ by EMMI is also kindly acknowledged. The authors acknowledge the funding of the OTKA FK 128429 grant.

<figure><img src="image_117.png" /></figure>

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

## Appendix A

The probability density function for the log-normal distribution is expressed as:

$$P(x) = \frac{1}{x\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\ln x - \mu}{\sigma}\right)^2}, \quad (A1)$$

where the parameters are $-\infty < \mu < \infty$ (scale parameter) and $\sigma > 0$ (shape parameter), can be obtained by the fitting of the function to the known points. The integration of the log-normal distribution curves results in the cumulative distribution function, which is:

$$F(x) = \Phi\left(\frac{\ln x - \mu}{\sigma}\right), \quad (A2)$$

where $\Phi$ is the Laplace integral [68,69].
<table><caption>**Table A1.** Measured morphodynamic data of field [17,64] (at the introduced reach of the Danube River in Hungary) and laboratory [59] measurements and of the referenced laboratory experiments [17,59,64].</caption><thead><tr><th rowspan="2"></th><th colspan="5">Main Channel, Bed Armor</th><th colspan="5">Sand Aggradation, Near-Bank Points</th><th colspan="5">Side Bar</th></tr><tr><th>u*</th><th>τ*</th><th>Re*</th><th>Re<sub>p</sub></th><th>d<sub>50</sub></th><th>u*</th><th>τ*</th><th>Re*</th><th>Re<sub>p</sub></th><th>d<sub>50</sub></th><th>u*</th><th>τ*</th><th>Re*</th><th>Re<sub>p</sub></th><th>d<sub>50</sub></th></tr></thead><tbody><tr><td rowspan="20">Danube, rkm 1798-1795</td><td>0.054</td><td>0.020</td><td>498</td><td>3546</td><td>9.2</td><td>0.025</td><td>0.005</td><td>185</td><td>2512</td><td>7.3</td><td>0.022</td><td>0.002</td><td>276</td><td>5705</td><td>12.6</td></tr><tr><td>0.053</td><td>0.009</td><td>1031</td><td>10,835</td><td>19.4</td><td>0.010</td><td>0.001</td><td>55</td><td>1604</td><td>5.4</td><td>0.035</td><td>0.007</td><td>364</td><td>4242</td><td>10.4</td></tr><tr><td>0.046</td><td>0.009</td><td>657</td><td>6886</td><td>14.3</td><td>0.022</td><td>0.062</td><td>1</td><td>1</td><td>0.0</td><td>0.028</td><td>0.002</td><td>689</td><td>15,567</td><td>24.6</td></tr><tr><td>0.042</td><td>0.009</td><td>500</td><td>5252</td><td>11.9</td><td>0.129</td><td>0.062</td><td>16</td><td>6</td><td>0.1</td><td>0.033</td><td>0.006</td><td>478</td><td>6995</td><td>14.5</td></tr><tr><td>0.058</td><td>0.019</td><td>633</td><td>4583</td><td>10.9</td><td>0.128</td><td>0.062</td><td>136</td><td>140</td><td>1.1</td><td>0.036</td><td>0.008</td><td>347</td><td>3802</td><td>9.6</td></tr><tr><td>0.061</td><td>0.020</td><td>701</td><td>5005</td><td>11.6</td><td>0.045</td><td>0.399</td><td>14</td><td>23</td><td>0.3</td><td>0.041</td><td>0.009</td><td>443</td><td>4577</td><td>10.9</td></tr><tr><td>0.066</td><td>0.030</td><td>598</td><td>3429</td><td>9.0</td><td>0.010</td><td>0.016</td><td>4</td><td>33</td><td>0.4</td><td>0.037</td><td>0.009</td><td>347</td><td>3676</td><td>9.4</td></tr><tr><td>0.044</td><td>0.010</td><td>540</td><td>5481</td><td>12.3</td><td></td><td></td><td></td><td></td><td></td><td>0.044</td><td>0.010</td><td>519</td><td>5231</td><td>11.9</td></tr><tr><td>0.069</td><td>0.014</td><td>1465</td><td>12,559</td><td>21.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.058</td><td>0.003</td><td>3745</td><td>65,995</td><td>64.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.100</td><td>0.033</td><td>1855</td><td>10,182</td><td>18.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.064</td><td>0.026</td><td>605</td><td>3737</td><td>9.5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.064</td><td>0.017</td><td>946</td><td>7251</td><td>14.8</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.069</td><td>0.043</td><td>475</td><td>2283</td><td>6.9</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.094</td><td>0.029</td><td>1766</td><td>10,314</td><td>18.7</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.084</td><td>0.022</td><td>1659</td><td>11,230</td><td>19.8</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.087</td><td>0.049</td><td>813</td><td>3658</td><td>9.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.056</td><td>0.021</td><td>535</td><td>3715</td><td>9.5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.048</td><td>0.019</td><td>365</td><td>2657</td><td>7.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.087</td><td>0.049</td><td>813</td><td>3658</td><td>9.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.053</td><td>0.019</td><td>499</td><td>3631</td><td>9.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.058</td><td>0.021</td><td>581</td><td>4026</td><td>10.0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.077</td><td>0.007</td><td>3955</td><td>46,694</td><td>51.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.049</td><td>0.004</td><td>1956</td><td>31,713</td><td>39.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.079</td><td>0.012</td><td>2543</td><td>23,129</td><td>32.1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.061</td><td>0.005</td><td>3042</td><td>44,346</td><td>49.5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.082</td><td>0.006</td><td>5374</td><td>67,130</td><td>65.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.060</td><td>0.005</td><td>2896</td><td>42,850</td><td>48.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.073</td><td>0.006</td><td>3714</td><td>46,211</td><td>50.9</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.048</td><td>0.005</td><td>1449</td><td>21,195</td><td>30.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.066</td><td>0.020</td><td>913</td><td>6486</td><td>13.7</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.079</td><td>0.013</td><td>2345</td><td>20,525</td><td>29.6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="13">Laboratory measurements</td><td>0.065</td><td>0.032</td><td>528</td><td>2951</td><td>8.1</td><td>0.040</td><td>0.022</td><td>186</td><td>1262</td><td>4.6</td><td>0.035</td><td>0.011</td><td>228</td><td>2155</td><td>6.6</td></tr><tr><td>0.078</td><td>0.032</td><td>915</td><td>5087</td><td>11.7</td><td>0.033</td><td>0.013</td><td>177</td><td>1570</td><td>5.3</td><td>0.055</td><td>0.026</td><td>398</td><td>2481</td><td>7.2</td></tr><tr><td>0.070</td><td>0.023</td><td>930</td><td>6118</td><td>13.2</td><td>0.036</td><td>0.014</td><td>207</td><td>1727</td><td>5.7</td><td>0.055</td><td>0.021</td><td>478</td><td>3287</td><td>8.7</td></tr><tr><td>0.060</td><td>0.026</td><td>518</td><td>3194</td><td>8.6</td><td>0.043</td><td>0.022</td><td>219</td><td>1474</td><td>5.1</td><td>0.045</td><td>0.017</td><td>329</td><td>2537</td><td>7.4</td></tr><tr><td>0.086</td><td>0.033</td><td>1197</td><td>6548</td><td>13.8</td><td>0.041</td><td>0.016</td><td>257</td><td>2015</td><td>6.3</td><td>0.039</td><td>0.011</td><td>320</td><td>3016</td><td>8.3</td></tr><tr><td>0.058</td><td>0.017</td><td>704</td><td>5374</td><td>12.1</td><td>0.035</td><td>0.011</td><td>237</td><td>2220</td><td>6.7</td><td>0.045</td><td>0.017</td><td>330</td><td>2506</td><td>7.3</td></tr><tr><td>0.061</td><td>0.017</td><td>805</td><td>6151</td><td>13.3</td><td></td><td></td><td></td><td></td><td></td><td>0.048</td><td>0.018</td><td>387</td><td>2914</td><td>8.1</td></tr><tr><td>0.070</td><td>0.021</td><td>1018</td><td>7012</td><td>14.5</td><td></td><td></td><td></td><td></td><td></td><td>0.040</td><td>0.016</td><td>253</td><td>2005</td><td>6.3</td></tr><tr><td>0.056</td><td>0.016</td><td>697</td><td>5558</td><td>12.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.054</td><td>0.018</td><td>549</td><td>4084</td><td>10.1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.087</td><td>0.023</td><td>1785</td><td>11,739</td><td>20.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.050</td><td>0.019</td><td>417</td><td>3066</td><td>8.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.052</td><td>0.019</td><td>459</td><td>3339</td><td>8.8</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

References

1. Dey, S. Fluvial Processes: Meandering and Braiding. In *Fluvial Hydrodynamics: Hydrodynamic and Sediment Transport Phenomena*; Springer: Berlin/Heidelberg, Germany, 2014; pp. 529–562.

2. Leopold, B.L.; Wolman, M.G. *River Channel Patterns: Braided, Meandering and Straight*; Geological Survey Professional Paper; U.S. Government Printing Office: Washington, DC, USA, 1957; Volume 282.

3. Lane, E.W. *A Study of the Shape of Channels Formed by Natural Streams Flowing in Erodible Material*; U.S. Army Corps of Engineers Report; The Division: Omaha, NE, USA, 1957; Volume 141.

4. Leopold, L.B.; Wolman, M.G. River meanders. *Geol. Soc. Am. Bull.* **1960**, *71*, 769–794. [CrossRef]

5. Schumm, S.A. Evolution and response of the fluvial system, sedimentologic implications. *SEPM Spec. Publ.* **1981**, *31*, 19–29.
6. Jaeggi, M.N.R. Effect of Engineering Solutions on Sediment Transport. In *Dynamics of Gravel-Bed Rivers*; Billi, P., Hey, R.D., Thorne, C.R., Tacconi, P., Eds.; John Wiley & Sons Ltd.: Hoboken, NJ, USA, 1992.

7. García, M.H. Sediment Transport and Morphodynamics. In *Sedimentation Engineering: Processes, Measurements, Modeling, and Practice, Manuals and Reports on Engineering Practice*; ASCE: Reston, VA, USA, 2008; pp. 21–163.

8. Parker, G. Transport of Gravel and Sediment Mixtures. In *Sedimentation Engineering*; Garcia, M., Ed.; American Society of Civil Engineers: Reston, VA, USA, 2008; ISBN 978-0-7844-0814-8.

9. Schumm, S.A.; Dumont, J.F.; Holbrook, J.M. *Active Tectonics and Alluvial Rivers*; Cambridge University Press: Cambridge, UK; New York, NY, USA; Melbourne, Australia, 2000; ISBN 0-521-66110-2.

10. Edward, J. Hickin. Chapter 4. In *Fluid Mechanics*; Online course note; Simon Fraser University and R.S. Graphics and Printing: Burnaby, BC, Canada, 2009.

11. Piton, G.; Recking, A. The concept of travelling bedload and its consequences for bedload computation in mountain streams. *Earth Surf. Process. Landf.* **2017**, *42*, 1505–1519. [CrossRef]

12. Farkas-Iványi, K.; Guti, G. The Effect of Hydromorphological Changes on Habitat Composition of the Szigetköz Floodplain. *ACTA Zool. Bulg.* **2014**, *7*, 117–121.

13. Hankó, Z.; Starosolszky, Ö.; Bakonyi, P. Megvalósíthatósági tanulmány a Duna környezetének és hajózhatóságának fejlesztésére (Danube Environmental and navigation Project, Feasibility Study). *Vízügyi Közlemények* **1996**, *78*, 291–315.

14. Rákóczi, L. A Duna-meder sorsa Szap és Szob között (Destiny of the Danube channel between Szap and Szob). *Vízügyi Közlemények* **2000**, *82*, 262–280.

15. Goda, L. A Duna gázlói Pozsony-Mohács között (Shallows of the River Danube between Pozsony, Bratislava and Mohács. *Vízügyi Közlemények* **1995**, *77*, 71–102.

16. Baranya, S.; Józsa, J.; Török, G.T.; Ficsor, J.; Mohácsiné Simon, G.; Habersack, H.; Haimann, M.; Riegler, A.; Liedermann, M.; Hengl, M. A Duna hordalékvizsgálatai a SEDDON osztrák-magyar együttműködési projekt keretében (Introduction of the joint Austro-Hungarian sediment research under the SEDDON ERFE-project). *Hidrologiai Közlöny* **2015**, *95*, 41–46.

17. Török, G.T.; Baranya, S. Morphological investigation of a critical reach of the upper Hungarian Danube. *Period. Polytech. Civ. Eng.* **2017**, *61*. [CrossRef]

18. Észak-Dunántúli Vízügyi Igazgatóság. *Nagyvizi mwe haederkezelési terv (High Water River Management Plan) 01.nmt.02. (egyeztetési terv)*—Duna 1809,76—1786,00 fkm; 2014; Professional report; Észak-Dunántúli Vízügyi Igazgatóság: Győr, Hungary, 2014.

19. Liedermann, M.; Gmeiner, P.; Pessenlehner, S.; Haimann, M.; Hohenblum, P.; Habersack, H. A Methodology for Measuring Microplastic Transport in Large or Medium Rivers. *Water* **2018**, *10*, 414. [CrossRef]

20. Van Rijn, L.C. Sediment Transport, Part I: Bed Load Transport. *J. Hydraul. Eng.* **1984**, *110*, 1431–1456. [CrossRef]

21. Wilcock, P.R.; Crowe, J.C. Surface-based transport model for mixed-size sediment. *J. Hydraul. Eng.* **2003**, *129*, 120–128. [CrossRef]

22. Shields, A. Application of Similarity Principles and Turbulence Research to Bed-Load Movement. *Mitt. Preuss. Versuchsanst. Wasserbau Schiffbau* **1936**, *26*, 47.

23. Bergillos, R.J.; Rodríguez-Delgado, C.; Millares, A.; Ortega-Sánchez, M.; Losada, M.A. Impact of river regulation on a Mediterranean delta: Assessment of managed versus unmanaged scenatios. *Water Resour. Res.* **2016**, *52*, 5132–5148. [CrossRef]

24. Bergillos, R.J.; Ortega-Sánchez, M.; Masselink, G.; Losada, M.A. Morpho-sedimentary dynamics of a micro-tidal mixed sand and gravel beach, Playa Granada, southern Spain. *Mar. Geol.* **2016**, *379*, 28–38. [CrossRef]

25. Fischer-Antze, T.; Rüther, N.; Olsen, N.; Gutknecht, D. Three-dimensional (3D) modeling of non-uniform sediment transport in a channel bend with unsteady flow. *J. Hydraul. Res.* **2009**, *47*, 670–675. [CrossRef]

26. Wilcock, P.R.; Kenworthy, S.T.; Crowe, J.C. Experimental Study of the Transport of Mixed Sand and Gravel. *Water Resour. Res.* **2001**, *37*, 3349–3358. [CrossRef]

27. Wilcock, P.R. A two-fraction model for the transport of sand/gravel mixtures. *Water Resour. Res.* **2002**, *38*, 1–12. [CrossRef]

28. Sziło, J.; Bialik, R.J. Grain Size Distribution of Bedload Transport in a Glaciated Catchment (Baranowski Glacier, King George Island, Western Antarctica). *Water* **2018**, *10*, 360. [CrossRef]
29. Bialik, R.J. Numerical study of saltation of non-uniform grains. *J. Hydraul. Res.* **2011**, *49*, 697–701. [CrossRef]

30. Einstein, H.A. *The Bed-Load Function for Sediment Transportation in Open Channel Flows*; Department of Agriculture, Soil Conservation Service: Washington, DC, USA, 1950.

31. Ashida, K.; Michiue, M. Study on hydraulic resistance and bedload transport rate in alluvial streams. *Trans. Jpn. Soc. Civ. Eng.* **1972**, *206*, 59–69. [CrossRef]

32. Meyer-Peter, E.; Müller, R. Formulas for Bed-Load Transport. In *Proceeding of the IAHSR 2nd Meeting*, Stockholm, Sweden, 7–9 June 1948.

33. Wiberg, P.L.; Smith, J.D. Model for calculating bed load transport of sediment. *J. Hydraul. Eng.* **1989**, *115*, 101–123. [CrossRef]

34. Török, G.T.; Baranya, S.; Rüther, N. 3D CFD Modeling of Local Scouring, Bed Armoring and Sediment Deposition. *Water* **2017**, *9*, 56. [CrossRef]

35. Gaeuman, D.; Andrews, E.D.; Krause, A.; Smith, W. Predicting fractional bed load transport rates: Application of the Wilcock-Crowe equations to a regulated gravel bed river. *Water Resour. Res.* **2009**, *45*, 1–15. [CrossRef]

36. Janssen, S.R. Testing Sediment Transport Models under Partial Transport Conditions. Master's Thesis, University of Twente, Enschede, The Netherlands, 2010.

37. Török, G.T.; Baranya, S.; Rüther, N. Three-dimensional numerical modeling of non-uniform sediment transport and bed armoring process. In *Proceedings of the 18th Congress of the Asia & Pacific Division of the International Association for Hydro-Environment Engineering and Research, Jeju, Korea*, 19–23 August 2012.

38. Rüther, N.; Olsen, N.R.B. Modelling free-forming meander evolution in a laboratory channel using three-dimensional computational fluid dynamics. *Geomorphology* **2007**, *89*, 308–319. [CrossRef]

39. Rüther, N.; Olsen, N.R.B. 3D modeling of transient bed deformation in a sine-generated laboratory channel with two different width to depth ratios. In *Proceedings of the International Conference on Fluvial Hydraulics, Lisbon, Portugal, 6–8 September 2006; Taylor & Francis: New York, NY, USA, 2006*.

40. Rüther, N.; Olsen, N.R.B. Three-dimensional modeling of sediment transport in a narrow 90° channel bend. *J. Hydraul. Eng.* **2005**, *131*, 917–920. [CrossRef]

41. Olsen, N.R.B.; Melaanen, M.C. Three-dimensional calculation of scour around cylinders. *J. Hydraul. Eng.* **1993**, *119*, 1048–1054. [CrossRef]

42. Bihs, H.; Olsen, N.R.B. Numerical Modeling of Abutment Scour with the Focus on the Incipient Motion on Sloping Beds. *J. Hydraul. Eng.* **2011**, *137*, 1287–1292. [CrossRef]

43. Deltares. 3D/2D modelling suite for integral water solutions. *Delft3D-FLOW, User Manual*; Deltares: Delft, The Netherlands, 2014; Volume 710.

44. Skinner, C.J.; Coulthard, T.J.; Schwanghart, W.; Van De Wiel, M.J.; Hancock, G. Global sensitivity analysis of parameter uncertainty in landscape evolution models. *Geosci. Model Dev.* **2018**, *11*, 4873–4888. [CrossRef]

45. Unal, N.E. Shear Stress-Based Analysis of Sediment Incipient Deposition in Rigid Boundary Open Channels. *Water* **2018**, *10*, 1399. [CrossRef]

46. Habersack, H.; Kreisler, A. Sediment Transport Processes. In *Dating Torrential Processes on Fans and Cones*; Beniston, M., Ed.; Springer: Dordrecht, The Netherlands; Heidelberg, Germany; New York, NY, USA; London, UK, 2013; ISBN 9789400743359.

47. Zanke, U.C.E. *Hydromechanik der Gerinne und Küstengewässer*; Blackwell; Vieweg+Teubner Verlag: Berlin, Germany, 2002; ISBN 978-3-322-80212-5.

48. Parker, G.; Klingeman, P.C. On why gravel rivers are paved. *Water Resour. Res.* **1982**, *18*, 1409–1423. [CrossRef]

49. Cushman-Roisin, B. Rivers and Streams, Chapter 15. In *Environnemental Fluid Mechanics*; John Wiley & Sons, Inc.: Hanover, NH, USA, 2019.

50. Parker, G. Surface-based bedload transport relation for gravel rivers. *J. Hydraul. Res.* **1990**, *28*, 417–436. [CrossRef]

51. Dey, B.S.; Papanicolaou, A. Sediment Threshold under Stream Flow: A State-of-the-Art Review. *KSCE J. Civ. Eng.* **2008**, *12*, 45–60. [CrossRef]

52. Dwivedi, A.; Melville, B.W.; Shamseldin, A.Y.; Guha, T.K. Analysis of hydrodynamic lift on a bed sediment particle. *J. Geophys. Res.* **2011**, *116*. [CrossRef]

53. Morsi, S.A.; Alexander, A.J. An investigation of particle trajectories in two phase flow systems. *J. Fluid Mech.* **1972**, *55*, 193–208. [CrossRef]
54. Rubinow, S.I.; Keller, J.B. The transverse force on a spinning sphere moving in a viscous fluid. *J. Fluid Mech.* **1961**, *11*, 447–459. [CrossRef]

55. Saffman, P.G. The lift on a small sphere in a slow shear flow. *J. Fluid Mech.* **1965**, *22*, 385–400. [CrossRef]

56. Dey, S. Sediment threshold. *Appl. Math. Model.* **1999**, *23*, 399–417. [CrossRef]

57. Fernandez Luque, R. Erosion and Transport of Bed-load Sediment. BSc Thesis, Delft Technical University, Delft, The Netherlands, 1974.

58. Fernandez Luque, R.; van Beek, R. Erosion and Transport of Bed-load Sediment. *J. Hydraul. Res.* **1976**, *14*, 127–144. [CrossRef]

59. Török, G.T.; Baranya, S.; Rüther, N.; Spiller, S. Laboratory analysis of armor layer development in a local scour around a groin. In Proceedings of the International Conference on Fluvial Hydraulics, RIVER FLOW 2014, Lausanne, Switzerland, 3–5 September 2014; Taylor and Francis Group: Lausanne, Switzerland, 2014; pp. 1455–1462.

60. Detert, M.; Weitbrecht, V. User guide to gravelometric image analysis by BASEGRAIN. In *Advances in River Sediment Research*; Fukuoka, S., Nakagawa, H., Sumi, T., Zhang, H., Eds.; CRC Press: Boca Raton, FL, USA, 2013; ISBN 9781138000629.

61. Kim, S.C.; Friedrichs, C.T.; Maa, J.P.Y.; Wright, L.D. Estimating bottom stress in tidal boundary layer from Acoustic Doppler velocimeter data. *J. Hydraul. Eng.* **2000**, *126*, 399–406. [CrossRef]

62. Cea, L.; Puertas, J.; Pena, L. Velocity measurements on highly turbulent free surface flow using ADV. *Exp. Fluids* **2007**, *42*, 333–348. [CrossRef]

63. Przyborowski, Ł.; Loboda, A.M.; Bialik, R.J. Experimental investigations of interactions between sand wave movements, flow structure, and individual aquatic plants in natural rivers: A case study of Potamogeton Pectinatus, *L. Water* **2018**, *10*, 1166. [CrossRef]

64. Török, G.T. Methodological Improvement of Morphodynamic Investigation Tools for Rivers with Non-Uniform Bed Material By. Ph.D. Thesis, Budapest University of Technology and Economics, Budapest, Hungary, 2018.

65. Koris, K. *Hidrológia II. (Hydrology II)*; 1. kiadás.; Class note (Egyetemi jegyzet): Budapest, Hungary, 2014; ISBN 978-963-12-0752-1.

66. Mueller, E.R.; Pitlick, J.; Nelson, J.M. Variation in the reference Shields stress for bed load transport in gravel-bed streams and rivers. *Water Resour. Res.* **2005**, *41*. [CrossRef]

67. Vowinckel, B.; Jain, R.; Kempe, T.; Fröhlich, J. Entrainment of single particles in a turbulent open-channel flow: A numerical study. *J. Hydraul. Res.* **2016**, *54*, 158–171. [CrossRef]

68. *EasyFit Professional, Version 5*; MathWave Technologies: Washington, DC, USA, 2015.

69. Aristizabal, R.J. Estimating the Parameters of the Three-Parameter Lognormal Distribution. FIU Electronic Theses and Dissertation, Florida International University, St. Miami, FL, USA, 2012.

<figure><img src="image_118.png" /></figure>

© 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/).
Article

# Modelling of Soil Erosion and Accumulation in an Agricultural Landscape—A Comparison of Selected Approaches Applied at the Small Stream Basin Level in the Czech Republic

Jiří Jakubínský ¹,* , Vilém Pechanec ², Jan Procházka ² and Pavel Cudlín ¹

¹ Global Change Research Institute CAS, Bělidla 986/4a, 603 00 Brno, Czech Republic; cudlin.p@czechglobe.cz

² Department of Geoinformatics, Palacký University Olomouc, 17. listopadu 50, 771 46 Olomouc, Czech Republic; vilem.pechanec@upol.cz (V.P.); jan.prochazka07@upol.cz (J.P.)

* Correspondence: jakubinsky.j@czechglobe.cz

Received: 28 December 2018; Accepted: 20 February 2019; Published: 26 February 2019

<figure><img src="image_119.png" /></figure>

check for
updates

**Abstract:** This article deals with the modelling of erosion and accumulation processes in the contemporary cultural landscape of Central Europe. The area of interest is the headwater part of the small stream catchment—the Kopaninský Stream in central Czech Republic. It is an agricultural and forest-agricultural landscape with a relatively rugged topography and riverbed slope, which makes the terrain very vulnerable to water erosion. The main aim of this article is to compare the results of four selected soil erosion and sediment delivery models, which are currently widely used to quantitate the soil erosion and sediment accumulation rates, respectively. The models WaTEM/SEDEM, USPED, InVEST and TerrSet work on several different algorithms. The model outputs are compared in terms of the total volume of eroded and accumulated sediment within the catchment per time unit, and further according to the spatial distribution of sites susceptible to soil loss or sediment accumulation. Although each model is based partly on a specific calculation algorithm and has different data pre-processing requirements, we have achieved relatively comparable results in calculating the average annual soil loss and accumulation. However, each model is distinct in identifying the spatial distribution of specific locations prone to soil loss or accumulation processes.

**Keywords:** soil loss; sediment delivery; erosion modelling; environmental change; agriculture; Czech Republic

## 1. Introduction

Soil erosion is a physical process occurring almost everywhere in the world, even in landscapes with low slopes. Its main cause is, in addition to natural terrain properties, the character of the soil substrate and local climatic conditions, especially the prevailing land-use type and intensity of the anthropogenic pressure to which the given landscape is subjected. According to Pimentel [1], about 10 million hectares of fertile land are eroded each year, and about 115 million hectares of soil are affected by erosion in Europe [2]. Generally, this topic has been dealt with in the long term by both the professional community and political decision makers – as evidenced by the fact that soil erosion, as the main soil degradation process, has been identified as a key priority within the Soil Thematic Strategy of the European Commission [3]. In Europe, a number of research projects have been carried out in the past focusing on soil erosion and its modelling, and the importance of the theme is also illustrated by the fact that the European Soil Data Centre (ESDAC) was established within the Joint Research Centre in Ispra, Italy [4].
It is obvious that humans—as the initiators of soil loss, increasing the susceptibility of the soil to erosion—act directly through specific anthropogenic influences most often connected with land-use changes, as well as indirectly, usually in connection with the impacts of climate change, which is also conditioned by human activities to a certain extent. Indirect impacts of human activities affecting the soil erosion rates are most often manifested in the form of extreme hydrometeorological phenomena (flash floods or drought episodes). The relationship between the impacts of climate change and the soil loss potential is being given increasing attention, e.g., [5–8]. This is also closely related to the evapotranspiration rates and plant biomass production, among other factors [9]—it is a mechanism to reduce biomass production that could lead to long-term increases in soil erosion due to sparser vegetation cover [10].

When the upper layer of the soil cover is dredged from the agricultural land, there is a rapid decrease in the fertility, which leads to an obvious decrease in the yields of the crops grown there. An easily erodable soil substrate on sloping areas also presents a certain risk during significant precipitation episodes, where specific flood events in the form of debris flow (a fast-moving landslide made up of liquefied, unconsolidated, and saturated sediment mass) may occur. De Vente et al. [11] consider such debris flow an off-site effect of soil erosion. This phenomenon is manifested, for example, in the clogging of riverbeds, increased flood risk, and reduced lifetime of water reservoirs [12]. This makes soil erosion a critical problem [13] that requires holistic solutions involving physical and socioeconomic approaches. For these reasons, accurate quantification of soil loss and sediment delivery rates is very important, and information on the spatial distribution of these phenomena in a given landscape is also crucial. In addition to processes characterizing the current state, the possibility of predicting soil erosion-accumulation rates in the near future also becomes significant. Modelling of these phenomena takes into account, in particular, the expected impacts of climate change on the current cultural landscape as well as locally specific conditions that usually reflect the decision-making of regional stakeholders—usually landowners or local governments. The outlined factors influencing the prediction of soil erosion are, on the practical level, reflected through the land-use changes.

The Universal Soil Loss Equation (USLE) [14] and its more recent version, the Revised Universal Soil Loss Equation (RUSLE), revised by Renard et al. [15], are applied for the purposes of determining the volume of soil erosion almost everywhere in the world. Algorithms of a large number of soil loss and sediment delivery models and procedures are based on these equations—the models applied in this study are no exception. In the long term, soil erosion issues are highly relevant, especially in the Mediterranean region [16], which is, according to de Vente et al. [17], particularly sensitive to erosion due to its climate characterized by dry summers followed by intense autumn rainfall and often a steep topography with fragile soils. These findings illustrate the soil erosion risk maps created, for example, in Italy [18], Greece [19] or Spain [20]. Similar conditions, when the erosion-accumulation phenomena in the landscape were additionally amplified by the impacts of climate change starting in the late 1980s, were analyzed and subsequently described also in North Africa, as exemplified by Hallouz et al. [21] in a case study from Algeria. In Central European conditions, the increased susceptibility of soil to erosion is still mainly due to the use of inappropriate farming methods on sloping land, but the influence of climatic conditions and extreme hydrometeorological situations is increasing, which greatly increases the area of cultivated land with potential for higher soil erosion risk (e.g., [22–24]). The choice of specific crops grown on those lands plays a key role in the erosion rate risk—in Central Europe, this issue was dealt with, e.g., by Jones et al. [25] or Lieskovský and Kenderessy [26].

In the Czech Republic, soil erosion is one of the most widespread types of soil degradation [27], with more than 50% of agricultural land at present being threatened by water erosion and more than 10% by wind erosion [28]. The modelling of erosion processes in the Czech landscape was dealt with, e.g., by Dostál et al. [29], Krása et al. [30], Van Rompaey et al. [31] or Konečná et al. [32]. Although the currently used approaches allow relatively accurate identification of sites subject to erosion risk only based on modelling techniques (i.e., without the need for field measurements), there is still considerable uncertainty associated with the determination of the volume of sediment eroded. This is
also evidenced by the relatively different results achieved by individual models applied to the same territory within this study. The aim of this article is to determine the main causes of the differences between four selected, commonly used soil loss and sediment delivery models. Specifically, it is a comparison of the models in terms of their ability to identify sites with the potential for soil erosion or accumulation, as well as the exact quantification of the volume of transported material within the analyzed area. Based on this information, it is possible to determine the range of the most suitable application tasks in which the models achieve the highest-quality results.

## 2. Materials and Methods

### 2.1. Case Study Area

Selected soil loss and sediment delivery models were applied to the study area located in central Czech Republic in the Bohemian-Moravian Highlands, approximately 90 km southeast of Prague. This is the headwater area of the Kopaninský Stream (IV. order stream) with an average altitude of 565 m, with an area of 7.1 km², diverting water to the Želivka River (the left-side tributary of the Sázava River). The exact position of the selected watershed is shown in the map in Figure 1. The area is characterized by rugged terrain and a stream from its headwater area to the final profile near Velký Rybník Village at a distance of 4.3 km, surpassing an elevation of 108 m. Long-term average discharge in the final profile (between 1991 and 2000) was 0.027 m³ s⁻¹ [33]. The basin of the Kopaninský Stream is classified as a slightly warm and damp climatic area, characterized by an average annual temperature of 7.1 °C. In terms of precipitation totals, the study area is located in a rather drier region (average annual rainfall is 665 mm per year) under the influence of the continental climate.

<figure><img src="image_120.png" /><figcaption>**Figure 1.** Kopaninský Stream catchment and its location within the Czech Republic.</figcaption></figure>
The soil type is cambic hyperskeletic Leptosol, with extremely permeable shallow loamy sands. Along the streambed, there are mainly Gleysols and Pseudogleysols in terms of soil cover. From the granulometric point of view, clay soils predominate, supplemented with sandy loam along the watercourse. The catchment area is approximately half covered with forests, and the rest is agricultural land, which is prone to erosion due to the relatively rugged relief (see Figure 2). The erosive potential is increased by the fact that these are compact and large-scale blocks of soil cover, without fragmentation made by line vegetation cover, etc. When evaluating the vegetation cover, the most common are the “distance to nature habitats” and “alien habitats,” which together cover approximately 660 hectares. Selected hydrometeorological and geomorphological characteristics of the area, important also because of their influence on the formation of rainfall-runoff processes in the landscape, are given in Table 1. For the catchment, an agricultural-forest landscape is typical, with approximately 50% of forest stands—these are mainly coniferous forest with an area of nearly 260 hectares. The rest of the territory comprises non-irrigated arable land and also land principally occupied by agriculture, with significant areas of natural vegetation. The composition of agricultural crops corresponds to natural conditions. Grains and oilseed rape predominate; there is also a significant proportion of maize and peas, and some root crops.

<figure><img src="image_121.png" /><figcaption>Figure 2. Typical agricultural landscape in the area of interest (photo: O. Cudlín, 2016).</figcaption></figure>

<table><caption>Table 1. Selected hydrometeorological and geomorphological characteristics of the study area.</caption><thead><tr><th rowspan="2">Precipitation (mm)</th><th colspan="3">Average Annual Values</th><th rowspan="2">Average Terrain Slope (°)</th><th rowspan="2">Average Stream Slope (°)</th></tr><tr><th>Discharge in Closing Profile (m³ s<sup>-1</sup>)</th><th>Specific Outflow (l.s<sup>-1</sup> km<sup>-2</sup>)</th><th>Air Temperature (°C)</th></tr></thead><tbody><tr><td>665.0</td><td>0.027</td><td>3.802</td><td>7.1</td><td>18.3</td><td>2.6</td></tr></tbody></table>

## 2.2. Data Sources

The quality of soil loss and sediment accumulation model outputs is directly dependent on the spatial resolution of the topography used in the form of the digital elevation model (DEM). In particular, the identification of sites with significant terrain slope and their extent, as well as the selection of various relief microforms that can influence the drainage processes at a given river basin by their morphometry, are very important. Often, these are anthropogenic forms of relief, such as road and rail embankments, flood levees, various water passes and other forms of relief, which, moreover, often occur in an agricultural landscape. However, the location and spatial extent of these forms are now relatively easy and reliably identifiable due to the high accuracy of commonly available terrain models based on the application of modern Earth remote sensing methods (e.g., LIDAR). This study used the digital elevation model of the Czech Republic of the 5th generation (DEM 5G) with a spatial resolution of 5 × 5 m and an average height error of 0.18 m in the exposed terrain and 0.30 m in forest terrain. The model was
created from data acquired by airborne laser scanning in 2009–2013 and the resulting data are provided by the Czech Office for Surveying, Mapping and Cadastre (COSMC). DEM is the basic input layer of all soil loss and sediment delivery models based on the USLE equation. Based on these altitudinal data, the slope length and steepness factor (LS factor) can be estimated, which is one of the most important input parameters affecting the results of erosion process calculations. All models used work with a DEM in the form of a raster file.

Another key input for erosion-accumulation models is the data on the nature of prevailing land-use types. In this study, data from the ZABAGED®dataset (COSMC), supplemented by data from the Land Parcel Identification System (LPIS), and subsequently verified on the basis of the current aerial images of the catchment and its field survey, were used for the purposes of characterizing the landscape coverage. The land cover was categorized based on the requirements of individual models applied in this study—because of different working requirements, each model worked with a different number of land-use categories. The values of the individual factors of the USLE equation are also crucial inputs of the models—they are commonly used in its basic (USLE) or revised form (RUSLE) for determining the vulnerability of agricultural soils to water erosion. This approach is given by the following equation [14]:

$$A = R \times K \times LS \times C \times P, \quad (1)$$

where A is the long-term average annual soil loss [t ha⁻<fnref n="1" /> year⁻<fnref n="1" />]; R is a rainfall erosivity factor expressed in terms of kinetic energy, total precipitations and intensity of high precipitations; K represents the soil erodibility factor expressed depending on the texture and structure of soil, organic matter content in the regolith layer and permeability of the soil profile; LS is the topographic factor expressing the effect of uninterrupted slopes and steepness on soil loss amount; C is cover and the management factor expressed depending on the state of vegetation cover and agricultural technologies used; and P represents the support practice factor.

These parameters are given in the form of raster files and in some cases as constant values for the entire area of interest. All models used in the study work on the basis of the RUSLE equation, but some parameters are not applicable for the conditions of the Czech Republic according to the RUSLE methodology because of a lack of data. The equations used are therefore a combination of RUSLE and USLE approaches. The specific values used for the individual equation factors and their sources are presented in Table 2. The spatial distribution of some factors (i.e., K, LS and C factor) within the basin is shown in the maps in Figure 3.

<table><caption>**Table 2.** Factors of the (R)USLE equation, its values and data sources used for the models applied.</caption><thead><tr><th>Factor</th><th>Factor description</th><th>Used Values</th><th>Unit</th><th>Data Source</th></tr></thead><tbody><tr><td>R</td><td>Rainfall erosivity factor</td><td>40.00</td><td>MJ ha<sup>-1</sup> cm h<sup>-1</sup></td><td>Janeček et al. [34]</td></tr><tr><td>K</td><td>Soil erodibility factor</td><td>0.32–0.48<sup>1</sup></td><td>t ha<sup>-1</sup></td><td>Own computation based on Desmet and Govers [35]</td></tr><tr><td>LS</td><td>Topographic factor (slope length and steepness)</td><td>0.03–50.47<sup>1</sup></td><td>-</td><td>Janeček et al. [34]</td></tr><tr><td>C</td><td>Cover and management factor</td><td>0.00–0.65<sup>1</sup></td><td>-</td><td>Janeček et al. [34]</td></tr><tr><td>P</td><td>Support practice factor</td><td>1.00</td><td>-</td><td>Janeček et al. [34]</td></tr></tbody></table>

<footnote n="1">¹ The values of the given factor are spatially variable and the actual distribution within the area of interest is shown in Figure 3.</footnote>
<figure><img src="image_122.png" /><figcaption>**Figure 3.** Spatial distribution of the K, LS and C factor values in the study area.</figcaption></figure>

## 2.3. Approaches Used to Soil Loss and Sediment Accumulation Modelling

The following paragraphs provide basic information on the individual models applied within this article. However, this is not a detailed description of the algorithms and equations on which the models are based; this information is contained in the original studies published by authors of the given models or in the cited articles, the authors of which have also applied the models in practice. For the purposes of this study we selected a total of four models, which can be categorized as soil loss and sediment delivery models according to their official descriptions. According to available documentation, each model uses the (R)USLE equation parameters to calculate the soil loss amounts and other different algorithms to calculate the soil accumulation. In all the models applied, this is the basic USLE equation (see Equation (1)), which determines the potential soil loss from the area of 1 ha per year (it means the long-term average value valid at each site of the river basin). The USLE equation is modified for use in the GIS environment by replacing the “L” and “S” factors by the combined topographic “LS” factor derived directly from the digital elevation model. Modification of the “R” factor is not applicable here because all models consistently expect a direct input of the rainfall erosivity value, not a sub-calculation of the rain erosion efficiency based on meteorological measurements. For ease of comparison, the support practice factor “P” is not considered, the value 1 is explicitly used. An overview of the models used, their selected parameters and input data requirements is provided in the Table 3.

<table><caption>**Table 3.** Selected parameters of the models applied.</caption><thead><tr><th>Model</th><th>InVEST</th><th>USPED</th><th>TerrSet</th><th>WaTEM/SEDEM</th></tr></thead><tbody><tr><td rowspan="6">Basic input parameters</td><td>DEM</td><td>Raster file</td><td>Raster file</td><td>Raster file</td></tr><tr><td>R factor</td><td>Raster file</td><td>Numerical value</td><td>Raster file</td><td>Numerical value</td></tr><tr><td>K factor</td><td>Raster file</td><td>Raster file</td><td>Raster file</td><td>Raster file or value</td></tr><tr><td>LS factor</td><td>Numerical value</td><td>Raster file</td><td>Numerical value</td><td>Numerical value</td></tr><tr><td>C factor</td><td>Numerical values (for each land-use category)</td><td>Raster file</td><td>Raster file</td><td>Raster file or value</td></tr><tr><td>P factor</td><td>Numerical value</td><td>Numerical value</td><td>Raster file</td><td>Not included</td></tr></tbody></table>
<table><caption>Table 3. *Cont.*</caption><thead><tr><th>Model</th><th>InVEST</th><th>USPED</th><th>TerrSet</th><th>WaTEM/SEDEM</th></tr></thead><tbody><tr><td>Other inputs (compulsory)</td><td>River basin layer (shapefile), Biophysical parameters, Land-use and land cover data, Flow accumulation threshold, Calibration parameters</td><td>Not included</td><td>Terrain properties to identify areas with similar erosion rates, SDR values</td><td>Watercourse layer (shapefile), Land-use and land cover data, Ptef, Parcel Connectivity Data, Calibration parameters</td></tr><tr><td>Optional inputs</td><td>Information on drainage systems</td><td>Calibration parameters</td><td>No data value, Units used, Land fragmentation</td><td>Algorithm of LS factor computation, Retention ponds, Ploughing direction, Soil roughness, Units used</td></tr><tr><td>Geodata format</td><td>Common raster file, Shapefile (.SHP), .CSV file (biophysical parameters only)</td><td>Common raster file</td><td>Idrisi file (.RST)</td><td>Idrisi file (.RST)</td></tr></tbody></table>

### 2.3.1. InVEST

InVEST is a large set of open-source tools (models) designed to map and evaluate ecosystems and natural processes. There is a wide range of tools such as models for carbon modelling, crop pollination modelling, sediment transport modelling, but also even for fisheries or recreation modelling. The individual tools are installed in the form of one package, but each is started separately. InVEST is developed within the Natural Capital Project, which is based on a partnership between several of the world's leading academic institutions—Stanford University, the Chinese Academy of Sciences, the University of Minnesota, The Nature Conservancy and the World Wildlife Fund. InVEST is one of the few in the range of used models to be further developed [36]. Earlier versions of the entire InVEST package are also implemented in TerrSet model as an Ecosystem Service Modeler tool.

Outputs from the sediment transport model include the sediment load delivered to the stream at an annual time scale (i.e., “sediment accumulation”), as well as the amount of sediment eroded in the catchment and retained by vegetation and topographic features. The sediment delivery module (used for computation of the sediment accumulation amount) is a spatially explicit model working at the spatial resolution of the input digital elevation model raster. For each pixel, the model first computes the amount of annual soil loss from that pixel, and then computes the sediment delivery ratio (SDR), which is the proportion of soil loss actually reaching the stream. Once sediment reaches the stream, we assume that it ends up at the catchment outlet, thus no in-stream processes are modelled. This approach was proposed by Borselli et al. [37] and has received increasing interest in recent years [38–40]. The amount of annual soil loss on pixel is given by the revised universal soil loss equation (RUSLE).

Based on the work by Borselli et al. [37], the model first computes the connectivity index (IC) for each pixel. The connectivity index describes the hydrological linkage between sources of sediment (from the landscape) and sinks (like streams). Higher values of IC indicate that source erosion is more likely to make it to a sink (i.e., is more connected), which happens, for example, when there is sparse vegetation or higher slope. Lower values of IC (i.e., lower connectivity) are associated with more vegetated areas and lower slopes. The sediment delivery ratio for a pixel is then derived based on the connectivity index IC, maximum theoretical SDR and calibration parameters, defining the shape of the SDR-IC relationship (which is an increasing function). The value of sediment transport from a pixel is the amount of sediment eroded from that pixel that actually reaches the stream [38].
2.3.2. USPED

USPED (Unit Stream Power-Based Erosion/Deposition Model) is a two-dimensional model of erosion and accumulation based on the RUSLE equation. This model is oriented to the spatial distribution of areas with topographical potential for soil erosion and accumulation. The model was developed in the USA in collaboration with experts from the U.S. Army Construction Engineering Laboratories, Illinois Natural History Survey, and the Faculty of Science of Comenius University in Bratislava, Slovakia. The first results were obtained by using this model in Central Illinois and Yakima Ridge, Washington, USA [41]. USPED has been applied to the basins of various scales in the USA [42] and in Italy [43,44]. Unlike the one-dimensional revised universal soil loss equation (RUSLE) models, which assume erosion mainly depends on rainfall detachment capacity, USPED model assumes that soil erosion and accumulation mainly depend on the sediment transport capacity of the surface runoff. If soil particles are already detached by rain, but there is not enough runoff to transport the soil particles because of the terrain shape or vegetation effect, the actual amount of erosion will be significantly reduced. The total surface runoff through a given pixel determines the sediment transport capacity of this area. USPED assumes that both surface runoff and sediment transport are continuous variables. The net erosion/accumulation of sediment can be positive, indicating soil deposition, or negative, indicating soil loss. The sediment transport capacity calculation is approximated using a static sediment flow rate.

In the Central European conditions, the model was used, for example, in a Slovak study by Dotterweich et al. [45] or in the Czech Republic by Bek [46] or more recently by Vysloužilová and Kliment [47]. The major advantage of the USPED model is that it does not require the installation of any special software and can be easily implemented in common GIS environments. Generally, there are a number of interrelated raster data operations that were implemented in ESRI ArcGIS 10.2 (using Raster Calculator and Spatial Analyst Tools). The result is a raster file (the size of the grid corresponding to the input data parameters) giving information about the nature of the phenomena prevailing in each cell of the given raster—positive values indicate predominant sediment delivery and negative values indicate soil loss. A detailed procedure for implementing the USPED model is described in [47] or [43].

2.3.3. TerrSet

TerrSet is commercial raster-based software that includes a large number of individual modelling tools. It also includes tools for modelling the impacts of climate change, and integrates a whole set of InVEST models (older versions of tools) or tools for processing of the remote sensing data. The software is developed and distributed by Clark Labs at Clark University in Worcester, MA, USA [48]. The calculation of soil loss and sediment accumulation was done using the RUSLE model combination, followed by the SEDIMENTATION tool from the IDRISI GIS Analysis tool kit. The main differences from other software are based on the fact that the actual calculation is performed on raster data, but all the analyses take place on the most homogeneous surfaces, the parameters of which are chosen when starting the RUSLE model and are viewed as mutually isolated patches.

According to Clark Labs [48] the SEDIMENTATION tool utilizes patch-level output from the RUSLE modelling tool to evaluate the net soil movement (erosion or deposition) within patches. The Sediment Delivery Ratio (SDR) is determined from the average annual sediment yield. Sediment yield (as in the case of InVEST model) is the amount of soil loss (erosion) that reaches a stream and is transported within the waterway. Determining net erosion or deposition begins with the total soil loss by patch produced from the RUSLE module. To determine net soil loss or deposition in each patch, SEDIMENTATION tool first determines the average elevation for each patch. Then the highest elevation in the river basin or the highest elevation in each field is located. The direction of movement of the soil is then established by the relative elevation differences between contiguous patches. Movement is always in the downslope direction. The amount of soil loss that moves into the surrounding lower patches is proportional to the length of the common boundary between the
higher patch and the lower patches. Next, the net soil loss or soil accumulation in all lower patches is calculated as follows. Using the output from RUSLE, the proportional soil loss for the higher patch is compared to the soil loss for the lower patch. The difference between the amounts of soil loss from the higher patch to the lower patch represents the net soil loss or accumulation in the lower patch.

#### 2.3.4. WaTEM/SEDEM

WaTEM/SEDEM is a specific tool developed at KU Leuven in Belgium, designed for modelling soil erosion and transport processes. The software is freely distributed and is based on the IDRISI GIS analysis tool and also uses its RST file format, which is related to the need to convert input data. The model can be used to determine soil loss and accumulation values due to water erosion, to identify areas of agricultural land susceptible to water erosion processes or to simulate the effects of different soil conservation scenarios [49,50]. The water erosion component of WaTEM/SEDEM uses an adapted version of the Revised Universal Soil loss equation (RUSLE) to calculate mean annual soil loss values. Runoff patterns are calculated with a flow algorithm that takes into account field borders, tillage direction and road infrastructure. Sediment is routed along these flow paths to the nearest river using a transport capacity term that is proportional to the potential rill erosion rate. The boundaries of the land and the topology of the river network are taken into account, which makes it suitable for modelling larger river basins.

The sediment production values do not correspond with erosion values as calculated with RUSLE, the latter being higher. RUSLE only predicts soil loss per pixel, not sediment deposition. On pixels where sediment deposition is lower than RUSLE erosion, netto erosion will be the case. It is the latter value which is used to calculate sediment production (soil accumulation). If sediment accumulation is higher than RUSLE soil loss, there will be netto sediment accumulation and this value will be used to calculate the catchment's total sediment accumulation value. This model uses the R, K, C and P factors of the RUSLE equation. The topographical L and S factors are automatically calculated from the DEM. In order to adapt the RUSLE to a two-dimensional landscape, the upslope length is replaced by the unit contributing area, i.e., the upslope drainage per unit of contour length. The use of contributing areas implies that the effects of flow convergence and divergence on topographically complex landscapes are explicitly accounted for [51].

The erosion rate is considered equal to the sum of the potential erosion unless the local transport capacity is exceeded. The transport capacity coefficient describes the proportionality between the potential for rill erosion and the transport capacity of the overland flow. The transport capacity is the maximum sediment mass that can be transported by the overland flow. If the sediment production is higher than this transport capacity, sediment will be deposited. Thus, the higher the transport capacity coefficient, the more sediment can be transported downslope. For each land use type, transport capacity can be different. The transport capacity on a given slope segment was considered to be directly proportional to the erosion potential. If the sediment inflow exceeds the transport capacity, deposition occurs, so that the amount of material equals the transport capacity. WaTEM/SEDEM version 2004, launched through the 2006 version, is used in this study, which brings together all available versions of the given software and corrects some of their errors. The algorithms used are listed in the program manual [52]. This software generally belongs to relatively frequently used tools—in Central Europe, its outputs were presented, for example, by Bezak et al. [53] or in the Czech Republic by Van Rompaey et al. [31] and by Krása [54], which describes in more detail the individual algorithms of the model.

## 3. Results

### 3.1. *Analysis of Soil Loss Caused by Water Erosion*

Selected models were applied within the area of interest to determine the soil loss amount caused by water erosion. Each model works on a slightly different principle of soil loss estimation.
The WaTEM/SEDEM (W/S) model calculates the "netto erosion" (the results range from negative values indicating soil erosion to positive values that represent the sediment accumulation). This principle is also applied within the USPED model. The InVEST model provides results in the form of "average erosion," which can be understood as a direct output from the USLE equation (the results range from zero to positive values that represent soil loss). In contrast, the TerrSet model calculates both of these options ("netto erosion" and "average erosion"). Comparison of the results of the four used models in terms of spatial distribution of erodible sites within the catchment is provided by the maps in Figures 4 and 5. From these outputs, a different principle of calculation is again evident—whereas the TerrSet model refers to the soil erosion amounts on individual patches corresponding to the sites with a homogeneous land-use type (the basic spatial unit is therefore a unitary land-use category), the W/S predicts the resulting soil erosion in a resolution corresponding to the extent of the input data, i.e., individual pixels with an area of 100 m². On the same principle as W/S, other applied models also work. For this reason, TerrSet's results may be somewhat less accurate in terms of identifying specific sites with the potential to soil loss. It can be concluded that in terms of spatial distribution, comparable results are achieved in particular by InVEST and TerrSet model, which identify almost identical sites, but differ in the estimated soil erosion volumes slightly (with total soil loss in the whole basin ranging from 43.8 to 55.4 tonnes per year). In the case of the USPED model, it is obvious that, as well as sediment accumulation modelling, soil loss volumes are significantly higher than in the case of other models, although significant soil erosion is expected in only a few localities within the catchment. When comparing the results of all four applied models, we can conclude that the models differ especially in terms of the total soil erosion volumes expected throughout the catchment—the values vary from 42.0 to 76.0 tonnes per year. This variance is largely influenced by the different area of the sites with soil erosion potential (see Table 4), which is also the result of a different approach to the areas where soil erosion and sediment accumulation were not identified. While W/S assigns a "0" value and includes these areas in following calculations, all other models assigns a "No data" attribute to that territory and do not work with it in any further steps. The difference in the number of pixels included in the calculations for each model is a key factor affecting the overall results.

<figure><img src="image_123.png" /><figcaption>**Figure 4.** Soil loss in the study area based on the results of the InVEST and USPED models.</figcaption></figure>
<figure><img src="image_124.png" /><figcaption>**Figure 5.** Soil loss in the study area based on the results of the TerrSet and WaTEM/SEDEM models.</figcaption></figure>

<table><caption>**Table 4.** Comparison of the basic soil loss characteristics of the four models applied.</caption><thead><tr><th>Model Used</th><th>InVEST</th><th>USPED</th><th>TerrSet</th><th>WaTEM/SEDEM</th></tr></thead><tbody><tr><td>Total soil loss in the basin (t year<sup>-1</sup>)</td><td>43.78</td><td>76.03</td><td>55.38</td><td>41.96</td></tr><tr><td>Total area with soil loss (ha)</td><td>349.38</td><td>186.2</td><td>354.77</td><td>305.10</td></tr><tr><td>Average soil loss (t ha<sup>-1</sup> year<sup>-1</sup>)</td><td>0.062</td><td>0.107</td><td>0.078</td><td>0.059</td></tr></tbody></table>

The calculated soil loss values and the respective areas of the sites where the erosion takes place are shown in the graph in Figure 6. In principle, this figure captures the dependence between the number of pixels and the corresponding soil erosion values of each output raster file. The resulting chart points to the fact that in terms of the nature of the relationship between the two variables studied, there are significant differences between the models applied. The W/S model achieves the lowest soil losses, with the values found in areas ranging from 0.01 ha to almost 10 hectares (according to the number of pixels with given values).

Relatively similar results are achieved by the InVEST model, which has slightly increased soil loss values and pixels with very low erosion values are not present in such amount as compared to W/S. The curve pattern of the InVEST model outputs indicates that there is almost the same number of pixels with heterogeneous values of soil loss. A similar situation with slightly higher values (with an area covering of about 3 ha) is also the case for the USPED model. Significantly different values that occur in the TerrSet model output are mainly due to the fact that the model works with larger patches characterized by a homogeneous and higher average value.
<figure><img src="image_125.png" /><figcaption>**Figure 6.** Relationship between the soil loss values and the corresponding number of pixels identified within the catchment where soil erosion takes place based on the results of the four models applied.</figcaption></figure>

## 3.2. Analysis of Sediment Accumulation as a Consequence of Soil Loss by Water Erosion

In addition to the soil erosion, all of the applied models also allow the computation of sediment accumulation, and since the study area is a coherent catchment (its upper part)—a hydrologically correctly defined spatial unit—the volume of sediments should correlate with the amount of eroded material to some extent (taking into account that part of the sediment volume is transported outside the area of interest through the closing profile).

Each model has achieved relatively different results not only in terms of total volumes of the sediment accumulation (Table 5), but also when comparing the spatial distribution of the observed phenomenon (Figures 7 and 8). The outputs of all the models have been unified for ease of interpretation in terms of the sediment accumulation amounts displayed, ranging from 0 to 20 tonnes per hectare per year. Although it is not obvious from map visualization, the highest values of sediment accumulation in the study area are predicted by the USPED model, which assumes soil accumulation of more than 0.11 tonnes per hectare per year. The main reason is the accumulation of extremely high amounts in the small area, which roughly corresponds to the results of the other models applied. In terms of total soil accumulation in the basin, the results of the InVEST and TerrSet models are relatively similar. However, since InVEST assumes sedimentation on a much larger area, it is clear from the comparison of average sedimentation (tonnes per hectare per year) that the outputs for InVEST and USPED are more comparable.

Based on the spatial pattern of soil accumulation values, a slightly different calculation algorithm is also apparent when comparing the three models—InVEST, USPED (in Figure 7) and WaTEM/SEDEM (on the right in Figure 8) with the TerrSet model outputs (on the left in Figure 8). While these three models counts soil accumulation (as well as soil loss) in each pixel alone and do not communicate
with their surroundings (“per pixel” method), the remaining TerrSet model counts using “per patch” method. This means that at the beginning the area is divided into patches, defined on the basis of the homogeneous terrain properties and landscape cover. This approach creates a group of adjacent pixels that together make up one area in which the calculation takes place. The maximum length of drainage path is given by the user’s value or by landscape elements interrupting the drainage path. Individual patches communicate with each other, which mean that they are receiving information about computed values in adjacent areas. This aspect has a significant effect on the L factor (slope length factor), since in the case of patches the uninterrupted length of drainage path is much larger compared to the size of the per pixel models (i.e., the length of the drainage line in each pixel is always the same).

<table><caption>**Table 5.** Comparison of the basic soil accumulation characteristics of the four models applied.</caption><thead><tr><th>Model Used</th><th>InVEST</th><th>USPED</th><th>TerrSet</th><th>WaTEM/SEDEM</th></tr></thead><tbody><tr><td>Total soil accumulation in the basin (t year<sup>-1</sup>)</td><td>21.24</td><td>80.84</td><td>22.42</td><td>39.86</td></tr><tr><td>Total area with accumulation (ha)</td><td>320.58</td><td>123.87</td><td>93.10</td><td>46.65</td></tr><tr><td>Average soil accumulation (t ha<sup>-1</sup> year<sup>-1</sup>)</td><td>0.030</td><td>0.114</td><td>0.032</td><td>0.056</td></tr></tbody></table>

<figure><img src="image_126.png" /><figcaption>**Figure 7.** Sediment accumulation as a result of soil loss caused by water erosion in the study area—comparison of InVEST and USPED outputs.</figcaption></figure>

When assessing the spatial distribution of the sediment accumulation capacity of the catchment, it is obvious that large sediment accumulation-prone areas are predicted, especially by the InVEST and TerrSet model, which expects the sediment accumulation to occur at the foot of practically all slopes (this is also reflected in the sediment accumulation area, achieving the highest extent from all the models in case of the InVEST). However, based on the data in Table 5, it is evident that InVEST identifies large areas susceptible to sediment accumulation, but at the same time the model assigns relatively low values to those sites, compared to other models (up to 269 tonnes/hectare/year). By comparing the graphical outputs of all models, the influence of a different calculation algorithm is obvious. While the
second pair of models (Figure 8) shows that the key factor for the prediction of the spatial distribution of sediment accumulation is the nature of land-use (the individual areas identified correspond to the extent and boundary position of the land-use units), for the first pair of models (Figure 7) the terrain characteristics play a crucial role—topography is projected onto the results, especially in the case of the InVEST model. Specific to the InVEST and USPED models, the calculation is also carried out on road surfaces, which the TerrSet and WaTEM/SEDEM models exclude from the calculation.

<figure><img src="image_127.png" /><figcaption>**Figure 8.** Sediment accumulation as a result of soil loss caused by water erosion in the study area—comparison of TerrSet and WaTEM/SEDEM outputs.</figcaption></figure>

The above-described differences in the results of the analyzed models are evident, especially in the detailed view of the graphical outputs, showing the spatial distribution of the observed phenomena (Figure 9). The W/S model attribute high sediment accumulation values to only a few relatively small areas and minimal accumulation values that are homogeneous on a relatively large area. In contrast, the InVEST model allows for much more extensive areas, characterized by increased, but not extreme, sediment accumulation. From the anthropogenically conditioned forms of relief, the road network plays a very important role in soil erosion and sediment accumulation modelling. This element influences the results of the models by modifying the morphological parameters of the terrain, which is then reflected in input data for modelling (DEM), but also in the manner in which a given model processes this type of landscape cover. The TerrSet and WaTEM/SEDEM models used exclude the road network from calculation (they gave the “No data” attribute to the pixels), which is reflected especially when comparing the average sediment accumulation or soil loss for the whole catchment area—each model therefore calculates with a little different area affected by the number of pixels with the “No data” attribute.
<figure><img src="image_128.png" /><figcaption>**Figure 9.** Detailed view on the selected site within the area of interest—comparison of the soil accumulation modelling results.</figcaption></figure>

For the efficient interpretation of differences in the results of the analyzed models, graph (Figure 10) was created. The graph shows the course of the dependencies between the amount of sediment accumulation (i.e., soil accumulation per an area of one hectare) and an area that is characterized by given specific soil accumulation capabilities. Based on Figure 10 it is obvious that the analyzed models have different results in terms of their quality. From the distribution of values within the correlation plot, it can be stated that the greatest differences is observed between the W/S and the USPED model, where W/S assigns a large number of pixels with very low soil accumulation values, the USPED assumes the occurrence of pixels with significantly higher values. The TerrSet and InVEST models have a very similar distribution of values within the correlation plot.
<figure><img src="image_129.png" /><figcaption>**Figure 10.** Relationship between the sediment accumulation values and the corresponding number of pixels (pixel area is 100 m²) identified within the area where the accumulation takes place based on the results of the four models applied.</figcaption></figure>

**4. Discussion and Conclusions**

Based on the application of the four soil loss and sediment accumulation models on a selected small stream catchment in a typical Central European landscape, it can be stated that there are obvious differences between their outputs, both in terms of the total expected soil erosion and sedimentation volumes and spatial distribution of these phenomena within the catchment. Different model results, however, cannot be considered as a symptom of a certain error in the processing, but rather as a result of different algorithms and approaches used to process the input data (which were otherwise identical for each model). One of the main causes of differences in the results is the fact that each model actually works with a different number of pixels within the used raster files. This fact is dependent on different approaches to the classification of the land cover, for example, the W/S model also includes areas where there is no soil erosion (assigns a “0” value), whereas the other applied models refer to these areas as “No data” and exclude them. Furthermore, the TerrSet and W/S models also earmark the road network from the final calculation. This different fragmentation of the study area and the continuity of individual patches are key factors influencing the resulting soil erosion and sediment deposition values, respectively. Another factor influencing soil loss values is the runoff algorithm (flow accumulation algorithm) used—while W/S and TerrSet use the “Multiple Flow Direction” algorithm, USPED and InVEST have only a “Single Flow Direction” mechanism integrated. The details about their differences and influences on the runoff values are described in Wolock and McCabe [55].

Apart from the above-described differences between the results of each modelling process and its probable causes, it is also necessary to take into account the data input quality, which can greatly influence the overall accuracy of the results. It is mainly the spatial and temporal uncertainty of
the input data that affects the quality of the result. An example is the "cover and management factor" (C factor), where the same value for each crop for the whole year was applied to all the models. These values are based on widely used tabular reports valid for the territory of the Czech Republic [34]. A similar principle is often used in other European countries. However, it is evident that the use of a single C factor value is largely inaccurate, especially in regions with alternating seasons, as it fluctuates during the course of the year, depending on the actual growth phase of vegetation. Much more accurate results that could subsequently be used for other decision-making processes in a given area can be achieved by determining the C factor value by means of a series of several soil erosion partial calculations, for different conditions, depending on a particular growing season (for more information on this issue, see Pechanec et al. [56]). However, to assess the quality of model outputs, the use of a single value is sufficient, especially with regard to the necessity of affecting the differences in the actual soil loss or sediment delivery calculation algorithm, not in the quality of input data.

The models used can be further compared according to their user requirements and input data requirements—in this regard, the InVEST model is probably the most user-friendly. It has a relatively purely processed user interface and has the lowest input and pre-processing requirements. The WaTEM/SEDEM model offers a fairly good user interface, but the need for high-quality input data and its pre-processing is considerable. Additionally, the need to work with data in IDRISI format can now be considered as a certain disadvantage. In contrast, the USPED model does not have any custom interface—it is essentially a series of calculations (mathematical operations with raster files) that can be done in any GIS software. Overall, this approach can be considered relatively easy; however, due to the large number of intermediate stages, it is a rather time-consuming process. The advantages and disadvantages of all analyzed models with respect to their practical applicability are summarized in Table 6.

<table><caption>**Table 6.** Rating of the models in terms of their demands on the input data and previous user experiences. 
"++" means the most positive rating – user-friendly model, low input data requirements, etc.; "-" indicates the most negative rating; and 0 means that the model does not have a user interface.</caption><thead><tr><th>Model Used</th><th>InVEST</th><th>USPED</th><th>TerrSet</th><th>WaTEM/SEDEM</th></tr></thead><tbody><tr><td>Data pre-processing requirements</td><td>++</td><td>-</td><td>+</td><td>-</td></tr><tr><td>User-friendly interface</td><td>++</td><td>0</td><td>+</td><td>+</td></tr><tr><td>Claims for prior user experiences</td><td>+</td><td>++</td><td>-</td><td>-</td></tr><tr><td>Failure sensitivity (including the difficulty of error detecting)</td><td>+</td><td>+</td><td>-</td><td>-</td></tr></tbody></table>

Another possibility to evaluate the quality of the outputs of individual models is offered based on their comparison with real soil erosion values, observed and measured in the field directly within the area of interest. Since the field measurement process of soil loss is a relatively demanding issue, it is not always possible to perform such a comparison and thus verify the results of model simulations. The advantage of the area studied in this article is the fact that it is a long-term observation area, which has been investigated in the past by a number of other Czech authors (i.e., Doležal et al. [57] or Fučík et al. [58]—among other measurements, soil erosion has been measured over a certain time period. Pavlík et al. [59] carried out field measurements in the Kopaninský Stream catchment in 2005 and 2006 and found out that the average annual soil erosion amounts vary in the range between 3.4 and 4.4 t ha⁻¹ year⁻¹. According to the Czech Hydrometeorological Institute (CHMI) in the area of interest the average annual precipitation amounts recorded in these years (2005 and 2006) were 798 resp. 826 mm, which is a significantly above-average state (long-term average annual value is 665 mm year⁻¹). The outputs can also be compared with the study realized by Panagos et al [8], which predicts a soil erosion of approximately 1.0–8.9 t ha⁻¹ year⁻¹ in the soil erosion map of Europe and soil database (ESDAC) for the given area of interest. Based on these results, it is evident that most models applied achieve highly underestimated values when compared with the outputs mentioned.
However, it is very important to emphasize that the comparison of the model outputs with the real values obtained by field measurements is only for the idea of how much the model outputs can differ from the real situation. But it is not possible to consider the models as inoperative because of the differences, as the main reason for these differences is the fact that the models do not capture the current conditions in the locality, but only the typical average values valid for a wider area. One of the key elements affecting significantly the outputs is the choice of the C factor (cover and management factor). In this study we were used year round constant value recommended for the Czech Republic by Janeček et al. [34]. Another reason for discrepancies between model outputs and the real situation may also be the unavailability of information on the existence of possible erosion control measures in the area of interest. For this reason, a constant value of "1" for the P factor (Support practice factor) is entered into the models, but for example, in the above-mentioned comparison, we do not have any detailed information on the particular state of the erosion control measures at the site where the measurements took place in the past. An important finding from the comparison of the four selected models is that all models achieve very similar results in assessing the overall average soil loss (the values found range from 0.125 to 0.156 tonnes per hectare per year). Thus, by applying these models, relatively similar results can be obtained, but in terms of the spatial distribution of the particular soil loss or accumulation areas, the outputs can vary considerably and they cannot be used for precise identification of specific locations. For the purposes of precise identification of such soil erosion-prone areas as is required for example in landscape planning or flood protection solutions, it is appropriate to use differently designed models, described, for example, in Pechanec et al. [60].

Overall, it can be summarized that for the purposes of application tasks the InVEST model is probably the most acceptable from the four models compared within this article. This model is less demanding in terms of claims on the input data and is characterized by a relatively high-quality user interface, which makes it usable even for less experienced users. The USPED model is also easy to apply, but the ease lies mainly in the fact that it is a sequence of operations feasible in any GIS software, therefore the steps are clear and the user can easily verify the quality of the partial results. The opposite is the W/S model, which works in the form of a black box (the user does not see the algorithms of the individual calculations) and the input data need to be entered in a specified format. However, all the models applied are characterized by the demands on a certain user experience in the field of environmental modelling and, in particular, the availability of the most accurate data, which completely determines the quality of the outputs.

**Author Contributions:** J.J. wrote the main paper and performed the USPED model computations. V.P. and J.P. performed the remaining model computations. V.P. also discussed the results and commented on the manuscript at all stages. P.C. discussed the results and commented on the manuscript in its initial phase.

**Funding:** This study was supported by the Ministry of Education, Youth and Sports of CR within the National Sustainability Program I (NPU I), grant number LO1415 and also by the project "LAND4FLOOD: Natural Flood Retention on Private Land", grant number LTC18025 and project "Small watercourses and riparian ecosystem management to mitigate the impacts of environmental change (SMART2Envi)", grant number LTC18069.

**Acknowledgments:** The authors would like to thank Ing. Petr Fučík, Ph.D. for his initial idea to carry out this type of research in the given area of interest and also to Ing. Ondřej Cudlín, Ph.D. for providing photos from the area of interest shown in the article.

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

## References

1. Pimentel, D. Soil erosion: A food and environmental threat. *Environ. Dev. Sustain.* **2006**, *8*, 119–137.
   [CrossRef]

2. Kibblewhite, M.G.; Miko, L.; Montanarella, L. Legal frameworks for soil protection: Current development and technical information requirements. *Curr. Opin. Environ. Sustain.* **2012**, *4*, 573–577. [CrossRef]

3. Communication from the Commission to the Council, the European Parliament, the European Economic and Social Committee and the Committee of the Regions—Thematic Strategy for Soil Protection. Brussels 2006.
Available online: https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:52006DC0231 (accessed on 17 November 2018).

4. Panagos, P.; Van Liedekerke, M.; Jones, A.; Montanarella, L. European Soil Data Centre (ESDAC): Response to European policy support and public data requirements. *Land Use Policy* **2012**, *29*, 329–338. [CrossRef]

5. Morgan, R.P.C.; Quinton, J.N. Erosion modeling. In *Landscape Erosion and Evolution Modeling*; Springer: Boston, MA, USA, 2001; pp. 117–143.

6. Nunes, J.P.; Seixas, J.; Pacheco, N.R. Vulnerability of water resources, vegetation productivity and soil erosion to climate change in Mediterranean watersheds. *Hydrol. Process*. Int. J. **2008**, *22*, 3115–3134. [CrossRef]

7. Lorencová, E.; Frélichová, J.; Nelson, E.; Vačkář, D. Past and future impacts of land use and climate change on agricultural ecosystem services in the Czech Republic. *Land Use Policy* **2013**, *33*, 183–194. [CrossRef]

8. Panagos, P.; Borrelli, P.; Poesen, J.; Ballabio, C.; Lugato, E.; Meusburger, K.; Montanarella, L.; Alewell, C. The new assessment of soil loss by water erosion in Europe. *Environ. Sci. Policy* **2015**, *54*, 438–447. [CrossRef]

9. Nearing, M.A.; Pruski, F.F.; O’Neal, M.R. Expected climate change impacts on soil erosion rates: A review. *J. Soil Water Conserv.* **2004**, *59*, 43–50.

10. Pruski, F.F.; Nearing, M.A. Runoff and soil-loss responses to changes in precipitation: A computer simulation study. *J. Soil Water Conserv.* **2002**, *57*, 7–16.

11. De Vente, J.; Poesen, J.; Verstraeten, G.; Van Rompaey, A.; Govers, G. Spatially distributed modelling of soil erosion and sediment yield at regional scales in Spain. *Glob. Planet. Chang.* **2008**, *60*, 393–415. [CrossRef]

12. Verstraeten, G.; Poesen, J. The nature of small-scale flooding, muddy floods and retention pond sedimentation in central Belgium. *Geomorphology* **1999**, *29*, 275–292. [CrossRef]

13. Montgomery, D.R. Soil erosion and agricultural sustainability. *Proc. Natl. Acad. Sci. USA* **2007**, *104*, 13268–13272. [CrossRef] [PubMed]

14. Wischmeier, W.H.; Smith, D.D. Predicting Rainfall Erosion Losses—A Guide to Conservation Planning; Science and Education Administration, USDA: Hyattsville, MD, USA, 1978; 62p.

15. Renard, K.G.; Foster, G.R.; Weesies, G.A.; McCool, D.K.; Yoder, D.C. Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE); United States Department of Agriculture: Washington, DC, USA, 1997.

16. Guerra, C.A.; Maes, J.; Geijzendorffer, I.; Metzger, M.J. An assessment of soil erosion prevention by vegetation in Mediterranean Europe: Current trends of ecosystem service provision. *Ecol. Indic.* **2016**, *60*, 213–222. [CrossRef]

17. Grimm, M.; Jones, R.J.; Rusco, E.; Montanarella, L. Soil erosion risk in Italy: A revised USLE approach. *Eur. Soil Bureau Res. Rep.* **2003**, *11*, 23.

18. Van Rompaey, A.; Bazzoffi, P.; Jones, R.J.; Montanarella, L. Modeling sediment yields in Italian catchments. *Geomorphology* **2005**, *65*, 157–169. [CrossRef]

19. Kouli, M.; Soupios, P.; Vallianatos, F. Soil erosion prediction using the revised universal soil loss equation (RUSLE) in a GIS framework, Chania, Northwestern Crete, Greece. *Environ. Geol.* **2009**, *57*, 483–497. [CrossRef]

20. García-Ruiz, J.M. The effects of land uses on soil erosion in Spain: A review. *Catena* **2010**, *81*, 1–11. [CrossRef]

21. Hallouz, F.; Meddi, M.; Mahé, G.; Toumi, S.; Rahmani, S.E.A. Erosion, Suspended Sediment Transport and Sedimentation on the Wadi Mina at the Sidi M'Hamed Ben Aouda Dam, Algeria. *Water* **2018**, *10*, 895. [CrossRef]

22. Boardman, J.; Poesen, J. Soil erosion in Europe: Major processes, causes and consequences. *Soil Erosion Eur.* **2006**, *4*, 477–487.

23. Scholz, G.; Quinton, J.N.; Strauss, P. Soil erosion from sugar beet in Central Europe in response to climate change induced seasonal precipitation variations. *Catena* **2008**, *72*, 91–105. [CrossRef]

24. Mullan, D. Soil erosion under the impacts of future climate change: Assessing the statistical significance of future changes and the potential on-site and off-site problems. *Catena* **2013**, *109*, 234–246. [CrossRef]

25. Jones, R.J.; Spoor, G.; Thomasson, A.J. Vulnerability of subsoils in Europe to compaction: A preliminary analysis. *Soil Tillage Res.* **2003**, *73*, 131–143. [CrossRef]

26. Lieskovský, J.; Kenderessy, P. Modelling the effect of vegetation cover and different tillage practices on soil erosion in vineyards: A case study in Vráble (Slovakia) using WATEM/SEDEM. *Land Degrad. Dev.* **2014**, *25*, 288–296. [CrossRef]
27. Novotný, I. *Příručka Ochrany Proti Vodní Erozi*; Ministerstvo Zemědělství: Praha, Czech Republic, 2014; 73p. (In Czech)

28. Novotný, I. *Příručka Ochrany Proti Erozi Zemědělské Půdy*; Ministerstvo Zemědělství a Výzkumný ÚSTAV MELIOrací a Ochrany Půdy, v.v.i.: Praha, Czech Republic, 2017; 86p. (In Czech)

29. Dostál, T.; Krása, J.; Váška, J.; Vrána, K. The map of soil erosion hazard and sediment transport in scale of the Czech Republic. *Vodní Hospodářství* **2002**, *52*, 46–48. (In Czech)

30. Krása, J.; Dostál, T.; Vrána, K.; Plocek, J. Predicting spatial patterns of sediment delivery and impacts of land-use scenarios on sediment transport in Czech catchments. *Land Degrad. Dev.* **2010**, *21*, 367–375. [CrossRef]

31. Van Rompaey, A.; Krása, J.; Dostál, T. Modelling the impact of land cover changes in the Czech Republic on sediment delivery. *Land Use Policy* **2007**, *24*, 576–583. [CrossRef]

32. Konečná, J.; Karásek, P.; Fučík, P.; Podhrázská, J.; Pochop, M.; Ryšavý, S.; Hanák, R. Integration of Soil and Water Conservation Measures in an Intensively Cultivated Watershed – A Case Study of Jihlava River Basin (Czech Republic). *Europ. Countrys.* **2017**, *1*, 17–28.

33. Doležal, F.; Kvítek, T.; Soukup, M.; Kulhavý, Z.; Tippl, M. Czech highlands and peneplains and their hydrological role, with special regards to the Bohemo-Moravian Highland. *IHP HWRP-Berichte* **2004**, *2*, 41–56.

34. Janeček, M. *Ochrana Zemědělské Půdy Před Erozí. Metodika*; Česká Zemědělská Univerzita, Fakulta Životního Prostředí: Praha, Czech Republic, 2012; 113p. (In Czech)

35. Desmet, P.J.J.; Govers, G. A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units. *J. Soil Water Conserv.* **1996**, *51*, 427–433.

36. Natural Capital Project. Available online: https://www.naturalcapitalproject.org/ (accessed on 5 May 2018).

37. Borselli, L.; Cassi, P.; Torri, D. Prolegomena to sediment and flow connectivity in the landscape: A GIS and field numerical assessment. *Catena* **2008**, *75*, 268–277. [CrossRef]

38. Cavalli, M.; Trevisani, S.; Comiti, F.; Marchi, L. Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments. *Geomorphology* **2013**, *188*, 31–41. [CrossRef]

39. Lopez-Vicente, M.; Poesen, J.; Navas, A.; Gaspar, L. Predicting runoff and sediment connectivity and soil erosion by water for different land use scenarios in the Spanish Pre-Pyrenees. *Catena* **2013**, *102*, 62–73. [CrossRef]

40. Sougnez, N.; Wesemael, B.; Van Vanacker, V. Low erosion rates measured for steep, sparsely vegetated catchments in southeast Spain. *Catena* **2011**, *84*, 1–11. [CrossRef]

41. Mitášová, H.; Hofierka, J.; Zlocha, M.; Iverson, L.R. Modelling topographic potential for erosion and deposition using GIS. *Int. J. Geogr. Inf. Syst.* **1996**, *10*, 629–641.

42. Warren, S.D.; Mitášová, H.; Hohmann, M.G.; Landsberger, S.; Iskander, F.Y.; Ruzycki, T.S.; Senseman, G.M. Validation of a 3-D enhancement of the Universal Soil Loss Equation for prediction of soil erosion and sediment deposition. *Catena* **2005**, *64*, 281–296. [CrossRef]

43. Pistocchi, A.; Cassani, G.; Zani, O. Use of the USPED model for mapping soil erosion and managing best land conservation practices. In Proceedings of the 1st International Congress on Environmental Modelling and Software, Lugano, Switzerland, 24–27 June 2002.

44. Pelacani, S.; Märker, M.; Rodolfi, G. Simulation of soil erosion and deposition in a changing land use: A modelling approach to implement the support practice factor. *Geomorphology* **2008**, *99*, 329–340. [CrossRef]

45. Dotterweich, M.; Stankoviansky, M.; Minár, J.; Koco, Š.; Papčo, P. Human induced soil erosion and gully system development in the Late Holocene and future perspectives on landscape evolution: The Myjava Hill Land, Slovakia. *Geomorphology* **2013**, *201*, 227–245. [CrossRef]

46. Bek, S.; Chuman, T.; Šefrna, L. The Usability of Contours in Erosion Modelling: A Case Study on ZABAGED, Czech Republic. *Acta Universitatis Carolinae* **2008**, *1–2*, 77–86.

47. Vysloužilová, B.; Kliment, Z. Modelování erozních a sedimentačních procesů v malém povodí. *Geografie* **2012**, *2*, 170–191. (In Czech)

48. CLARK LABS. Available online: https://clarklabs.org/about/ (accessed on 12 May 2018).

49. Verstraeten, G.; Van Oost, K.; Van Rompaey, A.; Poesen, J.; Govers, G. Evaluating an integrated approach to catchment management to reduce soil loss and sediment pollution through modelling. *Soil Use Manag.* **2002**, *18*, 386–394. [CrossRef]
50. WaTEM/SEDEM Homepage. Available online: http://geo.kuleuven.be/geography/modelling/erosion/watemsedem/index.htm (accessed on 5 May 2018).

51. Hamel, P.; Falinski, K.; Sharp, R.; Auerbach, D.A.; Sánchez-Canales, M.; Dennedy-Frank, P.J. Sediment delivery modeling in practice: Comparing the effects of watershed characteristics and data resolution across hydroclimatic regions. *Sci. Total Environ.* **2017**, *580*, 1381–1388. [CrossRef] [PubMed]

52. Notebaert, B.; Vaes, B.; Verstraeten, G.; Govers, G. WaTEM/SEDEM Version 2006 Manual; KU Leuven, Physical and Regional Geography Research Group: Leuven, Belgium, 2006.

53. Bezak, N.; Rusjan, S.; Petan, S.; Sodnik, J.; Mikoš, M. Estimation of soil loss by the WaTEM/SEDEM model using an automatic parameter estimation procedure. *Environ. Earth Sci.* **2015**, *74*, 5245–5261. [CrossRef]

54. Krása, J. Hodnocení Erozních Procesů Ve Velkých Povodních za Podpory GIS. Dizertační Práce; České Vysoké Učení Technické v Praze, Fakulta Stavební, Katedra Hydromeliorací a Krajinného Inženýrství: Praha, Czech Republic, 2004; 176p. (In Czech)

55. Wolock, D.M.; McCabe, G.J. Comparison of Single and Multiple Flow Direction Algorithms for Computing Topographic Parameters in TOPMODEL. *Water Resour. Res.* **1995**, *5*, 1315–1324. [CrossRef]

56. Pechanec, V.; Mráz, A.; Benc, A.; Cudlín, P. Analysis of spatiotemporal variability of C-factor derived from remote sensing data. *J. Appl. Remote Sens.* **2018**, *12*, 016022. [CrossRef]

57. Doležal, F.; Kulhavý, Z.; Kvítek, T.; Soukup, M.; Čmelík, P.; Fučík, P.; Novák, J.; Peterková, E.; Pilná, P.; Pražák, M.; et al. Hydrologický výzkum v malých zemědělských povodích. *J. Hydrol. Hydromech.* **2006**, *54*, 217–229. (In Czech)

58. Fučík, P.; Kaplická, M.; Kvítek, T.; Peterková, J. Dynamics of Stream Water Quality during Snowmelt and Rainfall—Runoff Events in a Small Agricultural Catchment. *Clean Soil Air Water* **2012**, *40*, 154–163. [CrossRef]

59. Pavlík, F.; Dumbrovský, M.; Podhrázská, J.; Konečná, J. The influence of water erosion processes on sediment and nutrient transport from a small agricultural catchment area. *Acta Univ. Agric. Silvic. Mendel. Brun.* **2012**, *60*, 155–165. [CrossRef]

60. Pechanec, V.; Burian, J.; Kiliánová, H.; Němcová, Z. Geospatial analysis of the spatial conflicts of flood hazard in urban planning. *Morav. Geogr. Rep.* **2011**, *19*, 41–49.

<figure><img src="image_130.png" /></figure>

© 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/).
<figure><img src="image_131.png" /></figure>
MDPI

St. Alban-Anlage 66

4052 Basel

Switzerland

Tel. +41 61 683 77 34

Fax +41 61 302 89 18

www.mdpi.com

Water Editorial Office

E-mail: water@mdpi.com

www.mdpi.com/journal/water

<figure><img src="image_132.png" /></figure>
<figure><img src="image_133.png" /></figure>
MDPI

St. Alban-Anlage 66
4052 Basel
Switzerland

Tel: +41 61 683 77 34

Fax: +41 61 302 89 18

www.mdpi.com

<figure><img src="image_134.png" /></figure>

ISBN 978-3-03921-432-7