Dune Contribution to Flow Resistance in Alluvial Rivers

: One of the most relevant features of alluvial rivers concerns ﬂow resistance, which depends on many factors including, mainly grain resistance and form drag. For natural sand-bed rivers, dunes furnish the most signiﬁcant contribution and this paper provides an insight on it. To achieve this aim, momentum balance equations and energy balance equations are applied to free ﬂow in alluvial channels, assuming hydrostatic pressure distribution over the cross sections conﬁning the control volume, which includes a reference bed form pattern. The resulting equation in terms of energy grade accounts for an empirical bed form drag coe ﬃ cient resulting from the actual ﬂow pattern and bed form geometry. The model has been validated using a large selection of ﬁeld data and it seems somewhat sensitive to the dune geometry and to the Nikuradse equivalent roughness, whereas it is shows greater sensitivity to the adopted grain surface resistance formula (e.g., Manning–Strickler formula).


Introduction
In alluvial channels, flow resistance is the result of several factors, including the presence of bends, vegetation, local acceleration due to flow unsteadiness, and channel geometry variations, in addition to bed surface roughness and bed forms drag. Morphological evolution of active channels also affects flow resistance, since it is responsible for the evolution of large-scale bed forms (i.e., dunes and bars). Since the 1950s great efforts have been devoted to river morphodynamic models. Most of them consist of coupled systems of depth-averaged flow mass and momentum equations, a bed evolution equation, as well as a sediment transport formula. Valuable state-of-the-art technologies for river sedimentation and morphology modeling may be found in [1,2]. These conventional morphodynamic models mainly refer to nonlinear shallow water equations, Exner equation, and an empirical formula for sediment transport. Considerable uncertainty derives from the large number of such formulae available in the literature. Moreover, sediment transport in open channel flows includes bedload and suspended load; but, although different mechanisms govern these two modes of transport, a reliable method to account for both of them has not yet been provided [3]. In conventional models further uncertainty arises from the lack of understanding of some fundamental mechanics related to sediment transport. Among them, the effect of bed slope, which has been considered by Maldonado and Borthwick, who recently proposed a simplified two-layer model [4] for bedload-dominated scenarios. Alternatively, morphodynamic models based on two-phase and two-layer approaches may also be considered [5][6][7]. Although they are scientifically more insightful than the conventional approaches, they prove to be mathematically more complex and computationally more demanding.
Even in the case of steady and quasi-uniform flow, confined within the main active channel and involving non-cohesive sediment, most of the existing resistance formulae (e.g., [8][9][10][11]) commonly provide inaccurate predictions of flow resistance. This leads to considerable error in stage-discharge prediction, both in terms of water depth (±20%), flow velocity (±15%) [12] or bed friction, for which the difference between measured and predicted values may be ±50% [11] and, where field data is used, discrepancy is even wider [13,14]). Since the middle of the twentieth century this topic has attracted the interest of many scientists and several criteria-based methods have been proposed, based on linear and non-linear approaches.
In the latter case (e.g., [12,[15][16][17]) the resistance coefficient is kept as a single factor and the method does not require any knowledge of bed form geometry (i.e., empirically-based approach). By contrast, the linear approach originally introduced by Meyer-Peter and Mueller [18] considers an overlapping of effects. In fact, flow resistance depends on forces acting on individual particles (i.e., skin friction) and on forces affected by bed form configuration (i.e., form drag). Among the others Azareh et al. [19] investigated the contribution of form friction to the total friction factor, providing evidence that form friction contributes up to 65% of the total friction factor of a gravel-bed river.
Since the bed form contribution to flow resistance mainly depends on dune geometry and flow conditions, recently big efforts have been devoted to researching bed form geometry. Yalin and Karahan [20] found that the spacing of dunes was dependent on the relative depth (i.e., d 50 /y with d 50 being bed material average grain size and y water depth). Julien and Klaassen [21] made a field investigation on the Meuse River and the Rhine River during large floods, and discovered that dune steepness remains relatively constant with discharge, and suggested a linear proportion between wavelength and water depth. The proposed linear coefficient differs from the one empirically proposed by Yalin [22], which either takes into account flume and river data, or theoretically derived [15]. Agarwal et al. [23] stressed that, for specific range of relative depth, dune spacing may greatly differ from what has been suggested by Yalin [15,24], or by Julien and Klaassen [21]. Aberle et al. [25] used a statistical approach to investigate bed forms during different flood conditions in the Elbe River showing how statistical parameters may be used to predict the flow-dependent bed roughness.
Despite the uncertainties regarding the correct representation of dune geometry, the linear approach to flow resistance in presence of bed forms still remains fascinating, in particular because it offers the possibility to explore each component (i.e., skin roughness and bed form drag) separating the independent parameters involved in the process. Recently, Yang [26] proposed an empirical approach identifying the effective toss length over the dune where the flow is in contact with the bed (i.e., associated with the skin roughness), and the separation zone behind the dune crest, which is responsible for the bed form contribution to flow resistance. An analytical model was introduced by Engelund [27] and Yalin [22]. These authors considered the form drag as a result of a sudden flow expansion on the bed form lee side. More recently, Ferreira Da Silva and Yalin [28] generalized two modes of bed forms drag (by virtue of additivity of losses), also taking into account the ripples superimposed on dunes. In any case, they applied momentum, energy and mass balance equations to a reference bed form, assuming that the effects of a bed form on the flow are comparable to a sudden expansion of a pipe flow.
In this paper, a semi-analytical model for the bed form drag is proposed, considering the effects of a sudden expansion of a free surface flow rather than of a pressure flow. The energy balance equation is applied to a 2D steady flow over a reference dune bed. In order to account for the actual flow and bed pattern, an empirical coefficient is introduced in the bed form drag formula and its dependence on dune geometry is analyzed. The skin resistance is calculated assuming Pandtl-Karman velocity distribution and, according to superposition of the effect, the hydraulic energy gradient is calculated. The model is validated by comparing observed and calculated energy gradients. Since this paper focuses on sand rivers in presence of dunes, field data related to sand streams and large sand rivers are considered, accounting for a large span of different hydraulic and sedimentary conditions.

Flow Resistance
In alluvial rivers with sediment transport, flow resistance is generated by surface roughness and bed forms, which act as a macroroughness. In presence of macroroughness the average velocity distribution and cross section geometry are not irrelevant. In fact, the velocity distribution and its gradient, as well as bed shear stress, vary along the channel reach, accounting for actual geometry and flow field.
In a two-dimensional two-phase steady flow in presence of dunes (see Figure 1), the velocity distribution relates to the following characteristic parameters [15,28]: where U is the mean flow velocity, y the mean flow depth, S the friction slope, ν and ρ the kinematic viscosity and the density of water, respectively, γ s the specific weight of immersed sediment, d s the representative sediment diameter, ∆ the dune height, Λ the dune length, and g the gravity acceleration. Therefore, it must hold a dimensionless general relationship of dependency on the following seven dimensionless parameters: Water 2019, 11, 2094 3 of 16

Flow Resistance
In alluvial rivers with sediment transport, flow resistance is generated by surface roughness and bed forms, which act as a macroroughness. In presence of macroroughness the average velocity distribution and cross section geometry are not irrelevant. In fact, the velocity distribution and its gradient, as well as bed shear stress, vary along the channel reach, accounting for actual geometry and flow field. In a two-dimensional two-phase steady flow in presence of dunes (see Figure 1), the velocity distribution relates to the following characteristic parameters [15,28]: where U is the mean flow velocity, y the mean flow depth, S the friction slope, ν and ρ the kinematic viscosity and the density of water, respectively, γs the specific weight of immersed sediment, ds the representative sediment diameter, Δ the dune height, Λ the dune length, and g the gravity acceleration. Therefore, it must hold a dimensionless general relationship of dependency on the following seven dimensionless parameters: s s y U y U S d g y g y y 2 , , , , , , 0 Assuming a constant value for the relative immersed weight of sediment in natural channels, and introducing the following dimensionless variables (i.e., relative submergence Z, Reynolds' number Re and Froude number F): ; the following relationship for the global friction slope is obtained: Assuming a constant value for the relative immersed weight of sediment in natural channels, and introducing the following dimensionless variables (i.e., relative submergence Z, Reynolds' number Re and Froude number F): the following relationship for the global friction slope is obtained: So far, no theory analytically derives such a function, even for the simple reference case of steady, uniform 2D flow in a steady sediment transport condition. A possible approach consists in separating the global flow resistance into two individual contributions: surface roughness and macroroughness (i.e., bed form). Meyer-Peter and Muller [18] and Einstein [29] proposed a linear superposition assuming that the total bed shear stress (τ) could be separated into two contributions-plane bed stress (τ ) and additional shear stress-associated with the bed forms (τ"): According to the assumed superposition of the effects, τ' is equal to the shear stress acting on a sand-bed plane having the same grain-size distribution and the same hydrodynamic condition (i.e., velocity and flow depth). Later, Engelund [27] proposed that the total mean energy gradient per unit length of the stream (S) could be considered as a contribution of a friction loss (S ) due to the stress acting along the dune surface and of a cumulative loss due to a sudden flow expansion just downstream from the dune crest (S"): and Equation (5) leads to: Therefore, the total head loss (∆H) over a stream reach of length L results from head loss due to the grain resistance (∆H ) and head loss related to the dune (∆H"):

Grain Contribution to Flow Resistance
The characteristics of the roughness elements located along the wetted perimeter affects the vertical profile of the downstream component and the resistance to flow. For fully developed turbulent flow over a rough sand-bed plane (∆/y = Λ/y = 0), Re and F become irrelevant and the stream velocity profile can be expressed by a logarithmic law: where u is the flow velocity at elevation z above the bed, k = 0.4 is the Von Karman's constant, u * = τ /ρ is the shear velocity related to the skin roughness and k s is the equivalent grain roughness. The vertically-averaged velocity U corresponds to the local velocity u at the relative depth z/y = e −1 = 0.368. By integration over the flow depth, it results: where χ is the Chezy's coefficient and C the dimensionless conveyance coefficient related to the skin roughness k s '. According to different Authors, k s is proportional to a characteristic sediment size d s , with subscript "s" being equal to 35, 50, 65, 84, 90 (i.e., the percentage of the finer particle size distribution by weight). Typically, k s ranges between 1.25 d 35 and 5.10 d 84 [1][2][3][4][5][6][7][8][9]16,[30][31][32][33], whereas Millar [34] concluded that in gravel streams there is no significant difference between using d 35 , d 50 , d 65 , d 84 or d 90 . In the present analysis, the mean grain size diameter d 50 is assumed as a characteristic sediment size. Different values for k s were considered and the most appropriate value resulted k s = 2·d 50 . In terms of slope friction, Equation (10) results: which is consistent with Equation (4).

Sand Dune Contribution to Flow Resistance
In a two-dimensional two-phase steady flow in presence of dunes (see Figure 1), the momentum balance equation is applied to the reference control volume of a portion of 2D dune bed. The control volume is bounded between cross section 1 upstream, in correspondence to the crest of the dune, and cross section 2, downstream, in the trough zone where the streamlines are assumed to be parallel to the bed: where P represents the pressure vector acting on the boundary of the fixed control between the two consecutive cross sections 1 and 2, G represents the mass force vector acting on the control volume and M represents the momentum flux vector through the boundary. The x-component of the previous equation gives: being P 1 , P 2 , P L the pressure force acting on the upstream cross section, downstream cross section and on the lee side of the dune, respectively. T L is the shear stress on the lee side of the dune, M 1 and M 2 the momentum flux across the upstream and downstream cross sections 1 and 2, respectively, ω the angle between the lee dune side and the horizontal, and σ the average bed slope.
In a sand-bed river channel, the active longitudinal component of mass force G·sin(σ) can usually be disregarded with respect to the other forces. Since streamlines are parallel to the bed at cross sections 1 and 2, and the velocity in the flow separation zone may be considered null, hydrostatic pressure holds at both cross sections 1 and 2. The surface force component related to the shear stress (T L cos(ω)) is already included in the grain roughness component of the total resistance to flow. Introducing the momentum coefficient β (accounting for non-uniform velocity distribution over the cross section), the momentum balance equation per unit width of the channel becomes: where y i is water depth at cross section (i = 1,2). The energy balance equation applied between cross sections 1 and 2 results: where z i (i = 1,2) is the bed level and α i the Coriolis coefficient (i = 1,2) referred to cross sections 1 and 2.
For the sake of simplicity from now on α 1 = α 2 = β 1 = β 2 = 1. The location of the stagnation point is not known a-priori and, taking into account for a mild slope of the toss face of the dune, it may be assumed z 1 −∆ = z 2 . Moreover, considering a reference dune bed pattern having a vertical dune lee side (i.e., ω = 90 • ), Equation (15) becomes: and momentum balance Equation (14) leads to: Disregarding the difference in water levels between cross sections 1 and 2, Equation (16) results: It is convenient to refer to the mean water depth: where y(x) is the local water depth (see Figure 1). For a dune bed-dominated river, the following approximation is considered [15]: Introducing the water discharge per unit width of the channel: and accounting for the mean water depth y (Equation (20)), results: After algebraic simplification: where the dune geometric correction function Γ ∆/y is: It is worth noting that Yalin [22] and Engelund [27], independently, applied a 1-D momentum conservation equation to the same dune bed reference pattern, but they assumed a constant pressure over the cross sections 1 and 2 (i.e., Borda-Carnot's theorem for pressure pipe flow) rather than a hydrostatic pressure distribution, and they obtained: where ∆H P " indicates the bed form energy loss according to the pressure pipe flow approach. In this case: A comparison of Equations (23) and (25) demonstrates the difference between the approach based on free surface flow instead of pipe flow, when the effects of sudden flow expansion in presence of dune is considered (Figure 2). Over the typical range of relative dune depth ∆/y = 0.1 − 0.3 the ratio between the two approaches, in terms of energy loss, varies from around 10 to 6. It supports the use of Equation (23) instead of Equation (25). Using ΔH'' = S''⋅Λ, and after multiplying and dividing by the mean water depth (y), the energy slope related to the dune form drag results:

Empirical Coefficient for Bed Form Drag
Tokyay and Altan-Sakarya [35] tested the local energy losses at negative steps in subcritical open channel flows and demonstrated that energy losses are higher for inclined steps than for abrupt vertical steps. Van der Mark [36], following experimental results of flow over a single bed form with a vertical lee face concluded that an empirical coefficient should be applied to Equation (25) in order to take into account that flow velocity distribution is not uniform over the cross-section at the crosssections 1 and 2. Engel [37] made experiments on the length of flow separation over artificial dunes in a flume, using fixed bed geometries. He showed that the separation length is independent of the Froude number, but it is dependent on the relative dune height Δ/y. Shen et al. [38] carried out laboratory experiments on rigid bed forms using both a rough and smooth surface. He showed that pressure drag coefficient relative to the bed form depends on the relative dune height, and on the dune slope steepness, whereas flow velocity and grain roughness does not influence it.
Therefore, the empirical correction coefficient κ, which is a function of the dune steepness δ = Δ/Λ, is introduced into Equation (27), to also take into account any differences between the theoretical and the actual flow field and bed form geometry, as a consequence of the assumptions so far introduced: Using the breakdown of the energy slope S into its two components related to skin roughness S' and bed form drag S'' (Equation (6)), Equation (28) gives the following expression for the empirical coefficient κ: Once the bed dune geometry is known, considering the measured energy slope S, and after S' is calculated by Equation (11) (assuming ks = 2·d50), it is possible to determine the empirical coefficient Using ∆H" = S"·Λ, and after multiplying and dividing by the mean water depth (y), the energy slope related to the dune form drag results:

Empirical Coefficient for Bed Form Drag
Tokyay and Altan-Sakarya [35] tested the local energy losses at negative steps in subcritical open channel flows and demonstrated that energy losses are higher for inclined steps than for abrupt vertical steps. Van der Mark [36], following experimental results of flow over a single bed form with a vertical lee face concluded that an empirical coefficient should be applied to Equation (25) in order to take into account that flow velocity distribution is not uniform over the cross-section at the cross-sections 1 and 2. Engel [37] made experiments on the length of flow separation over artificial dunes in a flume, using fixed bed geometries. He showed that the separation length is independent of the Froude number, but it is dependent on the relative dune height ∆/y. Shen et al. [38] carried out laboratory experiments on rigid bed forms using both a rough and smooth surface. He showed that pressure drag coefficient relative to the bed form depends on the relative dune height, and on the dune slope steepness, whereas flow velocity and grain roughness does not influence it.
Therefore, the empirical correction coefficient κ, which is a function of the dune steepness δ = ∆/Λ, is introduced into Equation (27), to also take into account any differences between the theoretical and the actual flow field and bed form geometry, as a consequence of the assumptions so far introduced: Using the breakdown of the energy slope S into its two components related to skin roughness S and bed form drag S" (Equation (6)), Equation (28) gives the following expression for the empirical coefficient κ: Once the bed dune geometry is known, considering the measured energy slope S, and after S is calculated by Equation (11) (assuming k s = 2·d 50 ), it is possible to determine the empirical coefficient κ. A selection of 132 field data collected on 7 different sand-bed rivers with dune bed forms is considered (see Table 1). It refers to the Savio and Fiumi Uniti rivers [39], Calamus River [40], Missouri (data from Shen 1978, as reported by Brownlie [41]), Jamuna and Parana [42], Bergsche Maas (data from Adriaanse 1986, as reported by Julien [42]), and Meuse [42]. The database reports for each data set the number of measurements N, and the range of the following observed parameters: water discharge Q, mean depth y, mean flow velocity U, measured energy gradient S, Froude number F, mean grain size diameter d 50 , dune height ∆ and dune length Λ. Despite the relative scatter of field data, Figure 3 shows how the drag coefficient decreases with increased dune steepness. Some outliers are present, corresponding to few data recorded on the Jamuna and Missouri rivers [42], characterized by very low values of dune steepness (i.e., ∆/Λ < 0.01), and filtered in the fitting procedures. The best fitting equation is: with m = 0.053 and n = −0.2. Table 1. Summary of field data used to determine empirical coefficient κ (Equation (29)). Each column reports maximum and minimum value.  κ. A selection of 132 field data collected on 7 different sand-bed rivers with dune bed forms is considered (see Table 1). It refers to the Savio and Fiumi Uniti rivers [39], Calamus River [40], Missouri (data from Shen 1978, as reported by Brownlie [41]), Jamuna and Parana [42], Bergsche Maas (data from Adriaanse 1986, as reported by Julien [42]), and Meuse [42]. The database reports for each data set the number of measurements N, and the range of the following observed parameters: water discharge Q, mean depth y, mean flow velocity U, measured energy gradient S, Froude number F, mean grain size diameter d50, dune height Δ and dune length Λ. Despite the relative scatter of field data, Figure 3 shows how the drag coefficient decreases with increased dune steepness. Some outliers are present, corresponding to few data recorded on the Jamuna and Missouri rivers [42], characterized by very low values of dune steepness (i.e., Δ/Λ < 0.01), and filtered in the fitting procedures. The best fitting equation is: with m = 0.053 and n = −0.2.  Table 1). Figure 3. Empirical coefficient κ as a function of dune steepness δ = ∆/Λ. (Data A-H see Table 1).
According to Equation (29), the estimated drag coefficient k depends on the skin energy slope (S ) which, in turn, depends on the equivalent roughness k s . Among the selected field data (river A-G, Table 1), Fiumi Uniti, Savio, Missouri and Meuse River reported information about sediment gradation or d 90 . Figure 4 shows a comparison between the drag coefficient κ obtained considering equivalent roughness k s = 2·d 50 or k s = 3·d 90 .
The general trend remains confirmed, and it may be concluded that coefficient κ is almost independent of the equivalent skin roughness adopted to estimate the grain resistance contribution S . gradation or d90. Figure 4 shows a comparison between the drag coefficient κ obtained considering equivalent roughness ks = 2·d50 or ks = 3·d90.
The general trend remains confirmed, and it may be concluded that coefficient κ is almost independent of the equivalent skin roughness adopted to estimate the grain resistance contribution S'.  Table 1).
The dune steepness may be expressed by means of relative dune height and relative dune length, hence Equation (30) becomes: The previous expression is particularly convenient, since according to several authors [10,15,24], the relative dune length may be considered as a constant. Thus, the drag coefficient remains as a function of the relative dune height and, after substituting Equation (31) in Equation (28), the bed form contribution to the energy slope is:  Table 1).
The dune steepness may be expressed by means of relative dune height and relative dune length, hence Equation (30) becomes: The previous expression is particularly convenient, since according to several authors [10,15,24], the relative dune length may be considered as a constant. Thus, the drag coefficient remains as a function of the relative dune height and, after substituting Equation (31) in Equation (28), the bed form contribution to the energy slope is: therefore, once the constant value for the relative dune length has been defined, the energy slope related to the dune drag results as a function only of the Froude number and of the relative dune height. Introducing a constant value for the relative dune length Λ/y, the best fitting coefficients m and n may change from the previous values (m = 0.053 and n = −0.20) as will be discussed here in the model validation section.

Model Validation and Sensitivity Analysis to the Bed Form Geometry and Skin Roughness
In order to validate the proposed model a large selection of 491 field datasets is considered (Table 2), including those already used for the preliminary calibration of the empirical coefficient κ ( Table 1). The Table 2    The model is validated comparing the observed and the estimated total energy slope (i.e., S = S' + S" Equation (6)), on the basis of the hydraulic parameters listed on Tables 1 and 2.
In order to assess energy slope S, S is calculated by Equation (11) assuming k s = 2·d 50 , and S" is calculated using Equations (32) and (24): It results a parametric function of Froude number, dune geometry, and parameters m and n referred to the dune drag coefficient function (i.e., Equation (30)): Therefore the energy slope related to the bed form drag (S") in Equation (33) involves calculation of few parameters on the right-hand side of Equation (33). There are many empirical relationships in literature related to the bed form geometry. In the present work, the following have been considered (Van Rjin [10] and Karim [44], Equation (35) and Equation (36) respectively): It is worth noting that Equation (36) gives relative dune height as a function of only hydraulic parameters, and equivalent skin roughness, because of introducing Equations (11) and (35) in Equation (36). Consequently, the best fitting parameters in Equation (32) result: m = 0.07 n = −0.19 (37) which are slightly different from those preliminary assessed (m = 0.053, n = 0.20). Figure 5 shows the comparison between calculated and measured energy slope. Accounting for the observed water depth y and Froude number F (see column 5 and 8 in Tables 1 and 2), the following equations are involved in energy slope calculation: S is calculated with Equation (11) assuming k s ' = 2·d 50 (see column 9 in Tables 1 and 2); S" is calculated with Equation (33) along with parameters m and n from Equation (37), and eventually bed form steepness using Equations (35) and (36).
It is worth noting that in the range of observed mean flow velocity U, ranging 0.5-2.0 m/s, the skin roughness contribution gives a shear velocity u * ' = (g S y) 0.5 in the range of about 0.02-0.06 m/s, and the dimensionless conveyance coefficient C = (g S y) 0.5 respectively correspond to approximately 7-14 and 22-27.
Among the field data, 93.5% lies within the ±30% discrepancy band error and 68.2% of the measured energy slope lies within the ±20% band error. Considering the large uncertainties in field measurements, this can be considered as a good agreement. and the dimensionless conveyance coefficient C = (g S y) 0.5 respectively correspond to approximately 7-14 and 22-27. Among the field data, 93.5% lies within the ±30% discrepancy band error and 68.2% of the measured energy slope lies within the ±20% band error. Considering the large uncertainties in field measurements, this can be considered as a good agreement. Skin roughness, as well as the adopted resistance formula, affects the grain flow resistance component and the related energy gradient calculation. Following the linear approach, the contribution due to the bed surface is equal to that of plane bed with identical hydraulic condition and sediment characteristic, without any bed form. The concepts of Karman-Prandtl logarithmic velocity distribution and Nikuradse equivalent grain roughness are considered (Equations (9) and (10)), where the latter is proportional to the characteristic grain size [45]. To test the sensitivity of the model, Nikuradse equivalent roughness k s = 1.0·d 50 [46] and a different approach, reflecting Manning-Strickler Equation (38), are considered: n = 0.0416·d 50 0.165 (38) where n is Manning-Strickler coefficient related to the grain roughness. Hence: and Equation (39) is used instead of Equation (11) to calculate the energy slope skin roughness contribution when Manning-Strickler formula is considered.
On the other hand, the assessment of bed form drag contribution involves dune geometry, in terms of dune steepness δ = ∆/Λ or, equivalently, the dimensionless dune height ∆/y and dimensionless dune length Λ/y (see Equation (33)). Different approaches proposed by Yalin [15,24], via empirical consideration based on field data, and reflecting theoretical consideration (Equations (40) and (41), respectively) were considered to explore the sensitivity of the proposed model to the dune geometry.
Λ /y = 6.28 (41) Table 3 reports the results of sensitivity analysis. Decreasing the value of Λ/y leads to a decreasing of the model accuracy. In fact, the field data included within the ±30% error band decrease from 93.5% to 91.3%, whereas the data included within the ±20% error band remains almost stable. By reducing the Nikuradse equivalent roughness (i.e., k s = 1 d 50 ), field data included within the ±30% error band reduce from 93.5% to 90.7%, and the data included within the ±20% error band decrease from 68.4% to 65.5%. On the contrary, using the Manning-Strickler resistance formula (Equation (38)) only 88.2% of the predicted energy slope are within the ±30% band error and the data included within the ±20% band error are reduced to 65.5% (see Figure 6).

Conclusions
The focus of the present paper is the bed form contribution to the flow resistance in natural sandbed rivers with dune bed forms, assuming a linear separation approach between skin roughness and dune bed form drag. To this aim, the momentum balance and the energy balance equations are applied to 2D flow in open channel, assuming hydrostatic pressure distribution over the cross sections bounding the control volume, which includes a reference bed form pattern. The related energy loss deviates from that derived by Borda-Carnot's applied theorem. The resulting equation in terms of energy grade takes into account an empirical correction factor due to the actual flow field and bed form pattern. The empirical coefficient results as a power function of the dune steepness and is slightly dependent on the Nikuradse equivalent roughness of the grains. Thus, the dune contribution to the energy slope remains an explicit function of dune geometry, in terms of relative length and height, as well as Froude number. The model has been validated using 491 field measurements referred to sand rivers in presence of dunes, showing a good agreement. A sensitivity analysis on estimated energy slope was carried out in terms of dune geometry and the skin roughness

Conclusions
The focus of the present paper is the bed form contribution to the flow resistance in natural sand-bed rivers with dune bed forms, assuming a linear separation approach between skin roughness and dune bed form drag. To this aim, the momentum balance and the energy balance equations are applied to 2D flow in open channel, assuming hydrostatic pressure distribution over the cross sections bounding the control volume, which includes a reference bed form pattern. The related energy loss deviates from that derived by Borda-Carnot's applied theorem. The resulting equation in terms of energy grade takes into account an empirical correction factor due to the actual flow field and bed form pattern. The empirical coefficient results as a power function of the dune steepness and is slightly dependent on the Nikuradse equivalent roughness of the grains. Thus, the dune contribution to the energy slope remains an explicit function of dune geometry, in terms of relative length and height, as well as Froude number. The model has been validated using 491 field measurements referred to sand rivers in presence of dunes, showing a good agreement. A sensitivity analysis on estimated energy slope was carried out in terms of dune geometry and the skin roughness adopted model. In particular, different relative dune length approaches were considered. Decreasing the relative dune length more than 30% of the reference value does not substantially degrade the model accuracy, whereas the model is relatively more sensitive to the Nikuradse equivalent roughness, and to the adopted resistance formula (e.g., Manning-Strickler formula).
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest. Manning coefficient related to grain roughness u local flow velocity u * shear velocity related to the skin roughness U mean flow velocity P 1 , P 2, P L pressure force acting on the upstream, downstream cross section, and on the lee side of the dune, with respect to the control volume of the streamflow P pressure vector acting on the boundary surface of the control volume of streamflow q water discharge per unit with of the channel Re

Abbreviations
Reynolds' number S friction slope S friction slope due to the grain resistance S" friction slope due to the dune drag T L shear stress on the lee side of the dune y mean flow depth y local water depth z vertical elevation above the river bed Z relative submergence α i Coriolis coefficient β momentum coefficient γ s specific weight of immersed sediment Γ ∆/y and Γ P ∆/y dune geometric correction function δ dune steepness ∆ dune height ∆H total head loss ∆H head loss due to the grain resistance ∆H" head loss due to the dune drag ∆H P " head loss due to the dune drag according to pressure pipe flow approach κ empirical correction coefficient Λ dune length ν kinematic viscosity of water ρ density of water σ average river bedslope τ total bed shear stress τ' grain shear stress contribution to the total bed shear stress τ" bedform shear stress contribution to the total bed shear stress ω angle between the lee dune side and the horizontal