2.3.4. Spatial Discretization

An overview of the spatial discretization is given for all numerical models and all case studies in Table 2. The area of C.1 is covered by 109 cross-sections in the 1D model as well as 7980/47,880 computation nodes in 2D/3D, which are based on a horizontal, rectangular discretization of 0.1 × 0.03 m. For the numerical calculations of C.2, 43 cross-sections (1D), 18,918/59,220 (2D/3D) computation nodes were employed. The 2D computation mesh consists of rectangular cells with a node distance of 1.0 × 1.0 m, while in the case of 3D a polyhedral mesh of hexagonal cells (node distance: 1.0 × 1.0 m) was employed. The area of C.3 is covered by 22 cross-sections in the 1D model, 28,396/170,376 computation nodes in 2D/3D. The 2D computation mesh is based on triangular cells with a node distance of 1.0 × 1.0 m, while the 3D polyhedral mesh consists of 1.0 × 1.0 m hexagonal cells. Independent of the case study the water column was subdivided into six vertical layers. In general, several spatial discretizations of varying node distances (double and half the selected distance) were tested to ensure grid-independent solutions.


**Table 2.** Overview of the spatial discretization of all the numerical models (h. = hexagonal, r. = rectangular, and t. = triangular mesh discretization).

#### 2.3.5. Boundary and Initial Conditions

The inflow and outflow boundary conditions are listed in Table 3 for all simulation scenarios. At the inflow the normal velocities are derived from the given discharge and prescribed as boundary conditions. The normal velocities are calculated from a division of the total discharge and the boundary area, resulting in the specific discharge for 1D and 2D and the cell-based velocities in 3D, which are

assigned orthogonal to the inflow cell faces. At the outflows energy slopes (1D and 2D models) or water surface levels (3D models) are used.


**Table 3.** Overview of boundary conditions for all simulation scenarios.

The 1D and 2D numerical simulations were performed without particular initial conditions (dry channel), while in the case of 3D numerical simulations the initial free surface elevations were set based on the results of the 1D numerical simulations.

#### *2.4. Comparison of Modelling Results*

In order to ensure comparability between the different model dimensions as well as case studies, result ratios *Mr*,*d*,*c*,*<sup>q</sup>* were considered,

$$M\_{r,d,c,q} = \left\{ \frac{r\_{d,c,q,s,n}}{R\_{c,q,s}} \right\} \tag{10}$$

where *rd*,*c*,*q*,*s*,*<sup>n</sup>* denotes the 2D or 3D simulated result (water depth, depth-averaged flow velocity or bed shear stress), while the 1D hydrodynamic result *Rc*,*q*,*s* serves as reference in the calculation. The resulting ratio *Mr*,*d*,*c*,*<sup>q</sup>* is the basis of the boxplot. The indices used point out the origin of individual sets distinguished by results *r* (water depth, depth-averaged flow velocity or bed shear stress), dimensions *d* (2D or 3D), case studies *c* (C.1 to C.3) and discharges *q* (*QL* to *HQ*1). Additionally used indices representing the cross-section *s* and computational node *n* sliced by the cross-section are necessary for the ratio calculations. The average value is calculated separately for each individual set *Mr*,*d*,*c*,*q*.

#### **3. Results**

#### *3.1. Calibration of Numerical Models*

In the course of the calibration, the flow resistances were adjusted in the numerical models until a match with measured water surface elevations was achieved. Good agreements with average differences below 1 cm between measurements and simulations were obtained for C.1 (Figure 2a—solid lines) for all numerical models excluding the inlet area between distance 10.5 m and 9.5 m, where maximum differences of 3 cm occurred, which are related to issues of the physical model setup. The calibrated roughness values are *kSt* = 52 m1/3 s−<sup>1</sup> for the 1D and 2D model and *kS* = 0.005 m for the 3D model.

