Model-Based Analysis of Macrophytes Role in the Flow Distribution in the Anastomosing River System

: The impact of vegetation on the hydrology and geomorphology of aquatic ecosystems has been studied intensively in recent years. Numerous hydraulic models developed to date help to understand and quantitatively assess the inﬂuence of in-stream macrophytes on a channel’s hydraulic conditions. However, special focus is placed on single-thread rivers, leaving anastomosing rivers practically uninvestigated. To ﬁll this gap, the objective of this study was to investigate the impact of vegetation on ﬂow distribution in a complex anastomosing river system situated in northeastern Poland. The newly designed, one-dimensional, steady-ﬂow model, dedicated for anastomosing rivers used in this study indicated high inﬂuence of vegetation on water ﬂow distribution during the whole year in general, but—as expected—signiﬁcantly higher in the summer season. Simulations of in-stream vegetation removal in selected channels reﬂected in Manning’s coefﬁcient alterations caused relatively high discharge transitions during the growing season. This proved the signiﬁcance of feedback between process of plants growth and distribution of ﬂow in anabranches. The results are unique and relevant and could be successfully considered for the protection of semi-natural anabranching rivers.


Introduction
In recent years, there has been an increasing recognition of the role of vegetation in hydrology and fluvial geomorphology, supported by numerous studies documenting its influence on river channel morphological processes and geomorphic forms [1][2][3]. Riparian and aquatic vegetation, depending on their structure, density, coverage, etc. can affect river morphology in different ways, by altering the flow and subsequently the sediment balance through hydraulic resistance [4][5][6]. Numerous studies have described the empirical relationships between biological and physical characteristics in aquatic ecosystems [7,8]; nonetheless, more quantitative understanding, linking current knowledge with predictive models is required [9]. One of the examples of model-based assessment of vegetation impact on river pattern and morphodynamics was presented by Oorschot et al. [10]. In their study of a single-channel meandering river, they used a 2D sediment routing model coupled with a dynamic vegetation model. The study results indicated that models with dynamic vegetation, compared to static, give a different morphological response of a meandering channel. The first approach indicated meandering maintenance and the second one a reduction of morphodynamics over time and transformation of the river system into being vegetation-dominated. It stresses the fact that vegetation affects but also is affected by morphological processes.
The complexity of the biological and physical links escalates for multi-channel rivers characterized by intricate flow conditions and hydraulics. Most previous research has used numerical models in multichannel rivers that focused on braided rivers exploring flow-vegetation interactions [11][12][13] as well as sediment transport as a control on channel pattern [14,15]. A broad review evaluating the modeling frameworks appropriate for simulating the morphodynamics of braided rivers was given by Williams et al. [16]. The review covers a whole range of models varying from 1D to 3D and from reach to catchment scale. According to the authors, a 2D physics-based approach offers the highest potential for simulating braided river morphodynamics. In the study by Jowett and Duncan [17], a comparison of 1D and 2D model predictions of water depth and velocity is presented. The results indicate that the correlation between predicted and measured depths and velocities was higher for 1D than for 2D models, but 2D models gave more accurate predictions outside the calibration range of 1D models. Therefore, either 1D or 2D models are suitable for braided rivers modeling, and if done well, only small differences between the results should be obtained [17].
Anastomosing rivers constitute different type of multi-channel river systems. Their stable channels, which are separated by vegetated (floodplain) islands, occur in a variety of climatic conditions: subarctic, temperate, tropical and arid regions [18,19]. They are reported to form via avulsion, in which new channels (anabranches) form when water breaks through erosion-resistant banks and begins to incise gradually into the floodplain [20][21][22][23][24]. Numerous hypotheses explain avulsion triggering, which is indispensable for the creation and maintenance of anabranches. Among them the most widely applicable are a highly variable flood-prone flow regime [19,23] and localized flow vegetation obstacles [20,22]. The role of vegetation in influencing anastomosing river hydrology and geomorphology is complex, and-as stated by Tooth and Nanson [25]-needs to be assessed more fully, based on the data from a variety of rivers, before any generalizations can be made. The most comprehensive studies describing the impact of vegetation on anastomosing system formation concerned the River Okavango (Botswana) [20,22,[26][27][28]. The authors, in their studies, created a conceptual model of the formation of the anastomosing system, which was assumed to be mainly controlled by the growth of bank and in-channel vegetation. Moreover, a theoretical description of the interaction of vegetation-channel geometry was provided. According to the authors, vegetation-induced flow velocity decrease and sediment aggradation initiates feedback mechanisms involving further invasion of aquatic plants in the channel, reducing flow rate and causing elevation of the channel water surface. This in turn increases the rate of water loss from the channel, accelerating blockage and aggradation. Similar mechanisms were described by Wende and Nanson [29], who examined the mechanisms of anabranch formation in Western Australia and emphasized the crucial role of vegetation.
However, besides qualitative, more quantitative, process-based studies recognizing the role of vegetation in influencing anabranching river hydrology and geomorphology are required. The hydraulic geometry of anabranching rivers remained unstudied until Tabata and Hickin [30] investigated a 10-km long reach of the upper Columbia River (Canada). They concluded in their study concerning the hydraulic efficiency of the anastomosing Columbia River that there is a great need for additional studies investigating inter-channel hydraulic properties and flow conditions in anastomosing rivers. Moreover, the little data available for anabranching rivers hints at formative processes not examined by hydraulic geometry and vegetation [31]. One of the examples of anastomosing river modeling was presented by Van et al. [32] who developed a 1D, steady-state HEC-RAS model for the anastomosing section of the Mekong river. However, they declared no field measurements on the flow distribution and river cross-sections, and limited their research only to one bifurcation without taking into account the detailed role of vegetation in flow distribution. Another attempt of anastomosing river modeling was made by Gibson and Pasternack [33], who developed a 1D and 2D model for three morphologically different reaches among which, anastomosing could be found. Although only one bifurcation was modeled, the authors found 1D model simulations underwhelming and split flow modeling to be an intermediate level of complexity. Therefore, and additional investigation improving split flow modeling is required.
Hydraulic modeling using 1D [34][35][36] and 2D [10,17,37] models is widely used to describe water flow in rivers. Modeling of hydraulic processes, flow and water-surface elevations in a fluvial system is influenced by vegetation contributing to roughness or overall flow resistance. Total roughness can be either partitioned to determine the relative contributions of vegetation type, distribution and morphology [38], or combined into a single parameter including the impact of vegetation density and morphology [39]. In hydraulic models, vegetation induced roughness is very often expressed in terms of the Manning's coefficient through the Manning's equation for uniform flow [40][41][42]. Although it was stated that Manning's equation might be inappropriate in cases where the drag force is dominant over the bed friction [43], at the reach scale, the results are reliable [44,45]. Therefore, it is a widely accepted methodology in stream resistance evaluation, especially in vegetated areas [46] and was successfully incorporated in numerous studies investigating impact of vegetation on water flow, depth and velocity [47,48]. The importance of Manning's coefficient assignment in the accuracy of hydraulic models was raised by Yang et al. [49] who found it highly difficult as roughness depends on flow circumstances, a stream's geomorphology and physical conditions. To reduce the uncertainties to the greatest possible extent, field measurements can be conducted following Frias et al. [50], who recommended accurate and dense measurements of the river bed using GPS instead of incorporating Digital Elevation Models. The Manning's coefficient is varied due to the changing discharge and diverse vegetational influence over the year [51][52][53]. There are numerous methods estimating flow resistance based on Manning's coefficient that describe open-channel flow in a simplified but accurate way [54]. Most models focus on the calibration of the roughness parameter which, together with the geometry, is considered to have significant impact on flow modeling [55].
A model's sensitivity analysis can be used to evaluate the magnitude of the dependency of model computations on certain modeling assumptions, e.g., parameter values. In other words, it is the study of how the different input variations of a model influence the variability of its output. It can be divided into two broad categories: local and global sensitivity analysis [56]. The local sensitivity is the sensitivity of the model output when the impact of small input parameter perturbations near its nominal values is considered. Global sensitivity determines the effect of the input parameters' variation (due to parameter uncertainty) on model output [57,58]. In this study, focused on investigation of vegetation influence on water distribution in anastomosing river system, local sensitivity analysis is evaluated.
Studies on the modeling of multi-channel river systems mainly use 2D models, solving the momentum equation and mass conservation in distributing flow between side channels, which does not require detailed discharge measurements [11][12][13]17]. In cases where 1D models were used and the energy conservation equation was solved, the modeling was limited to only one bifurcation, as for proper model validation discharge measurements in all anabranches are required [32]. In experiments investigating vegetation-flow interactions and more complex geomorphological responses [10][11][12][13], upstream boundary conditions of flow were fixed either as a hydrograph or a constant flow. For single-channel rivers, it is appropriate and it is usually derived from discharge measurements. In multi-channel rivers, flow is distributed between side-channels and the distribution is highly related to flow obstacles in each channel and the channel geometry. A decrease of flow in one anabranch causes an increase in another and the effect propagates to downstream and upstream channels. This means that upstream boundary conditions of river flow in each anabranch remain unknown. Therefore, to be able to model complex interactions between biological and physical characteristics in anastomosing rivers, a reliable model, reflecting accurate flow distribution, is required. In other words, we cannot investigate complex feedbacks in an anastomosing system unless we take "one step back" and make sure that the flow rate in each anabranch is accurately simulated. 2D models are usually evaluated against water depths and velocities, but a remaining hindrance to wide acceptance of 2D modeling is the lack of standards for 2D-model evaluation [59]. Their structure is based on Digital Elevation Model (DEM) [40]. Light Detection and Ranging (LIDAR) derived DEMs are suitable for reflecting the geometry of a floodplain, but, since red laser does not penetrate through water table, it does not represent channel's geometry accurately [60]. That the river reach investigated in this study consists of interconnected river channels with average width equal to 3 m convinced us that the use of 1D model supported by detailed field measurements of discharge distribution and channel's geometry is the most suitable.
Against this background, the objective of this study was to investigate the impact of vegetation on flow distribution in a complex anastomosing river system situated in northeastern Poland. To this end, a newly designed 1D, steady-flow model, dedicated for anastomosing rivers was used. By means of mathematical simulations, model sensitivity analysis on vegetation changes was conducted to reflect its impact on flow distribution between anabranches.

