Macroinvertebrate habitat requirements in rivers: overestimation of environmental flow calculations in incised rivers

Flow variability determines the conditions of river ecosystems and river ecological functioning. The variability of ecological processes in river ecosystems gradually decreases due to river channelization and incision. Prediction of the environmental flow allows us to keep biological diversity and river health developed as a response to the degradation of aquatic ecosystems overexploited by humans. The goal of the study was to test the influence of river incision on environmental flow estimation based on the Biological Monitoring Working Party (BMWP) macroinvertebrate index. A total of 240 macroinvertebrate assemblages of 12 waterbodies differing in bed substrate and amplitude of discharge were surveyed in southern Poland. Variations in the distribution of 151 466 macroinvertebrates belonging to 92 families were analysed. The similarity of benthic macroinvertebrates reflects the typological division of the rivers into three classes: Tatra mountain streams, mountain flysch rivers, and upland carbonate and silicate rivers. As a response variable reflecting the macroinvertebrate distribution in the river, the BMWP_PL index was chosen. The river incision significantly increased the values of e-flow calculations in relation to redeposited channels. The area of optimal habitat for macroinvertebrates decreased with the bed incision intensity. In highly incised rivers, the environmental flow values are close to the mean annual flow, suggesting that a high volume of water is needed to obtain good macroinvertebrate conditions. As a consequence, river downcutting processes and impoverishment of optimal habitats will proceed.

1 Introduction

Human water demand, including irrigation to increase crop productivity, dams, reservoirs to control the timing of stream flow, and water withdrawal from rivers, has increased dramatically over the last 100 years (Vörösmarty et al., 2010; Veldkamp et al., 2017). Maintenance of a suitable water flow in an active river channel should not only secure human needs, but above all also ensure the proper functioning of aquatic ecosystems (Anderson et al., 2006). This has become particularly important since riverbeds began to be perceived not only as channels filled with water, but also as complex ecological systems in which biological elements play a key role (Poff et al., 1997; Bunn and Arthington, 2002; White et al., 2016). The Water Framework Directive (WFD; European Community, 2000/60/EC, 2000) was introduced by European countries to protect and improve the state of aquatic ecosystems and formalize a water flow framework that would maintain this state (Chen and Olden, 2017).

Discharge intensity is one of the most important factors influencing communities of aquatic and water-dependent organisms (Tharme, 2003; Arthington et al., 2006; Higgisson et al., 2019). It is a parameter which shapes the morphology (Michalik and Książek, 2009) and hydraulic flow conditions (water depth and flow velocity), and it influences the diversity and quality of habitats for fauna and flora in the active channel and in the floodplain (Allan, 1995; Poff et al., 1997; Ward and Tockner, 2001; Skalski et al., 2016, 2020). Furthermore, flow significantly influences abiotic elements, such as water temperature and oxygenation, as well as nutrient cycles in the aquatic ecosystem (Monk et al., 2008; Laini et al., 2019). This applies in particular to rivers subjected to strong human impact (e.g. channel regulation and incision, dams, or retention reservoirs, as well as a continuous increase in water abstraction). Artificial restriction and control of a range of water flow values lead to substantial impoverishment of biological diversity (Pander et al., 2019). Environmental flow is the amount of water required to maintain biological diversity in the river ecosystem (Arthington et al., 2006). This definition requires quantification of the ecological response of aquatic elements to flow alteration, for which data are rather scare in the literature (Poff and Zimmerman, 2010). Therefore, it appears crucial to estimate empirical ranges of environmental flows that ensure optimal habitat conditions for living organisms (Bunn and Arthington, 2002; Acreman et al., 2014).

