**Kurt Glock 1,\*, Michael Tritthart 1, Helmut Habersack <sup>2</sup> and Christoph Hauer <sup>1</sup>**


Received: 20 December 2018; Accepted: 26 January 2019; Published: 29 January 2019

**Abstract:** For centuries, scientists have been attempting to map complex hydraulic processes to empirical formulas using different flow resistance definitions, which are further applied in numerical models. Now questions arise as to how consistent the simulated results are between the model dimensions and what influence different morphologies and flow conditions have. For this reason, 1D, 2D and 3D simulations were performed and compared with each other in three study areas with up to three different discharges. A standardized, relative comparison of the models shows that after successful calibration at measured water levels, the associated 2D/1D and 3D/1D ratios are almost unity, while bed shear stresses in the 3D models are only about 62–86% of the simulated 1D values and 90–100% in the case of 2D/1D. Reasons for this can be found in different roughness definitions, in simplified geometries, in different calculation approaches, as well as in influences of the turbulence closure. Moreover, decreasing 3D/1D ratios of shear stresses were found with increasing discharges and with increasing slopes, while the equivalent 2D/1D ratios remain almost unchanged. The findings of this study should be taken into account, particularly in subsequent sediment transport simulations, as these calculations are often based on shear stresses.

**Keywords:** hydraulic models; numerical simulations; flow resistance; bed shear stress

#### **1. Introduction**

The interaction between gravity and flow resistance is the mechanism underlying both fundamental hydraulic processes and river morphological changes [1,2].

Following Powell [2], flow resistance is composed of the boundary resistance, vegetation resistance [3–7], channel resistance [8–12], spill resistance [13–18] and sediment transport resistance [19–25] and affects bed shear stress and energy dissipation. The challenge in river hydraulics is to achieve the complexity of flow resistance through physical relationships and engineering approaches to gain knowledge about flow velocities, water depths, bed shear stresses and energy losses in channels [2].

Classical approaches for the calculation of flow variables under the influence of flow resistance are the equations of Darcy–Weisbach, Chezy and Manning–Strickler [26–28]. These are empirical methods, which were developed from observations in channels in which the flow resistance was determined almost exclusively from the boundary of the flow—much like a block sliding down a plane [29]. Attempts for modifications of these classical approaches are particularly found for relatively rough beds in numerous subsequent research works [7,25,30–37]. However, the Manning–Strickler equation

and the related bed roughness values of Manning's *n* and Strickler's *kSt* are still used in engineering practice and particularly in numerical modelling.

For the determination of the bed roughness value, *kSt*, Strickler [38] developed an equation considering the grain size distribution of the bed material. On the basis of this investigation Meyer-Peter and Müller [39] as well as Jäggi [40] proposed improved approaches by separately considering the surface layer of the bed material or the variable water depths in the channel, respectively. In addition, guidelines for the determination of *kSt*-values [41–43] have been developed for different channel types. However, the Strickler equation was developed for simplified channel geometries without considering lateral heterogeneities. To address these differences in compound and composite cross sections different approaches were elaborated in various research works [44–47]. Of particular interest in this context was vegetation in riverine landscape, which has a significant impact on flow resistance [48–50].

Despite some uncertainties, including the water depth dependencies of *kSt*-values, related to the use of Strickler's value, this approach is employed in a substantial percentage of all numerical simulation studies for the definition of the flow resistance. Most well-known 1D numerical models based on the Saint Venant equations (e.g., HEC-RAS, MIKE 11, etc.) and also 2D numerical models employing the depth-averaged Reynolds-equations (e.g., Hydro\_AS-2D, MIKE21, TELEMAC-2D, etc.) define the roughness by Strickler's or Manning's value. In contrast, most 3D numerical models (e.g., Delft3D, Fluent, RSim-3D, SSIIM) solving the fully three-dimensional Reynolds-equations with different turbulence closures, typically permit the application of the equivalent sand roughness *ks* [51], which is based on the sediment size and also allows accounting for the average height of bed forms. To date, numerous studies [52–61] have been published based on the mentioned numerical models to calculate hydrodynamic properties and their impacts on rivers.

Moreover, the modelled hydrodynamics are often used for predicting the sediment transport as well as morphological changes in rivers. In this context, the bed shear stress is of particular importance, due to the fact that well-known bed load transport equations [39,62–65] are based on this variable and are used both in analytical calculations [66,67] as well as numerical simulations [68–71].

These issues show the paramount importance of roughness values in the calculation of bed shear stresses in the context of numerical simulations. Independent of the dimensionality of the models this roughness value usually represents the main opportunity for parameterizing flow resistance in numerical simulations. Moreover, it is the main variable parameter of choice to achieve successful model calibration implying reasonable difference between measured and modelled water surface elevations as well as flow velocities.

Thus, the central research question of this paper is whether valid calculations of water surface elevations in different model dimensions (1D, 2D and 3D) imply a consistent depiction of flow velocities and near-bed flow forces. Of particular interest is a comparison of simulated bed shear stresses of different models, due to their substantial influence on the simulation of bed load transport. The influence of certain characteristic morphologies including artificial bathymetries, pool-riffle sections as well as river bends is investigated on the basis of three different case studies.