Site Description
The Upper River Narew is a lowland, low-gradient river situated in northeastern Poland ( Figure 1). The river valley, due to its precious environmental nature, has also been designated a Wetland of International Importance under the terms of the RAMSAR Convention and is protected as a part of Natura 2000 (Bird and Habitat directive) and the Narew National Park (NNP). Within the NNP borders, the river consists of a complex network of small interconnected, unconfined channels draining peat substrate. The riparian corridor vegetation structure is dominated by reed communities (Phragmition) occupying 35% and sedge communities (Magnocaricion) also covering 35% of the valley. Mature trees and scrubs (Alnetea glutinosae) cover 15% of the NNP area mostly at the edges of the valley and the remaining 15% is occupied by grass (Molinio-Arrhenatheretea) ( Figure 1). The uppermost reach of the Narew catchment was dammed and a reservoir with total capacity of 79.5 million m 3 was constructed ( Figure 1). The reservoir, which has operated since 1992, caused changes in extreme flow events, increasing low and decreasing high flows [61].
Water 2018, 10, x FOR PEER REVIEW 4 of 16 equal to 3 m convinced us that the use of 1D model supported by detailed field measurements of discharge distribution and channel's geometry is the most suitable. Against this background, the objective of this study was to investigate the impact of vegetation on flow distribution in a complex anastomosing river system situated in northeastern Poland. To this end, a newly designed 1D, steady-flow model, dedicated for anastomosing rivers was used. By means of mathematical simulations, model sensitivity analysis on vegetation changes was conducted to reflect its impact on flow distribution between anabranches.

