Flow Management to Control Excessive Growth of Macrophytes – An Assessment Based on Habitat Suitability Modeling

Mediterranean rivers in intensive agricultural watersheds usually display outgrowths of macrophytes – notably alien species – due to a combination of high concentrations of nutrients in the water runoff and low flows resulting from water abstraction for irrigation. Standard mechanical and chemical control is used to mitigate the problems associated with excessive growth of plant biomass: mainly less drainage capacity and higher flood risk. However, such control measures are cost and labor-intensive and do not present long-term efficiency. Although the high sensitivity of aquatic vegetation to instream hydraulic conditions is well known, management approaches based on flow management remain relatively unexplored. The aim of our study was therefore to apply physical habitat simulation techniques promoted by the Instream Flow Incremental Method (IFIM) to aquatic macrophytes – the first time it has been applied in this context – in order to model shifts in habitat suitability under different flow scenarios in the Sorraia river in central Portugal. We used this approach to test whether the risk of invasion and channel encroachment by nuisance species can be controlled by setting minimum annual flows. We used 960 randomly distributed survey points to analyze the habitat suitability for the most important aquatic species (including the invasive Brazilian milfoil Myriophyllum aquaticum, Sparganium erectum, and Potamogeton crispus) in regard to the physical parameters ‘flow velocity,’ ‘water depth,’ and ‘substrate size’. We chose the lowest discharge period of the year in order to assess the hydraulic conditions while disturbances were at a low-point, thus allowing aquatic vegetation establishment and subsistence. We then used the two-dimensional hydraulic River2D software to model the potential habitat availability for different flow conditions based on the site-specific habitat suitability index for each physical parameter and species. Our results show that the growth and distribution of macrophytes in the hydrologically stable vegetation period is primarily a function of the local physical instream condition. Using site-specific preference curves and a two-dimensional hydraulic model, it was possible to determine minimum annual flows that might prevent the excessive growth and channel encroachment caused by Myriophyllum aquaticum.

Mediterranean rivers in intensive agricultural watersheds usually display outgrowths of macrophytes -notably alien species -due to a combination of high concentrations of nutrients in the water runoff and low flows resulting from water abstraction for irrigation. Standard mechanical and chemical control is used to mitigate the problems associated with excessive growth of plant biomass: mainly less drainage capacity and higher flood risk. However, such control measures are cost and labor-intensive and do not present long-term efficiency. Although the high sensitivity of aquatic vegetation to instream hydraulic conditions is well known, management approaches based on flow management remain relatively unexplored. The aim of our study was therefore to apply physical habitat simulation techniques promoted by the Instream Flow Incremental Method (IFIM) to aquatic macrophytes -the first time it has been applied in this context -in order to model shifts in habitat suitability under different flow scenarios in the Sorraia river in central Portugal. We used this approach to test whether the risk of invasion and channel encroachment by nuisance species can be controlled by setting minimum annual flows. We used 960 randomly distributed survey points to analyze the habitat suitability for the most important aquatic species (including the invasive Brazilian milfoil Myriophyllum aquaticum, Sparganium erectum, and Potamogeton crispus) in regard to the physical parameters 'flow velocity,' 'water depth,' and 'substrate size'. We chose the lowest discharge period of the year in order to assess the hydraulic conditions while disturbances were at a low-point, thus allowing aquatic vegetation establishment and subsistence. We then used the two-dimensional hydraulic River2D software to model the potential habitat availability for different flow conditions based on the site-specific habitat suitability index for each physical parameter and species. Our results show that the growth and distribution of macrophytes in the hydrologically stable vegetation period is primarily a function of the local physical instream condition. Using site-specific preference curves and a two-dimensional hydraulic model, it was possible to determine minimum annual flows that might prevent the excessive growth and channel encroachment caused by Myriophyllum aquaticum.