Environmental flow has been studied by many researchers, resulting in numerous methods for its determination. The simpler ones include hydrological methods, which are based on historical hydrological data and mean annual discharge (Tennant, 1976; Jowett, 1997; Tharme, 2003; Rosenfeld, 2017). Analysis of such data makes it possible to specify a percentage of the mean annual flow as the critical value below which severe degradation of biotic elements occurs. Unfortunately, hydrological methods do not take into account the morphology of the riverbed, which is a key factor shaping the river habitat (Książek et al., 2020). Therefore, a number of hydraulic methods based on simple hydraulic variables such as critical riffle analysis and wetted area/wetted perimeter have been introduced (Gippel and Stewardson, 1998; Książek et al., 2019). Determination of discharge values (Q) for environmental flow involves defining the breaking point of the hydraulic variable discharge curves as the e-flow (Gippel and Stewardson, 1998; Vezza et al., 2012; Tare et al., 2017). Over time, hydraulic methods have developed in the direction of habitat simulation methods. They have additionally focused on the habitat requirements of selected groups of model organisms, most commonly water depth, flow velocity, and bed substrate (Jowett and Davey, 2007; Li et al., 2009; Muñoz-Mas et al., 2016). Based on the analysis of these environmental factors, habitat–discharge curves were drawn for organisms, and from these it was possible to read the optimal flows maintaining the normal ecological functions of aquatic ecosystems. Another type of method, which emphasizes the importance of the natural flow regime for the entire ecosystem, is holistic methods. They attempt to maintain the natural flow regime as well as flow variability. In this case, environmental flow is defined in the category of deviation from the natural flow regime (Yarnell et al., 2015).

The methods presented above focus on the fish distribution and rarely on diversity and availability of habitats for freshwater macroinvertebrates, which are the most important and sensitive indicators of the ecological state of the ecosystem (Jowett et al., 2008; Birk et al., 2012). The diversity and taxonomic composition of aquatic organisms living in freshwater streams and rivers are used as indicators in the evaluation of environmental flow (Pander et al., 2019). In many cases, macroinvertebrate assemblages are considered (Hayes et al., 2014; Laini et al., 2019), as numerous studies confirm that they are relatively good indicators of ecological water quality and integrity (Buss et al., 2015; Wyżga et al., 2016; Schneider and Petrin, 2017). Freshwater macroinvertebrates also play an important role in the processing of nutrients and organic energy in running water ecosystems, as well as in sustaining ecosystem integrity.

Another parameter, which is usually neglected in flow modelling, is associated with morphological channel modification and incision (Wyżga et al., 2012; Skalski et al., 2016). Incision and channel simplification constitute a global problem overwhelming most of the rivers in the mountain as well as in upland areas (Skarpich et al., 2020). During the last 100 years, anthropogenic processes related to river regulation (narrowing and straitening) have disturbed fluvial processes, leading to enormous river incision (Rinaldi et al., 2005; Wyżga, 2007). As a result, rivers have become vertically closed systems, losing the ability to store alluvial material. Moreover, incision up to the bedrock simplifies the microhabitat array of the river (Neachell, 2014) and leads to elimination of most of the habitats (Muñoz-Mas et al., 2016), as well as affecting ecosystem functioning (biodiversity loss and food web network simplification; Shields et al., 1998; Jeffres et al., 2008).

The goal of the study was to test the influence of river incision on environmental flow estimation based on the Biological Monitoring Working Party (BMWP) macroinvertebrate index. Specific aims of the study were as follows: (1) to establish the habitat preferences of macroinvertebrates communities (240 local assemblages) in mountain and upland rivers using generalized additive models, (2) to calculate the e-flow values combining the habitat requirements and hydraulic method of environmental flow calculation in relation to river hydromorphological parameters (redeposition and incision), (3) to identify reality of providing e-flow values for different hydromorphological modifications in relation to available amount of water (low low flow, mean low flow, and mean annual flow) and (4) to check and visualize the e-flow values in relation to available water volume on randomly chosen, incised, and redeposited rivers based on the CCED2D model. We expected that e-flow in incised rivers, allowing us to obtain the shelf zone level of the river, should be much higher than mean low flow. Such an assumption could determine the consecutive higher discharges and increase the bed degradation. Firstly, we should restore the sedimentation processes in incised rivers to obtain a hydrodynamic balance and then manage the proper volume of water. As a consequence, optimal habitats for invertebrates and fish will be enlarged.

2 Materials and methods

2.1 Study sites

The survey was conducted in 12 mountainous rivers assigned to three typological groups according to the Polish Water National Authority and the Water Framework Directive (Jusik et al., 2014): Tatra mountain rivers (Biały Dunajec, Dunajec, and Białka – Group 1), mountain flysch rivers (Raba, Brynica, Toszecki Potok, and Nysa Kłodzka – Group 2), and upland carbonate and silicate rivers (Sołokija, Warta, Ropa, Biała, and Odra – Group 3) (Fig. 1), varying in bed modification (incision intensity or redeposition).