Site Description
The Upper River Narew is a lowland, low-gradient river situated in northeastern Poland ( Figure  1). The river valley, due to its precious environmental nature, has also been designated a Wetland of International Importance under the terms of the RAMSAR Convention and is protected as a part of Natura 2000 (Bird and Habitat directive) and the Narew National Park (NNP). Within the NNP borders, the river consists of a complex network of small interconnected, unconfined channels draining peat substrate. The riparian corridor vegetation structure is dominated by reed communities (Phragmition) occupying 35% and sedge communities (Magnocaricion) also covering 35% of the valley. Mature trees and scrubs (Alnetea glutinosae) cover 15% of the NNP area mostly at the edges of the valley and the remaining 15% is occupied by grass (Molinio-Arrhenatheretea) ( Figure 1). The uppermost reach of the Narew catchment was dammed and a reservoir with total capacity of 79.5 million m 3 was constructed ( Figure 1). The reservoir, which has operated since 1992, caused changes in extreme flow events, increasing low and decreasing high flows [61].
The study focused on a 4-km long stretch of the anastomosing river section in the NNP ( Figure  2). The analyzed reach is characterized by low gradient (0.0002 m·m −1 ) and mean annual discharge (1951-2016) recorded at the Suraż gauging, located upstream of the analyzed section, equal to 15.5 m 3 ·s −1 . Average width/depth ratio is very low and equals 2. According to long-term records of sediment concentrations monitored at the Suraż water quality monitoring point, average concentration of sediment equals 19 mg·L −1 .  The study focused on a 4-km long stretch of the anastomosing river section in the NNP ( Figure 2). The analyzed reach is characterized by low gradient (0.0002 m·m −1 ) and mean annual discharge  recorded at the Suraż gauging, located upstream of the analyzed section, equal to 15.5 m 3 ·s −1 . Average width/depth ratio is very low and equals 2. According to long-term records of sediment concentrations monitored at the Suraż water quality monitoring point, average concentration of sediment equals 19 mg·L −1 . In the NNP, the abundance of aquatic vegetation in the river channels strongly depends on the water depth and current velocity. Active channels deeper than 2 m are usually free of vegetation (Figures 3C,D and 4D) but shallower reaches are extensively colonized by both emergent and submerged macrophytes ( Figures 3A,B and 4A-C). It is worth stating that in-stream vegetation in the analyzed area densely overgrows channels not only during the summer season, but also survives to some extent the dormancy season, influencing the hydraulics of channels for the whole year.
Studies show that in recent decades uncontrolled expansion of common reed (Phragmites australis) became a serious problem causing significant loss of biodiversity [62] and the loss of anabranches [61,63,64]. Therefore, a special focus on a detailed recognition of flow conditions, improving the mechanistic understanding of anastomosing in lowland river floodplains and facilitating the development of sustainable conservation strategies for these important habitats is required. In the NNP, the abundance of aquatic vegetation in the river channels strongly depends on the water depth and current velocity. Active channels deeper than 2 m are usually free of vegetation ( Figures 3C,D and 4D) but shallower reaches are extensively colonized by both emergent and submerged macrophytes ( Figures 3A,B and 4A-C). It is worth stating that in-stream vegetation in the analyzed area densely overgrows channels not only during the summer season, but also survives to some extent the dormancy season, influencing the hydraulics of channels for the whole year.
Studies show that in recent decades uncontrolled expansion of common reed (Phragmites australis) became a serious problem causing significant loss of biodiversity [62] and the loss of anabranches [61,63,64]. Therefore, a special focus on a detailed recognition of flow conditions, improving the mechanistic understanding of anastomosing in lowland river floodplains and facilitating the development of sustainable conservation strategies for these important habitats is required.