INTRODUCTION
Aquatic macrophytes play an important role in riverine ecosystems, providing habitats for many organisms and affecting the hydraulic and chemical instream condition (Carpenter and Lodge, 1986). Their distribution and abundance are primarily determined by the hydrologic regime (frequency, duration, and intensity of flood events) (Riis and Biggs, 2003;Franklin et al., 2008), which controls biomass loss and gain processes. Whereas loss processes are caused by increased drag forces during high flood events that cause stem breakage and uprooting of plants, biomass gain processes happen while disturbances are absent during medium to low flow conditions (Riis et al., 2008). In these stable interflood periods, macrophyte growth is controlled by several physical and chemical factors, including flow velocity and depth (Chambers et al., 1991;Riis and Biggs, 2003), light availability (Carr et al., 1997;Köhler et al., 2010), water temperature (Barko et al., 1986;Carr et al., 1997), and riverbed grain size (Baattrup-Pedersen and Riis, 1999), as well as the nutrient content of the riverbed and water (Barko et al., 1986;Demars and Edwards, 2009). Anthropogenic disturbances, such as high nutrient concentrations from water runoff (Jones et al., 2002;Mainstone and Parr, 2002), low suspended sediment concentrations and the resulting increase in light availability from river damming (Madsen et al., 2001;Köhler et al., 2010) and stabilization of the flow regime (less floods) (Riis and Biggs, 2003;Franklin et al., 2008) can alter the ecological equilibrium of the system and have been shown to stimulate excessive growth of aquatic vegetation, notably invasive alien species (Bunn and Arthington, 2002). This is known to cause various forms of ecological and economic damage (Brundu, 2014), including changes in species composition and richness (Bunn and Arthington, 2002;O'Hare et al., 2006), increased flood risk through higher flow resistance (Vereecken et al., 2006;Nikora et al., 2008), and interferences with human water uses such as water abstraction, hydropower, recreation and river navigation (Halstead et al., 2003;Gómez et al., 2013). Management of aquatic macrophytes by mechanical (cutting) or chemical (herbicides) means is therefore common practice in many rivers worldwide (Madsen, 2000;Hussner et al., 2017).
Especially in regulated Mediterranean rivers flowing through intensive agricultural watersheds and presenting prolonged spells of low flows the outgrowth of aquatic vegetation, and notably alien species, is a common phenomenon Aguiar and Ferreira, 2013). Despite their high costs, mechanical control measures are widely applied in Portugal .
Although the growth and distribution of aquatic macrophytes in unshaded streams is mainly influenced by local hydraulic conditions (depth/velocity/sediments) (Chambers et al., 1991;Riis and Biggs, 2003), whose impact overshadows that of hydrochemistry (Steffen et al., 2014), little attention has thus far been paid to the possibility that channel encroachment and invasion can be controlled by establishing minimum annual flows. One common way of exploring the effectiveness of such ecosystem-regulation measures is ecological modeling, because model-based testing is faster and requires less financial inputs than actual physical experiments (Perona et al., 2009;Schmolke et al., 2010). Modeling species distribution or habitat suitability as functions of environmental factors is frequently used to provide spatial decision support for environmental management, weed or pest species risk assessments and studies of climate-change impacts (Franklin, 2013). In the case of river ecosystems, the instream flow incremental method (IFIM) (Bovee, 1982;Raleigh et al., 1986) is probably still the most widely used and accepted methodology for predicting the response of aquatic biota to the instream physical condition (Jowett et al., 2008;Conallin et al., 2010). However, its concepts have never been directly applied to the management of aquatic macrophytes.
Against this background, the main aim of this study was, for the first time, to apply and validate the hydraulic habitat modeling techniques promoted by the IFIM for the assessment of annual minimum flows with the ability to reduce the risk of channel encroachment and invasion by the alien Myriophyllum aquaticum in a heavily regulated Mediterranean river. Our hypothesis was that summer low flows further intensified by water abstraction for irrigation create physical instream conditions that favor the excessive growth of M. aquaticum over the autochthonous Sparganium erectum and Potamogeton crispus, and that this situation can be mitigated by establishing minimum flows above a certain threshold.