Figure 1 Map of the studied mountainous rivers in the Carpathian–Sudetian region of Poland.

The first group comprises rivers located in an alpine granitoid region, characterized by calcareous and silicate bedrock. The second group consists of rivers flowing through much lower mountain ranges (up to the timber zone), where the bedrock contains sandstone rock formations. The third group represents rivers of upland landforms with various carbonate and silicate sediments and rocks. The typology of river channel modification was obtained from field observation and channel measurements (cross-sections, longitudinal profile and cover, height of the floodplain). Narrow channels with downcutting to the floodplain and simplified channel morphology were defined as incised.

All rivers are routinely monitored by the nearest monitoring station of the Environmental Agency (Environmental Agency Data, 2018), and all 12 rivers have consistently been assigned a similar average chemical status in recent years. Analysis of variance (ANOVA) showed no variation between the river groups in incision bed modification (F=1.56, p=0.26) as well as in the following physicochemical properties: dissolved oxygen, conductivity, hardness, pH_max, NH₃, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, total N, and ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{3}-}$. Only water temperature and pH_min significantly depended on the river group. All habitat variables (flow, depth, and substrate type) were significantly dependent on river group (Table 1); meanwhile the incision was not influenced by the parameters variation.

Table 1 Mean values ± standard deviation of the physicochemical and habitat variables of the three river groups, with results of one-way ANOVA. Items in bold indicate statistically significant differences.

2.2 Macroinvertebrate sampling

Benthic invertebrate samples were collected in two seasons: autumn (October 2017) and spring (April 2018). No flood waves occurred between these surveys, and the channel morphology remained the same throughout the sampling period. We collected 20 subsamples (1 m² each subsample) from each low-flow channel along a representative 100 m section of each river according to the sampling procedure for the BMWP_PL index (Bis and Mikulec, 2013). A total of 480 subsamples were taken from a wide range of water depths and flow velocity. Following Jowett et al. (1991) and Muñoz-Mas et al. (2016), the substrate types were converted to a single index by summing the weighted percentages of each type.

Macroinvertebrate samples were collected with a D-frame net according to the Environmental Agency's sampling protocol for biomonitoring assessment using a kicking motion for 3 min across all habitats (Bis and Mikulec, 2013). All collected material was preserved in the field with 4 % formaldehyde. Aquatic macroinvertebrates were separated from the rest of the material in the laboratory using a stereoscopic microscope, and then they were identified to the family level (Tachet et al., 2000), except Oligochaeta, Porifera, and Hydrozoa, which were recorded as such. Due to the varied preferences of macroinvertebrates to habitat conditions, the BMWP_PL index was adopted as the best qualitative index. The Biological Monitoring Working Party (BMWP) is one of the most commonly used biotic indices in various rivers and streams around the world (Roche et al., 2010; Wyżga et al., 2013). It has been adopted in many countries, including Poland (Dz.U. 2019 poz. 2149, 2019). The BMWP index was originally developed to represent water quality, but subsequent studies showed that it reflects ecological quality of the waterbodies and can also be related to hydromorphological impoverishment such as incision or straightening (Mutz et al., 2013; Wyżga et al., 2013; Mikuś et al., 2021). This index best considers the sensitivity of invertebrates to environmental variables because families with similar stress tolerances are grouped together (Armitage et al., 1983).

2.3 Data analysis

ANOVA was used to verify the statistical significance of the differences in environmental data between the three river groups (Statsoft, 2013). Non-metric multidimensional scaling (NMDS) was used to test the relationship between the macroinvertebrate taxonomic composition of the assemblages of the 12 rivers assigned to three groups (Group 1, Group 2, and Group 3) and hydromorphological variables (water velocity and depth) during the spring and autumn. Descriptive physical properties (water depth and velocity) were classified into two or three categories: low, medium, and high. We used minimum and maximum values of depth and velocity range in each river group and divided them into 33 percentile ranges of the total value variability. In the case when the ranges were less than 0.5 m depth, we have chosen two groups of 50 percentiles of the depth ranges. The significance of differences between depth and velocity classes was tested by analysis of similarities (ANOSIM) (p values of pairwise comparison with Bonferroni correction) on the Bray–Curtis dissimilarity matrix with 499 permutations of the data. PAST software (version 3.13) was used to analyse NMDS and ANOSIM (Hammer et al., 2001).

