Quantifying Hydraulic Roughness From Field Data: Can Dune Morphology Tell the Whole Story?

Hydraulic roughness is a fundamental property in river research, as it directly affects water levels, flow strength and the associated sediment transport rates. However, quantification of roughness is challenging, as it is not directly measurable in the field. In lowland rivers, bedforms are a major source of hydraulic roughness. Decades of research have focused on dunes to allow parameterization of roughness, with relatively little focus on field verification. This study aims to establish the predictive capacity of current roughness predictors, and to identify reasons for the unexplained part of the variance in roughness. We quantified hydraulic roughness based on the Darcy‐Weisbach friction factor (f) calculated from hydraulic field data of a 78 km‐long trajectory of the river Rhine and river Waal in the Netherlands. This is compared to predicted roughness values based on dune geometry, and to the spatial trends in the local topographic leeside angle, both inferred from bathymetric field data. Results from both approaches show the same general trend and magnitude of roughness values (0.019 < f < 0.069). Roughness inferred from dune geometry explains at best 31% of the variance. Efforts to explain the remaining variance from statistics of the local topographic leeside angles, which supposedly control flow separation, were unsuccessful. Unexpectedly, multi‐kilometer depth oscillations explain 34% of the total roughness variations. We suggest that flow divergence associated with depth increase causes energy loss, which is reflected in an elevated hydraulic roughness. Multi‐kilometer depth variations occur in many rivers worldwide, which implies a cause of flow resistance that needs further study.