In C.2 (Figure 2b—solid lines), comparisons of measurements and numerical calculations show average differences of 2 cm excluding maximum outliers of 12 cm (1D) and 9 cm (2D and 3D) at a distance of 90 m. The single occurrence of the outlier indicates a local measurement issue. Best agreement was achieved by employing roughness values in section 215–150 m of *kSt* = 9.1/9.0 m1/3 s−<sup>1</sup> (1D/2D-model) and *kS* = 1.5 m (3D model), while the areas up- and downstream of this section (300–215 m and 150–0 m) received roughness values of *kSt* = 18.2/18.0 m1/3 s−<sup>1</sup> (1D/2D) and *kS* = 0.25 m (3D model). The higher flow resistance in section 215–150 m reflects the stronger turbulent flow conditions in this area, which is characterized by a steeper bed slope and coarser bed sediments.

The average differences in C.3 (Figure 2c) are below 1 cm except the comparison at distance 220 m, where outliers of 9 cm (2D) and 8 cm (3D) occur. The single occurrence of the outlier again indicates a local measurement issue. The calibrated roughness values in C.3 are *kSt* = 25.0 m1/3 s−<sup>1</sup> for the 1D and 2D model and *kS* = 0.20 m for the 3D model.

In the natural river sites (C.2 and C.3), additional discharge scenarios, defined by the mean discharge of each area (Table 3), have been simulated with the same roughness values and are depicted in Figure 2b,c with dotted lines. Due to lack of data, an evaluation with measurements was not possible, but the results were compared among each other. In C.2 (Figure 2b—dotted lines) the maximum average difference of 3 cm occurred between 2D and 3D simulations, while the maximum average difference in C.3 (Figure 2c—dotted lines) was 1 cm, determined between 1D and 3D simulations.

**Figure 2.** Calibration—Longitudinal profile plots of (**a**) C.1, (**b**) C.2 and (**c**) C.3 including bed elevations (dashed black line) and comparisons of measured (black crosses) and modelled water surface elevations (W.S.E.) at low flow (*QL*) conditions (results of 1D simulations = blue solid line, 2D simulations = red solid line and 3D simulations = black solid line); Comparison of modelled W.S.E. at an additional discharge scenario (mean flow *QM* ) in C.2 and C.3 (results of 1D simulations = blue dotted line, 2D simulations = red dotted line and 3D simulations = black dotted line).

#### *3.2. Sensitivity Analysis of Roughness Values in the Case of High Flow Conditions*

The difficulties of using numerical simulations in cases of high flow conditions (*HQ*1) without monitoring data is assessed in the form of a sensitivity analysis in C.2, demonstrating the differences between various models and the influences of changing roughness values. The water surface elevations simulated by 1D (blue solid line), 2D (red solid line) and 3D (black solid line) models is shown in Figure 3 including a bandwidth (dotted lines with corresponding colors) calculated with roughness values, which were changed by +/−20% referred to calibration (Section 3.1). Good agreements were found between 1D and 2D simulations, using roughness values of calibration resulting in average difference of around 5 cm, while the comparison to 3D simulations yields increasing differences starting with values of 1 cm at the outflow boundary (river section 0 m) up to 60 cm at the inflow boundary (river section 300 m). Following a roughness reduction of 20% in the case of 1D and 2D models, the highest differences compared to 3D simulations at reference conditions reduce to 40 cm, while a roughness increase of 20% in the case of 3D modelling leads to the smallest differences of 15 cm compared to 1D using roughness values of the calibration.

**Figure 3.** Sensitivity analysis—Longitudinal profile plot of C.2 at high flow conditions (*HQ*1) including bed elevations (dashed black line) and modelled water surface elevations (W.S.E.) (results of 1D simulations = blue solid line, 2D simulations = red solid line and 3D simulations = black solid line) with a bandwidth considering +/−20% of roughness changes (dotted lines with corresponding colors).

#### *3.3. Analysis of Bed Shear Stresses*