To develop habitat suitability functions of macroinvertebrates, reflecting the optimal conditions in the river, generalized additive model (GAM) procedures were chosen. The advantage of the method described by Jowett and Davey (2007) is that it calculates the probability of relations between dependent biotic variables and independent flow parameters. To choose the best-fitting model, we have ranked the available models according to the Akaike information criterion procedure and Δ AICc values, which reflect the difference of AIC between a given model and the lowest AIC. The best-fitting model, describing the relationship between independent variables (depth and velocity and two-way interaction between them) and the macroinvertebrate BMWP_PL index, was the generalized additive model with Poisson error distribution and log link function. We have also measured the accuracy of the GAM procedures (Shearer et al., 2015). The total deviance explained calculated as the relative difference between the residual and the null deviances of the model ([null deviance-residual deviance] / null deviance) was adopted. The course of the regression line of the BMWP-PL and depth and velocity for each group of the bed material rivers was obtained using CurveExpert software, where the best-fitted line for the set of non-linear curves was applied and ranked. The BMWP_PL curve maximum values were regarded as the most optimal for invertebrates and the most preferred. We were interested in calculation of optimal condition for depth and velocity separately to obtain the optimal conditions, allowing us to calculate the discharge which is needed for hydraulic and CCHE2D modelling. The preferred depths and velocities for each season and riverbed material groups were used to calculate the hydraulic discharges which are the most optimal for BMWP_PL variables and recognized as environmental flow.

2.4 Hydraulic modelling

We used the hydraulic method for the assessment of the environmental flow of each river because of the relationship between the hydraulic parameters of watercourses (depth and velocity) and the quality of the aquatic environment (BMWP_PL–GAM relations). We used rating curves for each river describing the water depth–flow relations to obtain environmental flow for a given optimal depth. Detailed description of the applied hydraulic method of environmental flow calculation is given in Książek et al. (2019). To compare the environmental flow in relation to hydromorphological parameters (incision and redeposition), we used the proportion of environmental flow (Q_env) to mean hydraulic parameters of the minimum discharge: low low flow (LLF – the lowest low flow), mean low flow (MLF – average of the minimum annual flows), and mean annual flow (MAF – average of the annual flows). These metrics show the position of the calculated environmental flow in relation to available water volume (flow characteristics from hydrological year-to-year 1961 to 2017 observations).

2.5 Case study 2D modelling methodology

We provided the detailed CCHE2D modelling of randomly chosen incised and redeposited rivers using simple randomization procedure based on the single sequence of throwing a dice. The model is a depth-averaged two-dimensional numerical model for simulating unsteady, turbulent, free-surface flow in open channels with a moveable bed. The CCHE2D model solves depth-integrated shallow water equations for all hydraulic calculations (Wu et al., 2000; Duan et al., 2001). The CCHE2D package consists of two modules: a mesh generator (MG) and a graphical user interface (GUI). The main function of the MG is designing a complex mesh system. The mesh is generated based on the surveyed topography and/or a digital terrain model (DTM). The model was applied in two representative rivers, varying in riverbed morphology – from incised bed rock channels to a channel with natural sediment structures (with redeposition). The mesh for each sector of the river was generated by interpolating cross-sections. A total of 5112 observations were used, Raba 3033 (incision) and Ropa 2079 (redeposition). The shape of the channels was fairly regular along the reach under study, and its pattern presented little complexity (i.e. a single channel with no islands), but riffle–pool sequences were observed. The 153–200 m long meshes were composed of cells and nodes (length and number of modes, respectively, for Ropa 153 m and 49 715 and for Raba 200 m and 99 200). Data used for the initial conditions were extracted from field measurements. Special attention was devoted to bed roughness due to its importance for water surface level. Roughness values ranged from 0.01 in hydraulic smooth bed zones to 0.07 in rough areas. Finally, the model time step was defined at 0.1 or 0.25 s, depending on the model structure. The model was calibrated by comparing the measured and computed water surface levels for measured discharges in all cells and nodes (Fig. 2).

Figure 2 Comparison of calculated and measured water surface levels: the Ropa River for discharge of 6.71 m³ s⁻¹ and the Raba River for discharge of 10.29 m³ s⁻¹ (Δh – difference between measured and calculated water surface level, F – area of particular differences, percentage).

