Spatiotemporal Variation in Benthic-Invertebrates-Based Physical Habitat Modelling : Can We Use Generic Instead of Local and Season-Specific Habitat Suitability Criteria ?

Generic habitat suitability criteria (HC) are often developed from spatially and temporally variable hydroecological datasets to increase generality, cost-effectiveness, and time-efficiency of habitat models. For benthic macroinvertebrates (BMIs), however, there is no prior knowledge on the spatiotemporal variation in their habitat preferences and how this may be reflected in the final environmental flow (e-flow) predictions. In this study, we used a large, spatiotemporally variable BMI-hydroecological dataset and developed generic, local, and season-specific subsets of HC for three seasons and two river types within various data pre-treatment options. Each subset was used to train a fuzzy habitat model, predict the habitat suitability in two hydrodynamically-simulated river reaches, and develop/compare model-based e-flow scenarios. We found that BMIs shift their habitat preferences among seasons and river types; consequently, spatiotemporally variable e-flow predictions were developed, with the seasonal variation being greater than the typological one. Within this variation, however, we found that with proper data pre-treatment, the minimum-acceptable e-flows from the generic models mostly (65–90%) lay within the acceptable e-flows predicted by the local and season-specific models. We conclude that, within specific limitations, generic BMI-HC can be used for geographically extended, cost-effective e-flow assessments, compensating for the within-limits loss of predictive accuracy.


