Towards i5 Ecohydraulics: Field Determination of Manning’s Roughness Coefficient, Drag Force, and Macroinvertebrate Habitat Suitability for Various Stream Vegetation Types

Abstract

Ecohydraulic models have commonly used the flow velocity, water depth, and substrate type (i3 models) as the three fundamental determinants of the distribution of freshwater biota, but a fourth determinant has largely been neglected: stream vegetation. In this study, we provide the hydraulic and habitat information required to develop vegetation-adapted ecohydraulic models (i4 models) in streams. We calculated drag forces and Manning’s roughness coefficients (n_V) for nine types of submerged, emergent, and overhanging stream vegetation. In addition, we developed habitat suitability curves (HSCs) for benthic macroinvertebrates for these stream vegetation types. Hydraulic modules can now be upgraded to simulate stream vegetation by including the vegetation-adapted n_V values within an additive approach in which n_V is added to the n value of the inorganic substrate to which the vegetation is rooted. Habitat modules can also be upgraded to include macroinvertebrate HSCs for stream vegetation, again by adding the vegetation-adapted habitat suitability to that of the inorganic substrate to which the vegetation is rooted. In combination, i4 ecohydraulic models (including vegetation) can now be designed and applied, and we suggest that ecohydraulic research should further focus on including a fifth variable (water temperature) to ultimately advance to i5 ecohydraulic models that will optimally simulate the hydroecological reality.

1. Introduction

It has been more than 60 years since the development of the first computational fluid dynamics (CFD) model in the Los Alamos National Laboratory (1957; Los Alamos, NM, USA; [1]), which could solve the equations of fluid motion to simulate flow dynamics. Since then, CFD models have been routinely coupled with ecological/habitat models (and termed ‘ecohydraulic models’) to facilitate the management of freshwater ecosystems by, e.g., assessing the environmental flows downstream of dams [2,3], evaluating the ecological effectiveness of river restoration [4,5], or ecologically optimizing fish passage through fishways [6,7].

Ecohydraulic models have been founded on a long-proven theory that physical (hydraulic) variables (mainly flow velocity (V), water depth (D), and substrate type (S)) are major determinants of the distribution of freshwater biota [8,9]. Regardless of purpose, these models commonly include the following processes: A hydraulic model (CFD-based) calculates V and D in various water discharges across a computational mesh simulating a river reach by solving the equations of fluid motion (conservation of mass and x-wise, y-wise, and z-wise momentum, depending on the spatial resolution required). Then, a coupled habitat model combines the calculated V and D at each node of the mesh with S and compares all V-D-S combinations (node by node) with the a priori known V-D-S preferences of selected freshwater biota (mainly fish, benthic macroinvertebrates, and aquatic plants) to simulate/predict their instream distribution or to calculate the suitability of each V-D-S combination in the river reach in the various simulated discharges [10,11]. However, despite the progress that has been made, a fourth fundamental determinant of the distribution of freshwater biota has been neglected by most ecohydraulic simulations: stream vegetation [12,13].

Stream vegetation (i.e., vegetation on the bottom of a river or riparian vegetation that falls into the water) shapes freshwater communities either directly, by providing additional instream habitats for freshwater biota, thus increasing habitat and biological diversity [14,15], or indirectly, by altering local hydraulics (increasing flow resistance, mostly leading to slower and deeper waters [15,16]). Consequently, ecohydraulic models not accounting for the hydroecological influence of stream vegetation may produce less robust environmental flow recommendations [17] or less confident results on the effectiveness of river restoration [18,19,20]. To account for stream vegetation, ecohydraulic models require specific adjustments in both their hydraulic and ecological/habitat modules.

Hydraulic models/modules consider stream vegetation as a variable that increases water flow resistance [18,21], and thus hydraulic research has focused on calculating the appropriate roughness coefficients to account for vegetation-induced hydraulic effects on the flow field [22,23,24]. To this end, stream vegetation has been commonly categorized into emergent and submerged, rigid and flexible. The flow field has been subdivided into layers, and simple or complex equations have been formulated to appropriately calculate vegetation-adapted bed roughness coefficients (Manning’s hydraulic roughness coefficient (n) has been commonly used) or drag force terms for the above vegetation types [15,24,25,26,27], accounting for vegetation density, flexibility, submergence, shape, etc. [28,29]. However, despite the progress that has been made, hydraulic simulations in the presence of stream vegetation are associated with uncertainty regarding the accuracy of the hydrodynamic calculations [30] (similar to other types of hydraulic models, e.g., sea-based [31,32]) due to the lack of appropriate data for validation. Consequently, although there is a growing theoretical background for incorporating stream vegetation in hydraulic models, there is currently a lack of field data around real vegetation patches to enable the practical application of these equations in hydraulic (or ecohydraulic) models [33].