In the case of the Raba River, for 70 % of the calculated nodes, the difference between the calculated and measured water surface level (WSL) was in the range ±0.02 m. A total of 84 % of Ropa River nodes were in the range of −0.02 to 0.06 m. In all described models, Δh in the main channel does not cross ±0.02 m, but the visible differences are related to the horizontal layout of WSL in the cross-section. Evaluation of the compatibility measures of the numerical model showed very good agreement (Książek et al., 2010), and the prepared models did not need recalibration.

For each research section, we chose 20 points at each subsampled area differing in water velocity and water depth as the main environmental variables creating habitat heterogeneity for macroinvertebrates. Then, according to the relationship between hydromorphological habitat attributes (water depth and velocity) and the BMWP_PL index values (describing the ecological quality of the river), we constructed a GAM model as the best-fitted method to mark out the range of hydromorphological attributes (where the BMWP_PL suitability index obtained from the GAM model curve is the highest). Based on the optimal depth values, environmental flow was established using rating curves.

Two rivers (located in the same Carpathian region) representing opposite bed modifications (incision and redeposition) were chosen for the model as a case study. The modelled sectors of the river had channels with a pool–riffle sequence and fluvial deposits but varied in terms of degradation of the bed structure. The hydrological characteristics of the modelled river are presented in Fig. 3.

Figure 3 Changes in hydrological regime of the Raba and Ropa rivers. The horizontal line indicates the mean annual flow (MAF).

The Raba was selected to represent incised channel rivers (bottom material mainly gravel and small stones, substrate index 14.9). The Dobczyce retention reservoir, which influences the hydrology and morphology of the river, is located upstream of the examined sector of the river (12 km). Construction of the retention reservoir in 1986 led to a significant decline in average annual flow values (MAF values varied from 12.22 m³ s⁻¹ in 1951–1985 to 10.57 m³ s⁻¹ in 1986–2015, F=49.90, p<0.0001) and broke the continuity of the sediment transport. The reduction in flow, blockade of sediment supply, and longitudinal training work of the Raba led to incision of the riverbed and permanent compactness of the bed material. The Ropa River, chosen to represent the redeposition processes, was located among upland, carbonate, and silicate rivers, with the lowest human impact, i.e. agricultural land. The bottom material consists mainly of gravel and sand (substrate index 7.2), where bedload transport remains undisturbed.

We also wanted to estimate minimum flow values for two rivers which were modelled using CCHE2D. The values of depth and velocity corresponding to the highest BMWP_PL, obtained from the GAM model for each group of river and season, were plotted against the number of pixels having optimal values. Given those calculations, we were able to obtain the weighted usable area of macroinvertebrate communities (WUA) showing the most optimal habitat parameters (GAM depth and GAM velocity). WUA is often defined as an index to various ecological parameters at different organization levels: population (such as biomass, microhabitat area, size classes) (Muñoz-Mas et al., 2016) or other community level (diversity indices or ecological metrics) (Jowett, 1997, 2003; Theodoropoulos et al., 2015; Pander et al., 2019). Each pixel covered 0.25 m² of total river area, so the numbers of counted cells of the given values of velocity and depth were summarized and multiplied by the surface area. Based on these calculations using the CCHE2D model, we were able to find the relationship between usable area and flow values. To calculate the optimal environmental flow values, the curve between flow and optimal area was created. The low border of optimum of environmental flow was estimated as 50 % of WUA values (Jowett et al., 2008) for CCHE2D modelled rivers.

A hydraulic habitat 2D model of each river section was used for spring and autumn as an example to estimate habitat prediction in terms of calculated environmental flow during the season. Environmental flow that did not meet the conditions of 100 % habitat suitability for macroinvertebrates was expressed as the critical instream environmental flow value (Q_env critical), below which the parameters of aquatic macroinvertebrate communities dramatically declined.

3 Results

3.1 Environmental flow based on benthic invertebrate distribution in relation to river hydromorphology

A total of 151 466 individuals belonging to 92 benthic invertebrate families from 480 macroinvertebrate assemblages were identified. High variation was shown in the taxonomic composition of aquatic invertebrates, depending on the hydromorphological parameters (water depth and velocity) and the season (Fig. 4). In the case of rivers classified as Group 1, water velocity was found to significantly affect the taxonomic composition of the macroinvertebrates in both spring and autumn (Table 2).