Fieldwork and Modeling Tool
Broad field measurements conducted for the purpose of this research included discharge measurements at all 14 confluences and channel geometry as well as water level measurements using GPS (76 cross-sections). Cross-sections of all river reaches were measured at every inlet and outlet part and additionally every 500 m. Discharge measurements were conducted in two different periods: (1) during bankfull water stages in the dormancy season (16 m 3 ·s −1 ); and (2) during low water stages in growing season (7 m 3 ·s −1 ). This ensured the recognition of different in-stream vegetation conditions varying between the seasons. In both cases, only river flows below bankfull stages were considered.
In this study, a one-dimensional steady flow model developed in the MATLAB computing environment, hereafter referred to as Multi-Channel Flow Routing Model (MFRM), was used. The

Fieldwork and Modeling Tool
Broad field measurements conducted for the purpose of this research included discharge measurements at all 14 confluences and channel geometry as well as water level measurements using GPS (76 cross-sections). Cross-sections of all river reaches were measured at every inlet and outlet part and additionally every 500 m. Discharge measurements were conducted in two different periods: (1) during bankfull water stages in the dormancy season (16 m 3 ·s −1 ); and (2) during low water stages in growing season (7 m 3 ·s −1 ). This ensured the recognition of different in-stream vegetation conditions varying between the seasons. In both cases, only river flows below bankfull stages were considered.
In this study, a one-dimensional steady flow model developed in the MATLAB computing environment, hereafter referred to as Multi-Channel Flow Routing Model (MFRM), was used. The

Fieldwork and Modeling Tool
Broad field measurements conducted for the purpose of this research included discharge measurements at all 14 confluences and channel geometry as well as water level measurements using GPS (76 cross-sections). Cross-sections of all river reaches were measured at every inlet and outlet part and additionally every 500 m. Discharge measurements were conducted in two different periods: (1) during bankfull water stages in the dormancy season (16 m 3 ·s −1 ); and (2) during low water stages in growing season (7 m 3 ·s −1 ). This ensured the recognition of different in-stream vegetation conditions varying between the seasons. In both cases, only river flows below bankfull stages were considered.
In this study, a one-dimensional steady flow model developed in the MATLAB computing environment, hereafter referred to as Multi-Channel Flow Routing Model (MFRM), was used. The rationale for developing a new model instead of using existing ones was that none of the available ones were equipped with automatic optimization procedure of the roughness parameters. Whilst in single channel systems this procedure is relatively simple and manually feasible, in complex interconnected system of channels it becomes difficult and laborious. While the full description of model setup, calibration and validation was presented in the study of Marcinkowski et al. [65], here we provide a brief overview, important in the context of the main goal of this paper.
MFRM operates under the following assumptions: (1) steady flow conditions; (2) one-dimensional flow is considered, where only the velocity components in the flow direction are taken into account, and the elevation of the energy line is the same in the whole cross-section; (3) river flow can be expressed in terms of energy conservation equation; (4) discharge within each river branch is uniform; and (5) flow is subcritical. Manning's roughness coefficients were used to describe the flow resistance resulting from vegetation presence in the channel. It was assumed that each channel can be characterized with a single value of the coefficient. Floodplain roughness was not considered, as only channel flows were investigated (below bankfull stages). The stationarity of the model required elaborating two different sets of Manning's coefficients: (1) for highly vegetated growing season; and (2) for less intensively vegetated dormancy season.
The model was calibrated and validated to minimize the difference between calculated and observed water levels at every channel junction of the analyzed river system, ensuring that computed water energy heads at parallel anabranches are equal (at channel splitting junctions). Manning's coefficients were determined using the automatic optimization procedure minimizing the difference in water energy levels between branches and measurements.