Study Area
The study area is located along the Sorraia river in central Portugal (Figure 1). The river basin has an accumulated area of 7719 km 2 and a semi-arid Mediterranean climate in which most of the annual rainfall (600-800 mm) occurs between October and May and the mean annual temperature is 16-19 • C. The fieldwork was carried out along a naturally braided, unconfined segment of the river. The riparian corridor from the edge of the active channel to the adjacent agricultural areas consists mostly of willow shrubs, and willows (Salix alba) in higher areas, and extends an average of 60 m either side of the river. The active channel has an average width of 15 m and is mostly unshaded. The segment's substrate is dominated by sands, gravels and cobbles. Surrounding land is given over to intensive rice, maize, and tomato cultivation. We chose a calibration reach of approximately 1000 m in length for the model-building, and a model reach with a length of 320 m directly downstream for testing and application. Both reaches contain all the different mesohabitats (pool/run/riffle) found in the segment.
The Sorraia's hydrological regime presents a high intraand inter-annual discharge variability, which is characteristic of Mediterranean watersheds (Gasith and Resh, 1999). The mean annual discharge is 20.14 m 3 /s (available data for 1933-1980, "Ponte Coruche" Gauging station). The heaviest winter floods can attain 887 m 3 /s, while during the summer months (June-September) the mean discharge is 3.2 m 3 /s and low flow spells are common (Figure 2).  The area between the upper (0.9) and lower (0.1) quantiles is shaded gray; the black line represents the mean daily discharge; the gray line represents the median daily discharge.
The flow regime is heavily regulated by a system of reservoirs, weirs and canals that was implemented between 1933 and 1958. Water abstraction for agricultural irrigation is managed by a local farmers' association, which mechanically cleans the river channel of aquatic macrophytes and riparian vegetation every few years to reduce flood risk.

Aquatic Vegetation
The main aquatic macrophyte species occurring in study area are M. aquaticum, S. erectum, and P. crispus. Other species that presented less prevalence and were therefore not considered were Ceratophyllum demersum and Typha domingensis. Based on their growth form, M. aquaticum and S. erectum are classified as sediment-rooted plants with floating or emergent shoots/leaves, whereas P. crispus is a sediment-rooted submerged plant (Den Hartog and Van Der Velde, 1988). Following the definition of Pyšek et al. (2004) M. aquaticum is considered an invasive species in Portugal. It was first reported in 1936 (Aguiar and Ferreira, 2013), but massive spreading was only observed in the 1970s . M. aquaticum is displacing native aquatic species, including P. crispus and C. demersum, in many parts of the River Tagus (Ferreira and Moreira, 1995).

IFIM Overview
The instream flow incremental methodology (Bovee, 1982;Raleigh et al., 1986) is a framework which the United States Fish and Wildlife Service developed in the late 1970s to determine appropriate minimum annual flows by considering the effects of flow changes on instream habitat suitability of aquatic biota. It is probably still the most widely used and accepted methodology for predicting the response of aquatic biota to the instream physical condition (Jowett et al., 2008;Conallin et al., 2010). Its main feature is a hydraulic habitat suitability model that can be separated into a hydraulic component and a habitat component. The hydraulic model predicts water velocity, depth and other hydraulic variables. The habitat model is based on local habitat suitability curves (HSC) that describe the optimum range of a physical parameter affecting the species and are built on expert knowledge or field analyses of local species occurrence and habitat availability. Integrating the two components makes it possible to calculate a composite suitability index (CSI) that combines the suitability information for each physical parameter at a given flow. The weighted usable area (WUA) for the target species is quantified by multiplying the CSI by its area of influence. In order to assess an appropriate minimum annual flow, the hydraulic habitat suitability model is applied to a range of flows to produce a WUA-vs.-discharge graph.