Bed shear stresses simulated in C.1 are depicted in Figure 4a for all numerical models. In the case of the 1D simulation, the bed shear stresses are 4.9 Nm−<sup>2</sup> without any differences in the lateral or longitudinal direction. Along the thalweg, similar values are obtained in the higher dimensional models, but due to the two consecutive bends higher bed shear stresses up to 6.5 Nm−<sup>2</sup> and 10.0 Nm−<sup>2</sup> occur along the inner banks of the 2D and 3D simulations, respectively, while the outer banks are characterized by lower values of around 2.0/0.4 Nm−<sup>2</sup> (2D/3D).

In C.2 at mean flow conditions (*QM*; Figure 4b), the highest bed shear stresses arise in the area between cross-sections 21 and 30, which is a riffle section [72] with low water depths, high flow velocities as well as high roughness values, resulting in values of 79/108/93 Nm−<sup>2</sup> depending on the model dimension (1D/2D/3D). In the other sections, the bed shear stresses are substantially lower with values between 2 and 16 Nm−<sup>2</sup> upstream the riffle section and between 2 and 41 Nm−<sup>2</sup> downstream of this section in all models. Independent of the section, it is obvious that the 1D model is not able to predict any lateral variability, leading to higher values at the river banks and lower bed shear stresses in the main stream compared to higher dimensional models. In addition, it is shown that the values in the main stream of the 2D model tend to be larger in comparison to the 3D model, while at the banks a slight reverse trend occurs.

The modelling results of C.3 at mean flow conditions (*QM*; Figure 4c) also yield the highest bed shear stresses in the riffle sections, which are located between cross-sections 7 and 11 (21/33/44 Nm−<sup>2</sup> for 1D/2D/3D) as well as cross-section 16 and 19 (13/16/19 Nm−<sup>2</sup> for 1D/2D/3D), and are independent of the model dimension. While in the 2D simulations, the peaks occur in the main stream of the river, the highest values in the 3D simulation arise at the river banks. The reduced spatial discretization in the case of 1D simulations is again shown by constant values in the lateral direction.

**Figure 4.** Comparison of bed shear stresses calculated by 1D/2D/3D models depicted for (**a**) C.1, (**b**) C.2 at *QM* and (**c**) C.3 at *QM*.

#### *3.4. Generalized Comparison of Numerical Models*

A generalized representation of modelling results including water depths, depth-averaged flow velocities and bed shear stresses is depicted in Figure 5 for all case studies and discharge scenarios. The water depth ratio between 2D/1D and 3D/1D (Figure 5a) at low flow and mean flow conditions is almost unity, stating that no substantial differences in the calculated water depths occurred by numerical models of different dimensions. An exception is the ratio 3D/1D at high flow conditions with quartiles of 0.80 and 0.90. Figure 5b shows the depth-averaged flow velocity ratio between 2D/1D and 3D/1D for all scenarios. The quartiles in the case of 2D/1D are between 0.64 and 1.23 and 0.61 and 1.43 in the case of 3D/1D. The bed shear stress ratios between 2D/1D and 3D/1D (Figure 5c) are characterized by quartiles of 0.45 and 1.28 and 0.43 and 1.15, respectively. Additionally, a summarized analysis (Figure 5d) including the averaged ratio between 2D and 3D values referred to 1D values shows almost no differences in the case of water depths (blue) at low flow and mean flow conditions resulting in values close to unity. An exception is the ratio 3D/1D at high flow conditions with a value of 0.85. In the case of flow velocities (grey), averaged ratios of 0.88 to 0.99 (2D/1D) and 0.88 to 1.13 (3D/1D) were obtained, whereby the ratios increase with the discharges in each case study. The largest deviations were found again for bed shear stresses (brown) with values of 0.90 to 1.00 (2D/1D) and 0.62 to 0.86 (3D/1D). Independent of the scenario, the calculated ratios are below the line of perfect agreement, indicating on average lower bed shear stresses in 3D models compared to 2D models.