Introduction
Spatiotemporal variation has long been a subject of discussion in physical habitat modelling [1][2][3][4][5][6][7].From a management-oriented perspective, the use of multiple local and season-specific habitat suitability criteria to develop model-based environmental flow recommendations is considered unrealistic [2,3]; usually, the costs and time required to generate site-and season-specific habitat preferences cannot be covered by the limited funding available and, thus, the ability to develop appropriate local hydroecological datasets is also limited [3].Research, however, suggests that aquatic organisms often shift their habitat preferences among seasons and may have different requirements in different geographical locations.This has been repeatedly confirmed for fish, as there is a complete agreement between studies on different species and from different regions [1,[8][9][10][11], and thus, in fish-based studies, site-and season-specific habitat suitability criteria have been preferred over generic ones [12,13].In contrast, the fewer available studies on benthic macroinvertebrates (BMIs) reach contrasting conclusions; although most studies suggest increased performance of habitat models when local and season-specific criteria are applied [4,5,11,14], others conclude that the observed differences are not significant enough to inhibit the application of generic habitat suitability criteria [2,3,15].However as Kelly et al. indicate [14], 'some balance needs to be found between generality and specificity' to widen the applicability of model-based environmental flow assessments.
Typically, habitat preference datasets consist of microhabitat samples relating the flow velocity (V), the water depth (D) and the type of substrate (S) with a habitat suitability value, calculated either by using the abundance of target aquatic organisms or by applying traits-or metrics-based approaches.These preferences are often visualized in the two-dimensional space as habitat suitability curves [4,16,17] but other approaches have been also implemented [18][19][20].The habitat preference datasets are then used as training data in hydraulic habitat models, that is, the reference data which will be used by the model to predict the habitat suitability (K) in samples/microhabitats with known V, D, and S values and unknown K. Currently, there are various alternatives available for the development of habitat suitability criteria, reflecting the aforementioned challenging effort to balance the sources of error and variation towards cost-effectiveness and time-efficiency.In this effort, hydroecological data from geographically and hydrologically various river types and in different time periods, often collected within different projects, are either treated separately to develop site-and season-specific criteria or are aggregated to increase sample size and/or extend the geographical and typological applicability of a model [21][22][23].In BMI-based studies, typical aggregation schemes include spatiotemporal pooling-of samples from different sites and/or seasons-without pre-treatment [17,24]; spatial aggregation of samples only from specific, usually low-flow periods [25,26]; and spatiotemporal pooling of samples after proper pre-treatment [4,19,27].
Although it has been previously acknowledged that spatial and temporal variation may result in error and uncertainty in model-based environmental flow predictions [28,29], case studies including and quantifying this spatiotemporal component are currently missing.Do benthic macroinvertebrates shift their habitat preferences among seasons and in different geographical locations?Will a model-based environmental flow assessment reach similar prediction by applying local and season-specific-versus generic habitat suitability criteria?Are the possible differences in the predictions variable enough to inhibit the application of generic habitat suitability criteria?Can a generic dataset be treated appropriately to enable its use in multiple locations?These questions have not been properly addressed due to the lack of relevant comparative assessments, but within the context of 'trading-off between modelling accuracy and generality' [30], such assessments are necessary to properly guide future habitat modelling efforts.
The purpose of this study was (i) to explore the possible seasonal and typological variation in the habitat preferences of benthic macroinvertebrates and (ii) to investigate whether this spatiotemporal variation can be-at least partially-omitted to facilitate the development of generic BMI habitat suitability criteria.Based on a large BMI-hydroecological dataset collected from Greek streams and rivers, we developed local and generic habitat suitability criteria for benthic macroinvertebrates, including (i) three seasons, (ii) two river types, and (iii) five K-calculation alternatives.With the application of a two-dimensional hydrodynamic habitat model, we used these criteria to develop model-based environmental flow predictions in two river reaches in Greece.Our ultimate goal was to calculate the probability of agreement between the environmental flow predictions based on the different habitat suitability criteria.Low probability would suggest that local and season-specific habitat suitability criteria cannot be replaced by generic ones due to the high loss in the model's predictive accuracy.High probability of agreement would suggest that generic habitat suitability criteria can be used in model-based EFAs within the concept of balancing between a model's predictive performance, generality, cost-effectiveness, and time-efficiency.

Overview of the Analysis
The steps followed to implement our analysis are outlined below: 1.
Development-acquisition of a benchmark microhabitat preference dataset.

2.
Calculation of habitat suitability (κ) for each microhabitat sample using BMI metrics.

4.
Training of a habitat model using each of the five κ-normalization options within three seasons and two river types.

5.
Application of hydrodynamic simulations in two river reaches to acquire V and D values in multiple discharges.6.
Prediction of κ in the two reaches using the various training alternatives (subsets).7.
Development of spatially and temporally varying environmental flow scenarios, comparisons and discussion on the selection of the minimum acceptable and optimal environmental flows within the various subsets.

Calculation of Habitat Suitability Using BMI Metrics
As detailed in Theodoropoulos et al. [18], for the calculation of habitat suitability, each of the aforementioned metrics was weighted based on a combination of expert-judgment and previous literature [32][33][34][35][36] to reflect its relevant contribution-significance, and the final, unnormalized κ has been calculated as follows: where κ is the unnormalized habitat suitability of the ith microhabitat; n i denotes the number of BMI taxa (families) found at the ith microhabitat; H i denotes the Shannon's diversity index for the ith microhabitat; EPT i is the number of EPT taxa found at the ith microhabitat; α i is the abundance of benthic macroinvertebrates found at the ith microhabitat.
As the unnormalized κ values arbitrarily ranged from 0 to higher than 3000, and in accordance with previous literature on the calculation of K [16,[36][37][38][39], a normalization process was afterwards applied to scale the κ values to the 0-1 range.Correlations between the abiotic predictors and the selected BMI metrics are shown in the Appendix A (Table A1).

Normalization of κ in the 0-1 Scale, Using Five Alternatives
We applied five alternatives to normalize the unnormalized habitat suitability values: 1.
The unnormalized habitat suitability values are divided by the maximum κ observed at the whole dataset (seasonal and typological variation not accounted; no seasonal grouping, no typological grouping applied-hereafter called maxall).
K i is the normalized habitat suitability of the ith microhabitat ranging from 0 to 1; n i denotes the number of BMI taxa (families) of the ith microhabitat; H i denotes the Shannon's diversity index for the ith microhabitat; EPT i is the number of EPT taxa found at the ith microhabitat; α i is the abundance of benthic macroinvertebrates of the ith microhabitat; κ max denotes the maximum unnormalized habitat suitability value observed in the whole dataset.

2.
The unnormalized habitat suitability values are divided by the maximum κ observed at each season (seasonal variation; no typological grouping applied-hereafter called maxseason).
K n i is the normalized habitat suitability of the ith microhabitat at the nth season, ranging from 0 to 1; n n i denotes the number of BMI taxa (families) found at the ith microhabitat at the nth season; H n i denotes the Shannon's diversity index for the ith microhabitat of the nth season; EPT n i is the number of EPT taxa found at the ith microhabitat at the nth season; α n i is the abundance of benthic macroinvertebrates of the ith microhabitat at the nth season; κ n max denotes the maximum unnormalized habitat suitability value observed at the nth season.

3.
The unnormalized habitat suitability values are divided by the maximum κ observed at each river type (typological variation; no seasonal grouping applied-hereafter called maxtype).
is the normalized habitat suitability of the ith microhabitat of the pth river type, ranging from 0 to 1; n p i denotes the number of BMI taxa (families) found at the ith microhabitat of the pth type; H p i denotes the Shannon's diversity index for the ith microhabitat of the pth type; EPT p i is the number of EPT taxa found at the ith microhabitat of the pth type; α p i is the abundance of benthic macroinvertebrates found at the ith microhabitat of the pth type; κ p max denotes the maximum unnormalized habitat suitability value observed of the pth type.

4.
The unnormalized habitat suitability values are divided by the maximum κ observed at each site at each season (seasonal and typological grouping-hereafter called maxsite).
i is the normalized habitat suitability of the ith microhabitat at the jth site at the nth season ranging from 0 to 1; n j,n i denotes the number of BMI taxa (families) of the ith microhabitat at the jth site at the nth season; H j,n i denotes the Shannon's diversity index for the ith microhabitat of the jth site at the nth season; EPT j,n i is the number of EPT taxa found at the ith microhabitat of the jth site at the nth season; α j,n i is the abundance of benthic macroinvertebrates of the ith microhabitat at the jth site at the nth season; κ j,n max denotes the maximum unnormalized habitat suitability value observed at the jth site at the nth season.

5.
The value of each metric is divided by the maximum value of this metric observed at each site at each season (seasonal and typological grouping-hereafter called maxmetric).
i is the normalized habitat suitability of the ith microhabitat at the jth site at the nth season ranging from 0 to 1; n j,n max denotes the maximum number of BMI taxa (families) of the jth site at the nth season; H j,n max denotes the maximum Shannon's diversity index of the jth site at the nth season; EPT j,n max is the maximum number of EPT taxa found at the jth site at the nth season; α j,n max is the maximum abundance of benthic macroinvertebrates of the jth site at the nth season.

Training of a Habitat Model Using Each of the Five κ-Normalization Options within Three Seasons and Two River Types
A fuzzy rule-based Bayesian algorithm (FRB) based on Brookes et al. [40], detailed in Theodoropoulos et al. [18] and implemented in the HABFUZZ software [41] was trained and cross-validated (see Supplementary Material-Table S2 for details).In brief, the FRB algorithm applies a fuzzification process to convert the numerical input values of V and D to membership degrees (MDs) in overlapping, trapezoidal-shaped fuzzy sets [19] (Table 1).By this process, each numerical input is assigned to one or more fuzzy sets with a MD ranging from zero to one; K values are also classified in five classes (0 ≤ bad ≤ 0.2; 0.2 < poor ≤ 0.4; 0.4 < moderate ≤ 0.6; 0.6 < good ≤ 0.8; 0.8 < high ≤ 1).The training dataset (in this case the multiple alternatives of the benthos-GR dataset), with a priori calculated K values, is used to develop sets of data-driven IF-THEN rules (Table S2), relating the simultaneous occurrence of a series of fuzzy sets to a specific K class.The fuzzy MD of each input is afterwards accounted as the probability of occurrence of the particular fuzzy set, such as 'IF V is low with a MD of 1 AND D is moderate with a MD of 1 AND S is gravel with a MD of 1 THEN K is good with a MD of 0.3 and moderate with a MD of 0.7'.The IF-THEN rules are then combined using the Bayesian joint probability, so that (referring to the previous example) the probability of the specific microhabitat's K being good is the joint probability that V is low AND D is moderate AND S is gravel AND K is good (1 × 1 × 1 × 0.3 = 0.3), while the probability of K being moderate is the joint probability that V is low AND D is moderate AND S is gravel AND K is moderate (1 × 1 × 1 × 0.7 = 0.7).Based on a utility function [18,40], a score is assigned at each K class (bad: 0.1; poor: 0.3; moderate: 0.5; good: 0.7; high: 0.9) and the habitat suitability for each microhabitat is predicted using the following equation: where, K p is the predicted habitat suitability; M ij denotes the joint probability of occurrence of each κ class; S ij denotes the score of each κ class; For the previous example, K p equals to 0.7 × 0.5 + 0.3 × 0.7 = 0.56.The predicted K value is then compared with the K observed in the training dataset for the specific microhabitat and the predictive accuracy of the algorithm is estimated within a three-K-class system (low: 0-0.2; moderate: 0.2-0.6;acceptable: 0.6-1) (Table S2), using a 10-fold cross-validation process [42]: The initial dataset is randomly partitioned in ten equal-sized subsamples.Nine subsamples are used as the training dataset and the remaining subsample is used for model validation.This process is repeated ten times (folds), using a different subsample for validation at each iteration.The performance of each model is evaluated as the average percentage of the correctly classified instances (%CCI) between each iteration of the ten-fold cross-validation process (Appendix A-Figures A6 and A7).

Hydrodynamic Simulation of Two River Reaches to Acquire V, D and S Values in Multiple Discharges (Q)
Topographic data (longitude, latitude, bottom elevation) and hydrometric data (Q, V and D at randomly selected points in two different discharges), were recorded in two river reaches in Greece (Parapeiros River and Oinoi Stream-Figure 2) at multiple sampling campaigns.The Blue Kenue software [43] was used to construct a triangulated computational mesh for each river reach based on the relevant topographic information (Parapeiros-reach mesh properties: 277-m long; 9875 elements; 5170 nodes, Oinoi-reach grid properties: 370-m long; 7140 elements; 3938 nodes).The TELEMAC 2D hydrodynamic model (Electricité de France, Paris, France) [44] was afterwards used to simulate V and D values at each node of each computational mesh in multiple discharges (11 Q scenarios for the Parapeiros reach and 16 Q scenarios for the Oinoi reach).The models were calibrated-validated as follows: initial values for the Manning's roughness coefficient (n) were acquired from the tables given in Chow [45] through an on-site visual estimation of the major types of substrate at each reach.The n values were afterwards properly adjusted within multiple model runs until an acceptable agreement was observed between the field-recorded and the simulated values of V and D at each of the randomly selected points at each river reach.The models were calibrated in the lower discharges and validated in the higher discharges.As this study is focused on habitat modelling, for further details on the hydrodynamic simulations of the Parapeiros and Oinoi reaches, please refer to [23,46], respectively.

Prediction of K in the Two Test Reaches Using the Various Training Alternatives
Within the process mentioned in Section 2.4, we developed 20 training alternatives (subsets) to explore the seasonal and typological variation in the BMI habitat preferences and the relevant variation in the selection of the minimum acceptable and optimal environmental flow scenario within each subset (Figure 3).As it is obvious, not all normalization alternatives can be used in all cases.For example, since in the maxall option, the κ values are divided by the maximum observed κ at the whole benthos-GR dataset, the maxall option cannot be applied only for the spring (or only for the RM4) samples; in contrast, the maxseason, maxsite, and maxmetric options can be applied in a seasonal subset, and likewise, the maxtype, maxsite, and maxmetric options can be applied in a typological subset.Each of the aforementioned 20 subsets was used as a training dataset in the HABFUZZ habitat model, based on which the integrated TELEMAC 2D + HABFUZZ model predicted K at each node of each computational mesh according to the V, D, and S values simulated by TELEMAC 2D at each discharge.In total, the K values in 540 discharge scenarios were simulated; 20 subsets × 11 Q scenarios for the Parapeiros reach and 20 subsets × 16 Q scenarios for the Oinoi reach.

Development of Environmental Flow Scenarios
Similarly to the fish-based 'Weighted Usable Area' approach of Bovee [16], to advance from the microhabitat scale to the reach scale, the K multimetric index (microhabitat scale) was combined with abiotic metrics describing habitat connectivity, habitat availability, and the percentage of wetted nodes, within the Optima Flow Scenario (OFS) index (reach scale) as follows: 1.
Normalized OSI (nOSI): nOSI = OSI w where K i (from 0 to 1) denotes the habitat suitability at each node of the computational grid; w denotes the total No. of wetted nodes in the computational grid at each Q scenario.

Certainty of prediction (COP):
The ratio of the No. of microhabitat combinations actually found in the training dataset to the total No. of nodes in the computational mesh; instead of requiring the user's interference to manually adjust the missing fuzzy rules, HABFUZZ is completely data-driven and when a microhabitat combination is not found in the training dataset, instead of returning some arbitrary K value for a particular node (e.g., −1), it uses the K value of its neighboring node in the mesh, with a simultaneous assessment of the relevant prediction error.4.
Percentage of wetted nodes in the computational mesh at each Q scenario (w).

Habitat connectivity (C):
The ratio of connected (neighboring) nodes with K > 0.6 to the total number of wetted nodes with K > 0.6.

6.
Habitat availability (A): The ratio of connected (neighboring) nodes with K > 0.6 to the total number of nodes in the study reach (wetted and dry).
All OFS i values were normalized in a 0-1 scale by dividing each OFS i with the maximum OFS observed.Based on the status classification system introduced in the Water Framework Directive 2000/60/EC [47], all Q scenarios with OFS values higher than 0.6 were considered acceptable as environmental flows; Seasonal and typological comparisons were applied based on the minimum and optimal environmental flow predictions to ultimately explore the variation of the different environmental flow scenarios developed by the various subsets, corresponding to spatially and temporally different habitat suitability criteria.

Results
The seasonal and typological OFS means and standard deviations (SD) within the various κ-normalization alternatives (see Tables A2-A5 for the actual mean and SD values), along with the OFS values for each discharge simulated in the various seasons and river types, are depicted in Figures 4-6 for the Parapeiros reach and in Figures 7-9 for the Oinoi reach.Two trends are evidenced in both reaches: i.
The seasonal differences in the OFS values for the same Q within the various subsets were greater and more variable than the relevant typological differences.For example, in the Parapeiros reach, the OFS value for Q = 0.3 m 3 /s was 0.18 in spring, 0.77, in summer and 0.04 in autumn (mean: 0.33; SD: 0.39; maxseason normalization).For the same Q, the OFS value for the RM1-2 type was 0.46 and for the RM4 type was 0.86 (mean: 0.66; SD: 0.29; maxtype normalization).In the Oinoi reach, the OFS value for Q = 0.05 m 3 /s was 0.91 in spring, 0.61 in summer, and 0.26 in autumn (mean: 0.59; SD: 0.33; maxseason).For the RM1-2 type, the OFS for Q = 0.05 m 3 /s was 0.89, and for the RM4 type it was 0.65 (mean: 0.77; SD: 0.17; maxtype).ii.
The observed seasonal and typological variation decreased when site-based or metric-based normalizations were applied (maxsite and maxmetric, respectively) with the maxmetric option mostly showing the lowest variation among seasons and river types.In the seasonal comparisons of the Parapeiros reach, the maxmetric option had the lowest OFS-SD for all Q values (100% lower SD) when compared with the maxseason option and for 7 out of 11 discharges (64%) when compared with the maxsite option.In the relevant typological comparisons, the maxmetric option showed 81% lower SD, in comparison with the maxtype option and 64% lower SD when compared with the maxsite option.The seasonal comparisons for the Oinoi reach showed that the SD for the maxmetric option was 81% lower (13 out of 16 Q values) when compared with both the maxseason and maxsite alternatives.The relevant typological comparisons however, showed only 31% lower SD values (5 out of 16 Q values) in comparison with the maxtype option.The comparison between the maxsite and maxmetric options, once again showed 68% lower SD for the maxmetric normalization.When comparing the seasonal and typological differences in the prediction of the minimum acceptable environmental flows (OFS > 0.6) within the various κ-normalization options, it can be seen that the predictions were highly variable; for the Parapeiros reach (Figure 10), the minimum acceptable environmental flow varied from 0.01 m 3 /s (maxseason-spring) to 2 m 3 /s (maxsite-summer).For the Oinoi reach (Figure 11) the minimum environmental flow varied from 0.01 m 3 /s (maxseason-spring; maxtype-RM1-2) to 0.9 m 3 /s (maxsite-RM4).It must be noted that, within this variation, when environmental flows are predicted by pooling all seasons and types under the various κ-normalization options, the minimum environmental flows predicted within the maxsite and maxmetric alternatives in both reaches, mostly lay within the acceptable environmental flows predicted by all other options (hereafter called between-models' degree of agreement-DoA), in contrast to the other training alternatives (maxseason, maxtype, maxall).In the Parapeiros reach, maximum DoA was reached by pooling the dataset within the maxsite option (90%; within acceptable limits in 18 out of the 20 predictions).The maxmetric option showed 65% DoA and all other alternatives had 35% agreement.In the Oinoi reach, the maxsite option showed again maximum agreement (90%), with the maxmetric, maxseason, maxall, and maxtype options showing 75%, 75%, 30%, and 15% DoA, respectively.

Seasonal and Temporal Variation in the Habitat Preferences of Benthic Invertebrates
The results of the study suggest that benthic macroinvertebrates shift their habitat preferences among seasons and in different geographical locations, with the seasonal shifts being greater than the typological ones (the between-seasons SD values were constantly higher than the between-types SD).During summer, benthic invertebrates preferred faster-flowing and deeper habitats (the relevant habitat suitability peaked at higher discharges, reflecting higher V and D values).During spring and autumn, their preferences were more variable, but in general, a preference for shallower and slow-flowing habitats was observed, with the habitat suitability peaking at the lowest discharges mostly in spring.Although seasonal differences in macroinvertebrate-community abundance and composition have long been reported in literature [48][49][50][51][52], community-shifts in their habitat preferences have not been studied; we assume that such 'shifting behavior' may be associated with the seasonal fluctuation of environmental variables, which are known primary drivers of BMI-community changes, that is, shading, water temperature and dissolved oxygen [53].As summer approaches, shading from the surrounding riparian vegetation is reduced, the water temperature increases, the dissolved oxygen concentration decreases [54] and consequently, faster flowing, deeper, better oxygenated habitats provide shelter against the changing environment.In contrast to the under-studied seasonal shifts, typological differences in the habitat preferences of BMIs have been previously documented, suggesting a preference for shallower water depths and lower velocities in smaller rivers [15,55].Our results also confirm this trend; the habitat suitability in the RM1-2 type, which includes mid-and lowland sites draining larger catchments, peaked at higher discharges, reflecting a BMI preference for higher velocities and water depths, in comparison to the RM4 type consisting of highland sites with small-sized catchments.

Seasonal and Typological Variation in the Environmental Flow Predictions within the Various κ-Normalization Options
As it was expected, the seasonal and typological variation in the habitat suitability of benthic macroinvertebrates was highly reflected in the minimum acceptable and optimal model-based environmental flow (e-flow) predictions.E-flow recommendations based on habitat suitability criteria developed from spring samples only, were slightly or highly different than e-flow recommendations based on habitat suitability criteria developed from summer or autumn samples of the same set of sites (maxseason option-Figures 6 and 9).The same trend, but with smaller differences, was observed with the use of habitat suitability criteria from different river types (maxtype option-Figures 6 and 9).This suggests that local and season-specific habitat criteria should be used to increase predictive accuracy-it also confirms most previous BMI-based studies (but see [56])-accounting for the observed spatial and temporal variation [4,5,11,14].
However, despite the observed and acknowledged spatiotemporal variation, two key-findings of our study support the use of generic habitat suitability criteria after proper data pre-treatment: (i) among the various κ-normalization options, the maxmetric option significantly reduced the observed seasonal and typological variation (maxmetric-Figures 6 and 9), thus being a potential candidate for generic applications; (ii) when the habitat suitability values in the reference dataset were normalized per site or per metric (maxsite and maxmetric options, respectively) the final environmental flow predictions mostly lay within the acceptable environmental flows range predicted by all other options.This means that such a generic-criteria-based assessment will have a high probability, ranging from 65% (maxmetric) to 90% (maxsite) to develop environmental flow recommendations, which will finally benefit the benthic-invertebrate community.Similarly to the majority-vote approach commonly followed in machine-learning algorithms, all models, either season-or type-specific, agreed-voted that the minimum acceptable e-flow predictions developed by pooling all samples, pre-treated using the maxsite and maxmetric options were also acceptable (though not minimum) by the other models (Figures 10 and 11).The fact that this trend was observed twice in this study (for both river reaches studied) leads us to assume that it was not due to chance but future studies are of course required to either further confirm or disprove the aforementioned trend.The use of generic habitat suitability criteria has also been suggested by the results of Komínková et al. [56], which showed only 10% seasonal variation in their BMI-based e-flow predictions for an urban creek, with the seasonal, minimum-acceptable e-flows being almost equal.

Issues to Be Considered in Spatiotemporal Pooling of Hydroecological Data
Habitat models predict complex distributional patterns based on a reduced set of predictor variables [28].Due to the stochastic nature of hydroecological data [57], mismatches between the models' prediction and the actual observed are inevitable; spatiotemporal pooling of hydroecological data has been considered as an additional source of prediction error [28].The benefits and limitations of spatiotemporal pooling have long been discussed in literature and despite the multitude of data-mining methods, research on the development of robust data-mining applications is still ongoing [58].Our study, inter alia, confirmed the aforementioned, showing varying results in many internal aspects of the habitat modelling process, which should be taken into account when spatiotemporal pooling of hydroecological data is attempted: 1.
The fuzzy rules developed from the various κ-normalization options were very different; the habitat suitability class for the same microhabitat combination varied from bad to high (Table S2), reflecting the importance of the normalization process (seasonal, typological, siteor metric-based).Theoretically, normalizing a spatially variable dataset, that includes multiple river types, only by season, does not take into account the geographical-typological variation.Normalizing a temporally variable dataset, that includes multiple seasons, only by type, does not take into account the seasonal variation.The very low environmental flow values (0.01 m 3 /s) predicted by the seasonal models are possibly indicative of a relevant inadequacy of these κ-normalization options.We assume that a per-site or per-metric normalization could partially account for both seasonal and typological variation, but in the absence of a field-validation for our models (they were cross-validated but not field-validated), we can currently only discuss the within-models-agreement trends on the final environmental flow prediction.

2.
Low sample sizes often resulted in decreased cross-validation accuracy; this was evident by the %CCI values of the autumn samples (n = 60; Figures A6 and A7), which were the lowest observed within the various models developed, varying from 46% to 52%.The current FRB algorithm requires 5 × 5 × 8 = 200 fuzzy rules to adequately predict the habitat suitability (although not all rules may be necessary, for example, a microhabitat combination of very high V and very deep D will rarely be observed in a river reach).Thus, a small dataset will be inherently incapable of providing the observations required to develop an effective rules-database [28].As the sample size increases however, more microhabitat combinations will be included, the fuzzy-rules database will increase, and the model's performance is also expected to increase [59,60].

Can We Use Generic Instead of Local and Season-Specific Habitat Suitability Criteria?
Before answering, it must be noted that the comparisons applied in our study included typologically different sites, but of similar hydromorphological and hydraulic properties, that is, flow velocities ranging from 0 m/s to 1.2 m/s, water depths ranging from 0.01 m to 1 m and substrate types consisting mainly of boulders, cobbles, pebbles, and gravel.It is obvious that the application of generic habitat suitability criteria from these sites in rivers with completely different hydromorphology (e.g., deep, slow-flowing habitats composed of sand and silt) is prohibitive.
The results of our study clearly showed spatial and temporal shifts in the habitat preferences of benthic invertebrates.In combination, however, with the key-findings from the comparison between the various κ-normalization options, and within the aforementioned limitations, the answer to the question is the following: In ideal applications, local and season-specific habitat suitability criteria should be used to maximize predictive performance and develop seasonally and typologically variable environmental flow scenarios.In real-life applications, however, where a balance between predictive accuracy, generality, cost-effectiveness, and time-efficiency is usually sought, generic habitat suitability models, developed within proper data-treatment options to compensate between loss in performance and generality, can be applied.

Conclusions
This study showed that:

•
Benthic macroinvertebrates shift their habitat preferences among seasons and in different geographical locations; this resulted in highly variable model-based environmental flow predictions between the local, season-specific and generic habitat models.

•
Local and season-specific habitat suitability criteria should be used to maximize predictive accuracy, accounting for the observed spatiotemporal variation.

•
Spatiotemporal data pooling increases sample size and, possibly, predictive accuracy, but has inherent limitations, mainly associated with the normalization process, which should be carefully selected within pooling attempts.

•
With proper pre-treatment (per-site or per-metric κ-normalization), spatiotemporally variable datasets can be aggregated to develop generic habitat suitability criteria that can be used to implement model-based environmental flow assessments in multiple locations of similar hydromorphological and hydraulic properties; the loss of predictive accuracy from the data-aggregation process has a high probability ranging from 65% to 90% to lie within the acceptable range of environmental flow predictions that would be made by local and season-specific habitat models.

Figure 1 .
Figure 1.Three-dimensional representation of the 380 microhabitat samples of the benthos-GR dataset spanning three seasons (a) and two river types (b).V: flow velocity, D: water depth, S: substrate type.

Figure 2 .
Figure 2. The Parapeiros River and Oinoi Stream where the relevant reaches were hydrodynamically simulated.

Figure 3 .
Figure 3. Seasonal and typological subsets of the benthos-GR dataset developed within the various κ-normalization methods.In total, 20 subsets were developed and 20 environmental flow assessments were applied at each test case, using each subset as a training dataset for the habitat model.

Figure 4 .
Figure 4. Graphical representation of the Seasonal optimal flow scenario (OFS) means (grey shaded) and standard deviations within the various κ-normalization alternatives for the Parapeiros reach.The OFS values of spring, summer and autumn for each discharge (Q) have been averaged.

Figure 5 .
Figure 5. Graphical representation of the typological OFS means (gray shaded) and standard deviations within the various κ-normalization alternatives for the Parapeiros reach.The OFS values of RM1-2 and RM4 types for each discharge (Q) have been averaged.

Figure 6 .
Figure 6.Seasonal and typological variation in the selection of the OFS for the Parapeiros reach within the various habitat-suitability calculation-normalization options (Sp: Spring; Su: Summer; Au: Autumn).

Figure 7 .
Figure 7. Graphical representation of the seasonal OFS means (gray shaded) and standard deviations within the various κ-normalization alternatives for the Oinoi reach.The OFS values of spring, summer and autumn for each discharge (Q) have been averaged.

Figure 8 .
Figure 8. Graphical representation of the typological OFS means (gray shaded) and standard deviations within the various κ-normalization alternatives for the Oinoi reach.The OFS values of RM1-2 and RM4 types for each discharge (Q) have been averaged.

Figure 9 .
Figure 9. Seasonal and typological variation in the selection of the optimal flow scenario (OFS) for the Oinoi reach within the various habitat-suitability calculation-normalization options (Sp: Spring; Su: Summer; Au: Autumn).

Figure 10 .
Figure 10.Overall comparison on the selection of the minimum acceptable and optimal environmental flow in the Parapeiros River between the various habitat-suitability calculation-normalization options including seasonal and typological variation (Sp: Spring; Su: Summer; Au: Autumn).Predictions within the green area are considered acceptable (lying between the minimum and optimal environmental flows).

Figure 11 .
Figure 11.Overall comparison on the selection of the minimum acceptable and optimal environmental flow in the Oinoi Stream between the various habitat-suitability calculation-normalization options including seasonal and typological variation (Sp: Spring; Su: Summer; Au: Autumn).Predictions within the green area are considered acceptable (lying between the minimum and optimal environmental flows).

Figure A6 .
Figure A6.Mean (gray shaded columns) and standard deviation values from the 10-fold cross validation process for each κ-normalization alternative (model) applied for the Parapeiros reach.CCI: Correctly Classified Instances.

Figure A7 .
Figure A7.Mean (gray shaded columns) and standard deviation values from the 10-fold cross validation process for each κ-normalization alternative (model) applied for the Oinoi reach.CCI: Correctly Classified Instances.

Table 1 .
Class names and relevant class parameters of the hydraulic variables.Substrate type (S) and habitat suitability (K) were treated as crisp (not fuzzy) sets.V: flow velocity; D: water depth.

Table A2 .
Seasonal OFS means and standard deviations (SD) within the various κ-normalization alternatives for the Parapeiros reach.The OFS values of spring, summer and autumn for each discharge (Q) have been averaged.

Table A3 .
Typological OFS means and standard deviations (SD) within the various κ-normalization alternatives for the Parapeiros reach.The OFS values of RM1-2 and RM4 types for each discharge (Q) have been averaged.

Table A4 .
Seasonal OFS means and standard deviations (SD) within the various κ-normalization alternatives for the Oinoi reach.The OFS values of spring, summer and autumn for each discharge (Q) have been averaged.

Table A5 .
Typological OFS means and standard deviations (SD) within the various κ-normalization alternatives for the Oinoi reach.The OFS values of RM1-2 and RM4 types for each discharge (Q) have been averaged.