Hydraulic Habitat Suitability Modeling
In order to calibrate (train) the habitat suitability model, a total of 961 sample points were distributed systematically (2 m × 2 m), with a randomly chosen starting point along each mesohabitat (pool, run, and riffle) found in the calibration reach. The mesohabitats were visually delimited in the field.
The occurrence of the main macrophyte species and physical habitat characteristics -flow velocity, water depth and grain size of the bed material -were analyzed at each sample point. The fieldwork was done in August 2016 and July 2017, during measured discharges of around 0.3 m 3 /s. We chose the lowest discharge period of the year in order to assess the hydraulic conditions during the period of least disturbance, which allows aquatic vegetation establishment and subsistence. Locations shaded by riparian vegetation (less than 5% of the analyzed reach) were excluded, since in this situation aquatic plant growth is mainly constrained by insufficient light (Carr et al., 1997). Depths were measured with a simple meter ruler and classified in intervals of 20 cm. Flow velocities were measured with a water flow probe (model FP101, Global Water Instrumentation, United States) positioned in the flow direction at 60% of the flow depth and using 0.05 m/s intervals. The bed grain size was assessed visually and classified according to the Wentworth scale (sand: 0.62-2 mm; gravel: 2-64 mm; cobble: 64-256 mm). The habitat preferences for M. aquaticum, S. erectum, and P. crispus were then calculated by dividing habitat-utilization (amount of species occurrences in each class of the physical parameters) by habitat-availability (total amount of each class of the physical parameters). The final preference values were normalized, from a minimum value of 0 for unsuitable to 1.0 for optimal habitats (the class of the physical parameter with the highest amount of species occurrences), and expressed as a HSC for each physical parameter.
In order to apply and test the hydraulic habitat suitability model, we selected a 320 m-long reach directly downstream from the calibration reach. We chose a two-dimensional approach for the hydraulic simulation: the River2D model (Steffler and Blackburn, 2002). Two-dimensional hydraulic models predict depth and velocity laterally and longitudinally along the whole length of the river channel. They are therefore better able to simulate the complex flow patterns found in braided rivers than the more conventional (with regard to the IFIM) onedimensional models that only predict depth and velocity across channel transects (Benjankar et al., 2015). The topography of the riverbed of the model reach, which is the main input into the hydraulic model, was measured in July 2016 with a Leica TCR703 Total Station (angle accuracy 3 ) along 970 points. The initial bed roughness values were estimated based on substrate size and vegetation distribution. To determine the boundary condition and calibrate the model, water depth and velocity were assessed along six transects including the down-and upstream cross-section, with measurements taken every 20 cm along the cross-section. The hydraulic model was calibrated by adjusting bed roughness until simulated water surface elevations matched measured water surface elevations.
The model was then used to simulate the physical instream conditions for a series of potential annual minimum flows of between 0.3 and 10 m 3 /s, representing a common flow range during the vegetation period. The WUA concept was used to evaluate the shift in habitat suitability for each discharge (Bovee, 1982). The WUA computation is based on the habitat suitability evaluated at every node of the topographic mesh and the "tributary area" of that node. We also calculated the Hydraulic Habit Suitability (HHS) for each discharge by dividing the WUA by the inundated area. The HHS can be understood as the percentage of the WUA from the inundated area at a given discharge. A value of 1 would mean that the whole of the wetted area classifies as usable area for a certain species or species group.
We used two different methods to calculate the habitat suitability. The classical, deterministic approach of the IFIM calculates a CSI as the geometric mean of the separate suitability indices for depth, velocity, and substrate size. It is directly integrated into the River2D Model on the basis of the HSC for each species.
VSI -Velocity Suitability Index DSI -Depth Suitability Index SSI -Substrate Suitability Index In addition to the deterministic approach, we computed the habitat suitability for each species based on the random forest algorithm (RF) for classification (Breiman, 2001). We used the R package "randomForest" (Liaw and Wiener, 2002) to grow 1000 trees based on bootstrap samples of the same training data as that used to build the HSC, and incorporated 50% class weights into the classifier to account for the low prevalence of P. crispus and S. erectum.

Model Validation
We mapped the true presence and absence of the main macrophyte species (M. aquaticum, S. erectum, and P. crispus) in the model reach with a Global Positioning System unit (Ashtech, model Mobile Mapper 100; accuracy < 50 cm) during the same period (July/August) and with the same discharge (0.3 m 3 /s) as those when the data for the model calibration was collected. We then modeled the macrophyte distribution using the deterministic and the random forest approach based on the hydraulic simulation for the same discharge, and tested the agreement between observed and predicted distribution by assessing the area under the receiver operating characteristic curve (AUC) (Fielding and Bell, 1997). The AUC of a model is equivalent to the probability that the model will rank a randomly chosen species-presence site higher than a randomly chosen absence site. In addition, we transformed the predicted occurrence probabilities of both models to a binary presence/absence format for each species using the threshold of occurrence that maximizes the sum of sensitivity and specificity (Cantor et al., 1999;Liu et al., 2005). In order to assess the accuracy of the binary classification, we used the "True Skill Statistic" (TSS; sensitivity + specificity -1), because it accounts for the effect of the species prevalence (Allouche et al., 2006). All accuracy measurements were carried out using the R package "SDMtools" (VanDerWal et al., 2014).
In order to investigate whether our models accounted for all the factors causing the species' distributional pattern, we checked the observed species occurrence in the model reach for spatial autocorrelation using the Ripley's K function, and tested the error between observed and predicted species occurrence for clustering with the Moran's I index. The spatial analyses were done with the spatial statistics toolbox from ArcGIS for desktop (version 10.4.1).