Figure 4 Non-metric multidimensional scaling (NMDS) of macroinvertebrate taxonomic composition of three groups of rivers in the spring and autumn season according to velocity and depth ranges.

Table 2 Results of ANOSIM analysis comparing macroinvertebrate assemblages between classes of velocity and depth measured for the three river groups in the spring and autumn season. p values marked with bold indicate statistically significant differences.

Significance level (p with Bonferroni correction): ^∗ p<0.05. ${}^{\ast \ast }$ p<0.01. ${}^{\ast \ast \ast }$ p<0.001.

In spring, there were significant differences between velocity classes (low and high and medium and high), while in autumn, before overwintering, significant differences were only noted for medium and high classes. In neither season were the differences noted in taxonomic composition depending on the range of depth statistically significant in the case of rivers of the second abiotic group (Group 2); more significant differences were observed between velocity and depth classes (three depth classes were adopted due to the greater amplitude of these parameters). In the spring, significant differences were visible in all velocity classes, while in the case of depth they were noted only in the comparison of the low and middle depth classes. In autumn, differences were found for all classes in the case of variation in both velocity and depth. In the case of Group 3 rivers (carbonate and silicate fine sediments and rocks), the velocity parameter only taxonomically differentiated macroinvertebrate communities in the spring between the high and medium velocity classes. In the case of depth, differences were observed in both seasons – in spring between the deepest and shallowest environments and those with medium depth and in autumn only between the deepest and the shallowest zones (Table 2).

Each of the hydromorphological parameters was evaluated by the GAM model, which provided the best fit to the data (Table 3). There were significant effects of depth and velocity and its combination on variation of BMWP_PL index. Generally, the percentage of the total deviance was the highest for the combination of both hydrological parameters; however the depth parameter alone described a similar level of the total deviance. Velocity explained 38.1 % and 44.5 % of the total deviance of BMWP_PL variation in the mountain rivers (Group 1) for spring and autumn respectively. In other river groups the total deviance described for velocity varied between 6 % and 29 %. Bringing into consideration that both hydrological parameters alone described more of the total deviance, we regarded them in further analyses separately. The curves of the generalized additive models for the biotic index BMWP_PL in spring and autumn are presented in Fig. 5.

Table 3 Summary of the generalized additive models for BMWP_PL according to velocity and water depth parameters in the three river groups for spring and autumn season. Res. dev. – residual deviance, % deviance – percentage of total deviance, Res. d.f. – residual degrees of freedom, p – significance value.

Figure 5 Optimal habitat curves using generalized additive models of the BMWP_PL index for water velocity and depth in spring and autumn season for the three river groups.

These models were made for each of the three river groups: calcareous and silica bedrock alpine rivers (Group 1), sandstone mountain rivers (Group 2), and carbonate and silicate upland rivers (Group 3). In the first group, with a gravel bottom, the BMWP_PL index reached its highest values at high water velocity and in shallower zones (by the shores). In the second group of river, the BMWP_PL index was highest at medium velocities in spring and at high velocities in autumn. In both seasons, higher values for the biotic index were associated with shelf environments, as in the case of Group 1. Similar relationships with depth were noted in the Group 3 rivers, where BMWP_PL values were highest in the shallow environments at low velocity in both spring and autumn (Fig. 5).

Using the optimal depth characteristics reflecting the habitat suitability (Fig. 5), the environmental flow based on hydraulic method (rating curve) was defined. The results are shown in Table 4.

Table 4 Environmental flow and flow proportion (S) in different abiotic and bed modification types (I – incision, R – redeposition) of 12 mountainous rivers.

LLF – low low flow; MLF – mean low flow; MAF – mean annual flow; SLLF – the proportion of environmental low and the lowest low flow; SMLF – the proportion of environmental flow and the average of the minimum annual flows; SMAF – the proportion of environmental flow and the average of the annual flows.

There is a high variation of the Q_env, related to its own channel properties and volume of water. To obtain the relation to hydraulic river parameters, the mean Q_env relative similarity to MAF, MLF, and LLF was measured. There was no relation to the abiotic group of river (Table 5). The only significant relation was linked to channel modification (Fig. 6). In all cases, the relative similarity of flow was significantly higher in incised channels than redeposited ones.

Table 5 General linear modelling results for hydrological flow similarity (S) in relation to bed modification (incision and redeposition), season, and abiotic river group. SS – sum of squares, d.f. – degrees of freedom, MS – mean square. p values marked with bold indicate statistically significant relations.