Habitat models/modules consider stream vegetation as another habitat-defining variable, along with V, D, and S. The preference of freshwater biota for V, D, and S is commonly calculated and visualized either with habitat suitability curves (HSCs) [34] or (fuzzy) habitat suitability rules [35,36] in which a zero value reflects a complete absence of aquatic organisms (totally unsuitable habitat) and a value of one represents the maximum presence (totally suitable habitat). However, although there are several V-, D-, and S-based suitability curves available worldwide (and thus, V-D-S-based ecohydraulic models can be easily formulated), as far as we are aware, there is a lack of vegetation-based HSCs, which, however, are a prerequisite to practically apply vegetation-adapted ecohydraulic models.

The purpose of this study was to provide the hydraulic and habitat information required to develop and apply vegetation-adapted ecohydraulic models in streams. Using hydraulic and biological field measurements from two stream reaches in Greece, we calculated drag forces and Manning’s roughness coefficients for various types of submerged, emergent, and overhanging stream vegetation. We also developed HSCs for benthic macroinvertebrates for these stream vegetation types. Ecohydraulic scientists and practitioners can now use the calculated vegetation-adapted Manning’s roughness coefficients to calibrate and validate a robust vegetation-including hydraulic module and can use the vegetation-based HSCs to properly include stream vegetation as a fourth habitat-defining variable, thus increasing the accuracy of the habitat module and ultimately producing ecohydraulic models that optimally simulate hydroecological reality.

2. Materials and Methods

2.1. Presampling

Prior to sampling, stream vegetation was categorized into nine types: (1) algae, (2) bryophytes, (3) emergent fine-leaved plants (EFL), (4) emergent broad-leaved plants (EBL), (5) submerged fine-leaved plants (SFL), (6) submerged broad-leaved plants (SBL), (7) overhanging (tree trunks, bush branches, and stems that fall from the river bank to the water surface), (8) free-floating vegetation (e.g., water lilies), and (9) detritus (snags, dead branches, leaves, and roots).

2.2. Sampling

We sampled 124 vegetated microhabitats (rectangular-shaped; 0.0625 m²) of the above types from two stream reaches in the region of Attica, Greece: (i) the Oinoi Stream and (ii) the upper Great Rafina River (Lykorema Stream) (Figure 1). A total of 74 microhabitats were sampled in spring 2021 (high-flow period) and 50 microhabitats were sampled in summer 2021 (low-flow period) within a stratified sampling framework to evenly allocate vegetation types among the samples.

At each microhabitat/sample, we (1) recorded the type and density (percent cover) of vegetation (visual estimation; photos were taken to further justify our selection), (2) measured the flow velocity and water depth before, after, and on each side of the microhabitat (using the Swoffer 2100 current velocity meter and a water depth measurement rod attached to the velocity meter), and (3) collected a sample of benthic macroinvertebrates (Figure 2). Macroinvertebrates were sampled using a 0.25 × 0.25 = 0.0625 m² Surber sampler (or a hand net of the same dimensions) with a 500 μm mesh size. Samples were preserved in plastic bottles containing 70% ethanol. In the laboratory, all specimens were sorted and identified to the family level using a stereo microscope (magnification 6.39–639) and macroinvertebrate identification guides for the Mediterranean region [37,38,39].

2.3. Postsampling Analysis

2.3.1. Manning’s Roughness Coefficients and Drag Force

Vegetation-adapted Manning’s roughness coefficients (n_v) and drag force terms (F_dx; F_dy; and F_D) were calculated for each vegetation type (excluding algae, free-floating vegetation, and detritus) for three density (d) classes: (1) sparse—0% < d ≤ 50%, (2) moderate—50% < d ≤ 80%, and (3) dense—80% < d ≤ 100%. Manning’s n for emergent vegetation (EFL, EBL, and overhanging) was calculated using the following equation by Petryk and Bosmajian [25] for emergent rigid vegetation:

$${\mathrm{n}}_{\mathrm{v}}={\mathrm{n}}_{\mathrm{o}}\sqrt{1+\left(\frac{{\mathrm{C}}_{\mathrm{d}}\sum {\mathrm{A}}_{\mathrm{i}}}{2\mathrm{gAL}}\right){\left(\frac{1}{{\mathrm{n}}_{\mathrm{o}}}\right)}^2{\mathrm{R}}^{\frac{4}{3}}}$$