Model Performance
Manning's coefficients in MFRM for all river branches ranged from 0.027 m −1/3 ·s to 0.054 m −1/3 ·s (0.040 m −1/3 ·s on average) and from 0.027 m −1/3 ·s to 0.089 m −1/3 ·s (0.061 m −1/3 ·s on average) for dormancy and growing season, respectively. It appears to be reasonable comparing those values with the literature [66][67][68]. The overall assessment of model accuracy can be expressed by average differences in: (1) simulated water levels, which for bot, growing and dormancy seasons did not exceed 6.6 cm comparing to measured water levels; and (2) energy level between branches which reached accuracy at 10 −5 m level (Table 1). Figure 5A,B presents the spatial variability of MFRM Manning's coefficients for two separate vegetation conditions. As expected, on the left-hand side of the valley, Manning's coefficients are greater, which is caused by channels' shallow geometry and the vegetation blockage rate, i.e., height of vegetation (H) either equal to water depth (h) (h/H = 1) (cf. Figure 3A,B) or higher than water depth (h/H < 1) (cf. Figure 3A-C), overgrowing the whole channel width. For the main channel (cf. Figures 3C,D and 4D), on the right-hand side, this value is significantly lower, as the blockage rate is significantly lower. Although the Manning's coefficients differ between seasons (greater values in growing season), similar spatial variation is maintained in both seasons. Water flow distribution determined by vegetation resistance significantly differs between seasons. During the winter season, when mean water flow equals 16 m 3 ·s −1 nearly 20% (3.1 m 3 ·s −1 ) is contributed into side anabranches ( Figure 5C), whereas, in the summer season, with mean flow 7 m 3 ·s −1 , the contribution rate decreases to only 2% (0.17 m 3 ·s −1 ) ( Figure 5D). Table 1. Differences in measured vs. modeled water levels in MFRM for calibration and validation.  Figure 5A,B presents the spatial variability of MFRM Manning's coefficients for two separate vegetation conditions. As expected, on the left-hand side of the valley, Manning's coefficients are greater, which is caused by channels' shallow geometry and the vegetation blockage rate, i.e., height of vegetation (H) either equal to water depth (h) (h/H = 1) (cf. Figure 3A,B) or higher than water depth (h/H < 1) (cf. Figure 3A-C), overgrowing the whole channel width. For the main channel (cf. Figures  3C,D and 4D), on the right-hand side, this value is significantly lower, as the blockage rate is significantly lower. Although the Manning's coefficients differ between seasons (greater values in growing season), similar spatial variation is maintained in both seasons. Water flow distribution determined by vegetation resistance significantly differs between seasons. During the winter season, when mean water flow equals 16 m 3 ·s −1 nearly 20% (3.1 m 3 ·s −1 ) is contributed into side anabranches ( Figure 5C), whereas, in the summer season, with mean flow 7 m 3 ·s −1 , the contribution rate decreases to only 2% (0.17 m 3 ·s −1 ) ( Figure 5D).

Model Sensitivity Analysis
The sensitivity analysis of MFRM on vegetation changes was calculated based on the central difference formula (Equation (1)) [56].

Model Sensitivity Analysis
The sensitivity analysis of MFRM on vegetation changes was calculated based on the central difference formula (Equation (1)) [56].
where θ j stands for the Manning coefficient for j-th channel branch, ∆θ j is the change of the Manning's coefficient value, and Q is the discharge for all river branches as model output. It was assessed iteratively by altering the Manning's coefficient (by 0.005 m −1/3 ·s) in one branch at a time and launching the model computations. The relative coefficient change was tested for both growing season and dormancy season. Finally, sensitivity results were presented in relation to maximum discharge change, and hence, they ranged between 0 and 1.