#### **2. Materials and Methods**

#### *2.1. Case Studies*

The locations of the different case studies are depicted in Figure 1a. Case study 1 (C.1) "S-Shaped Trapezoidal Channel" (Figure 1b) represents controlled conditions within a laboratory experiment and was selected to focus on the influence of river bends without any morphological structures. C.2 "Ybbs" (Figure 1c) was chosen to investigate the effects of a pool-riffle section in a stretched natural river site [72]. The complex morphology of the second natural river site C.3 "Sulm" (Figure 1d) combines the different characteristics representing a river meander with pool-riffle sections [72].

**Figure 1.** Location of case studies in Austria (**a**); Overview of case studies including cross-sections used for a generalized comparison (grey and black lines) and the indicated cross-section for detailed analysis (encircled): Design of laboratory experiment representing C.1 (**b**), bed elevations of natural river site "Ybbs" representing C.2 (**c**) and bed elevations of natural river site "Sulm" representing C.3 (**d**).

### 2.1.1. S-Shaped Trapezoidal Channel (C.1)

The laboratory experiment was developed at the Hydraulic Engineering Institute of the University of Innsbruck, consisting of a physical laboratory channel 3.0 m in length followed by two consecutive bends with a mean radius of 4.0 m and apex angle of 60 each, and finally an exit channel being 2.8 m long with a bed slope of 0.005. The banks are sloped by a ratio of 2:3. The bed sediments are characterized by a mean diameter *dm* of around 5 mm. One set-up of the physical experiment, with a discharge of 0.09 m<sup>3</sup> s<sup>−</sup>1, was part of the validation study of the numerical model RSim-3D [73], which is also used within this research. Measurements of the water surface elevation gained during the experiment form the basis for calibrating the numerical models.

#### 2.1.2. Ybbs (C.2)

The river site C.2 "Ybbs" (Figure 1c) has a total length of around 300 m with an average bed slope of 0.0055 and is classified as a straight river due to natural morphological constraints [74]. The bed sediments are characterized by a mean diameter *dm* of around 70 mm [72], excluding the riffle section (located between cross-section 21 and 31) with a *dm* of around 93 mm. The hydrology is defined by conditions during low flow (*QL* = 1.6 m3 s<sup>−</sup>1), mean flow (*QM* = 4.5 m<sup>3</sup> s<sup>−</sup>1), and a flood with a return period of one year (*HQ*<sup>1</sup> = 64.0 m<sup>3</sup> s<sup>−</sup>1). In the course of a monitoring campaign in 2006, the bathymetry as well as the water surface elevation were measured at low flow conditions in a local reference system and are the basis for calibrating the numerical models within this study. The calibrated models are further used in the course of the comparison study including the simulation

of *QM* and *HQ*1. Due to the fact that the valley margins of C.2 are close to the river bed, the inundation areas are greatly reduced, thus this river site is particularly suitable for a sensitivity analysis of flow resistances during floods.

#### 2.1.3. Sulm (C.3)

The river site C.3 "Sulm" (Figure 1d) has a total length of around 500 m with an average bed slope of 0.0010 and is classified as a meandering river [74] with a mean diameter *dm* of the surface layer of around 42 mm [72]. The hydrology is characterized by conditions during low flow (*QL* = 2.7 m<sup>3</sup> s<sup>−</sup>1), mean flow (*QM* = 9.0 m<sup>3</sup> s<sup>−</sup>1) and a flood with a return period of one year (*HQ*<sup>1</sup> = 95.0 m<sup>3</sup> s<sup>−</sup>1). During a monitoring campaign in 2002, the bathymetry as well as the water surface elevation were measured at low flow conditions and are the basis for calibrating numerical models within this study. Additionally, *QM* is simulated with the calibrated models in the course of the comparison study.

#### *2.2. Basic Principles*

2.2.1. Empirical Relationship between the Equivalent Sand Roughness *ks* and the Strickler Value *kst*

Strickler [38] established the eponymous roughness value *kSt* considering the grain size distribution of the bed material, which is given by

$$k\_{\rm St} = \frac{a}{\sqrt[6]{k\_S}}\tag{1}$$

where *a* is an empirical parameter and *kS* is the equivalent sand roughness, which is a characteristic value of the grain size distribution. Following this approach Meyer-Peter and Müller [39] as well as USACE [43] adjusted *a* and *kS* (Table 1). According to Table 1, *kS* ranges between grain size diameters of *d*<sup>50</sup> and *d*90, for which 50% and 90% of the cumulative mass are finer, respectively. Jäggi [40] retained the fundamental concept of Strickler, whereby he developed a new approach for calculating *a*

$$a = (2.5 \, g^{0.5} \, \ln \left( 6.1 \, Z \right) / Z^{0.16}) \left( 1 - e^{-0.02 Z / I^{0.5}} \right)^{0.5} \tag{2}$$