Habitat Suitability Curves
The habitat sampling resulted in 224 M. aquaticum, 135 P. crispus, and 85 S. erectum presences in a total of 961 habitat samples.
Myriophyllum aquaticum displayed a substantial liking for low flow conditions, only having colonized areas with relatively slow velocities and low depth. It was already nearly absent at velocities over 0.1 m/s. The most suitable depths were 0-20 cm. In addition, it was found almost exclusively on sandy substrate. On the contrary, P. crispus seemed to prefer higher-flow areas. Its greatest presence occurred in medium velocities of 0.08-0.2 m/s and it clearly favored depths of more than 80 cm. Its preferred substrate was gravel. S. erectum displayed a preference profile similar to that of M. aquaticum, but was more tolerant of greater depth. The results show a distinct preference profile of the exotic M. aquaticum with regard to flow velocity and water depth (Figure 3).

Model Validation
In overall terms, the hydraulic habitat model based on the deterministic approach displayed a good discriminatory ability. In the case of M. aquaticum, accuracy was even in the excellent range (AUC = 0.9), while for P. crispus it was good (AUC = 0.87), and for S. erectum fair (AUC = 0.79). The performance of the binary classification differed more drastically between the species. Considering a threshold of occurrence for M. aquaticum of 0.24, the TSS score of the model was 0.66. It correctly predicted 86% of the actual presences (sensitivity) and 80% of the actual absences (specificity). The occurrence threshold for Potamogeton was set to 0.24. The TSS score was 0.62. Its occurrence was correctly predicted in 88% of cases, and its absence in 70%. The model's worst performance was for S. erectum, with an occurrence threshold of 0.08 (TSS = 0.44; Sensitivity = 0.7; Specificity = 0.66).
The random forest model did not perform as well as the deterministic approach. On the contrary, only the prediction of M. aquaticum achieved a similar accuracy (AUC = 0.85), whereas the predictions for P. crispus (AUC = 0.7) and S. erectum (AUC = 0.65) were less accurate. This was also visible in the binary prediction. Considering a threshold of occurrence of 0.6 for M. aquaticum, the model's TSS score was 0.66 (sensitivity = 0.8; specificity = 0.86). The prediction of Potamogeton based on a threshold of 0.5 returned a TSS score of 0.38 (sensitivity = 0.66; specificity = 0.72). Once again, the model performed worst for S. erectum (threshold = 0.2; TSS = 0.28; sensitivity = 0.66; specificity = 0.62).
The species occurrence as well as the errors between the observed and predicted distributions presented a similar degree of positive spatial autocorrelation (clustered pattern), indicating that although our models have a medium to high degree of accuracy, they do not account for all the factors explaining the species distribution.