(1)

where ${\mathrm{n}}_{\mathrm{v}}$ is the Manning’s n for each emergent vegetation type, ${\mathrm{n}}_{\mathrm{o}}$ is the boundary n (the Manning’s coefficient for the substrate without vegetation, which was set to 0.03 for clean natural streams based on the previous literature [40]), ${\mathrm{C}}_{\mathrm{d}}$ is the drag coefficient (set to 1 for vegetated channels [41]), ${\mathrm{A}}_{\mathrm{i}}$ is the cross-sectional flow area (m²), g is the acceleration of gravity (m/s²), and R is the hydraulic radius (m).

Manning’s n for submerged vegetation (SFL, SBL, and bryophytes) was calculated using the following equation by Freeman [42] for flexible submerged vegetation:

$${\mathrm{n}}_{\mathrm{v}}={\mathrm{K}}_{\mathrm{n}}0.183{\left(\frac{{\mathrm{E}}_{\mathrm{s}}{\mathrm{A}}_{\mathrm{s}}}{{\mathsf{\rho }\mathrm{A}}_{\mathrm{i}}{\mathrm{V}}_{*}^2}\right)}^{0.183}{\left(\frac{\mathrm{H}}{{\mathsf{{\rm Y}}}_{\mathsf{o}}}\right)}^{0.243}{\left({\mathsf{{\rm M}}\mathrm{A}}_{\mathrm{i}}\right)}^{0.273}{\left(\frac{\mathsf{\nu }}{{\mathrm{V}}_{*}\mathrm{R}}\right)}^{0.115}\left(\frac{1}{{\mathrm{V}}_{*}}\right){\mathrm{R}}^{\frac{2}{3}}{\mathrm{S}}^{\frac{1}{2}}$$

(2)

where ${\mathrm{n}}_{\mathrm{v}}$ is the Manning’s n for each submerged vegetation type, ${\mathrm{K}}_{\mathrm{n}}$ is a unit conversion factor (1.0 m^1/3/s), ${\mathrm{E}}_{\mathrm{s}}$ is the modulus of plant stiffness (Nm⁻²), ${\mathrm{A}}_{\mathrm{s}}$ is the total area covered by vegetation (m²), ρ is the fluid density (kg/m³), ${\mathrm{A}}_{\mathrm{i}}$ is the cross-sectional flow area (m²), V is the mean flow velocity (m/s), H is the average undeflected plant height (m), ${\mathsf{{\rm Y}}}_{\mathsf{o}}$ is the water depth (m), M is the relative plant density, V* is the shear velocity (m/s), R is the hydraulic radius, and S is the bed slope.

To facilitate the use of the vegetation-adapted coefficients in ecohydraulic modeling, we further averaged the coefficients of EFL, EBL, and overhanging to derive a single n for emergent vegetation and the coefficients of SFL, SBL, and bryophytes to derive a single n for submerged vegetation.

Additionally, the drag force terms (i.e., F_x and F_y; N/m²) were estimated according to the work of Sun et al. [15] as follows:

$${\mathrm{F}}_{\mathrm{x}}={\mathrm{S}}_{\mathrm{x}}+{\mathrm{F}}_{\mathrm{dx}} |{ \mathrm{F}}_{\mathrm{y}}={\mathrm{S}}_{\mathrm{y}}+{\mathrm{F}}_{\mathrm{dy}}$$

(3)

where ${\mathrm{S}}_{\mathrm{x}}$ and ${\mathrm{S}}_{\mathrm{y}}$ express the bed friction terms derived from the Manning–Strickler formula, while ${\mathrm{F}}_{\mathrm{dx}}$ and ${\mathrm{F}}_{\mathrm{dy}}$ are the drag force terms resolved in terms of the Chézy’s coefficient (x: longitudinal and y: transverse axis). The drag forces (i.e., F_dx and F_dy; N/m²) were estimated according to the work of Sun et al. [15] as follows:

$${\mathrm{S}}_{\mathrm{x}}=-\frac{\mathrm{a}}{\cos \left(\mathsf{\theta }\right)}\frac{{\mathrm{gn}}^2}{{\mathrm{h}}^{\frac{4}{3}}}\mathrm{U}\sqrt{{\mathrm{U}}^2+{\mathrm{V}}^2}|{\mathrm{S}}_{\mathrm{y}}=-\frac{\mathrm{a}}{\cos \left(\mathsf{\theta }\right)}\frac{{\mathrm{gn}}^2}{{\mathrm{h}}^{\frac{4}{3}}}\mathrm{V}\sqrt{{\mathrm{U}}^2+{\mathrm{V}}^2}$$

