Feedback Mechanism in Bifurcating River Systems: the Effect on Water-Level Sensitivity

Accurate and reliable estimates of water levels are essential to assess flood risk in river systems. In current practice, uncertainties involved and the sensitivity of water levels to these uncertainties are studied in single-branch rivers, while many rivers in deltas consist of multiple distributaries. In a bifurcating river, a feedback mechanism exists between the downstream water levels and the discharge distribution at the bifurcation. This paper aims to quantify the sensitivity of water levels to main channel roughness in a bifurcating river system. Water levels are modelled for various roughness scenarios under a wide range of discharge conditions using a one-dimensional hydraulic model. The results show that the feedback mechanism reduces the sensitivity of water levels to local changes of roughness in comparison to the single-branch river. However, in the smaller branches of the system, water-level variations induced by the changes in discharge distribution can exceed the water-level variations of the single-branch river. Therefore, water levels throughout the entire system are dominated by the conditions in the largest branch. As the feedback mechanism is important, the river system should be considered as one interconnected system in river maintenance of rivers, flood-risk analyses, and future planning of river engineering works.

. Map of the study area in the Netherlands. The three triangles indicate the representative locations for the water levels in the three downstream branches under consideration in this study. The six circles indicate the locations at which measurements of bedform dimensions are available.

Available Bedform Measurements
Bedform measurements for various periods of time and under different flow conditions are available ( Table 2). The majority of these data are obtained from [25], who have used single-beam and multibeam (mostly since 1997) measurements of the river bed at the Pannerdensche Kop to derive bedform dimensions during large flood events. This data set includes records taken during the discharge wave with the highest recorded discharge in the river Rhine in 1995 as well as from flood waves in 1997 and 1998. Additionally, [26] have analysed the bedforms around the IJsselkop during a moderate flood event in 2004. Finally, [27] provides data under moderate discharges along the Waal branch. Data points from every source at least contain a bedform height and a bedform length estimate as well as the local water depth, depth-averaged flow velocity and discharge estimate. Table 2. Characteristics of the available bedform measurements of [25][26][27]. The measurement locations are shown in Figure 1. For location WA2, measurements were taken for the northern (WA2a) and southern (WA2b) half of the main channel.