, the predictor of Van Rijn (1984) is also calibrated on field data, and the predictor by  is based on numerical experiments. Unfortunately, the predictive value of those equations is often limited (Warmink et al., 2013), and there is a large variation in predicted roughness among alternative formulations (Warmink, 2011).
The limited predictive value can be attributed to two factors. First, hydraulic roughness is, just like bed shear stress, scale dependent (Vermeulen et al., 2013). In other words, roughness from a point measurement will be different from roughness integrated over a dune field, over a cross-section, or over a longitudinal transect Hidayat et al., 2011;Hoitink et al., 2009;Sassi et al., 2011). Roughness predictors are often empirical equations derived from laboratory flume studies and are based on an integration over 2D dunes in a dune field, with limited consideration of the scale-dependence. Second, relations drawn under laboratory conditions are translated to field situations by nondimensionalizing the results. However, recent research has shown that dunes in laboratory settings have a steeper leeside angle than dunes observed in field conditions (Cisneros et al., 2020;Kostaschuk & Venditti, 2019;van der Mark et al., 2008). The slip face angle, the steepest part of the leeside, and the relative dune height determine the presence and strength of the flow separation zone (Best & Kostaschuk, 2002;, and hence determine the total form roughness (Lefebvre et al., 2013). As a result, low-angle dunes produce less form roughness than angle-of-repose dunes (Kwoll et al., 2016;. Relations drawn under laboratory conditions are therefore not directly applicable to field conditions and may lead to considerable uncertainties in model outcome (Warmink, 2011).
Roughness coefficients in numerical models are often calibrated based on measured water levels and discharges, instead of roughness predictors based on dune geometry. Unfortunately, those calibrated values are only valid for the conditions used for calibration (Klemes, 1986). To improve of operational models, it is essential to identify and quantify the spatiotemporal roughness variation, yet geographical insight into roughness is limited. Beyond using dune roughness predictors, bathymetric field data can be a potential source of knowledge about roughness induced by the river bed. In many lowland rivers, the bathymetry is measured regularly to control the navigation depth. Regular bathymetric measurements can identify spatiotemporal dynamics of bedforms, and disclose information about dune geometry. For example, recent research of de Ruijsscher et al. (2020) shows that geometrical properties of dunes change while they migrate over river banks, thereby exerting a different degree of resistance to the flow. The aim of this study is to establish the predictive capacity of well-established roughness predictors, and to identify reasons for the unexplained part of the variance in estimated roughness. We quantified the Darcy-Weisbach friction factor based on water surface slope measurements, over a trajectory of 78 kilometer of the main branch of the river Rhine in the Netherlands. We compare this to the predicted roughness value based on dune geometry, and to the spatial distribution of various measures of the leeside angle, both inferred from bathymetric field data. Ultimately, the insights could be used to improve roughness quantification in models, such as the state-of-the-art operational numerical flow model of the Rhine branches in the Netherlands (RWS-WVL & Deltares, 2017).

Field Site
The study area is the lowland sand bedded part of the river Rhine with its main distributary, the river Waal. The river Rhine (also called Bovenrijn in Dutch) enters the Netherlands at Spijk at river kilometer (RK) 857, 5 km upstream of Lobith, and bifurcates into two branches (Waal and Pannerdensch Kanaal) at the Pannerdense Kop (RK 867). Our study area was the 78 km-long reach (RK 857-935; Figure 1) of the Dutch part of the river Rhine, including its main branch, the river Waal. We focused on the reach between the Dutch border, and the location where the tidal motion starts to influence the water levels, near the city of Zaltbommel (RK 935).
The discharge entering the Netherlands at Lobith (RK 862) is on average 2,300 m 3 s −1 . It fluctuates significantly between 800 m 3 s −1 during low flow conditions, and reached up to 12,000 m 3 s −1 during a high discharge event in 1995 (Schielen et al., 2007). If the discharge exceeds 4,000 m 3 s −1 , floodplains convey part of the discharge. The main branch, the river Waal, receives about two thirds of the water discharge measured at Lobith (Schielen et al., 2007).
In this study area, the river changes from a river with relatively coarse sand and gravel to a fine sand bedded river. The width of the conducting section of the river channel varies between 220 and 350 m, and generally increases   Ruijsscher et al., 2020), replacing the groynes in the inner bends of the river. They split the river in a main channel and two bank-connected side channels of approximately 90 m width. Most field measurements were taken during the construction of the LTDs, and care should be taken when comparing data in this region from different months. Moreover, around Nijmegen (RK 883-885) and St. Andries (RK 925-928) fixed non-erodible layers in the outer bend (with a width of approximately 150 m) have been constructed (Sloff et al., 2006), and at Erlecom (RK 873-876) bendway weirs with a wavelength of 50 m fix the bed. Downstream of these layers, a long scour hole has developed, on average 3 m deep and 500 m long, and backwater raises water levels upstream. Therefore, the regions with a fixed bed were excluded from our current analysis.

Methods
We describe the field data used for calculation of the Darcy-Weisbach friction factor, and for prediction of dune roughness. Then, we discuss a method to define local leeside angles. To assess how much of the variance in hydraulic roughness was explained by dune roughness predictors and leeside angle statistics, we introduce the coefficient of determination.

Data Availability and Preprocessing
Field data were collected between 2014 and 2016 during periods of low river discharge, varying between 781 and 1,353 m 3 s −1 at Tiel (RK 913; Table 2).

Water Level and Discharge
Gauging stations (at each indicated place name in Figure 1)  Note. Water surface slope (WSS) measurements, with additional SBES bed-level measurements (z), were done at the same time as the three MBES campaigns. From these measurements, pressure slope (S p ) was derived. MBES campaigns measured x, y and z, and two of the surveys measured an additional local topographic leeside angle γ.

Table 2
Overview of Field Data in Cartesian Coordinates roll and yaw (measured 50 times per second by a motion sensor) was used to correct the position. Subsequently, the measurements were related to the water level measured at gauging stations, which was used to correct the water surface slope profiles with a 95% confidence level. The average deviation between the measured water level with the sensor and the gauging stations during a survey was 0.11 m, 0.08 and 0.02 m respectively. Fluctuations of water level during a measurement campaign due to slowly changing river discharge during the measurement period (about 24 hr), were in the order of a few centimeters, and were disregarded. Water surface slope profiles were measured during relatively low discharge conditions. This assured that only the main river channel was active, groynes were not flooded, and roughness imposed by the floodplains did not influence the water surface slope. Simultaneously with the water surface slope measurements, a singlebeam echosounder (SBES) scanned the underlying bathymetry along the same line. Corresponding multibeam echosounding (MBES) measurements in the same period were conducted separately on part of this transect, from the city of Dodewaard (RK 895) to Zaltbommel (RK 935).

Multibeam Echosounding
Data from a multibeam echosounder (MBES) were gridded onto a 1 × 1 m 2 grid by Rijkswaterstaat. Only grid cells with a minimum 10 hits per m 2 were analyzed, but in general a much larger number of data points were collected per cell. The resulting five MBES datasets contain x, y and z Cartasian coordinates. For two datasets an additional processing step was performed, where a surface was fitted through the individual points within the grid cell, resulting in an additional value for slope and orientation per grid cell. These two special campaigns contain data from the Dutch border (RK 857) until Zaltbommel (RK 935), measured in September 2014 and October and November 2014, over the whole river width. The average discharge at Tiel (RK 913) during the field measurements was 1249 and 1353 m 3 s −1 respectively ( Table 2). The surveys took approximately two weeks, wherein the river discharge was relatively constant. The discharge differences at the start and end of the surveys were 157 and 104 m 3 s −1 , respectively.
The remaining three MBES datasets were taken around the same date as the water surface slope measurements (Section 2.1.1) and only comprise of x, y and z coordinates without the additional processing step for slope and orientation. They were limited to river kilometers 895 through 936 (Table 2). Next, all bed-level data were converted from Cartesian (x, y) coordinates to curvilinear coordinates (s, n) maintaining the same spatial resolution (Vermeulen et al., 2014). Herein, s is the longitudinal direction, parallel to the river, and corresponds with river kilometers. n is the cross-sectional direction, where n = 0 m is defined as the river axis, which roughly coincides with the thalweg. Besides transformation of the x, y-coordinates, the vector rotation of the cells was calculated to transform the orientation of the fitted surface to the s-direction.

Grain Size
Grain size samples were taken in 2020 with a Hamon Sampler, in which the upper 25 cm of the river bed was taken. The samples were taken at every 500 m at the center line, and subsequently analyzed with sieve sizes between 63 μm and 90 mm (Reneerkens, 2020). From this, the 50 th and 90 th percentile (D 50 and D 90 ) were determined.

Determining River Geometry From Field Data
The river geometry was parameterized by river width, cross-sectional area, curvature and transverse bed slope. Man-made structures fix the river width and curvature to near-constant values, however, cross-sectional area is dependent on water level. River width W (m) was determined from a polygon following a longitudinal river line through the groyne heads, and was taken constant over time. Under low discharge conditions, such as in this study, this measure was considered to be the discharge carrying section of the river. Assuming a trapezoidal shaped channel with a top width W, where the measured water depth represents the width averaged water level, and noting that the slope of the groynes equals 1/3, the cross-sectional area A (m 2 ) can be calculated. Curvature r (km −1 ) was defined as the inverse of bend radius, following the approach of de Ruijsscher et al. (2020). Finally, transverse bed slope ξ (-) was defined as the slope between the two sides of the main river channel, longitudinally discretized in parts of 50 m.

Determining Hydraulic Parameters From Field Data
The water depth d (m) was calculated by subtracting the bed level z from the water surface level h obtained from the water surface surveys (Section 2.1.1). The bed level was derived from the corresponding bed-level measurements, which were taken simultaneously with the water-surface measurements. We chose to use bed-level measurements from the simultaneously taken SBES measurements, since the corresponding MBES measurements were not available over the full length of the study area. The validity of this procedure was checked by constructing a filtered width-averaged bed level from these three MBES surveys and comparing this with the filtered SBES profiles in the part of the study area where there was data available from both datasets. The datasets showed a very similar large-scale bed-level profile (correlation 95%), hence the use of the SBES data set is acceptable.
River discharge Q (m 3 s −1 ) at Lobith (RK 862), Pannerdense Kop (RK 867) and Tiel (RK 913) was calculated via a multistation rating curve, and was subsequently corrected with a correction factor derived from Acoustic Doppler Current Profiler measurements (a fifth degree polynomial fit; see Figure S1 in Supporting Information S1). Before correction, the discharge at Tiel was generally underestimated (up to 15%) while at Lobith, the rating relations could overestimate the discharge up to 10%. Discharge was set constant between Lobith and the Pannerdense Kop since no major confluences or bifurcations occur. From there until Tiel, discharge was assumed to vary linearly due to several small supply and drainage channels. From Tiel until Zaltbommel, discharge was considered to be constant. At the section of the LTDs, approximately 12% of the discharge was conveyed by the bank connected side channels (Ruijsscher et al., 2019;Sieben, 2020). Width-averaged flow velocity u (m s −1 ) was determined by dividing the discharge in the main channel by the cross-sectional area between the groynes.

Roughness Inferred From Water Surface Slopes
We estimated the Chézy coefficient directly from field measurements according to: where u = depth-width averaged flow velocity (m s −1 ), s = longitudinal river distance (m), d = water depth (m), S 0 = bottom slope (=∂z/∂s, in which z = bed level relative to Amsterdam Ordnance Datum, m), S p = pressure slope (=∂d/∂s), Q is discharge (m 3 s −1 ) and A is cross-sectional area (m 2 ). Appendix A offers a derivation of Equation 1. For uniform flow, the expression reduces to the Chézy equation.
All input parameters were discretized per kilometer, after smoothing with an 8 km LOESS filter (Cleveland, 1979;Cleveland & Devlin, 1988;de Ruijsscher et al., 2018). This span was considered to be the best trade-off between accuracy and resolution of the Chézy coefficient (see Figure S2 in Supporting Information S1). From the Chézy coefficient, we calculated the dimensionless Darcy-Weisbach friction factor according to Silberman et al. (1963): where g = gravitational acceleration (m s −2 ).

Roughness Predicted From Dune Characteristics
Existing roughness predictors based on primary dune characteristics including length (λ), height (Δ) and leeside angle (ϕ) were used to determine form roughness. Those characteristics were determined using a well-established Bedform Tracking Tool (BTT) (Van der Mark & Blom, 2007), following the methodology described by de Ruijsscher et al. (2020).
The MBES bed-elevation data were initially detrended by subtracting a reference surface from the Dutch national water authority based on the minimum depth of the fairway, established for dredging. A longitudinal river profile was constructed at the central river axis, and from this, bedform characteristics were determined. A filter span constant c = 1/6 was chosen to filter out small features, which means that the length on which filtering takes place is 1/6 of the main bedform length. Two bedform lengths of interest were defined (Zomer et al., 2021): 5 m ± 5 7 of 22 (hereafter referred to as superimposed bedforms) and 100 m ± 30 (hereafter referred to as dunes), and the corresponding span values (P0) were used as input for detrending the profile (Figure 2). Bedforms smaller than 5 m could not be detected due to the grid cell size of 1 × 1 m 2 . The span values were based on a spectral analysis to determine the dominant wave lengths in each section. The bedform lengths of interest (estimated bedform length) was also an input parameter for the smoothing of the profile. Based on the zero-crossings profile (Figure 2), dune characteristics were described every kilometer. We only considered the three characteristics which are used in the roughness prediction equations: dune height Δ (m), that is, the vertical distance between dune top and downstream trough, dune length λ (m), that is, the horizontal distance between two subsequent crests, and leeside angle ϕ (°), obtained from a linear fit of the leeside of the dune, excluding the upper and lower 1/6 of the dune height (Van der Mark & Blom, 2007). Dune steepness Δ/λ can readily be inferred from this. Form roughness (f f ) imposed by dunes can be predicted employing several previously developed equations, which commonly require primary dune height and length as input. The total predicted hydraulic roughness (̂ ) consists of form friction (f f ), and grain friction (f g ) (Einstein, 1950). Under the assumption that dunes are the fundamental structures causing form resistance, we calculated the total predicted hydraulic roughness using various approaches described below.
Van Rijn (1984) developed an equation based on calibration of field and lab data.
in which k s consists of roughness height k s f and grain roughness height k s g : in which D 90 is the 90 th percentile of the grain size distribution, and γ d is taken as 0.7 in field conditions. Soulsby (1997) used a different formulation for k s f : Bartholdy et al. (2010) defined k s f as: Lefebvre and Winter (2016) in which: To calculate the total hydraulic roughness, they summed grain roughness and form roughness as in ̂= + . Herein, grain roughness was assumed to only depend on grain size distribution and water depth, and was calculated following the approach of Van Rijn (1984).

Analysis of Leeside Angle Statistics
We explored if hydraulic roughness variability was correlated to leeside angles, by inferring hydraulic roughness information from MBES data avoiding assumptions about bedforms. To do so, we calculated the local topographic leeside angle γ (°) for each square meter in the fairway, based on the two MBES surveys that allow to do so. If the orientation of a slope in the 1 × 1 m tile was directed downstream, defined as within ±30° of the central axis -approximately parallel to the flow direction-, the slope was defined as the leeside.
The MBES data were influenced by side effects such as groynes. To purely focus on the fairway of the river, the groyne influenced part of the river bed was disregarded (Figure 3 and Figure S3 in Supporting Information S1).
To identify the part influenced by groynes, the river was subdivided in sections of 1 km in the stream-wise direction, and 1 m in the transverse direction. In these sections, the gradient of the adjacent mean leeside angles was taken and smoothed with a 20 point moving average filter. If the gradient at a certain river width was larger than an absolute value of 0.19, data at a larger width was removed from the analysis. This threshold of 0.19 coincides with 3 times the standard deviation of the gradient of the mean leeside angle in the central 100 m of the river, which was relatively undisturbed by side effects.
The mean local topographic leeside angle γ, calculated from 1 × 1 m 2 tiles of the river bed, was averaged over the river width and over 1 km along the river. It therefore includes the characteristics of 3D variations in bed geometry over the full width of the non-groyne influenced river bed.

Coefficient of Determination as a Tool for Explaining Variance
The coefficient of determination (R 2 ) shows how much of the variance in hydraulic roughness f is explained by the predictors. in which SS tot is the total sum of squares, SS res the residual sum of squares, f i the i th observation of f, ̄ the mean of the observations f and m i the i th model prediction output. R 2 will have a value below 1. A value between 0 and 1 explains how much of the variance of the observation is explained by the model, and values below 0 indicate that the model performs worse than simply taking the average of the observations. The coefficient of determination, R 2 , is only equal to the square of the correlation, if a linear regression with an offset is applied. In this manuscript, we chose to use the coefficient of determination which shows the correlation also for non-linear relations. Uncertainty is quantified as the standard deviation of the residuals. The residuals were calculated as the difference between the modeled values (the predicted roughness, ̂ ) and the observations (the roughness from water surface slope, f). Figure 4a shows the water surface profiles during the three surveys, indicating a relatively constant water surface slope over distance, for all three discharge conditions. Derived values of f fluctuated between 0.019 and 0.069 (Figure 4b). With a mean of 0.035, the river Waal can be characterized as a natural, winding stream (Fetter, 2001). A low discharge, and correspondingly low flow velocities, generally resulted in a high roughness (f averages 0.037, 0.034, 0.033 at a discharge of 1,271, 880, 781 m 3 s −1 respectively). This is especially visible upstream of the Pannerdense Kop (RK 867), while downstream this apparent relation between roughness and discharge becomes less clear.

Grain Size and Dune Geometry Observations
Grain size changed significantly over the study reach, with a D 50 ranging from 10 mm upstream to 0.7 mm downstream (Figure 5a). This decrease in grain size was reflected in a decrease in grain-related roughness.
Spatial variation in bedform geometry was substantial. Bedform height averaged per kilometer varied between 0.1 and 1.5 m, with an average of 0.7 m (Figure 5c). Between Lobith and Pannerdense Kop (RK 857-867), an almost flat bed was observed. Dune height was almost zero in the first 8 kilometers, and increased downstream, but did not exceed 0.7 m (excluding the bendway weir at Erlecom (RK 873-876). Relatively constant dune heights, lengths and leeside angles were observed between river kilometers 885 and 915. Leeside and stoss side angles followed the same spatial trends as dune height and length. Throughout our whole research transect, smaller dunes were imposed on the primary dunes. Those superimposed dunes, or secondary dunes, were on average 0.1 m high and 10 m long, and were clearly distinguishable from the primary dunes being on average 0.7 m high and 55 m long. Additionally, in certain locations the bed geometry was imposed by the fixed layers. For example, in Erlecom (RK 873-876) the sine shaped fixed layer was clearly visible by a deviating "dune" length of 50 m and a height of 1.5 m (Figures 5d and 5e).
Temporal variation in dune geometry was smaller than spatial variation. A higher discharge generally resulted in higher dunes (on average 0.76 m, 0.82 m, 0.78 m, 0.67 and 0.58 m between river kilometer 895-936 for Q = 1353, 1271, 1249, 880 and 781 m 3 s −1 , respectively; Figure 5). No obvious relation between dune height and water depth was found, neither in space nor for varying discharges.

Dune Roughness Prediction
Predicted values of the roughness coefficient (̂ ) followed the same pattern and were in the same order of magnitude as the values of f inferred from the water surface slope ( Figure 6). The general trend shown by all predictors reflected a relatively low roughness between Lobith and Pannerdense Kop (RK 857-867), and an increase in roughness between river km 885-925, followed by a slight decrease until the city of Zaltbommel (RK 935). Between Lobith and the Pannerdense Kop, almost all hydraulic roughness was caused by grain roughness (Figure 6), and grain size seemed to provide an upper limit for the predicted dune roughness values. More downstream, dunes increased in size ( Figure 5) and the difference between total roughness and grain roughness also increased.
The alternative predictors performed differently, and not all variations in roughness inferred from the water surface slopes were captured ( Figure 6). The differences between alternative predictors reached 0.021. At certain locations, the predicted dune roughness by Bartholdy et al. (2010) was 1.6 times as high as the predicted roughness by Van Rijn (1984). The equation by , accounting for the low leeside angle, strongly under-predicts the roughness. Variations in roughness due to changing discharge conditions were strikingly smaller than differences related to the choice of predictor. The coefficient of determination (R 2 ) shows how much of the variance in f is explained from ̂ for various predictors. For example, R 2 for ̂ predicted with Van Rijn (1984) is 0.31, meaning that 31% of the variability of f has been explained. The other predictors exhibit similar trends, yet they all have an R 2 value of less than 0. This means they perform worse than simply taking the average of the roughness inferred from water surface slope data.

Relation Between Roughness and Leeside Angle Statistics
The mean local topographic leeside angle was on average 3.5°, which is slightly higher than the mean dune leeside angle (being 2.3°). The distribution of the low topographic leeside angle was highly positively skewed, meaning low angles dominated and higher angles were less frequently occurring (see Figure S4 in Supporting Information S1). A very low mean local topographic leeside angle was observed until the Pannerdense Kop (RK 867; mean 1.7°; Figure 8). Further downstream, the topographic leeside angle slowly increases until Tiel (RK 913). The trends were similar for both discharge conditions. Other statistical measures of the local topographic leeside angle follow the same pattern.
Mean local topographic leeside angle seems to be unrelated to flow velocity or water depth, and grain size sets an upper limit by limiting sediment mobility (Figure 8b, 8c and 8d). The coefficients of determination (R 2 ) between γ and f were negative, which contradicts our suspicion that mean leeside angles can better explain effective roughness than roughness predictors based on dune height and length (Figure 7).

Influence of Depth Variation
Since dune predictors explain less than one third of the variance of f, it is likely that other roughness imposing elements cause a significant contribution to hydraulic roughness. Figure 9 shows variations in river geometry: the  Figure 9d). To the authors' knowledge, such relation has never been established before. The downstream stretch is relatively straight, with persistent dune fields and relatively small fluctuations in the regional bed level. It is characterized by large dunes (mean Δ = 0.81 m) with a relatively high leeside angle (mean ϕ = 3.0°, mean γ = 4.4°), a comparatively small difference between hydraulic roughness and the roughness predicted based on dune geometry (f -̂ = 0.005), small fluctuations in the smoothed detrended bed level (standard deviation of z = 0.1) and a negative coefficient of correlation between f and dz det /ds (corr = −0.67, see Figure 10). The negative coefficient of correlation indicates an inverse relation between smoothed detrended bed level and roughness. Upstream of river kilometer 893, the out-of-phase relation between dz det /ds and f is lost. This is a highly diverse reach, with coarse sediment, large fluctuations in grain size ( Figure 5), impacts of fixed layers, and strong curvature ( Figure 9). The reach is characterized by low dunes (mean Δ = 0.3 m) with a low leeside angle (mean ϕ = 1.2°, mean γ = 2.4°), a large difference between the hydraulic roughness and roughness predicted based on dune geometry (f -̂ = 0.01), and large fluctuations in the smoothed, detrended bed level (standard deviation of z = 0.2, Figure 10).

Discussion
We will discuss our findings related to the counter-phase fluctuation of bed-elevation gradient and friction factor, the choice of dune roughness predictor, and the impact of various measures of the (local topographic) leeside angle. Since the predictor of Van Rijn (1984) captured the effect of bedforms on mean roughness relatively accu- rately, other sources of roughness hardly contributed to the mean roughness. However, the variance in roughness was not fully explained. Potential sources for the remaining variability in roughness will be discussed, including the presence superimposed dunes, lateral variation in river characteristics, and man-made structures.

Counter-Phase Fluctuation of Bed-Elevation Gradient and Friction Factor
The most remarkable outcome of our analysis focused on the factors controlling hydraulic roughness, is the counter-phase relation between bed-elevation gradient and the friction factor in the downstream reach of the study area (RK 893). The fluctuation in roughness f seems to be imposed by the large scale bed elevation. They occur on a slightly larger spatial scale than alternating bars related to channel curvature, which are indicated by the alternating transverse bed slope in Figure 9c. Large-scale fluctuations of the river depth are found all around the world in studies demonstrating river geometrical variation (Gallo & Vinzon, 2005;Leuven et al., 2021;Trigg et al., 2009;Venditti et al., 2019), and are often related to stratigraphy or river-floodplain interaction.
The observed counter-phase relation between bed-elevation gradient and the friction factor suggests that when depth increases, and the accommodation space for the flow expands, flow energy is lost. We propose an explanation for the observed relation in Figure 11. When the flow runs into a deeper part of the river, the increase in water depth will cause the flow to diverge and decelerate. The resulting energy losses manifest as an increased value of f. This increase in roughness is not directly related to bed roughness from dunes or grains, but rather to internal energy loss where the flow diverges. The normalized values of the gradient in bed level and roughness have a coefficient of determination of 0.34, which means that the inverse gradient in bed elevation explains 34% of the variability in roughness.
Although the counter-phase relation between bed-elevation gradient and roughness can be considered a new finding at a scale of kilometers, it is known to occur at smaller scales. At smaller scales, flow divergence enhances turbulence, which leads to energy loss (Pope, 2000;Tennekes & Lumley, 1972). For example, at the stoss side of a river dune, the flow converges, turbulence is suppressed, and little energy loss takes place. The majority of the energy loss occurs at the lee side of the dune, where the depth increases and flow diverges, enhancing turbulence in the wake of the dune, and dissipating energy. Similarly, in a Venturi meter, the bulk of the energy loss occurs Figure 9. River geometry, curvature and hydraulic roughness, from datasets with comparable discharges (Q = 1271 and 1249 m 3 s −1 ). Shaded gray bars indicate manmade structures in the main channel (fixed layers, bendway weirs and longitudinal training dams at Tiel). (a) bed level, detrended with a third degree polynomial and smoothed with a 8 km LOESS filter, (b) river width, (c) river curvature r (black) and transverse bed slope ξ (green), (d) roughness calculated from the water surface slope (green) and gradient of detrended bed level (z det ; black).
in the expanding section of the flow meter. Hydraulic roughness is closely related to energy loss, it is scale-dependent, and it has an intimate relation with the cascade of turbulent vortices that occur at multiple scales. Our finding suggests additional flow energy dissipates as a result of flow divergence at the scale of kilometers, which is typically viewed as a scale beyond what is relevant for hydraulic roughness.
The bed-elevation gradient depends on the level of data smoothing ( Figure 12). In this study, the choice of filtering was motivated by a spectral gap between bed-level variations by bedforms, including dunes, and regional scale variations imposed by the geological setting, human-made constraints such as groynes, and alternating bars. We expected that part of the variance of f could be explained by each of the latter factors. The counter-phase fluctuation was only visible in the less disturbed middle reach of the study area. The absence of this relation in the more diverse upper reach indicates that the observed relation is not simply a result of bed level being an input parameter for the Chézy equation (Equation 1), or depends on the choice of a smoothing filter.

Roughness Prediction
When simply looking at the explaining variance R 2 , hydraulic roughness from bedforms and friction at the grain scale seemed to have limited predictive capacity. However, the estimation based on the equation by Van Rijn (1984) did predict the mean roughness successfully, as it underestimated the total roughness only by 3.1%.
It is not surprising that the predictor of Van Rijn (1984) performs best, possibly because this predictor has been calibrated using data of river Rhine.
Grain size directly relates to grain roughness. The grain size distributions in the river Rhine and river Waal showed a clear pattern of downstream fining. Up until the Pannerdense Kop (RK 867), average roughness was almost fully explained by grain roughness. In that part of the river, the grain size distribution impacts the potential for dune formation and dune growth (van den Berg & van Gelder, 1993), in such a way that in the Rhine dunes only form under high discharges (Wilbers & Ten Brinke, 2003). Dunes found under lower discharges can be remnants of former high flow conditions. Grain roughness decreased in downstream direction, where the contribution of grain roughness to the total roughness was lower.

The Impact of Leeside Angle on Roughness Prediction
The leeside angles observed in our study domain were low, which is common in large rivers (Cisneros et al., 2020;de Ruijsscher et al., 2020;Galeazzi et al., 2018;Hendershot et al., 2016). The mean leeside angle per kilometer did not exceed 5°. Mean dune leeside angles were slightly lower than local topographical leeside angles. This might indicate that the central axis of the river featured lower leeside angles than the sides, which were also included in the measure of local topographic leeside angle. Alternatively, the Bedform Tracking Tool ( Van der Mark & Blom, 2007) used herein can cause an underestimation in leeside angle, as it applies a linear fit of the middle 2/3 th of the leeside of the dune (Zomer et al., 2021).
The slip face of low-angled dunes, that is, the steepest part of the leeside angle, is considered to be important for flow separation (Paarlberg et al., 2007) and form roughness. Less strong and less frequent flow separation occurs at the leeside of low-angled dunes than at steeper dunes with permanent flow separation. Kornman (1995) found that form roughness is better quantified using the height of the slip face than using the bedform height.  directly included the leeside angle in their predictive equations, introducing a correction factor for low leeside angles. In our study, the predictor by Lefebvre and Winter (2016) (Equation 11) underestimated the total hydraulic roughness, which did not improve when the slip face angle instead of leeside angle was used. This does not necessarily prove the equation wrong, as a significant part of the discrepancy may be explained by alternative causes of roughness.
Against expectations, the various measures of leeside angle, including the mean local topographic leeside angles of low-angled dunes, explained little of the observed roughness variations. Attempts to explain the observed roughness variation from alternative statistics of local topographic leeside angles were unsuccessful ( Figure  S4 in Supporting Information S1). Furthermore, implementing various measures of leeside angle ( Figure S5 in Supporting Information S1) in Equation 11 does not improve roughness prediction. Contrary to the findings of , we found no hard evidence that the detailed shape of dunes matters, beyond properties that are included in classical roughness predictors such as the one from Van Rijn (1984).

Influence of Superimposed Bedforms on Roughness Prediction
Complexity in natural bed morphology at scales smaller than the primary dunes considered in this study may contribute substantially to roughness fluctuations. Secondary dunes are sometimes considered to be a separate roughness element (Julien et al., 2002;Yalin, 1992). In the approach of Julien et al. (2002), roughness heights are added as in k s = k sp + k ss + k sg (roughness height from primary dunes, secondary dunes and grains, respectively), and subsequently the roughness predictor from Van Rijn (1984) is used to quantify hydraulic roughness. Following the approach of Yalin (1992), the equation reads: in which C g represents the grain roughness, and the subscripts p and s denote primary and secondary, respectively. Neither of these approaches improve roughness predictions. Applying the equation of Yalin (1992) to our data resulted in an overestimation of the roughness by more than 100%. This overestimation may be caused the by interaction between primary and secondary dunes, which have a joint effect on the flow dynamics (Fernandez et al., 2006). Alternatively, the secondary dunes could be more important for roughness than the less steep primary dunes. We expect that the secondary dunes are an important cause of hydraulic roughness variation, but new tools are needed to quantify this effect.

Influence of Lateral Variations on Roughness Prediction
Lateral variation in river geometry, dune geometry and grain size, can influence hydraulic roughness parametrization and prediction. Regarding river geometry, the observed alternating bars in the river Waal can exert an indirect control on hydraulic roughness, which extends beyond the impact on dune properties (de Ruijsscher et al., 2020). Flow concentrates in the deepest part of the river, opposite to bars (Zomer et al., 2021). The resulting variation in flow velocity over the cross-section is significant, which means that for a uniform bedform field the effective roughness can vary over the cross-section. Consequently, considering the non-linearity of roughness relations, the 1D approach for quantification of bed geometry may not be entirely representative. Furthermore, width variations as visible in Figure 9b can further complicate the planimetric flow structure, causing the flow to converge or diverge. In our study area, width variations are small and gradual (except for the bifurcation at the Pannerdense Kop (RK 867) where the river Rhine bifurcates in the river Waal and Pannerdensch Kanaal) and are therefore not likely to strongly influence roughness.
Grain size can also show lateral variations. The passage of heavy ships induces erosion of the fine-grained beaches between the groynes, resulting in transport of fine-grained sand from the beaches to the river bed (Wilbers & Ten Brinke, 2003). Currents induced by ships are stronger during low discharge conditions (as present when measurements for this study were taken), and when ships are moving upstream (Bhowmik et al., 1995). Therefore, the southern side of the river Waal will be subject to more fine sand input than the northern side. The implications of variation in grain size over the cross-section on dune geometry remains to be explored. Roughness prediction may be improved by using a width-integrated approach, in which the spatial variability in both grain size and dune geometrical properties are both accounted for. A robust methodology for this is still to be developed.

Influence of Human Impact on Roughness Prediction
Finally, human made structures and dredging can contribute to the total roughness and its variability. For this reason, the fixed layers, bendway weirs, and their region of influence are excluded from the analysis (Section 2). Resistance by groynes contributes to the total hydraulic roughness of the channel. Groynes are known to cause turbulence and scour behind the groyne heads. Locally, the river bed can show undulations with the same wave length as the groyne spacing (Wilbers & Ten Brinke, 2003). Those undulations are fixed and do not interact with bedforms (Ouillon & Dartus, 1997;Wilbers, 1999), but may statically influence bedform geometry (de Ruijsscher et al., 2020). The influence of groynes changes with high discharge conditions, especially when groynes are submerged and become part of the conducting section of the bed (Möws & Koll, 2019;Yossef, 2004). The low flow conditions studied minimize the impact of groynes, but the ubiquitous presence of groynes along the opposite banks complicates quantifying groyne roughness. The expected undulations in the water surface profile are too subtle to be observed based on ship positioning. High resolution measurements of the water surface topography, which cover the full width of the river instead of merely a single track along the center line of river, may further elucidate the causes of hydraulic roughness variation. This will be to the benefit of models simulating flow, sediment transport and bed morphological change.

Conclusions
We quantified hydraulic roughness based on the Darcy-Weisbach friction factor f calculated from water surface slope measurements of a 78 km-long trajectory of the river Rhine and river Waal in the Netherlands. This was compared to predicted roughness values based on dune geometry, and to the spatial distribution of the local topographic leeside angle, both inferred from bathymetric field data.
In the upstream part of the study area , where dunes are likely inactive under low flow conditions and human-made hard layers cause nonuniform flow, total hydraulic roughness cannot easily be predicted from bedform dimensions. In the downstream part (RK 867-935), predictions of the Darcy-Weisbach friction factor f agree with estimates inferred from longitudinal surface level profiles. The best performing predictor explains 31% of the variance in f, which indicates that even in this part of the river, dune morphology is likely not the only factor explaining hydraulic roughness variation. We expected to explain part of the variance in f from statistics of local topographic leeside angle, which control flow separation and energy loss, but did not find confirmation of this expectation.
Serendipitously, we found that the longitudinal gradient of the smoothed river bed level oscillates in counter phase with f. Toward a topographic high, the friction factor decreases, and toward a topographic low, the friction factor increases. A deepening of the river thus corresponds with a higher hydraulic roughness, which may relate to flow divergence in the decelerating flow, and the corresponding energy loss. The depth variations explain 34% of the variance in hydraulic roughness. This effect is clearly visible in the downstream region, where the grain size is relatively constant, dunes are comparatively large, and dune predictors explain a large part of the variance in f.
Unresolved influences on hydraulic roughness include the effect of man-made structures such as groynes, secondary dunes and topographical steering of the river flow. Future work to further elucidate the effect of various scales of complex bed geometry on hydraulic roughness is warranted, as primary dunes explain only a limited part of the flow resistance.

Appendix A: Derivation of St. Venant Equation
The momentum balance of the St.-Venant equations expressed in water level and velocity reads: where u = depth and cross-sectional averaged flow velocity (Q/A, m s −1 ), Q = discharge (m 3 s −1 ) and A = wetted area (m 2 ), S p = pressure slope (= , where d is hydraulic depth), S 0 = bottom slope (= , where z = bed level relative to Amsterdam Ordnance Datum, m), and in which C = Chézy coefficient (m 1/2 s −1 ), and R = hydraulic radius (m). Assuming R ≈ d, which applies to lowland rivers for which W ≫ d, Equation A2 can be simplified as: = 2 2 (A3) For the lowland river subject to study the time derivative term is two orders of magnitude smaller than other terms in the momentum equation, so we set ≈ 0 . After rewriting, Chézy coefficient can be obtained as: Although commonly the simplified form of the Chézy equation is used ( = √ 0 ), this assumption does not hold when there are significant variations in either time or space (non-uniform flow). The bed-level profile (Figure 4a) does not have a constant slope, which causes significant variations in pressure slope.