(4)

$${\mathrm{F}}_{\mathrm{dx}}=-\frac{\mathrm{N}}{\mathrm{A}}\frac{\mathrm{D}}{2}{\mathrm{C}}_{\mathrm{d}}\mathrm{U}\sqrt{{\mathrm{U}}^2+{\mathrm{V}}^2} |{ \mathrm{F}}_{\mathrm{dy}}=-\frac{\mathrm{N}}{\mathrm{A}}\frac{\mathrm{D}}{2}{\mathrm{C}}_{\mathrm{d}}\mathrm{V}\sqrt{{\mathrm{U}}^2+{\mathrm{V}}^2}$$

(5)

where α denotes the vegetation porosity (%), θ is the channel slope, g is the acceleration of gravity (m/s²), n is Manning’s coefficient, h is the water depth between the bottom and the free surface (m), U and V are the longitudinal and lateral velocities (m/s), N is the total vegetation number per area A (m²), and D is the vegetation diameter (m) [43]. Subsequently, the divergence force (F_D; N/m²) of the two force components was calculated as:

$${\mathrm{F}}_{\mathrm{D}}=\sqrt{{\mathrm{F}}_{\mathrm{x}}{}^2+{\mathrm{F}}_{\mathrm{y}}{}^2}$$

(6)

2.3.2. Habitat Suitability

To develop habitat suitability curves (HSCs) for benthic macroinvertebrates, we applied a previously formulated macroinvertebrate-community-based habitat suitability index [2,36] that integrates four metrics that are commonly applied to assess macroinvertebrate habitat suitability [44,45,46,47]. These metrics are macroinvertebrate abundance, taxonomic richness, diversity (Shannon’s index), and Ephemeroptera–Plecoptera–Trichoptera richness (EPT richness), and they were calculated using the ASTERICS software version 3.1.1 (IRV Software, Vienna, AT, USA). Habitat suitability was calculated as follows:

$$\mathrm{K}=0.4{\mathrm{n}}_{\mathrm{m}}+0.3\mathrm{H}+0.2\mathrm{EPT}+0.1\mathrm{a}$$

(7)

where K represents the habitat suitability ranging from 0 (lowest) to 1 (highest), n_m denotes macroinvertebrate taxa richness (No. of families), H is the Shannon’s diversity index, EPT is the number of EPT taxa, and a denotes the abundance of benthic macroinvertebrates.

Each microhabitat’s K was normalized to a 0–1 scale by dividing by the maximum K value observed at the site of the sampled microhabitat. All K values (for each microhabitat) were further standardized by multiplying by vegetation density, and the overall K for each vegetation type was expressed as the average microhabitat K. An abundance-based habitat suitability was also calculated for each vegetation type and for each macroinvertebrate taxon by summing the abundance of each taxon among the samples of the same vegetation type. Thus, two habitat suitability indices were ultimately developed and applied:

$${\mathrm{K}}_{\mathrm{VC}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\frac{{\mathrm{K}}_{\mathrm{i}}{\mathrm{d}}_{\mathrm{i}}}{\mathrm{N}}$$

(8)

where ${\mathrm{K}}_{\mathrm{VC}}$ is the community-based habitat suitability for each vegetation type (ranging from 0 to 1), N is the number of samples of each vegetation type, ${\mathrm{K}}_{\mathrm{i}}$ is the habitat suitability of each sample, and di is the density of the vegetation at each sample (microhabitat) and:

$${\mathrm{K}}_{\mathrm{Va}}={{\displaystyle \sum }}_1^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}}$$

(9)

where ${\mathrm{K}}_{\mathrm{Va}}$ is the taxon-based habitat suitability for each vegetation type (ranging from 0 to 1), N is the number of samples of each vegetation type, and ${\mathrm{a}}_{\mathrm{i}}$ is the abundance of the taxon at each sample.

The calculated habitat suitability for each vegetation type was added to the habitat suitability of the inorganic substrate to which the vegetation was rooted based on a previous habitat suitability assessment for inorganic substrate types [36]. To ensure that no issues of spatiotemporal autocorrelation existed among the samples (for the calculated habitat suitability indices), we applied Moran’s autocorrelation coefficient (Moran’s I) in R version 4.0.3 [48] using the ape package [49].

3. Results