Source # Data Points
Location Period Flow Regime (s) [

Available Bedform Measurements
Bedform measurements for various periods of time and under different flow conditions are available ( Table 2). The majority of these data are obtained from [25], who have used single-beam and multibeam (mostly since 1997) measurements of the river bed at the Pannerdensche Kop to derive bedform dimensions during large flood events. This data set includes records taken during the discharge wave with the highest recorded discharge in the river Rhine in 1995 as well as from flood waves in 1997 and 1998. Additionally, [26] have analysed the bedforms around the IJsselkop during a moderate flood event in 2004. Finally, [27] provides data under moderate discharges along the Waal branch. Data points from every source at least contain a bedform height and a bedform length estimate as well as the local water depth, depth-averaged flow velocity and discharge estimate. Table 2. Characteristics of the available bedform measurements of [25][26][27]. The measurement locations are shown in Figure 1. For location WA2, measurements were taken for the northern (WA2a) and southern (WA2b) half of the main channel.

Source # Data Points Location Period
Flow Regime (s) [

Methodology
The sensitivity of water levels to variations in roughness is quantified for a single-branch river and for a bifurcating river system by modelling the water levels under various roughness scenarios ( Figure 2). Therefore, in each of the four Rhine branches, a realistic upper and lower value of the roughness is estimated (Section 3.1). This results in 16 possible combinations of system roughness consisting of a high or a low roughness value for each of the branches. The sensitivity analysis is performed using the one-dimensional hydraulic model of the Rhine branches (Section 3.2). A one-dimensional model allows for fast computations, while it is sufficiently accurate for the purpose of this study. Sensitivity is defined as the difference in water levels between the scenarios in which the local branch experiences a high and a low roughness. The single-branch river and the bifurcating river are differentiated by comparing the water levels at equal local branch discharge or at equal upstream discharge, respectively. Finally, the simulated water levels are linked to return periods to obtain design water levels (Section 3.3). This shows the impact of the water-level variations in the context of flood management.

Methodology
The sensitivity of water levels to variations in roughness is quantified for a single-branch river and for a bifurcating river system by modelling the water levels under various roughness scenarios ( Figure 2). Therefore, in each of the four Rhine branches, a realistic upper and lower value of the roughness is estimated (Section 3.1). This results in 16 possible combinations of system roughness consisting of a high or a low roughness value for each of the branches. The sensitivity analysis is performed using the one-dimensional hydraulic model of the Rhine branches (Section 3.2). A onedimensional model allows for fast computations, while it is sufficiently accurate for the purpose of this study. Sensitivity is defined as the difference in water levels between the scenarios in which the local branch experiences a high and a low roughness. The single-branch river and the bifurcating river are differentiated by comparing the water levels at equal local branch discharge or at equal upstream discharge, respectively. Finally, the simulated water levels are linked to return periods to obtain design water levels (Section 3.3). This shows the impact of the water-level variations in the context of flood management.

Deriving Roughness Scenarios from the Available Bedform Dimension Data
Sixteen roughness scenarios are defined that consist of combinations of an upper and a lower roughness value for each branch. These 16 scenarios represent the variability in main channel roughness in the river system. The high and low roughness value should be seen as an extreme, but realistic spread in roughness that could occur in each branch. It is not intended to give accurate predictions of roughness, given the large spread in the bedform observations, but merely to show the sensitivity of the water levels to variations in roughness. The assumptions that underlie the derivation of these roughness scenarios are discussed in Section 5.1.
The available bedform measurements are transformed into roughness values using empirical roughness predictors. The predictors of Van Rijn [15] and Vanoni-Hwang [28] are selected for this study as those showed the best performance for the Upper Rhine reach [3]. For every data sample, two roughness values are obtained, one for each of the two predictors.
The predictor of Van Rijn [15] is a function of the grain size and the bedform height and length: In addition to the bedform dimensions, the Vanoni-Hwang [28] predictor requires the magnitude of mean water levels and flow velocities:

Deriving Roughness Scenarios from the Available Bedform Dimension Data
Sixteen roughness scenarios are defined that consist of combinations of an upper and a lower roughness value for each branch. These 16 scenarios represent the variability in main channel roughness in the river system. The high and low roughness value should be seen as an extreme, but realistic spread in roughness that could occur in each branch. It is not intended to give accurate predictions of roughness, given the large spread in the bedform observations, but merely to show the sensitivity of the water levels to variations in roughness. The assumptions that underlie the derivation of these roughness scenarios are discussed in Section 5.1.
The available bedform measurements are transformed into roughness values using empirical roughness predictors. The predictors of Van Rijn [15] and Vanoni-Hwang [28] are selected for this study as those showed the best performance for the Upper Rhine reach [3]. For every data sample, two roughness values are obtained, one for each of the two predictors.
The predictor of Van Rijn [15] is a function of the grain size and the bedform height and length: Water 2020, 12, 1915 5 of 17 In addition to the bedform dimensions, the Vanoni-Hwang [28] predictor requires the magnitude of mean water levels and flow velocities: The friction coefficient c f [-] can be transformed into a Nikuradse roughness height as follows: In these predictors, k N is the Nikuradse roughness height (m), D 90 is the 90th percentile grain diameter (m), ∆ and Λ are the mean dune height and mean dune length (m), R is the hydraulic radius (m), h is the water depth (m), u is the depth-averaged velocity (m/s) and υ is the kinematic viscosity (m 2 /s).
The predicted roughness values that correspond to the data of the Waal branch show a large variation ( Figure 3). This scatter is a result of the variation in underlying bedform dimensions and variations in the flow parameters. Generally, the Van Rijn predictor results in higher roughness values than the Vanoni-Hwang predictor. This is coherent with earlier findings of [3]. No distinct discharge-dependency of the roughness is observed in the data at Waal discharges higher than 2000 m 3 /s. This is not consistent with the literature e.g., [29,30] as a larger discharge is generally accompanied by larger bedform dimensions. The predicted roughness values that correspond to the data of the Waal branch show a large variation ( Figure 3). This scatter is a result of the variation in underlying bedform dimensions and variations in the flow parameters. Generally, the Van Rijn predictor results in higher roughness values than the Vanoni-Hwang predictor. This is coherent with earlier findings of [3]. No distinct discharge-dependency of the roughness is observed in the data at Waal discharges higher than 2000 m 3 /s. This is not consistent with the literature e.g., [29,30] as a larger discharge is generally accompanied by larger bedform dimensions. The same analysis is carried out for the three other Rhine river distributaries (not shown here). The predicted roughness values of the Pannerdensch Kanaal (PK) observations show a slightly lower spread in roughness values (kN = 0.05 m up to kN = 0.35 m). However, discharge-dependency is more pronounced with the highest roughness values occurring at the highest observed discharges. This indicates that for higher discharges, higher roughness induced by bedforms may also occur, because of the larger grain size in the PK (Table 1)  The same analysis is carried out for the three other Rhine river distributaries (not shown here). The predicted roughness values of the Pannerdensch Kanaal (PK) observations show a slightly lower spread in roughness values (k N = 0.05 m up to k N = 0.35 m). However, discharge-dependency is more pronounced with the highest roughness values occurring at the highest observed discharges. This indicates that for higher discharges, higher roughness induced by bedforms may also occur, because of the larger grain size in the PK (Table 1) For the other branches, the lower and upper value are set to k N = 0.05 m and k N = 0.55 m as well, for consistency with the Waal branch. For these branches, fewer observations are available. As the bed material and flow conditions are fairly similar in all branches (see Table 1), there is no reason to expect very different roughness values in the other branches. Additionally, the literature indicates that the bedform variability in these branches is also significant [26,27]. Therefore, the same roughness values as for the Waal and Pannerdensch Kanaal are applied.

Sensitivity Analysis Using the Hydraulic SOBEK Model of the Rhine Branches
The SOBEK modelling environment is selected to perform the sensitivity analysis. In this environment, a detailed schematization is available for the study area. This official schematization 'Rijn-j16_5_v1 is used by the Dutch Ministry of Infrastructure and Water Management for operational uses. SOBEK is a modelling environment [31] for 1D hydraulic computations, numerically solving the 1D Saint-Venant equations. The river is schematized by a 1D network in which the other dimensions are assigned by imposing cross-sections. Cross-sections are divided into different parts, e.g., left and right floodplains and the main channel. Furthermore, the presence of minor embankments on the floodplains is accounted for by the addition of flow area in the floodplain area ( Figure 4). The governing equations are solved for the different parts separately and, after iteration, the division of discharge between the different parts of the cross-section is found. This allows for different flow velocities in the floodplains and the main channel.  Table 1), there is no reason to expect very different roughness values in the other branches. Additionally, the literature indicates that the bedform variability in these branches is also significant [26,27]. Therefore, the same roughness values as for the Waal and Pannerdensch Kanaal are applied.

Sensitivity Analysis Using the Hydraulic SOBEK Model of the Rhine Branches
The SOBEK modelling environment is selected to perform the sensitivity analysis. In this environment, a detailed schematization is available for the study area. This official schematization 'Rijn-j16_5_v1′ is used by the Dutch Ministry of Infrastructure and Water Management for operational uses. SOBEK is a modelling environment [31] for 1D hydraulic computations, numerically solving the 1D Saint-Venant equations. The river is schematized by a 1D network in which the other dimensions are assigned by imposing cross-sections. Cross-sections are divided into different parts, e.g., left and right floodplains and the main channel. Furthermore, the presence of minor embankments on the floodplains is accounted for by the addition of flow area in the floodplain area ( Figure 4). The governing equations are solved for the different parts separately and, after iteration, the division of discharge between the different parts of the cross-section is found. This allows for different flow velocities in the floodplains and the main channel. The model schematization has cross-sections with a longitudinal spacing of approximately 500 m. A schematisation of a typical cross-section is shown in Figure 4. In an earlier study [32], the crosssections were extracted from a 2D model, which is based on digital elevation maps constructed from echo-sounds in the main channel and Lidar measurements in the floodplains. After extracting the cross-sections from the 2D model, the cross-sectional dimensions were modified, such that such that modelled stage-discharge relationships match those of the 2D model [32]. Each cross-section is associated with discharge-dependent Chezy roughness values for the floodplains. The floodplains of the Rhine distributaries are vegetated, mostly with grass, but also with higher and denser vegetation. They are mainly used for agriculture and nature conservation and are regularly managed for those purposes. The roughness values are based on vegetation maps, vegetation roughness descriptions and two-dimensional model simulations [3,33]. In an earlier study, the 2D model simulations were used as reference case for the calibration of the floodplain roughness in the 1D SOBEK model [32]. Generally, Chezy values roughness are around 40 m 1/2 /s at lower discharges and 50 m 1/2 /s at higher discharges.
The upstream boundary of the model is located just upstream of Lobith in the Upper Rhine ( Figure 1). The downstream boundary of each branch consists of a stage-discharge relation based on measurements. These boundaries are located more than 60 km downstream of the area of interest, which is sufficiently far away to not influence the water levels in the area of interest (Figure 1). The model schematization has cross-sections with a longitudinal spacing of approximately 500 m. A schematisation of a typical cross-section is shown in Figure 4. In an earlier study [32], the cross-sections were extracted from a 2D model, which is based on digital elevation maps constructed from echo-sounds in the main channel and Lidar measurements in the floodplains. After extracting the cross-sections from the 2D model, the cross-sectional dimensions were modified, such that such that modelled stage-discharge relationships match those of the 2D model [32]. Each cross-section is associated with discharge-dependent Chezy roughness values for the floodplains. The floodplains of the Rhine distributaries are vegetated, mostly with grass, but also with higher and denser vegetation. They are mainly used for agriculture and nature conservation and are regularly managed for those purposes. The roughness values are based on vegetation maps, vegetation roughness descriptions and two-dimensional model simulations [3,33]. In an earlier study, the 2D model simulations were used as reference case for the calibration of the floodplain roughness in the 1D SOBEK model [32]. Generally, Chezy values roughness are around 40 m 1/2 /s at lower discharges and 50 m 1/2 /s at higher discharges. The upstream boundary of the model is located just upstream of Lobith in the Upper Rhine ( Figure 1). The downstream boundary of each branch consists of a stage-discharge relation based on measurements. These boundaries are located more than 60 km downstream of the area of interest, which is sufficiently far away to not influence the water levels in the area of interest ( Figure 1).
The roughness scenarios are implemented in the original schematization of the Rhine branches. In the original schematization, the main channel roughness is calibrated to attain measured water levels for a measured discharge distribution over the branches. The defined main channel roughness scenarios replace these calibrated values of main channel roughness. The roughness of the floodplains are kept at their calibrated values. Because the main channel roughness values are adapted in the scenario runs, they no longer correspond to their calibrated values. This caused shifts in the stage-discharge relationships. To attain values for water levels and bankfull discharges that are in the same order as the values from the calibrated model, a constant additional roughness wad added in all branches. This additional roughness accounts for all roughness inducing processes that are not accounted for in the model schematization, such as the flow resistance caused by groynes. A value of k N,rest = 0.3 m for all branches gives the best match with results from the original, calibrated model schematization (i.e., "Rijn-j16_5_v1"). To find this best match, the water levels averaged over all scenarios were compared with the results from the calibrated model. This value is added to the roughness of the main channel on top of the adopted roughness scenarios. This gives a constant low roughness value of k N = 0.35 m and a constant high roughness value of k N = 0.85 m for all branches. In summary, the floodplain roughness values remain at their calibrated values, while the main channel roughness values under the roughness scenarios are realistic compared to the fully calibrated model by the addition of a constant roughness.
The model is run for every of the 16 roughness scenarios. To assess the sensitivity of the water levels to discharge, the model is run for every scenario with constant discharges ranging from 3000 m 3 /s (above yearly-averaged discharge) up to 18,000 m 3 /s (maximum discharge at Lobith; see [2]) with steps of 500 m 3 /s. This gives a total amount of model runs of 496 (31 discharges and 16 roughness scenarios). The sensitivity to roughness at a certain discharge is defined as the mean difference in water levels between the scenarios in which the local branch experiences a high roughness (H) and a low roughness (L), respectively: Here, the differentiation is made between a single-branch river and a bifurcating river system. In the single-branch river system, the sensitivity to roughness is assessed at equal discharge inside the branch. In contrast, in the bifurcating river system, the sensitivity is assessed at equal upstream Lobith discharge, thereby allowing for variations in the discharge distribution under the various roughness scenarios.

Obtaining Design Water Levels (DWLs) Using Model Averaging
The return periods of modelled water levels are determined for the three locations of interest ( Figure 1) for the single-branch river and for the bifurcating river system. Model averaging (MA) is applied to combine the model results under the various scenarios into a single new model result. MA has been widely used in the field of hydrology [34] and recently MA has been used in river flood modelling as well [35]. Furthermore, it is common practice under the Dutch flood risk framework to apply MA to account for uncertainties, while still having a single model prediction [36]. MA is used to determine the average probability over all 16 scenarios of water level h occurring p(h): Water 2020, 12,1915 8 of 17 Here, p(m i ) is the probability of occurrence of a roughness scenario and p(h m i ) is the probability of water level h occurring under roughness scenario m i . It is assumed that every roughness scenario is equally likely to occur, so: For the sin gle − branch river : p(m i ) = 0.5; for i = 1 : 2 For the bifurcating river : p(m i ) = 0.0625; for i = 1 : 16 (6) The return periods from the modelling tool 'Generator of Rainfall and Discharge Extremes' (GRADE) are used for this analysis [37]. In the Netherlands, this tool is used in operational practices by the Dutch Ministry of Infrastructure and Water Management. Using these data, the return periods of input upstream discharges ranging from 3000 to 18,000 m 3 /s are obtained. For the single-branch river system, branch discharges are coupled to the upstream discharge and their return periods by taking the mean discharge distribution over all 16 scenarios.

Sensitivity to Roughness in a Single-Branch River
The results show that the imposed roughness scenarios induce large changes to the local stage-discharge relationships at the two locations in the Waal and IJssel branches ( Figure 5), representing the sensitivity of single-branch rivers for roughness variations. Averaged over all discharges, the difference between the water levels under the high roughness estimates and the water levels under the low roughness estimate shows water levels of 0.4 m higher at Nijmegen (Waal) and 0.3 m higher at De Steeg (IJssel). Similar sensitivities to the roughness are found for the location of Driel (Nederrijn). These values are in the same order of magnitude as values of water-level uncertainties caused by roughness uncertainty that were previously found in literature, e.g., a 95% confidence interval of 0.53 m at the design discharge in the Waal river [3] and 95% confidence intervals of 0.3-1.0 m in the Po and Garonne river [10]. by the Dutch Ministry of Infrastructure and Water Management. Using these data, the return periods of input upstream discharges ranging from 3000 to 18,000 m 3 /s are obtained. For the single-branch river system, branch discharges are coupled to the upstream discharge and their return periods by taking the mean discharge distribution over all 16 scenarios.

Sensitivity to Roughness in a Single-Branch River
The results show that the imposed roughness scenarios induce large changes to the local stagedischarge relationships at the two locations in the Waal and IJssel branches ( Figure 5), representing the sensitivity of single-branch rivers for roughness variations. Averaged over all discharges, the difference between the water levels under the high roughness estimates and the water levels under the low roughness estimate shows water levels of 0.4 m higher at Nijmegen (Waal) and 0.3 m higher at De Steeg (IJssel). Similar sensitivities to the roughness are found for the location of Driel (Nederrijn). These values are in the same order of magnitude as values of water-level uncertainties caused by roughness uncertainty that were previously found in literature, e.g., a 95% confidence interval of 0.53 m at the design discharge in the Waal river [3] and 95% confidence intervals of 0.3-1.0 m in the Po and Garonne river [10].
The modelled shapes of the stage-discharge relationships are typical for a lowland river with floodplains (see schematized cross-section in Figure 4). Under bankfull stages, marked in Figure 5, the stage-discharge relationship is steep. Once the water level exceeds the height of the minor embankments, the floodplains start accommodating discharge, causing flattening of the stagedischarge relationship. The flattening of the stage-discharge relationship is more pronounced for the IJssel than for the Waal branch, as the floodplains of the IJssel branches are relatively wider than those of the Waal branch (Table 1).  The modelled shapes of the stage-discharge relationships are typical for a lowland river with floodplains (see schematized cross-section in Figure 4). Under bankfull stages, marked in Figure 5, the stage-discharge relationship is steep. Once the water level exceeds the height of the minor embankments, the floodplains start accommodating discharge, causing flattening of the stage-discharge relationship. The flattening of the stage-discharge relationship is more pronounced for the IJssel than for the Waal branch, as the floodplains of the IJssel branches are relatively wider than those of the Waal branch (Table 1).
For single-branch rivers, the sensitivity of water levels to the roughness varies as a function of discharge. The highest sensitivity is around the bankfull stage with a maximum of 0.

Sensitivity to Roughness in a Bifurcating River System
In a bifurcating river system, the model results show that changes in the discharge distribution at the Pannerdensche Kop, indirectly induced by the roughness variations in all branches, have a strong effect on the water levels at Nijmegen ( Figure 6). Averaged over all discharges, the sensitivity of the water levels at Nijmegen to Waal roughness is 0.2 m in this bifurcating river system, while it was 0.4 m in the single-branch river. Minimum and maximum sensitivities are 0.09 m and 0.40 m, which occur at discharges that correspond to the same branch discharges as in the single-branch river, but in this single-branch river the minimum and maximum values were 0.22 m and 0.61 m. At higher discharges, the sensitivity is relatively constant at 0.21 m, while this was 0.45 m in the single-branch river. So, at all discharges, the sensitivity of the water levels to Waal roughness has been reduced, in comparison to the single-branch river, with approximately 0.2 m.
The reduced sensitivity in the bifurcating system compared to the single-branch river is associated with approximately 3% differences in Lobith discharge diverted towards the Waal (right panel of Figure 6). The effect of roughness variations on the ratio of the discharge distribution between the branches is nearly equal for all Lobith discharges. This implicates that absolute discharge variations linearly increase with upstream discharge, e.g., on average a difference of 155 m 3 /s and 440 m 3 /s at upstream discharges of 5000 m 3 /s and 16,000 m 3 /s, respectively.
At De Steeg in the IJssel branch, the difference in response of water levels to roughness variations between the bifurcating river and the single-branch river are even more pronounced. In the bifurcating system, water levels at De Steeg are insensitive to the local (IJssel) roughness (Figure 7). The sensitivity to IJssel roughness varies between 0.03 m (at Q Lob = 18,000 m 3 /s) and 0.16 m (at Q Lob = 5000 m 3 /s), while the sensitivity varied between 0.10 m and 0.58 m in the single-branch river. Still, looking at the overall water-level variations, these exceed the variations in the single-branch river. The difference between the best and worst case roughness scenarios ranges reaches a maximum value of 0.88 m (at Q Lob = 4500 m 3 /s), while it also exceeds 0.5 m at extremely high discharges (>16,000 m 3 /s). These large water-level variations are induced by the deviations in the discharge distribution (right panel of Figure 7).
comparison to the single-branch river, with approximately 0.2 m.
The reduced sensitivity in the bifurcating system compared to the single-branch river is associated with approximately 3% differences in Lobith discharge diverted towards the Waal (right panel of Figure 6). The effect of roughness variations on the ratio of the discharge distribution between the branches is nearly equal for all Lobith discharges. This implicates that absolute discharge variations linearly increase with upstream discharge, e.g., on average a difference of 155 m 3 /s and 440 m 3 /s at upstream discharges of 5000 m 3 /s and 16,000 m 3 /s, respectively.
(a) (b) Figure 6. Modelled water levels at Nijmegen in the Waal branch (a) and the fraction of upstream discharge towards the Waal (b) under the 16 roughness scenarios as a function of (upstream) Lobith discharge. Scenarios with low Waal roughness (blue shaded) cause relatively low water levels and high Waal discharges, while scenarios with high Waal roughness (red shaded) cause relatively high water levels and low Waal discharges.
At De Steeg in the IJssel branch, the difference in response of water levels to roughness variations between the bifurcating river and the single-branch river are even more pronounced. In the bifurcating system, water levels at De Steeg are insensitive to the local (IJssel) roughness (Figure 7). The sensitivity to IJssel roughness varies between 0.03 m (at QLob = 18,000 m 3 /s) and 0.16 m (at QLob = Figure 6. Modelled water levels at Nijmegen in the Waal branch (a) and the fraction of upstream discharge towards the Waal (b) under the 16 roughness scenarios as a function of (upstream) Lobith discharge. Scenarios with low Waal roughness (blue shaded) cause relatively low water levels and high Waal discharges, while scenarios with high Waal roughness (red shaded) cause relatively high water levels and low Waal discharges. The response of water levels to changes in roughness in a bifurcating river system is thus different from a single-branch river ( Figure 8). As seen in Figures 5-7, sensitivity to roughness is higher around the bankfull discharge (5000 m 3 /s) than at larger discharges, at which floodplains accommodate a portion of the discharge. In the bifurcating river, water levels are a function of the roughness of all branches. However, the contribution of the branches to the total water-level variation is not equal, illustrated by the unequal height of the bars. The water levels at every location are most The response of water levels to changes in roughness in a bifurcating river system is thus different from a single-branch river ( Figure 8). As seen in Figures 5-7, sensitivity to roughness is higher around the bankfull discharge (5000 m 3 /s) than at larger discharges, at which floodplains accommodate a portion of the discharge. In the bifurcating river, water levels are a function of the roughness of all branches. However, the contribution of the branches to the total water-level variation is not equal, illustrated by the unequal height of the bars. The water levels at every location are most sensitive to the roughness in the Waal branch, illustrated by the red bars being higher than the other bars of the bifurcating system. Oppositely, the roughness of the smaller branches has little influence on the water levels in the larger Waal branch, as expected. Notably, the maximum sensitivity in the bifurcating river system (total height of the bars) even exceeds the single-branch values (blue bars) for the smaller IJssel and Nederrijn branches. However, this is only reached under specific roughness conditions, e.g., a very high discharge in the IJssel occurs if the Waal has a high roughness, while the Pannerdensch Kanaal has a low roughness. A change in roughness in the Waal branch can almost offset the water-level variations caused by all other branches (equally high bars).
Water 2020, 12, x FOR PEER REVIEW 11 of 17 Figure 8. The sensitivity of water levels to main channel roughness in a single-branch river and in a bifurcating river system (Bif.) at upstream discharges at Lobith (Q) of 5000 and 16,000 m 3 /s. The shown locations are Nijmegen (Waal), De Steeg (IJssel) and Driel (Nederrijn). For every location, the singlebranch sensitivity is shown, as well as the sensitivity to the roughness of every branch in the bifurcating river system. Concluding, the model results showed that the sensitivity of water levels to roughness is very different in a single-branch river than in a bifurcating river system. The sensitivity to local roughness is smaller in a bifurcating river system due to changes in the discharge distribution over the branches. However, the smaller branches of the system experience large water-level variations that are related to the variations in discharge distribution. These discharge effects can exceed the local effects of roughness. These differences in sensitivities between a single-branch river and a bifurcating river system are important for flood-risk assessments. Next, the impact of these differences on design water levels is assessed.

Impact on Design Water Levels
Using model averaging, DWLs are determined at the three locations (Figure 1) for the singlebranch river case and for the bifurcating river case. Due to different degrees of water-level uncertainty, the DWLs are different in a bifurcating river system than in a single-branch river (Table  3).  Concluding, the model results showed that the sensitivity of water levels to roughness is very different in a single-branch river than in a bifurcating river system. The sensitivity to local roughness is smaller in a bifurcating river system due to changes in the discharge distribution over the branches. However, the smaller branches of the system experience large water-level variations that are related to the variations in discharge distribution. These discharge effects can exceed the local effects of roughness. These differences in sensitivities between a single-branch river and a bifurcating river system are important for flood-risk assessments. Next, the impact of these differences on design water levels is assessed.

Impact on Design Water Levels
Using model averaging, DWLs are determined at the three locations (Figure 1) for the single-branch river case and for the bifurcating river case. Due to different degrees of water-level uncertainty, the DWLs are different in a bifurcating river system than in a single-branch river (Table 3). Table 3. Change in design water levels (DWL) at different return periods (T) at the three locations if comparing the bifurcating river system (bif.) to the single-branch river cases. In the Waal branch, design water levels are lower in a bifurcating river system than in the single-branch Waal river (Figure 9). Water-level uncertainties in this branch, are smaller in the bifurcating system than in the single-branch river ( Figure 8). This can also be seen in Figure 9 as all blue dashed lines lie between the dashed red lines. Contrarily, for the IJssel branch, the DWLs are generally higher in the bifurcating river system than in the single-branch IJssel river ( Figure 10). This is caused by the higher water-level uncertainties in the bifurcating river system than in the single-branch river (Figure 8). The small IJssel branch is more sensitive to the discharge variations than the larger Waal branch. In Figure 10, many of the roughness scenarios (blue dashed lines) now fall outside of the single-branch scenarios (red dashed lines). Contrarily, for the IJssel branch, the DWLs are generally higher in the bifurcating river system than in the single-branch IJssel river ( Figure 10). This is caused by the higher water-level uncertainties in the bifurcating river system than in the single-branch river (Figure 8). The small IJssel branch is more sensitive to the discharge variations than the larger Waal branch. In Figure 10, many of the roughness scenarios (blue dashed lines) now fall outside of the single-branch scenarios (red dashed lines).
Contrarily, for the IJssel branch, the DWLs are generally higher in the bifurcating river system than in the single-branch IJssel river ( Figure 10). This is caused by the higher water-level uncertainties in the bifurcating river system than in the single-branch river (Figure 8). The small IJssel branch is more sensitive to the discharge variations than the larger Waal branch. In Figure 10, many of the roughness scenarios (blue dashed lines) now fall outside of the single-branch scenarios (red dashed lines).

Roughness Scenarios
In this study, roughness scenarios are defined on the basis of available data on observed bedform characteristics in the four Dutch Rhine branches. These scenarios represent a realistic estimate of variations in roughness, as roughness predictions from observations and values derived in the literature [26,29] fall within the defined roughness values. Due to limited data availability in the Pannerdensch Kanaal, IJssel and Nederrijn branches, it is assumed that the bedform dynamics in the Waal river are also a good representation of the bedform dynamics in the other branches. Additional measurements of bedform dimensions in mainly the Nederrijn and IJssel branch can show the validity of this assumption and may allow for an improvement of the roughness scenarios. These additional measurements could also show the presence of discharge-dependency, which is now not observed for the Waal branch data ( Figure 3). Generally, bedforms grow with increasing discharge [16,30]. However, the relationship of form roughness with discharge is even more complex. Amongst other factors, the shape of the bedform and its lee-side angle [38] as well as the water depth [26] play an important role. While in the upper reach of the Dutch Rhine discharge-dependency is clearly present [3,29], this was not as pronounced not found for downstream locations in the Waal branch [26,29].
Consequently, a higher than assumed spread of main channel roughness may occur at high discharges. Firstly, if the main channel roughness is more discharge dependent, higher roughness values at extremely high, unobserved discharges may be expected. However, secondly, flattening of bedforms in a condition of Upper Stage Plane Bed could also occur locally [19,39], leading to a significant decrease in roughness [18]. An increased spread in roughness at higher discharges is, therefore, possible. If the spread in roughness is larger than assumed, the bandwidth of modelled water levels will become somewhat larger. However, the difference in water-level response to roughness variations between a single-branch river case and a bifurcating river can still be expected. Likely, this difference grows with a larger spread in roughness, which would also lead to larger differences in design water levels.
Another assumption in the roughness scenarios was that the main channel roughness between the branches could change independently of the other branches. However, correlation of the roughness between the branches may be expected. The primary source of possible correlation between the roughness of the branches is the history of the upstream discharge by which the bedforms have developed [25,38]. Other factors that can cause correlation between the roughness of the branches can be grain characteristics [25,26], similar cross-sectional shape [25] and the shape of the discharge wave [40]. However, differences in these factors are present in the Rhine system. For example, the Waal branch has much finer sediment than the IJssel branch and also in the Waal branch a smaller portion of the discharge is accommodated by the floodplains than in the IJssel branch. Strong correlation is, thus, not apparent.
To test the effect of correlation of roughness between the branches, results are shown for the case if full correlation of the roughness is assumed (Figure 11). This leaves two roughness scenarios: either a high roughness in every branch or a low roughness in every branch. Sensitivities of water levels to roughness would generally decrease compared to the uncorrelated case. In the bifurcating river system, correlation between the roughness of the branches means that relative roughness differences between the branches diminish. Therefore, the discharge distribution will show less variation. Due to less variation in the discharge distribution, sensitivity is much closer to the single-branch river ( Figure 11). Furthermore, correlation causes the maximum water-level variation to decrease in the smaller branches of the system, which is dominated by the variations in the discharge distribution.
Water 2020, 12, x FOR PEER REVIEW 14 of 17 Figure 11. Sensitivities of water levels in a bifurcating river system if the main channel roughness were fully uncorrelated or fully correlated and in a single-branch river. The bar for the uncorrelated case is the sum of the individual components for each branch (total height in Figure 8).

Uncertainty Analyses in Bifurcating River Systems
In this study, the sensitivities of water levels to roughness have been quantified using a onedimensional SOBEK model of the Rhine branches. The 1D model results are expected to be reasonably accurate in support of the aim of this study. A two-dimensional model is preferred over a onedimensional model if more accurate predictions of absolute water levels are required [20]. However, in this study, the focus is on differences in water levels between the scenarios. The 1D approach is expected to be able to predict these differences accurately, as the various roughness scenarios are treated similarly in the model schematization. A comparison with an earlier study of water-level sensitivity in the Waal branch [3] showed similar values of water-level sensitivity as were found in this study. This paper has considered variability in discharge and uncertainty in main channel roughness, being the most important sources of uncertainty. Besides these sources, several other sources of uncertainty, such as floodplain roughness [33] and geometry [12], result in uncertain water-level estimates. If taking these sources of uncertainty into account, slightly larger water-level variability can be expected. Specifically at extreme discharges, uncertainty related to floodplain roughness and geometry will influence the water-level estimates. It can be expected that these sources of uncertainty show less correlation between the branches, as within a branch little correlation is observed e.g., [33] thereby inducing larger changes to the discharge distribution. Water levels modelled in a singlebranch river will, therefore, increasingly be different from the water levels modelled in a bifurcating river system. This will especially be true for sources of uncertainty that may occur very locally and in the vicinity of the bifurcating point, e.g., breaching of a flood defense [2], Upper Stage Plane Bed Figure 11. Sensitivities of water levels in a bifurcating river system if the main channel roughness were fully uncorrelated or fully correlated and in a single-branch river. The bar for the uncorrelated case is the sum of the individual components for each branch (total height in Figure 8).

Uncertainty Analyses in Bifurcating River Systems
In this study, the sensitivities of water levels to roughness have been quantified using a one-dimensional SOBEK model of the Rhine branches. The 1D model results are expected to be reasonably accurate in support of the aim of this study. A two-dimensional model is preferred over a one-dimensional model if more accurate predictions of absolute water levels are required [20]. However, in this study, the focus is on differences in water levels between the scenarios. The 1D approach is expected to be able to predict these differences accurately, as the various roughness scenarios are treated similarly in the model schematization. A comparison with an earlier study of water-level sensitivity in the Waal branch [3] showed similar values of water-level sensitivity as were found in this study.
This paper has considered variability in discharge and uncertainty in main channel roughness, being the most important sources of uncertainty. Besides these sources, several other sources of uncertainty, such as floodplain roughness [33] and geometry [12], result in uncertain water-level estimates. If taking these sources of uncertainty into account, slightly larger water-level variability can be expected. Specifically at extreme discharges, uncertainty related to floodplain roughness and geometry will influence the water-level estimates. It can be expected that these sources of uncertainty show less correlation between the branches, as within a branch little correlation is observed e.g., [33] thereby inducing larger changes to the discharge distribution. Water levels modelled in a single-branch river will, therefore, increasingly be different from the water levels modelled in a bifurcating river system. This will especially be true for sources of uncertainty that may occur very locally and in the vicinity of the bifurcating point, e.g., breaching of a flood defense [2], Upper Stage Plane Bed [19] or erosion of the river bed [41]. Accounting for the feedback mechanism by modelling the full bifurcating river system is thus essential in uncertainty analyses of water levels for flood risk management.

Conclusions
In this paper the sensitivity of water levels to main channel roughness is quantified for the bifurcating Dutch river Rhine system. In earlier studies, the effect of roughness uncertainty was evaluated for single-branch rivers. This study shows that the feedback mechanism between downstream water levels and discharge distribution at a river bifurcation strongly affects the sensitivity of water levels to roughness variations.
The results show that in a bifurcating river system water levels are less sensitive to changes of the local roughness than in a single-branch river. Changes in the discharge distribution counteract the effects of roughness on water levels. Comparing the bifurcating river with the single-branch rivers, sensitivity to local roughness at a very high upstream discharge of 16,000 m 3 /s reduced from 0.45 m to 0.22 m, from 0.17 m to 0.04 m and from 0.23 m to 0.07 m, in the Waal, IJssel and Nederrijn branches, respectively. In addition to sensitivity to local roughness, water levels in the bifurcating river system are also sensitive to the roughness of the other branches. Variations in roughness in the largest Waal branch lead to large variations in water levels throughout the entire system. In contrast, variations in roughness in the smallest branch lead to little change in water levels locally, and elsewhere in the system. In the smaller branches, the discharge variations that are governed by the larger branches, can exceed the sensitivity to roughness, thereby increasing the total variability of water levels in a bifurcating river system compared to a single-branch river. At an upstream discharge of 16,000 m 3 /s, the total variability of water levels in the IJssel and Nederrijn branches is 0.45 m and 0.38 m respectively, clearly exceeding the single-branch values of 0.17 m and 0.23 m. These principles are generally valid for lowland bifurcating river systems. Therefore, it implies that in these systems the water levels are dominated by water-level uncertainties in the largest downstream branch, if those uncertainties occur in the vicinity of the bifurcation point.
Design water levels are different in a bifurcating river system than in single-branch rivers. In the largest branch, design water levels are lower due to the lower water-level uncertainties, while in the smallest branch design, water levels increase. Therefore, it is essential that, in river management, the river system should be modelled as one interconnected system in the analysis of the water levels along the branches. This is important for the maintenance of rivers, the assessment of flood risks as well as for the future planning of river engineering works.