**Figure 5.** Relative comparison of the modeling results of (**a**) water depths, (**b**) depth-averaged flow velocities and (**c**) bed shear stresses consisting of ratios between 3D and 2D values, referring to 1D values depicted in box-plots for each case study and discharge scenario; (**d**) Average ratio between 3D and 2D values, referring to 1D values for water depths (blue), flow velocities (grey) and bed shear stresses (brown) depicted for each case study and discharge scenario (C.1 (\*), C.2-*QL* (-), C.2-*QM* (), C.2-*HQ*<sup>1</sup> (-), C.3-*QL* (), C.3-*QM* (•)).

Although models were calibrated to water surface elevations, the results exhibit clear differences in bed shear stress ratios. While in the case of 2D/1D, minor differences are noticed between the different river sites and discharges, the ratio decreased in the case of 3D/1D from lower to higher discharges (*QL*→*QM* and *QL*→*HQ*1) as well as from small-sloped to steeper river sites (C.3→C.2).

#### *3.5. Detailed Anaylsis of Simulated Hydrodynamics in Cross-Sections*

Figure 6 depicts the simulated hydrodynamics of all model dimensions including bed shear stresses, flow velocities, turbulent kinetic energy (3D) and water surface elevations in three representative cross-sections of all case studies (CS 49 of C.1, CS 28 of C.2, CS 17 of C.3—encircled in Figure 1). The different model dimensions are based on various definitions of flow velocities. The area-averaged flow velocity (light grey solid line) in the case of 1D, the depth-averaged flow velocity (2D—dark grey solid line) and the near-wall tangential flow velocity (3D—black solid line), which equals the near-wall flow velocity in the case of small bed cell slopes.

**Figure 6.** Detailed analysis of the hydrodynamics simulated by the 1D, 2D and 3D models in (**a**) CS 49 of C.1, (**b**) CS 28 of C.2 at *QM*, (**c**) CS 17 of C.3 at *QL* and (**d**) CS 17 of C.3 at *QM* including water surface elevations (blue dotted lines), bed shear stresses (brown large-spaced dashed lines), various flow velocities (area-averaged flow velocities (grey lines), depth-averaged flow velocities (dark grey lines), near wall tangential flow velocities (black lines)) and turbulent kinetic energies (green lines with rectangles) as well as bed elevations (black small-spaced dashed lines).

The area-averaged flow velocity in CS 49 of C.1 (Figure 6a) is constant with a value of 0.79 m s<sup>−</sup>1. The profile of the 2D depth-averaged flow velocities is characterized by a maximum of 0.84 m s−<sup>1</sup> close to the left bank and a minimum of 0.41 m s−<sup>1</sup> on the right bank. A similar profile was calculated for the near-wall tangential flow velocity in the 3D model, whereby the values are substantially lower, characterized by a peak of 0.72 m s−<sup>1</sup> and a minimum of 0.10 m s−1. The near-wall turbulent kinetic energy follows a falling trend starting with the highest value (0.04 m2 s<sup>−</sup>2) at the left bank and ending with the lowest value at the right bank (0.01 m2 s<sup>−</sup>2). The water depths are almost equal in all models with an average value of 0.43 m, whereby slight increasing trends from the left to the right bank occur in the 2D and 3D models. Due to the trapezoidal channel geometry, the lowest water depths arise at the banks. The comparison of bed shear stresses displays the dependency of this variable on flow velocities, turbulent kinetic energies and water depths. In the case of 2D, the profile of the bed shear stresses is similar to the profile of the flow velocities resulting in a maximum of 5.8 Nm−<sup>2</sup> in the area of the left bank and a minimum of 2.0 Nm−<sup>2</sup> close to the right bank, but due to the influence of the water depths an offset between the peaks of bed shear stress and flow velocity occurs. Comparable results were found in the case of 3D models, where similarities occur between the profiles of the near-wall tangential flow velocities and bed shear stresses, with a maximum (10.2 Nm−2) on the left bank and a minimum (0.4 Nm−2) on the right bank. However, the strong influence of the turbulent kinetic energy is obvious particularly at the left bank resulting in the highest bed shear stresses of the entire cross-section. In contrast, the one-dimensional approach results in a constant value of 4.8 Nm−<sup>2</sup> in the whole cross-section.