3.1. Vegetation-Adapted Manning’s Roughness Coefficients and Drag Force

In both streams that were examined in this study, the average flow velocity (V) was equal to 0.22 m/s (number of samples (N) = 124; standard error (SE) = ± 0.014; min = 0.033; max = 2.038), and the average water depth was 0.12 m (N = 124; SE = ± 0.004; min = 0.01; max = 0.64). All vegetation types had lower average V downstream of the vegetation patch and higher average V upstream of the vegetation patch (mean difference of 0.093 m/s). The highest V change (Vc) was observed for detritus (N = 23; Vc = −0.18 m/s), and the lowest Vc was observed for overhanging vegetation (N = 16; Vc = −0.03 m/s) (Figure 3). In contrast, the average water depths (D) were similar or slightly higher after the vegetation patches (mean difference of 0.00013 m). The highest D change (Dc) was observed for algae (N = 18; Dc = +0.025 m), and the lowest Dc was observed for emergent broad-leaved vegetation (N = 10; Dc = +0.001 m) (Figure 4).

Manning’s n values were higher for all vegetation types compared to n_o, gradually increasing from sparse to dense patches for both emergent and submerged vegetation. Manning’s n change (compared to n_o) was +31.70% for emergent vegetation types (N = 34) and +30.32% for submerged vegetation types (N = 45). The highest n change values were observed for dense emergent vegetation (N = 12; n change: 57.07%) and dense submerged vegetation (N = 22; n change = +48.01%). The lowest n change values were observed for sparse submerged vegetation (N = 16; n change = +23.5%) and sparse emergent vegetation (N = 17; n change = +16.49%) (

Table 1. Manning’s roughness coefficients (n) for emergent and submerged vegetation types in varying densities. Numbers in parentheses indicate the number of samples of each vegetation type.

Vegetation-Adapted Manning’s Roughness Coefficients (n)
	Veg. Type|Veg. Density	Dense (80% < d ≤ 100%)	Moderate (50% < d ≤ 80%)	Sparse (0% < d ≤ 50%)	Average
	Algae	n/a	n/a	n/a	n/a
Group 1: Emergent vegetation	Emergent fine-leaved	0.042 (6)	0.042 (6)	0.034 (6)	0.038 (12)
	Emergent broad-leaved	0.052 (1)	0.048 (2)	0.034 (4)	0.039 (6)
	Overhanging	0.047 (5)	0.045 (9)	0.036 (7)	0.041 (16)
	Group 1 average	0.047 (12)	0.045 (17)	0.035 (17)	0.039 (34)
Group 2: Submerged vegetation	Submerged fine-leaved	0.052 (7)	0.049 (9)	0.046 (6)	0.048 (15)
	Submerged broad-leaved	0.049 (1)	0.049 (1)	0.038 (4)	0.040 (5)
	Bryophytes	0.032 (5)	0.031 (6)	0.028 (4)	0.029 (10)
	Group 2 average	0.044 (13)	0.043 (16)	0.037 (14)	0.039 (30)
	Free-floating	n/a	n/a	n/a	n/a
	Detritus	n/a	n/a	n/a	n/a

).

Vegetation-adapted drag force (F_D) was highest in dense vegetation patches, lower in moderate patches, and lowest in sparse vegetation patches for both emergent and submerged vegetation (

Table 2. Drag force (F_D; N/m²) for emergent and submerged vegetation types in varying densities. Numbers in parentheses indicate the number of samples of each vegetation type.

Vegetation-Adapted Drag Force (FD)
	Veg. Type | Veg. Density	Dense (80% < d ≤ 100%)	Moderate (50% < d ≤ 80%)	Sparse (0% < d ≤ 50%)	Average
	Algae	n/a	n/a	n/a	n/a
Group 1: Emergent vegetation	Emergent fine-leaved	0.707 (6)	0.707 (6)	0.412 (6)	0.559 (12)
	Emergent broad-leaved	n/a	0.144 (2)	0.141 (4)	0.142 (6)
	Overhanging	0.252 (5)	0.252 (9)	0.352 (7)	0.302 (16)
	Group 1 average	0.478 (11)	0.367 (17)	0.301 (17)	0.334 (34)
Group 2: Submerged vegetation	Submerged fine-leaved	0.783 (7)	0.760 (9)	0.251 (6)	0.505 (15)
	Submerged broad-leaved	0.016 (1)	0.016 (1)	0.045 (4)	0.030 (5)
	Bryophytes	0.592 (5)	0.339 (6)	0.385 (4)	0.029 (10)
	Group 2 average	0.463 (13)	0.371 (16)	0.227 (14)	0.362 (30)
	Free-floating	n/a	n/a	n/a	n/a
	Detritus	n/a	n/a	n/a	n/a

). In more detail, F_D was highest in dense SFL (F_D = 0.783; N = 11), followed by dense EFL (F_D = 0.707; N = 6) and dense bryophytes (F_D = 0.592; N = 5), and was lowest in moderate SFL (F_D = 0.016; N = 1), followed by sparse EBL (F_D = 0.141; N = 4) and moderate overhanging vegetation (F_D = 0.252; N = 9).