Weighted Usable Area and Hydraulic Habitat Suitability
We only used the deterministic modeling approach to analyze the shifts in habitat suitability for incremental flows because of its better predictive performance.
The preference of M. aquaticum for low flow conditions is also reflected in the development of the WUA. From 1167 m 2 at Q = 0.3 m 3 /s, it rapidly increases until it reaches its maximum of 3085 m 2 at Q = 1.4 m 3 /s. The WUA drops steadily after that, although the inundated and therefore potentially invadable area continues to increase with rising flows. The WUA decreases more slowly from Q = 5 m 3 /s to Q = 8 m 3 /s, after which it remains nearly constant. At Q = 0.3 m 3 /s P. crispus has a WUA of 1017 m 2 , slightly lower than that of M. aquaticum and S. erectum. However, this then sharply increases, so that at Q = 3 m 3 /s the P. crispus WUA of 8004 m 2 is already three times higher than that of M. aquaticum. After that, the upward trend continues more slowly, but steadily. At Q = 10 m 3 /s, the P. crispus WUA of 10569 m 2 is over 10 times that of the invaders. The development of the WUA of S. erectum initially appears to be similar to that of M. aquaticum. However, it continues to gain area until Q = 3.5 m 3 /s, after which the WUA stays relatively constant at around 3900 m 2 , whereas the M. aquaticum WUA experiences a steady decline over the same range ( Figure 4A).
In the case of M. aquaticum, the HHS trends continuously downward as discharge increases. Whereas 36% of the wetted area is potentially suitable at Q = 0.3 m 3 /s, only about 10% remains suitable at Q = 4 m 3 /s. P. crispus experiences an increase in HHS with rising flows. The HHS only decreases slightly at around Q = 1 m 3 /s, due to a large increase in wetted area. From Q = 3.5 m 3 /s onward, the rate of change in HHS decreases. S. erectum also experiences a decline in HHS, sharply at first, to levels below even those of M. aquaticum, but remains nearly constant from Q = 2.5 m 3 /s onward ( Figure 4B).

DISCUSSION
In this study we wanted to explore setting minimum annual flows as an alternative management approach for controlling excessive  growth of macrophytes and invasion by M. aquaticum during the vegetation period in the Sorraia river. Following IFIM principles, we built a hydraulic habitat suitability model for M. aquaticum, S. erectum, and P. crispus, applied it to a range of discharges, and analyzed the changes in WUA and HHS. Our hypothesis was that low summer flows intensified by water abstraction for irrigation create physical instream conditions that stimulate excessive growth of M. aquaticum, and that this situation can be mitigated by establishing minimum flows above a certain threshold.
The modeling results support our hypothesis that the growth and distribution of macrophytes in interflood periods is primarily a function of the local physical instream condition, which is especially favorable to an invasion of M. aquaticum during the low flow range. Habitat suitable for M. aquaticum already declines above flows of 1.4 m 3 /s, while the autochthonous species, and especially P. crispus, continue to gain ground. It would therefore seem possible to reduce the risk of invasion and favor a more natural species composition by setting annual minimum flows. The combination of the artificial approximation of the habitat availability for both the exotic and the autochthonous species caused by stable periods of flows under 1.4 m 3 /s and the greater competitive ability of M. aquaticum may be the reason for the latter's successful expansion. Given that the mean annual flow during the vegetation period is 3.2 m 3 /s, it may well be that water managers can establish minimum annual flows above the 1.4 m 3 /s threshold and thereby avert this situation. This is an important result that can improve river restoration projects by preventing the degradation of natural aquatic vegetation communities.
However, we also observed that for the low flow range (0.3-1.4 m 3 /s), the WUA actually increases for M. aquaticum and that the rate of change in habitat suitability for all species is lower with high flows than with low flows. The explanation for this is that the suitable areas are concentrated in shallow waters along the banks of the stream, and these shallow areas initially increase when the river enters the floodplain and then remain relatively constant in size. In the case of M. aquaticum, this means that the WUA remains relatively constant above a discharge of 7 m 3 /s. Setting minimum annual flows will therefore not completely prevent an invasion; but it can contribute to an environmental flow regime that privileges autochthonous aquatic species and strengthens their competitive performance.
One major criticism of the IFIM habitat simulation to keep in mind when interpreting the results is the usage of the term WUA (Mathur et al., 1985), because it suggests a spatial extension of usable habitat when in fact it only actually describes the overall probability of use. So when we assess the effects of flow changes on aquatic biota, it is the shape of the WUA response curve that is more important than the magnitude (Jowett et al., 2008).
In addition, as with all modeling approaches, there are a number of different uncertainties that should be considered when interpreting the results.