In CS 28 of C.2 (Figure 6b), the area-averaged flow velocity is constant with a value of 0.54 m s<sup>−</sup>1. In the case of 2D, the profile of the depth-averaged flow velocities is characterized by a maximum of 0.69 m s−<sup>1</sup> in the main stream of the river located at distance 11 m of the cross-section. Similarities are present for the profile of the near-wall tangential flow velocity in the 3D model with a peak of 0.64 m s<sup>−</sup>1. The near-wall turbulent kinetic energies are characterized by a relative constant plateau in the cross-section center (11–19 m) with values between 0.70 and 0.80 m<sup>2</sup> s<sup>−</sup>1, a decreasing tendency to the banks and peaks up to 0.83 m<sup>2</sup> s−<sup>1</sup> on the banks. The water surface elevations are almost equal in all models with an average value of 100.68 m. Due to the straight river course in this cross-section, the water surface elevations do not change between the banks, hence the water depths only depend on the channel geometry. The comparison of bed shear stresses again indicates the dependency of this variable of flow velocities, turbulent kinetic energies and water depths. In the case of 2D as well as 3D, the profile of the bed shear stresses is similar to the profile of the depth-averaged and tangential flow velocities, resulting in maxima of 65.9 and 63.3 Nm−<sup>2</sup> in the main stream of the river. However, the influence of the turbulent kinetic energy is visible at the right bank, resulting in a peak of 11.0 Nm<sup>−</sup>2. In contrast, the one-dimensional approach leads to a constant value of 44.0 Nm−<sup>2</sup> in the whole cross-section.

The hydrodynamic variables for low flow *QL* (Figure 6c) and mean flow conditions *QM* (Figure 6d) are depicted for CS 17 of C.3. The area-averaged flow velocity is constant over the whole cross-section with values of 0.65 m s−<sup>1</sup> (*QL*) and 0.66 m s−<sup>1</sup> (*QM*). The profile of the 2D depth-averaged flow velocities is characterized by a maximum of 0.68 m s−<sup>1</sup> (*QL*) and 0.96 m s−<sup>1</sup> (*QM*), located in the main stream of the river. A similar profile was calculated for the near-wall tangential flow velocity in the 3D model, whereby the values are substantially lower with a peak of 0.41 m s−<sup>1</sup> (*QL*) and 0.59 m s−<sup>1</sup> (*QM*).

The near-wall turbulent kinetic energies are characterized by values of 0.03–0.04 m<sup>2</sup> s−<sup>2</sup> (*QL*), and 0.03–0.045 m<sup>2</sup> s−<sup>2</sup> (*QM*) in the main stream and peaks up to 0.14 m2 s−<sup>2</sup> (*QL*), and 0.10 m2 s−<sup>2</sup> (*QM*) at the banks. The water surface elevations are almost equal in all models with an average value of 270.95 m a.s.l. at low flow conditions and 271.36/271.34/271.33 m a.s.l. (1D/2D/3D) at mean flow conditions. Due to the straight river course in this cross-section, the water surface elevations do not vary between the banks, thus the water depths only depend on the present geometry. The comparison of bed shear stresses again displays the dependency of this variable on flow velocities, turbulent kinetic energies and water depths. Independent of the discharge, the profiles of the bed shear stresses