3.2. Macroinvertebrate Habitat Suitability of Various Stream Vegetation Types

In total, 17,248 macroinvertebrate individuals from 72 families were sorted and identified from all stream vegetation types. The taxonomic richness was highest in free-floating vegetation (N = 2; R = 10.2; SE = ± 0.26) and lowest in overhanging vegetation (N = 16; R = 6.44; SE = ± 0.16) (Figure 5). Shannon’s diversity varied slightly between vegetation types and was highest in free-floating vegetation (N = 2; H = 1.89; SE = ± 1.76) and lowest in emergent broad-leaved vegetation (N = 10; H = 1.12; SE = ± 0.36). The EPT richness was highest in algae (N = 18; EPT = 3.1; SE = ± 0.06) and submerged fine-leaved vegetation (N = 18; EPT = 3.1; SE = ± 0.11) and lowest in free-floating vegetation (N = 2; EPT = 1.5; SE = ± 0.35).

The macroinvertebrate abundance was highest in free-floating vegetation (N = 2; a = 241.5; SE = ± 20) and lowest in overhanging vegetation (N = 16; a = 28.94; SE = ± 1.6). The macroinvertebrate-community-based habitat suitability was highest for submerged broad-leaved vegetation (N = 6; K = 0.68; SE = ± 0.06) and lowest for bryophytes (N = 13; K = 0.2; SE = ± 0.01) (Figure 6). Emergent and submerged plants had higher habitat suitability compared to nonvascular vegetation (algae and bryophytes) and detritus.

Macroinvertebrate taxa had varying optimal preferences for different vegetation types. The abundance of Athericidae, Hydroptilidae, and Tipulidae was highest in algae, with 42%, 72%, and 50% of individuals, respectively, being found in algae (Figure S1; Supplementary material). The abundance of Elmidae (68%), Empididae (46%), Ephemeridae (100%), Glossiphoniidae (100%), Lepidostomatidae (92%), Leptoceridae (59%), Leptophlebiidae (100%), Nemouridae (100%), Rhyacophilidae (100%), and Scirtidae (80%) was highest in bryophytes. Calamoceratidae (100%), Dolichopodidae (100%), Dryopidae (57%), Leuctridae (61%), Ostracoda (100%), Physidae (65%), and Polycentropodidae (82%) were more abundant in detritus. The abundance of Calopterygidae (44%), Gyrinidae (80%), Hydrobiidae (68%), Neritidae (70%), and Sphaeriidae (84%) was highest in emergent broad-leaved plants. Ancylidae (60%), Corixidae (86%), Hydrachnidae (60%), and Psychodidae (100%) were more abundant in emergent fine-leaved plants. Hydrometridae (62.5), Libellulidae (43%), and Stratiomyiidae (67%) were more abundant in free-floating plants, and Aeshnidae (55%), Caenidae (53%), and Limoniidae (77%) were more abundant in submerged fine-leaved plants.

The habitat suitability values of each vegetation type (added to the inorganic substrate on which the vegetation was rooted) are shown in

Table 3. Habitat suitability (0: totally unsuitable; 1: totally suitable) for emergent and submerged vegetation types in combination with the inorganic bottom substrate type. Unvegetated means the habitat suitability of the inorganic substrate type without any vegetation. For each vegetation type, the habitat suitability was added to that of the inorganic substrate, e.g., SFL in silt (0.71) means that the habitat suitability for SFL is 0.71 − 0.30 = 0.41.

Vegetation-Adapted Macroinvertebrate Habit Suitability
	Veg. Type	Unvegetated	Algae	Bryophytes	SFL	SBL	EFL	EBL	Overhanging	Free-Floating	Detritus