which depends on the relative roughness *Z* = *h*/*kS*—including *kS* and the water depth *h*—and the bed slope *I*, resulting in a range of possible *a* values (Table 1).


**Table 1.** Overview of parameters, *a* and *kS*.

Independent of the different approaches, the equations are applicable only to medium-range values for the Strickler parameter. Naudascher [42] specified the applicability of these equations by defining a range of typical Reynolds numbers in rivers between 20 < 4*R*/*kS* < 1000 where *R* denotes the hydraulic radius [73].

#### 2.2.2. Calculation of Bed Shear Stress

Before going into the details of any results, the numerical principles of bed shear stress calculations in different models are demonstrated. The classical approach in hydraulic engineering for determining the average bed shear stress *τ* in a cross-section is known as the Du Boys equation and is given by

$$
\pi = \rho \text{g} \, R I\_f \tag{3}
$$

where *ρ* is the density of water, *g* is the gravitational acceleration, *R* is the hydraulic radius in a cross-section and *If* is the friction (or energy) slope. Due to the fact that in very broad channels, characterized by a channel width *B* higher than 30 times the water depth *h*, the hydraulic radius *R* can be replaced by the *h*, the Du Boys equation is simplified to

$$
\pi = \rho g h I\_f, \text{ for } B \ge \text{30h} \tag{4}
$$

1D numerical models usually employ the approach of Du Boys (Equation (3)).

This equation, however, is also applied in the case of 2D numerical models for determining bed shear stresses. Nevertheless, the calculation is improved due to the fact that hydraulic information (water depth *h*, flow velocity *v*) is available in each computational node leading to a higher resolution of results. The friction slope *If* is expressed by the Strickler-equation

$$I\_f = \frac{v^2}{k\_{St}^2 \mathcal{R}^{\frac{4}{3}}} \tag{5}$$

Considering the simplification of replacing the hydraulic radius *R* by the water depth *h* and taking into account Equation (5), the approach for calculating bed shear stresses in the case of 2D numerical models is given by

$$
\pi = \frac{\rho g v^2}{k\_{St} ^2 h^{\frac{1}{5}}} \tag{6}
$$

A different approach to 1D and 2D modelling based on the law of the wall is usually employed in the case of 3D numerical models. A well-established method is given by [75]

$$
\pi = \rho \mathbb{C}\_{\mu}{}^{\frac{1}{4}} k\_{\rho}{}^{\frac{1}{2}} \frac{\mathcal{V}\_{p,t}}{\upsilon^{+}} \tag{7}
$$

where *C<sup>μ</sup>* is an empirical constant of the k−ε model, *kP* denotes the turbulent kinetic energy (which can also be interpreted as turbulent flow velocity *vkP* <sup>=</sup> <sup>√</sup>*kP* at the near-wall node), *vp*,*<sup>t</sup>* is the near-wall tangential velocity and *v*<sup>+</sup> is representative of the velocity profile at turbulent flow conditions, according to the law of the wall [76] given by

$$
v^{+} = \frac{1}{\kappa} \ln \left(\frac{y}{y\_k}\right) \tag{8}$$

where *κ* is the Von Karman constant which takes the value of 0.41, *y* is the wall distance and *yk* is a parameter dependent of the equivalent sand roughness of Nikuradse (*ks*):

$$y\_k = k\_s^{-8.0\kappa} \tag{9}$$

Equations (8) and (9) describe the dependence of *v*<sup>+</sup> on the geometry as well as on the roughness. Considering that the square root of *kP* results in a velocity depending on the turbulence and *vp*,*<sup>t</sup>* as near-wall tangential velocity, the similarity between Equation (6) and (7) becomes visible. However, the bed shear stress in the case of the applied 3D numerical model is dependent on the turbulence as well as on the tangential velocity, whereas only the tangential velocity is taken into account in the case of the 2D numerical model.

### *2.3. Numerical Setup*

#### 2.3.1. 1.D Numerical Model

As a representative of 1D step-backwater models, the software HEC-RAS Version 4.1 [77] was used within this study. The spatial discretization in this software is provided by cross-sections including information about bed elevations and distances between cross-sections. Manning's value is used for the definition of roughness. The computational procedure for steady flow is based on the solution of the one-dimensional energy equation coupled with the momentum equation in situations where the water surface profile varies rapidly. For unsteady flow, the fully dynamic 1D Saint Vernant equation is solved using an implicit finite difference method.

#### 2.3.2. 2D Numerical Model

The 2D numerical model Hydro\_AS-2d [78], employing the depth-averaged Reynolds-equations, was applied in this study. The roughness is defined by Strickler's or Manning's value and the hydraulic conditions are calculated on a linear grid by a finite volume approach. The Surface-Water Modelling System (SMS) is used as a pre- and post-processing tool of the used 2D numerical model. The convective flow of the two-dimensional model is based on the Upwind-scheme by Pironneau [79] and the discretization of time is done by a second order explicit Runge Kutta method. The simulations were performed with the parabolic eddy viscosity model and a constant turbulent viscosity coefficient of 0.6.