calculated in the 2D model follow the profiles of the depth-averaged flow velocities, resulting in peak values of 9.4 Nm−<sup>2</sup> (*QL*) and 15.7 Nm−<sup>2</sup> (*QM*) in the main stream of the river. In the 3D model, comparable results were found in the main stream, where similarities between the profiles of the near-wall tangential flow velocities and bed shear stresses occur, with a maximum of 8.4 Nm−<sup>2</sup> (*QL*) and 12.2 Nm−<sup>2</sup> (*QM*). However, the strong influence of the turbulent kinetic energies close to the banks is obvious leading to the highest bed shear stresses of 9.4/15.8 Nm−<sup>2</sup> (*QL*/*QM*) in the entire cross-section. In contrast, the one-dimensional approach results in a constant value of 9.9/8.6 Nm−<sup>2</sup> (*QL*/*QM*) in the whole cross-section.

#### **4. Discussion**

The comparison of simulation results calculated in different model dimensions in general involves the risk of potential result dependencies on the applied spatial discretization. In order to overcome these issues, a sensitivity analysis including a variation of horizontal node distances was performed to ensure grid-independent solutions. However, in the case of 3D simulations the vertical discretization must be considered as well. In engineering practice, this discretization is usually a justified compromise between computational time and numerical accuracy. Hence, in sediment transport simulations, for which a correct prediction of the shear stress parameter is paramount, which depends on the vertical discretization, literature reports the usage of typically five to eight layers; good agreement with measurements has shown to justify this choice [53,58,59,69]. In line with this practice, the simulations in this study were performed with six layers.

Based on grid independency, the numerical models were successfully calibrated at low flow conditions. Although all site-specific heterogeneities (e.g., bed sediments, channel units and vegetation) have to be aggregated in roughness values, good agreement with measured water surface elevations was achieved in all areas investigated. The study shows that the roughness values are constant along the river course in all case studies except C.2, where a substantial coarsening of sediments in a riffle section led to a lower Strickler value and a higher equivalent sand roughness, respectively. The determined values are in line with guidelines [41–43].

However, the comparison of simulation results calculated by different model dimensions shows that the largest differences were found for bed shear stresses, due to the following reasons: (i) The averaged comparison ratio is lower than unity because of the area-averaged approach in the case of 1D models resulting in overestimations in bank areas. (ii) Moreover, the values are lower in general in the case of 3D compared to 2D simulations, which has also been reported by Lane et al. [80] and is based on the different calculation approaches. The bed shear stress calculations in 2D models are based on the depth-averaged flow velocities, which are in general higher than the near-wall flow velocities used in 3D models. (iii) Additionally, it was shown that the averaged comparison ratio between 3D/2D and 3D/1D is decreasing with increasing discharges, independent of the case study. This fact is again related to the different calculation approaches based on the corresponding roughness definitions.

A closer look into the different roughness definitions reveals the known fact that the Strickler value is a variable depending on the water depth and thus actually has to be adjusted according to the discharge (Ferguson [1]), while the equivalent sand roughness can be considered as constant as long as the river bed is stable. Jäggi [40] suggested an improved calculation of the Strickler value depending on the water depth, which leads to increasing *kSt* values—exemplarily calculated for case study C.2 to be 16/19/27 m1/3 s−<sup>1</sup> (*QL*/*QM*/*HQ*1)—while in this study, in line with predominant hydraulic engineering practice, the roughness is kept constant at 18 m1/3 s<sup>−</sup>1. In particular, the large differences in the case of high flow conditions have to be highlighted, which result in substantially higher bed shear stresses as well as water surface elevations calculated by the 1D and 2D models compared to the 3D model.

According to Jäggi [40], Lamb et al. [81] and Parker et al. [82], the Strickler value is additionally increasing with decreasing bed slopes. Considering this and the calculation of bed shear stresses in

1D and 2D models, where high *kSt* values result in low bed shear stresses, the decreasing 3D/1D ratios from small-sloped to steeper river sites (C.3→C.2) become obvious. Due to the fact that several parameters (e.g., water depth, equivalent sand roughness, etc.) vary between the different case studies, only this overall tendency can be identified.