Environmental Factors
Our study is based on the assumption that in hydrologically stable periods, physical habitat characteristics are the main limiting factor for aquatic species in streams. Indeed, several studies argue that flow velocity is the main environmental factor controlling the abundance and distribution of aquatic macrophytes (Chambers et al., 1991;Baattrup-Pedersen and Riis, 1999;Madsen et al., 2001;Janauer et al., 2010). Most studies relate the limiting effect of higher flow velocities on plant growth to increased drag forces on the plants and their anchoring ground, causing uprooting, or less frequently, stem breakage (Chambers et al., 1991;Riis and Biggs, 2003). However, a more recent study (Pollen-Bankhead et al., 2011) indicates that the preference of macrophytes for low velocities is less related to the drag forces on the plants and more to the conditions controlling erosion and deposition of fine substrate materials. The effect of substrate size has mainly been studied with regard to the distribution patterns of macrophytes, and not in terms of changes in biomass (Baattrup-Pedersen and Riis, 1999;Riis and Biggs, 2001;O'Hare et al., 2006). The findings indicate a niche separation between macrophytes based on different substrate size preferences. Apparently, submerged species favor coarser substrates (gravel and boulder), whereas species that grow both submerged and emergent, and species that only grow emergent, were associated with finer substrates (sand) typical of low flow conditions. This is coherent with our results. The influence of flow depth has been related to light availability, which decreases with greater depth (Koch, 2001). In situations of high turbidity or direct shading, for example through overhanging vegetation, light availability can also become the main limiting factor, which is why we excluded sample sites with these characteristics (Köhler et al., 2010). Temperature is also known to influence the growth rate of aquatic plants (Koch, 2001). It can, however, be assumed that temperature alterations in the analyzed flow range are marginal and are indirectly covered by the effects of velocity and depth (Gu et al., 1998). Besides the physical factors, geochemical properties of the stream and especially nutrient availability are known to have an influence on aquatic biota (Koch, 2001). Unnatural high concentrations of phosphorus, as often occur in agricultural watersheds, can stimulate excessive macrophyte growth (Mainstone and Parr, 2002). However, these factors are still most probably overshadowed by the hydraulic conditions (Barendregt and Bio, 2003;Steffen et al., 2014), as is also indicated by the high accuracy of our model.

Data Collection/Model Calibration
Different forms of data analysis for generating the HSC for each environmental factor are distinguished for the IFIM (Bovee, 1986): (a) expert knowledge; (b) analyses of actual habitat conditions used by the species (or presence only data); and (c) in situ species occurrence and habitat availability data (or presence/absence data). We based our model calibration solely on actual presence/absence data (c). It is the most highly recommended of the three methods (Jowett et al., 2008), and the only one that permits an estimation of the true probability of observing a species at a site (Guillera-Arroita et al., 2015). We kept geographical sampling bias to a minimum by selecting a calibration (training) reach and a model reach from the same river segment, and by applying a stratified, systematic sampling design with a random starting point. The detection error, which is crucial to the performance of many habitat suitability models (Lahoz-Monfort et al., 2014), can be considered negligible because of the sampling design, the small number of different species and their sessility.
Model calibration errors can also affect the two-dimensional hydraulic modeling, which can be compromised due to the collection of insufficient or erroneous bed topography data, insufficiently detailed substrate distribution mapping, erroneous model calibration, or failure to include effects of the bed topography upstream of the study site in the model (Jowett and Duncan, 2012).

Model Algorithm
The IFIM commonly uses a univariate algorithm to relate the abiotic characteristics to actual habitat suitability (Conallin et al., 2010). The univariate derivation of the CSI is criticized for being based on the assumption that organisms select each habitat variable independently, ignoring interactions and cumulative effects between them (Ahmadi-Nedushan et al., 2006), such as the influence of velocity on substrate stability and composition (Shields, 1936). Multivariate statistical models, such as Generalized Additive Models (Milner et al., 2001) and Artificial Neural Networks (Gozlan et al., 1999), are alternative means of fitting the suitability data that are able to account for interactions between the variables and overcome the problem of independence (Ahmadi-Nedushan et al., 2006). Another, increasingly popular, approach is the use of "fuzzy logic" to define a set of rules that classifies suitability according to a combination of different environmental factors. It allows consideration of uncertain measurements and vague expert knowledge, as well as multivariate effects, without requiring the input parameters to be independent (Noack et al., 2013). With random forests we also applied a distribution modeling technique that is capable of modeling complex interactions among predictor variables and is considered to have one of the greatest discriminatory capacities (Elith et al., 2006;Cutler et al., 2007).
However, random forest and all other approaches are static and ignore more complex processes that are known to shape the distribution patterns of macrophytes, such as interspecific competition and feedbacks between the plants and the physical environment known as niche construction (Corenblit et al., 2009). The latter has become very evident in the complex relationship between macrophytes and fine sediment, where macrophytes have been observed to create positive growth conditions through retention and stabilization of fine sediments, thereby also interacting with geomorphological processes (Schoelynck et al., 2012).

Model Validation
Ecological modeling is of little value if the prediction is not tested against independent data (Olden et al., 2002). We therefore separated the study reach from the calibration reach and collected field data in two different years. The overall model prediction capacity at Q = 0.3 m 3 /s was assessed as good using the thresholdindependent AUC statistic. The binary prediction, and especially the rate of observed absences of the species that fall in pixels of predicted presences (the commission error rate, which equals 1 minus specificity), was less convincing, but can in part be explained by the low prevalence of the species. A distinction must be made between two different types of commission error: real commission errors, in which combinations of environmental conditions that are not within the species' niche are falsely interpreted as suitable; and apparent commission errors, where absence represents a real feature of the species' distributional ecology due to interspecific interactions and historical factors (Peterson, 1999). A high commission error is therefore common among species that show a low prevalence, and can be an indicator that the species has not yet conquered the whole of its potential niche. If this interpretation is correct, it would support the use of our model as a screening tool for identifying areas that are at higher risk of invasion.
We can only speculate about the causes of the spatial autocorrelation in the errors between observed and predicted species distribution: disregard of interactions between the predictor variables, omission of important predictors (temperature and nutrients), or ecological processes (dispersal, competition, and niche construction) (Guisan and Thuiller, 2005). However, the model's good predictive performance against independent data nonetheless proves the usefulness of the IFIM approach for predicting macrophyte distribution.

Other Management Options and Conclusion
Mechanical methods are the most widely used measures for controlling aquatic macrophytes in both Portugal  and Europe as a whole (Hussner et al., 2017). They allow for containment or eradication, depending on the specific technique and frequency of application (Madsen, 2000). Although often regarded as environmentally less harmful, the most common and effective measures like mowing are not species-specific and can both harm non-target aquatic biota and cause sediment resuspension (e.g., Habib and Yousuf, 2014). Worldwide, chemical control is also applied. While proven very effective, even for eradicating nuisance weeds (Champion and Wells, 2014), herbicides will physiologically affect similar native aquatic plants and potentially also indirectly harm fish and invertebrates (Getsinger, 1998). The use of herbicides to control aquatic nuisance weeds is therefore severely restricted in various countries (especially in the EU). Biological measures also present a risk of off-target impacts, both directly and indirectly through alteration of the food web. Physical management methods are distinguished from mechanical techniques, because instead of the plants directly, it is their environment that is manipulated. Several physical techniques can be distinguished: dredging, drawdown, benthic barriers, shading or light attenuation, and nutrient inactivation (Madsen, 2000;Wersal et al., 2013). The control of nuisance weeds through flow regulation fits into the latter category, but has so far received little attention. Flushing flows have been successfully used to eradicate weeds in the Ebro river (Tena et al., 2013). However, frequency and magnitude of discharges (in the range of a 2-year flood) are not a viable option for intensive agricultural watersheds like the Sorraia, where both the side effects of the floodings and the competing water uses have to be considered.
Although most management techniques have some negative side effects on the ecosystem, so do the invasion and extreme growth of alien species. Maintaining minimum discharges in order to prevent channel encroachment may be an ecologically and financially advantageous addition to the range of commonly practiced control measures. We tested this approach by applying habitat suitability modeling techniques that are widely used to evaluate environmental flows and restoration measures aimed at fishes and invertebrates. Based on the specific habitat preferences of M. aquaticum, it seems possible to set minimum flows that reduce the invader's habitat while simultaneously promoting that of autochthonous and less invasive aquatic species. This measure can be recommended with a high level of confidence, given that when the model was checked against independent data, it displayed a good level of accuracy in predicting species distribution.

AUTHOR CONTRIBUTIONS
KO, RR, TF, and GE made substantial contributions to the conception and design of the study and the acquisition, analysis, and interpretation of data. They participated in drafting the article and revising it critically for important intellectual content. Gave final approval for the version of the manuscript to be submitted. Took public responsibility for the content.