LLF_sim – low low flow similarity, MLF_sim – mean low flow similarity, and MAF_sim – mean annual flow similarity.

Figure 6 The distribution of mean values ± SE (box) and whisker length (1σ) with distribution of jitter of e-flow proportion to low low flow (LLF), mean low flow (MLF), and mean annual flow (MAF) in relation to riverbed modification (I – incision, R – redeposition).

In each type of flow (MAF, MLF, and LLF), the relative similarity was higher in incised rivers than redeposited, showing that the incised rivers needed much more volume of water to sustain appropriate conditions for macroinvertebrates compared with the redeposited ones. More detailed analysis and visualization of spatial modelling were predicted by 2D modelling of randomly chosen rivers, presented below as a case study.

3.2 Case study

We calculated the detailed 2D modelling for two randomly chosen incised and redeposited rivers. According to the GAM macroinvertebrate habitat suitability model, WUA–flow curves were calculated for rivers with varying intensity of bed modification, Raba (incised) and Ropa (redeposited), as shown in Fig. 7.

Figure 7 Weighted usable area (WUA)–flow relation curves (spring and autumn season) of the rivers varying in bed modification: Raba River with incision and Ropa River with redeposition.

The environmental flow was defined as the lowest flow corresponding to 50 % of the value of the usable area, which ensures minimum optimal conditions for the development and functioning of aquatic macroinvertebrates (Jowett et al., 2008). Analysis of the curves for the Raba River shows a 50 % reduction in the usable area at the flow of about 10 m³ s⁻¹ for both spring and autumn. In the case of the Ropa River, the WUA–flow curves show a 50 % reduction in the usable area at the flow of about 2 m³ s⁻¹ in spring and 3 m³ s⁻¹ in autumn (Table 6).

Table 6 Environmental optimal and critical flow based on macroinvertebrate habitat suitability models of two mountainous rivers with mean MAF, MLF, and LLF in relation to the seasons.

A spatial visualization of macroinvertebrate habitat suitability for Q_env optimal conditions is presented in Fig. 8. In the case of the strongly incised Raba River, a very small optimal habitat area was observed, covering only the shelf zone. In the case of the Ropa River, where sediment transportation occurs, the usable areas constitutes more than 20 % of the environmental flow area. The modelling was also used to determine Q_env critical, at which the most valuable areas in terms of habitat (over 80 % suitability) disappear (Fig. 8 and Table 6). Below this Q_env critical value, a dramatic decline in macroinvertebrate diversity should be expected.

Figure 8 Probability of habitat suitability calculated as a percentage of the occurrence of optimal conditions of macroinvertebrate habitat suitability for calculated Q_env in spring and autumn seasons in incised (Raba) and redeposited (Ropa) rivers.

A comparison of the Q_env values (optimal and critical) and means, annual flow (MAF) and low flow (MLF), for the two types of rivers is presented in Table 7. In the highly incised river (Raba River), the Q_env optimal requirement for spring was lower but for autumn was higher than MAF, and Q_env critical was always higher than MLF. In the redeposited Ropa River, in spring as well as in the autumn season, Q_env optimal requirements were much lower than MAF, and MLF was higher than Q_env critical. Both findings are congruent with the former hydraulic calculations for all rivers.

4 Discussion

The present study showed that riverbed transformation, disturbing sedimentation processes and increasing the incision of the riverbed, vastly increases the environmental flow values for macroinvertebrates habitat suitability. This is important because incision processes are common in most European rivers (Gore, 1996). Channel incision decreases the area of optimal habitat for macroinvertebrates and increases the potential environmental flow to an extremely high level to obtain the minimum beneficial habitat capacity for macroinvertebrates (Bravard et al., 1997; Skalski et al., 2020). In incised channels, the degree of lateral connectivity between the river and floodplain is reduced, and the degree of modification of the substrate material is higher (Wyżga et al., 2012). As a consequence of channelization and incision, the continuity of the floodplain and shelf zone along the river is disrupted (Walther and Whiles, 2008; Kędzior et al., 2016; Anim et al., 2016; dos Reis Oliveira et al., 2019). Moreover, incision results in a concomitant decrease in sediment supply to the channels, reducing the microhabitat diversity and the quality of macroinvertebrate habitats (Wyżga, 2007; McKenzie et al., 2020). During the incision process, morphological changes in the channel, especially in the case of highly incised rivers, decrease the area of shelf habitat, and fluvial deposits are drastically reduced. Thus, to keep areas wet, flow requirements must be much higher than the mean annual flow and associated with inundation hazards.