Vegetation Changes Impact on Water Flow
To analyze the impact of in-channel vegetation changes on water flow distribution, plants, along with the degree of their abundance in the channel, were grouped based on their physical parameters (height and width), which determine the blockage factor (the fraction of the channel cross-section blocked by vegetation) [69]. As proven by Luhar and Nepf [70], the Manning roughness coefficient due to vegetation depends primarily on the blockage factor. Nepf [66] in her studies showed also that resistance to flow in channels with vegetation is influenced primarily by plants' submergence. The empirical nature of the Manning's equation has often been criticized, but it remains in wide use [40,42,49,[51][52][53]. Even though the most fundamental treatments link resistance coefficients to channel boundary roughness, representing many sources of flow resistance (i.e., vegetation) by such coefficients has been practiced for decades [67]. To each group of vegetation, Manning's coefficient value was assigned after Arcement and Schneider [68] and Rhee et al. [71] (Table 2). Afterwards, for the reaches with highest flow resistance (>0.05 m −1/3 ·s) (cf. Figure 5B), two different variants reflecting potential man-made interventions (i.e., mowing) in channel vegetation structure were tested in MFRM. After changing uniformly the Manning's coefficient values in aforementioned reaches to 0.05 m −1/3 ·s for summer season (which means absence of emergent vegetation) and 0.025 m −1/3 ·s for both seasons (which means absence of submerged vegetation along the channel bottom), discharge distribution was recalculated and compared with the calibrated variant.

Model Sensitivity Analysis
The results of MFRM sensitivity analysis induced by Manning's coefficient alteration are diverse in both magnitude and spatial distribution. There is also a considerable seasonal difference noted, with general higher discharge changes in the summer season (18% on average) compared to winter season (6% on average). In Figures 6 and 7, model sensitivity to Manning's coefficient changes for selected anabranches and their impact on water flow are presented. Although Manning's coefficient changes were tested for all 22 river reaches, only selected cases are shown, which represent characteristic groups of sensitivity (high, moderate and low). As shown, the model response is varied and strongly depends on location of change. The highest sensitivity is noted for reaches at the first confluence ( Figures 6A,B and 7A,B), where flow is distributed between main channel and side anabranches. Manning's coefficient changes in left-hand side reach (Figures 6A and 7A) significantly alters flow in further downstream side branches, but also in the main channel. Similar response is noted for change of Manning's coefficient in the main channel ( Figures 6B and 7B). Much lower sensitivity is noted for downstream side anabranches (Figures 6C and 7C) in which Manning's coefficient changes moderately affected river flow in part of the upstream reaches excluding main channel. There were also situations, especially for short (less than 500 m) reaches, with no impact on water flow noted ( Figures 6D and 7D).  Figures 6B and 7B). Much lower sensitivity is noted for downstream side anabranches ( Figures 6C and 7C) in which Manning's coefficient changes moderately affected river flow in part of the upstream reaches excluding main channel. There were also situations, especially for short (less than 500 m) reaches, with no impact on water flow noted (Figures 6D and 7D).

Impact of Vegetation Changes on Water Flow Distribution
Simulation of vegetation from Group 2 (cf. Table 2) indicated moderate impact on water flow distribution ( Figure 8A). Following sensitivity analysis results ( Figure 6A,C), Manning's coefficient changes in all left-hand side anabranches affected flow distribution in all reaches, including the main channel. Generally, it resulted in redirection of flow by 0.2 m 3 ·s −1 to side anabranches at the expense of the main channel. Although the change might be assumed as low, compared to average flow conditions in summer season (0.17 m 3 ·s −1 ), the contribution was more than doubled. Going further, simulation of vegetation from group (cf. Table 1) intensified the response. The simulation results indicate an overall increase of water flow contribution by 0.5 m 3 ·s −1 , which is three-fold higher than observed during the summer season ( Figure 8B).
The results observed for winter conditions are similar in terms of spatial distribution but slightly higher considering the magnitude of change. An overall increase by 0.6 m 3 ·s −1 was observed which was unequally redistributed between the side anabranches ( Figure 8C). Although the total  Figures 6B and 7B). Much lower sensitivity is noted for downstream side anabranches ( Figures 6C and 7C) in which Manning's coefficient changes moderately affected river flow in part of the upstream reaches excluding main channel. There were also situations, especially for short (less than 500 m) reaches, with no impact on water flow noted (Figures 6D and 7D).

Impact of Vegetation Changes on Water Flow Distribution
Simulation of vegetation from Group 2 (cf. Table 2) indicated moderate impact on water flow distribution ( Figure 8A). Following sensitivity analysis results ( Figure 6A,C), Manning's coefficient changes in all left-hand side anabranches affected flow distribution in all reaches, including the main channel. Generally, it resulted in redirection of flow by 0.2 m 3 ·s −1 to side anabranches at the expense of the main channel. Although the change might be assumed as low, compared to average flow conditions in summer season (0.17 m 3 ·s −1 ), the contribution was more than doubled. Going further, simulation of vegetation from group (cf. Table 1) intensified the response. The simulation results indicate an overall increase of water flow contribution by 0.5 m 3 ·s −1 , which is three-fold higher than observed during the summer season ( Figure 8B).
The results observed for winter conditions are similar in terms of spatial distribution but slightly higher considering the magnitude of change. An overall increase by 0.6 m 3 ·s −1 was observed which was unequally redistributed between the side anabranches ( Figure 8C). Although the total