Substrate
Silt		0.30	0.64	0.50	0.71	0.98	0.65	0.65	0.53	0.52	0.55
Sand		0.45	0.79	0.64	0.85	1.00	0.80	0.80	0.67	0.66	0.70
Fine gravel		0.37	0.71	0.57	0.78	1.00	0.72	0.72	0.60	0.59	0.62
Medium gravel		0.54	0.88	0.74	0.95	1.00	0.90	0.90	0.77	0.76	0.80
Large gravel		0.66	1.00	0.85	1.00	1.00	1.00	1.00	0.88	0.87	0.91
Small stones		0.65	0.99	0.84	1.00	1.00	1.00	1.00	0.88	0.87	0.90
Large stones		0.62	0.96	0.82	1.00	1.00	0.98	0.98	0.85	0.84	0.88
Boulders		0.63	0.97	0.82	1.00	1.00	0.98	0.98	0.86	0.85	0.88

. Considering that the habitat suitability was highest in SBL, all SBL added values yielded high suitability, reaching close to 1. Vegetated silt and sand had the lowest values, depending on the vegetation type, but they were highly increased compared to the unvegetated inorganic substrate. The habitat suitability was highest for vegetated small stones, large stones, and boulders.

4. Discussion

4.1. Parsimony and the Importance of Including Stream Vegetation in Ecohydraulic Models

A model is a simplification/approximation of a complex reality. Thus, a computer model that will simulate reality in its full complexity is elusive [50]. This is also the case for hydroecological models; these models can be useful when they include as many physical-hydrological-hydraulic variables as necessary to explain ecological processes, but they may become practically useless when more and more variables are included, especially when their explanatory power is low. Consequently, the challenge in hydroecological modeling is to formulate a parsimonious model that will accurately approximate reality with as few predictor variables as possible [51]. For ecohydraulic models, these ‘few predictors’ that primarily determine the distribution of freshwater biota had long been (i) flow velocity, (ii) water depth (both reflecting the microhabitat effect of water discharge fluctuations), and (iii) the type of substrate (boulders, pebbles, cobbles, etc.) [2,35,52]. Are these three variables enough to assume a parsimonious ecohydraulic model?

Previous research suggests that two more physical variables are major determinants of hydroecological processes in freshwaters: stream vegetation and water temperature. In this study, we focused on stream vegetation, the presence of which modifies flow and sediment transport patterns [12,16] and increases habitat and biological diversity [14,15]. Thus, V-D-S-based ecohydraulic models (i3 models) may oversimplify hydroecological processes; we argue that they are useful as a start and can indeed provide confident predictions [36] but should be enhanced to include stream vegetation and water temperature (i5 models) to increase their accuracy and optimally approximate hydroecological reality in a way that can be characterized as parsimonious. While the data to develop i3 ecohydraulic models are becoming more and more available worldwide, in this study we further collected the field data necessary (hydraulic and habitat) to enable advancing to i4 ecohydraulic models (including vegetation) and further on, when water temperature data become widely available, to i5 ecohydraulic models.

4.2. How to Include Stream Vegetation in the Hydraulic Module

What mostly matters in the hydraulic modules of ecohydraulic models is to accurately and cost-effectively calibrate and validate them [53,54]. At least two sets of V and D field observations (i.e., in two different discharges) are required, one for calibration and one for validation. Calibration is commonly applied by adjusting the Manning’s n values across the river reach (in different sections) until, for the first discharge dataset, there is a match (often > 90% is required) between the simulated and observed V and D values [55]. Then, the calibrated model is used to simulate V and D values for the second discharge dataset (no further n adjustment is applied), and the model is considered validated if there is another a > 90% match between the simulated and observed V and D data. In vegetated streams, this process may become very challenging, especially if the user starts calibrating with arbitrary n values. We suggest that prior to calibrating/validating a vegetation-including model the user should have at least a rough estimation of the vegetation distribution across the stream reach (either using photographic material or satellite imagery). Then, the user could apply our vegetation-adapted n values (n_v) to start calibrating. It is possible that the additive approach of Cowan [56] should be applied. In this approach, the n that should be used is the adjusted sum of several coefficients [18,21], as follows:

$$\mathrm{n}=\left({\mathrm{n}}_{\mathrm{o}}+{\mathrm{n}}_1+{\mathrm{n}}_2+{\mathrm{n}}_3+{\mathrm{n}}_4\right){\mathrm{m}}_5$$

(10)