A comparison to measure bed shear stresses could be a remedy to overcome the issue of considerable differences between the models. However, monitoring this parameter itself [83] is an ongoing challenge. Besides issues of handling the devices (installation, operation and service), concerns about the validity and explanatory power of time-averaged bed shear stresses were raised. Gmeiner et al. [84] and Liedermann et al. [85] showed that bed load transport takes place below expected discharges considering the concepts of Shields [86] and Zanke [87] for critical shear stresses. It is thus recommended to consider time series of bed shear stresses including fluctuations for which probabilistic approaches [88–91] are suitable. A part of these fluctuations is expected to be covered by considering the turbulent kinetic energies for calculating bed shear stresses in the case of 3D models. In this study, the dependency of bed shear stresses on turbulent kinetic energies are highlighted by peaks at the banks, which correlate with peaks of the turbulent kinetic energy. A proportional relationship between turbulent kinetic energy and bed shear stress for high strains has already been described by Menter [92]. Additionally, Rodi [93] stated an increasing influence of turbulence on shear stresses upon larger roughness values.

The aforementioned uncertainties of using deterministic bed shear stress calculations, as well as the additional limitations when employing these concepts in numerical models should always be considered by addressing the issues related to dependent processes including bed load transport or morphological changes. Particularly in the case of man-made interventions [94], including engineering measures [95] in rivers, reservoirs or oceans as well as management procedures [96] such as flushing, dredging or depositing sediments, it is recommended to focus on the evaluation of differences rather than on absolute values. In general, improving the results of 1D and 2D models by applying discharge-dependent Strickler values and preferring higher dimensional models over lower ones has the potential to enhance the validity of numerical simulations.

#### **5. Conclusions**

The application of numerical models in hydraulic engineering has become routine over recent decades, resulting in numerous published studies. Standard practice is the calibration of a single numerical model by comparing simulation results to monitoring data (e.g., water surface elevations, flow velocities, etc.) followed by an application of the models for any kind of work. Despite known uncertainties (e.g., model simplifications, mesh generation, roughness definition, etc.), it is generally assumed that the simulated results are valid for all hydrodynamic parameters, including properties for which no measurements are available. The present study addressed the differences between model dimensions, focusing on bed shear stresses considering different discharges as well as study sites.

The results show that, on average, the simulated bed shear stresses are 10% and 14–38% lower in the applied 2D and 3D model, respectively, as compared to the 1D model. At the same time, almost no differences were present in water surface elevations, which had been used for calibration. When comparing different river sites and discharges, minor differences in the ratio of bed shear stresses were found in 2D/1D models, while in the case of 3D/1D the differences became larger from lower to higher discharges as well as from small-sloped to steeper river sites. The major influence of different roughness definitions, as well as the various calculation approaches used in the numerical models, were identified as a cause of these differences. Moreover, reasons why, in general, the highest bed shear stresses in the main stream of the river occur in 2D models (e.g., depth-averaged versus near-bed flow velocities), and for bed shear stress peaks at the banks in 3D models were elaborated (e.g., consideration of turbulent kinetic energy).

The considerable differences between the numerical models present an issue in hydraulic engineering, due to the fact that bed shear stresses form the basis of many approaches including the calculation of sediment transport or evaluation of habitats. One of the key tasks of future works will be the comparison with measurement data considering the associated monitoring challenges. The question of whether the deterministic approach of calculating bed shear stresses, as well as the definition of bed shear stress, is adequately representing the occurring near-bed forces in an applied river context should also be addressed in future studies.

**Author Contributions:** K.G., M.T. and C.H. designed, set up and performed the numerical simulations and analyzed the results including the comparison with measurements as well as between the models. The draft text was prepared by K.G. and edited by M.T., C.H. and H.H.

**Funding:** This paper was written as a contribution to the Christian Doppler Laboratory for Sediment Research and Management. The financial support by the Austrian Federal Ministry for Digital and Economic Affairs and the National Foundation for Research, Technology and Development is gratefully acknowledged.

**Acknowledgments:** The authors thank the Hydraulic Engineering Institute of the University of Innsbruck for providing measurement data of the laboratory experiment.

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