Impact of Vegetation Changes on Water Flow Distribution
Simulation of vegetation from Group 2 (cf. Table 2) indicated moderate impact on water flow distribution ( Figure 8A). Following sensitivity analysis results ( Figure 6A,C), Manning's coefficient changes in all left-hand side anabranches affected flow distribution in all reaches, including the main channel. Generally, it resulted in redirection of flow by 0.2 m 3 ·s −1 to side anabranches at the expense of the main channel. Although the change might be assumed as low, compared to average flow conditions in summer season (0.17 m 3 ·s −1 ), the contribution was more than doubled. Going further, simulation of vegetation from group (cf. Table 1) intensified the response. The simulation results indicate an overall increase of water flow contribution by 0.5 m 3 ·s −1 , which is three-fold higher than observed during the summer season ( Figure 8B).
The results observed for winter conditions are similar in terms of spatial distribution but slightly higher considering the magnitude of change. An overall increase by 0.6 m 3 ·s −1 was observed which was unequally redistributed between the side anabranches ( Figure 8C). Although the total redistribution rate for winter conditions is higher compared to summer season, the relative change is significantly lower. For summer conditions, the flow increase reached 400%, whereas for winter conditions only 15%. It means that the impact of vegetation on flow distribution in analyzed anastomosing section is significantly higher in growing season.
It is worth stating that uniform change of Manning's coefficient in all side anabranches does not ensure the homogenous response in terms of water flow. As can be noted in Figure 8 for some reaches Manning's coefficient decrease resulted in flow decrease. Moreover, the response varies not only between certain reaches but also between seasons. It proves the complex nature of the system of interconnected channels that is strongly determined by resistance phenomena caused by in-channel vegetation in the first instance, and also by channels geometry (width and depth). It is worth stating that uniform change of Manning's coefficient in all side anabranches does not ensure the homogenous response in terms of water flow. As can be noted in Figure 8 for some reaches Manning's coefficient decrease resulted in flow decrease. Moreover, the response varies not only between certain reaches but also between seasons. It proves the complex nature of the system of interconnected channels that is strongly determined by resistance phenomena caused by in-channel vegetation in the first instance, and also by channels geometry (width and depth).

Discussion
Among studies concerning multi-channel rivers, braided patterns substantially prevail over anastomosing, leaving the latter poorly investigated. This study aimed at assessing the impact of instream vegetation on water flow distribution in the anastomosing river system, using a 1D steadyflow model. As stated by Kale et al. [72], the splitting discharge ratio at a bifurcation might be highly asymmetric due to changes in the planform and section of a river network. Following this assumption, Van et al. [32] developed a 1D steady-flow model for anastomosing section of the River Mekong, where measured discharge data entering each individual channel was not available. They used simple symmetry ratio to estimate the division of discharge at the bifurcations, based on the cross-sectional wetted areas of the splitting channels downstream of the bifurcation node. During the optimization procedure, they adjusted Manning's coefficient to obtain compatibility in water energy levels. It has to be noted that splitting rule based exclusively on the geometry (wetted area) is highly uncertain as it does not take into account vegetation drag force, which also controls the volume of water flowing into the channel. Here, the MFRM is optimized using actual measured entry discharge into each single channel of the complex anastomosing system, which is a significant advantage of this study. Gibson and Pasternack [33] developed a 1D model for anastomosing reach of the River Yuba where only one bifurcation was modeled. With no measurements on discharge division between the channels, the authors found 1D model simulations underwhelming and split flow modeling to be complex. In both studies, Manning's coefficient was used to describe channel roughness. Moreover, Van et al. [32] used Manning's coefficient to estimate the impact of different land cover change scenarios on water flow. Similar to this study, they used fixed Manning's values for different types