where ${\mathrm{n}}_{\mathrm{o}}$ is a basic n value for a straight, uniform, smooth channel; ${\mathrm{n}}_1$ is the adjustment factor for the effect of surface irregularity; ${\mathrm{n}}_2$ is the adjustment factor for the effect of variation in the shape and size of the channel cross-section; ${\mathrm{n}}_3$ is the adjustment factor for obstruction; ${\mathrm{n}}_4$ is the adjustment factor for vegetation (our n_v); and ${\mathrm{m}}_5$ is a correction factor for meandering channels. For example, in the case of emergent broad-leaved plants rooted in boulders, the user would likely add n_o = 0.07 + n₄ = 0.039 ⇔ n ≈ 0.11 or even apply the m₅ factor if there such information was available.

4.3. How to Include Stream Vegetation in the Habitat Module

Stream vegetation (B) in habitat modules/models can be included by expanding the two options for using habitat suitability: habitat suitability curves (HSCs) and fuzzy habitat suitability rules (HSRs). Within the HSC approach, at each node of the computational mesh representing the stream reach, V-suitability is often multiplied by D-suitability and S-suitability to calculate the overall node-specific suitability (or a geometric mean is often calculated). Macroinvertebrate B-suitability (based on Table 3) can thus be included by multiplying the V-, D-, and S-suitabilities with the B-suitability or by calculating a geometric mean of the four variables. If the HSR approach is to be applied, an advanced fuzzy rule-based algorithm should be formulated, which also includes vegetation-specific habitat suitability rules. For example, IF V is moderate AND D is shallow AND S is boulders AND B is EBL THEN habitat suitability is high (and so on) [36,57].

4.4. i5 Ecohydraulics: Implications and Future Research

Typically, stream vegetation has been introduced in (eco)hydraulic models by appropriately modifying various bed roughness coefficients [58]. These roughness coefficients have commonly been calculated based on equations including flow velocity, water depth, the type of substrate, and the density of the vegetation patch [40,59,60,61]. Morphological features of stream vegetation were further included in more advanced equations [28,42]. Our results showed that the effect of stream vegetation on Manning’s n was type-dependent (varying between emergent rigid and submerged flexible vegetation) and increased by increasing the vegetation density. Consequently, future studies may focus on this morphological variability of stream vegetation, potentially developing indices, such as the leaf area index [30], that would more accurately introduce the morphological variability of stream vegetation into ecohydraulic models.

As discussed earlier, i4 models (those that simulate the hydroecological effects of stream vegetation) may be considered an intermediate step to ultimately advance to the optimal (parsimonious) i5 ecohydraulic models (those that further simulate the hydroecological effects of water temperature), considering that all five physical/hydraulic variables are critical determinants of the distribution of freshwater biota. However, i5 models may require increased computational power, especially their habitat modules. The HSC approach is less demanding; stream vegetation only needs to be included in the equation that computes the overall habitat suitability (either a product or a geometric mean of V, D, S, and B habitat suitability). The HSRs approach, which takes into account the interactions between V, D, S, and B is closer to how natural/ecological processes operate and may dramatically increase computational times. For example, an HSR-based model including five V classes, five D classes, eight S classes, and our nine-class stream vegetation would require 1800 rules to operate. If another five classes of water temperature are included, the number becomes 9000 rules that need to be checked at each node. Still, the rapid increase in computational power will probably soon tackle this issue. We suggest that future research should focus on collecting field hydroecological data (i) from stream/river reaches around the world to compare and enable the application of localized i4 ecohydraulic models and (ii) on the influence of water temperature on freshwater biota that could be used in i5 ecohydraulic models.

5. Conclusions

Current ecohydraulic applications commonly use three predictor variables (i3 models) to simulate the habitat suitability of aquatic biota in various discharges. These i3 models have often resulted in accurate simulations, mostly depending on the purpose of each application, but they may also oversimplify the hydroecological reality (in real life, aquatic biota respond to a multitude of interacting environmental variables). In this study, we argue that i3 models should not be considered parsimonious and suggest that ecohydraulic research should focus on developing advanced ecohydraulic models (both conceptual models and the relevant computer software) that will include two additional predictors: stream vegetation (i4 models) and water temperature (i5 models).

We acknowledge, however, that the frequent lack of field data may inhibit such efforts towards i5 ecohydraulics. To this end, we hereby provide field-calculated Manning’s roughness coefficients for nine types of submerged, emergent, and overhanging stream vegetation, together with guidance on how to include these data in the hydraulic and habitat modules of i4 ecohydraulic models. Although locally restricted, the data of this study may be appropriately calibrated/adjusted for use in ecohydraulic applications in river reaches with similar hydroecological properties. We further argue that parsimony in ecohydraulic models should be sought in i5 models and suggest that future research should focus on collecting the data necessary to develop such models.