Linkage between mean annual flow and environmental flow estimation has been the subject of consideration for many years (Tennant, 1976), based on the assumption that to obtain good stream environment conditions, some percentage of the average flow is required (Richter et al., 2012; Van Niekerk et al., 2019). According to Tennant (1976), 10 % of the average flow is the minimum flow recommended to sustain short-term survival habitat for most aquatic life forms. A proportion of 30 % was recommended as a base flow to sustain good survival biota conditions, and 60 % provides an excellent to outstanding habitat for most aquatic life forms during their primary periods of growth and for most recreational uses. However, what about strongly channelized and incised rivers, which are the most common channel types in Europe? Our survey indicated that to obtain high macroinvertebrate diversity, we need a much higher volume of water than 10 % of MAF. In the case of incision, a high volume of water is needed to cover the shelves and sediment storage, which are the principal elements of macroinvertebrate habitats and refuges in a dynamic river system (Duan et al., 2009; Anim et al., 2018).

It is obvious that macroinvertebrates are closely linked to the substrate, which is highly variable in terms of particle size (Bravard et al., 1997; Merz and Ochikubo Chan, 2005; Duan et al., 2009). Alluvial processes are strongly disturbed in an incised river, leading to deepening of the channel and bed degradation (Wyżga, 2007). The areas shown in Fig. 7, which are 100 % optimal for macroinvertebrates, are extremely narrow in incised rivers throughout the spring and autumn. In most rivers with an augmented bed, the sedimentation process is disturbed, and thus only habitats located closer to the surface, where lateral erosion occurs, provide an optimal habitat for macroinvertebrates. Modern restoration efforts often involve the artificial addition of sediments to sand (dos Reis Oliveira et al., 2019) or modification of channel morphology to restore the sedimentation process (Violin et al., 2011; Anim et al., 2018).

The biotic integrity of rivers is primarily restricted by downstream transport of sediments controlling the integrity of fluvial ecosystems (Katano et al., 2009; White et al., 2016). Substrate characteristics such as size, stability, compactness, quality, and dynamics are key parameters determining the occurrence and variation in macroinvertebrate communities. High substrate stability, substrate heterogeneity, and low compactness determine high macroinvertebrate diversity (Beisel et al., 2000; Duan et al., 2009). On the other hand, fine sediments can be regarded as a potential stressor for macroinvertebrates (Meißner et al., 2019). In highly incised sectors of the river, a deficiency of sediment and its compactness as well as a lack of food sources (Shields et al., 1994; Jowett, 2003) lead to impoverishment of the taxonomic composition of macroinvertebrates and favour taxa adapted to high flow only (Wyżga et al., 2013). Our results indicates that prevention of optimal conditions requires more volume of water which exceeds the mean annual flow. This conclusion seems paradoxical and rather dangerous because increased discharge augments incision processes. We can thus fall into a kind of ecological trap. A solution may be to pay careful attention to the bed morphology, especially in the case of incised channels. There is still a problem of gathering information on the flow ecological response of any organisms, and an extended survey in an international context should be done (Poff and Zimmermann, 2010; Fornaroli et al., 2015). We then have two options to preserve the high biodiversity of invertebrates according to the EU Water Directive: to vastly increase the water volume or to restore sedimentation processes to obtain a hydrodynamic balance. As a consequence, optimal habitats for invertebrates and fish will be enlarged. The second option seems much more realistic. Only then will we be able to successfully maintain the diversity of aquatic biota.

5 Conclusions

In habitat modelling, careful attention should be paid to the morphology of the modelled river, its geometry, and the fluvial processes in the active channel. In incised channels where sedimentation processes are altered, for example, by dam reservoirs or bedrock downcutting, the area of optimal habitat is limited. Macroinvertebrate habitat preferences are strongly linked to shelf habitats, where sediment storage and redeposition of bed material are the highest. In that case, the recolonization pattern of invertebrates requires much higher flows, even higher than the mean annual flow. As a consequence, the river is endangered by downcutting processes and impoverishment of optimal habitats.