Discussion
Among studies concerning multi-channel rivers, braided patterns substantially prevail over anastomosing, leaving the latter poorly investigated. This study aimed at assessing the impact of in-stream vegetation on water flow distribution in the anastomosing river system, using a 1D steady-flow model. As stated by Kale et al. [72], the splitting discharge ratio at a bifurcation might be highly asymmetric due to changes in the planform and section of a river network. Following this assumption, Van et al. [32] developed a 1D steady-flow model for anastomosing section of the River Mekong, where measured discharge data entering each individual channel was not available. They used simple symmetry ratio to estimate the division of discharge at the bifurcations, based on the cross-sectional wetted areas of the splitting channels downstream of the bifurcation node. During the optimization procedure, they adjusted Manning's coefficient to obtain compatibility in water energy levels. It has to be noted that splitting rule based exclusively on the geometry (wetted area) is highly uncertain as it does not take into account vegetation drag force, which also controls the volume of water flowing into the channel. Here, the MFRM is optimized using actual measured entry discharge into each single channel of the complex anastomosing system, which is a significant advantage of this study. Gibson and Pasternack [33] developed a 1D model for anastomosing reach of the River Yuba where only one bifurcation was modeled. With no measurements on discharge division between the channels, the authors found 1D model simulations underwhelming and split flow modeling to be complex. In both studies, Manning's coefficient was used to describe channel roughness. Moreover, Van et al. [32] used Manning's coefficient to estimate the impact of different land cover change scenarios on water flow. Similar to this study, they used fixed Manning's values for different types of vegetation based on the literature and examined its influence on water flow. To our best knowledge, this study was the first that optimized the water flow model in multi-channel river system, based on the broad field measurements, i.e., channel cross-sections, water levels, entry discharges and in-situ vegetation mapping. Detailed fieldwork reduced the uncertainty of the results [50] and made them more reliable for further works. Although it was reported by Merwade et al. [73] that 1D hydraulic models are not suitable to model the detailed hydraulics of multi-channel rivers, the MFRM applied in this study helped to identify impact of vegetation on flow distribution, from which a more comprehensive hydraulic model could be constructed.
Manning's equation used in this study is widely accepted methodology in stream resistance evaluation, especially in vegetated areas [46], and was successfully incorporated in numerous researches investigating the impact of vegetation on water flow, depth and velocity [32,33,47,48]. A quantitative comparison of calibrated Manning's coefficient values with those reported in the literature can be made, linking plant traits (height and flexibility) to roughness. Shields et al. [67] gathered in their study various methods of Manning's coefficient estimation for rigid [74][75][76] and flexible plants [77,78] differentiating their values based on plant height (H) and flow depth (h). In the MFRM, reaches where rigid plants were mapped in the field ( Figure 4A-C), characterized by height ratio 0.5 < h/H < 1, calibrated Manning's coefficient varied between 0.04 m −1/3 ·s and 0.06 m −1/3 ·s, which is in agreement with the literature (0.03 m −1/3 ·s-0.06 m −1/3 ·s) [74][75][76]. Reaches with flexible vegetation mapped in the field ( Figure 3A,B), characterized by height ratio h/H = 1 varied between 0.07 m −1/3 ·s and 0.09 m −1/3 ·s, which is comparable with the results presented by Jarvela (0.07 m −1/3 ·s) [78] and Whittaker et al. (0.08 m −1/3 ·s) [77]. Methods compared by Shields et al. [67] differ in methodology but in general are all supported by field measurements, and therefore their results in most cases make physical sense [67]. Thus, optimized values of Manning's coefficient presented in this study can be assumed reasonable for different plant traits mapped in the field. Comparing the seasonal changes of Manning's coefficient, there is a notable difference between winter and summer season. Similar to study by Song et al. [46], roughness condition is significantly higher in summer compared to winter season. Therefore, the development of models for two different vegetational conditions presented in this study was justified, to capture the seasonal variations in flow conditions determined by the in-channel vegetation.
The impact of vegetation on anastomosing system flow characteristics is to date recognized mainly in theory [20,22,[26][27][28]. According to theoretical understanding, vegetation encroachment causes flow velocity decrease and sediment aggradation, which initiates feedback mechanisms involving further invasion of aquatic plants in the channel and water loss [26][27][28]. Marcinkowski et al. [63] recognized that anabranches loss in the anastomosing section of the River Narew is periodically triggered during the summer low flow conditions, when water discharge in the anabranches promotes sediment vertical accretion. Other studies on anastomosing rivers functioning indicate that localized flow blockages caused by vegetation are one of the key factors controlling their formation [20,22]. Therefore, it can be assumed that the most important factors controlling anabranches maintenance and persistence are related to water and sediment supply. This study, although not incorporating sediment routing, to our best knowledge, is the first that described quantitatively flow distribution under diverse vegetational conditions, using mathematical modeling. Complex system of splitting and rejoining channels of anastomosing river system responses differently to vegetation obstacles than single-channel river. In single-channel rivers flow blockage caused by in-channel vegetation decreases flow velocity and increases water levels, however, water volume flowing through the channel remains constant. In anastomosing rivers, any flow blockage increasing water levels causes impoundment effect upstream to bifurcation node where the flow splits. Subsequently, according to Bernoulli's conservation equation (for details of use in MFRM, see Marcinkowski et al. [65]), a certain volume of water is redirected from one channel (highly vegetated) to another (less vegetated). Therefore, flow distribution at each bifurcation in the anastomosing river depends on the channels geometry [72] and in-channel vegetation characteristics (e.g., height and density) determining their drag force [66]. Model results proved that the impact of vegetation changes under winter flow/vegetation conditions is significantly lower than under summer conditions during growing season. A great advantage of this study is that the MFRM allows to calculate flow redistribution in the entire river system as a response to vegetation changes in particular channels.