Numerical Modelling as a Support Tool for River Habitat Studies : An Italian Case Study

Numerical modelling is becoming a major tool for supporting environmental studies at different scales, thanks to the ability of up-to-date codes to reproduce the complex mechanisms of the natural environment in quite a reliable manner. In evaluating the habitat diversity of anthropized rivers, however, many issues are rising because of the intrinsic complexity of the physical processes involved and the limitations associated with numerical models. Using a reach of the Po River in Italy as a case study, the present works aims to provide a qualitative description of the changes of the Eco-Environmental Diversity index as a response to different constant flow discharges typically observed along this reach. The goals are achieved by means of two solvers of the freeware iRIC suite, applied in cascade to first simulate the 2D fluvial hydrodynamics and subsequently provide a qualitative estimate of the habitat conditions. Despite the several simplifications intrinsically present in the modelling cascade and the ones introduced for practical purposes, the results show that an extremely strong and long-lasting reduction of the flow discharge, like the one very recently observed, can ultimately threaten the overall biological status of the river. Because of the modelling uncertainties, these preliminary outcomes are only qualitative and show the need for more research, both in terms of data acquisition and numerical schematization, to adequately and quantitatively evaluate the effects of transient hydrology on the river ecosystems. Moreover, additional field surveys are necessary to calibrate and validate the used biological parameters, aiming to obtain sufficiently reliable estimates.


Introduction
The response of fluvial ecosystems to external and internal drivers has been a matter of study for more than a century [1,2], and the scientists involved have had very different backgrounds and applied manifold techniques.Research on the responses of river systems to external forcing has considered the influence of single [3], sequential [4,5] and multiple [6] stressors, showing that the cumulative effect produced by the latter is far from being completely understood, especially in large rivers.Given that large watercourses are among the most productive ecosystems worldwide [7], they are also highly affected by anthropogenic activities such as dredging, navigation, river training, bank stabilization or dam construction (see, among others, [2,[8][9][10][11][12][13]), involving significant morphological modifications, changes of riparian vegetation and species invasions [14][15][16] acting at manifold scales.Nowadays, the use of numerical models for evaluating the effects of multiple natural and anthropogenic stressors in impacted ecosystems is becoming paramount and widely diffused also because of the improvements in computational capacity.In fact, several modelling tools, both commercial and freeware, are now available and can be applied to evaluate the dynamics of river ecosystems at the appropriate spatial and temporal scales [2,3,17], as well as to assess the benefits deriving from small-and large-scale river restoration projects [18].
To provide additional insights on the use of numerical models in assessing riverine hydrodynamics and habitat conditions, the freeware iRIC (International River Interface Cooperative) suite (i-ric.org/en)was chosen because of its flexibility and the availability of different solvers, which can be combined and were developed for adequately simulating the eco-hydraulics of small and intermediate streams, but also large rivers [19].As one of the most anthropized watercourses in Italy, the Po River represents a perfect case study for testing the suite and evaluating its ability to numerically reproduce the effects of different hydrological conditions on the Eco-Environmental Diversity (EDD) index.Considering different flow discharges, it is possible to evaluate if the forecasted increase of drought periods [20] can affect the habitat conditions of this Italian watercourse.The EDD index, computed on the basis of Simpson's diversity index [21], is usually adopted to evaluate fish diversity along watercourses, especially in Japan [22], being one of a series of metrics frequently used to assess the hydraulic and biological characteristics of in-stream habitats that need to be considered for successful stream habitat restoration [23].Despite using a well-known and tested parametrization, because of the lack of spatially distributed information on fish habitat and the assumptions intrinsically correlated to the numerical codes, the present application represents only a preliminary estimate on EDD index, and is not pretending to propose any quantitative information, but rather to provide a qualitative indication of the habitat conditions along one of the most impacted Italian rivers.It is worth mentioning that, because of the introduced simplifications, the outcomes reported here are site-specific, and can hardly be extended to other watercourses.
After a short description of the study site, the schematization used in the two numerical models will be addressed, pointing out the main simplifications introduced and their impact on the modelling outcomes.Even if based on several schematizations and therefore intrinsically affected by uncertainties, as discussed in the following, the presented results highlight the ability of the iRIC suite to work in a combined way, simulating first the measured fluvial hydrodynamics in a rather good manner when compared with field information, and then the changes of the habitat conditions in terms of the EDD index.Because of a lack of data and the numerous assumptions, however, only qualitative indications on the possible changes of the Po River habitat can be retrieved from the modelling runs.The relationship between the EDD index and the flow discharge can provide additional insights into the possible consequences of the combined impact of climate change and anthropogenic pressure on the future habitat of the Po River.

Case Study
The Po River, the longest watercourse in Italy, flows eastward across northern Italy for about 660 km, from the Alps to the Adriatic Sea (Figure 1a).With a drainage area of about 74,000 km 2 , the watershed can be divided into three sectors based on the lithology and the maximum elevation: the Alpine sector of crystalline and carbonate rocks (maximum relief: around 4500 m asl), the Apennine sector, mostly composed of sedimentary rocks with high clay content (maximum relief: around 2000 m asl) and the central alluvial area that includes the Po plain and the Po River delta [24].The middle and lower portions of the Po River are subjected to high flood hazard and, consequently, have been heavily embanked for decades [25] to protect the highly anthropized surroundings.
The annual hydrograph shows two peaks in discharge, normally in autumn and spring, generated by rainfall and snowmelt, respectively.For the period 1924-2009, the long-term mean annual discharge recorded at the Piacenza gauging station was around 1000 m 3 /s [26], while the total annual sediment and freshwater discharges to the Northern Adriatic Sea were about 13 × 10 9 kg and 40-50 km 3 , respectively [27,28].
The iRIC models were applied to a reach of about 70 km between Borgoforte (Mn) and Sermide (Mn).As visible from the 2 m Digital Terrain Model (DTM) provided by the local water authority AIPO (Agenzia Interregionale per il fiume Po, Italy) in 2005 (Figure 1b), this reach has a mean width of 600-800 m, including the floodplains and the emerging bars, characterizing the Po as an intermediate stream.This DTM combines information retrieved by means of LiDAR (data collected using two different laser scanners: 3033 Optech ALTM and Toposys Falcon II), multi-beam sonar survey for the navigable portion and a traditional ground survey of some river cross sections.The bottom slope of the studied reach is in the order of 0.05%.For the present simulations, the hydrological conditions were derived from the analysis of the water discharges measured at the Borgoforte gauging station between 1986 and 2017 (Figure 1c), located at the upstream end of the studied reach.Even if a slight decrease in yearly flow discharge was observed during the studied period, no statistically meaningful trends can be recognized.The present river course has been heavily affected by human interventions during the last century, with evident modifications of the fluvial pattern.Gravel and sand mining from the riverbed, particularly intense during the 1960s, and the creation of river embankments to prevent flooding of the surrounding floodplains, as well as the construction of hydropower dams along the Alpine tributaries, promoted a severe decrease in the sediment load and, consequently, an important degradation of the overall river bed, with an oversimplification of its path [29][30][31].This long-term anthropogenic pressure altered the overall sediment dynamics, and nowadays the middle reach results in a morphological quasi-equilibrium [9], where the morphology is altered only by extreme events.At present, indeed, the lowland part of the Po River is like a single sandy-bedded canal flowing smoothly between fixed levees and vegetated banks, with a recognizable decrease in the habitat conditions fostered by an overstabilization of banks and bars.

Numerical Modelling
For the present simulations, the hydrological conditions were derived from the analysis of the water discharges measured at the Borgoforte gauging station between 1986 and 2017 (Figure 1c), located at the upstream end of the studied reach.Even if a slight decrease in yearly flow discharge was observed during the studied period, no statistically meaningful trends can be recognized.

Numerical Modelling
Several hydraulic models have been used to simulate instream habitat suitability at different temporal and spatial scale [23].Among others, PHABSIM (Physical Habitat Simulation System) is a well-established hydro-ecological model developed by the United States Geological Survey (www.usgs.gov/software/physical-habitat-simulation-phabsim-software-windows)and provides a suite of tools for the 1D numerical modelling of hydraulic habitat suitability for fish and invertebrate species [32,33].However, such a schematization can be too crude, given that complex environments require a 2D or even a 3D representation of the flow field [31,34,35].Therefore, the present application will involve a complete suite composed by a 2D solver for evaluating flow velocities and depths along the Po River, coupled with a 2D tool for estimating the instream habitat suitability at the reach scale.
In the iRIC (International River Interface Cooperative) structure, solvers are kept separate from the interface, communicating with it through a definition file that instructs the interface on how to configure itself for the solver to be used, and a binary data structure that exchanges information [19].Thus, iRIC is only an interface, which permits to investigate the same case study at different spatial and temporal scales by means of many interacting tools [36].In particular, the modelling of the habitat diversity along the Po River was performed applying two different solvers of the freeware iRIC suite (i-ric.org/en) in a cascade manner: River2D and DHABSIM.Indeed, the latter solver needs river geometry, flow depth and velocity field as input parameters, which can be derived from the output of the first solver.

River 2D Solver for the River Hydraulics
The original code River2D is a 2D (depth-averaged) finite-element solver developed at the University of Alberta, Canada [37,38].A thorough description of its mathematical background can be found in the solver manual [38], so it is not described here for the sake of simplicity, apart from a few characteristics.The original modelling structure was maintained in implementing it in the freeware iRIC solver (i-ric.org/en/solvers/river2d),which adopts a finite-element for solving the depth-averaged de St. Venant equations expressed in a conservative form.The code assumes that the pressure distribution along the vertical is hydrostatic, therefore it is applicable only for lowland reaches where the bed slope changes gently (namely, longitudinal variations of about 10% cannot be correctly modelled), like the presented case study.Moreover, Coriolis and wind forces are assumed to be negligible, as usually made for reproducing relatively short river reaches.The friction slope depends on the bed shear stresses, assumed related to the magnitude and direction of the depth-averaged velocity, and the bed roughness is computed using the White-Colebrook formulation of the Chézy coefficient.As for the depth-averaged transverse turbulent shear stresses, they are modelled following a Boussinesq type eddy viscosity formulation.
The River2D model distributes the flow across the inflow boundary proportional to the depth, resulting in the fastest velocity close to the thalweg.The river region included between the main levees was discretized by means of a triangular irregular network (TIN) grid, with grid elements having a maximum area of 250 m 2 .In this manner, the whole 70-km domain was covered using around 8 × 10 5 nodes, having a similar size regardless of being within the main active channel or along the floodplains.The default upwind coefficient (0.5) was used in River2D's Petrov-Galerkin finite element scheme [39,40], while the default groundwater parameters (maximum depth = 1 m, transmissivity = 0.1 m 2 /s, storativity = 1) were accounted for in the wet/dry algorithm, which changes surface flow equations to groundwater flow equations [38].Such a low value of the transmissivity means that the groundwater discharge can be neglected.The bed roughness was imposed equal to 0.035 s/m 1/3 in the main channel and 0.055 s/m 1/3 along the floodplains, as calibrated in previous applications [25,38].This solver uses a non-dimensional Chézy coefficient to close the stress terms, computed as a function of the effective roughness height.Therefore, the flow resistance along the stream is primarily due to the bed material roughness.
Provided that the River2D code simulates only constant discharges as input boundary conditions, six sets of scenarios were simulated for evaluating the possible flow fields associated with typical discharges observed along the stream.Because of these modelling limitations, yearly discharges were assumed to be more representative of the general river behaviour rather than the daily ones.Therefore, three scenarios correspond to the long-term statistics as measured at the gauging station of Borgoforte in the period 1986-2017 (Figure 1c): minimum, mean and maximum yearly discharges of 760 m 3 /s, 1300 m 3 /s and 2200 m 3 /s, respectively.Other scenarios are assumed for covering the whole hydrograph, accounting for extreme events, assuming that they last for a long period: 200 m 3 /s, 3000 m 3 /s and 5000 m 3 /s.As for the downstream conditions, a calibrated depth-unit discharge relationship was adopted, changing the depth-unit discharge parameter K for each simulation to obtain an equilibrium between the discharges entering and going out from the domain.This parameter spans from 0.12 for the case of 200 m 3 /s to 3.00 in the case of the flooding discharge 5000 m 3 /s.The exponent m of the depth-unit discharge was imposed equal to 0.5 for all the simulations.The fixed upstream (14.55 m asl) and downstream (5.51 m asl) water elevations were extracted from information retrieved from the local water authority.
The computational time step (maximum ∆t = 10 s) is reduced thanks to an adaptive time step, and the overall simulations are quicker with respect to other unstructured grid models [19].The total simulated time was 1 × 10 5 s, necessary to stabilize the code avoiding instabilities due to the initial conditions.For the sake of sensitivity, longer and shorter simulations were performed, pointing out that the adopted parameters are corrected, given that the modelling outcomes are indicative of the hydrodynamic behaviour of the reach.

DHABSIM Solver for the Ecosystem Dynamics
In Japan, where it was first developed, the DHABSIM (Diversity Habitat Simulation) code is frequently adopted to compute the EDD index, in an attempt to quantify the diversity of habitat within a home-range area [39], following an increasing trend of eco-hydrological modelling [41,42].Recent studies [19] found that this index was positively correlated with fish species richness, and can, therefore, provide a good estimate of the river status.
This 2D solver (i-ric.org/en/solvers/dhabsim)imports geometry, velocity, and depth from other flow calculation solvers implemented in iRIC, like the River2D, and then uses this information to compute the habitat conditions for small to medium-sized rivers without using Habitat Suitability Indices for specific fish species.Such a model cannot, therefore, be applied in large lowland rivers, characterized by high discharges and a very wide (kilometres) main channel.
Based on the diversity-based algorithm and the physical flow inputs, the solver outputs a habitat distribution map.The fish home-range derives from the Minns formula [43], and is computed based on the average fish length: where A hr represents the home range (m 2 ) and L the fish length (mm).
The EED index of each point is a combination of geometry, fish preferences on flow field (depth and velocity), sediment composition of the bed and vegetation cover (Equation ( 2)).This index ranges between 0 and 1, and an increment of 0.1 point means that one species of fish may increase in the target river.In the case of EDD index > 0.8 points, the habitat has good conditions for fish habitat.
where S is the total number of categories i (depth, velocity, bed, vegetation) considered, N is the number of nodes in the home range and n i is the node number of category i.Indeed, the code permits to spatially specify the categories by drawing local regions.
For the present simulation, the grid was directly exported from River2D and has a maximum area of 25 m 2 , equal to the 1/10 of the home-range and directly computed by the solver [22].As for the values representing the fish preferences (depth and velocity), they are assumed constant in the whole reach, due to the lack of more detailed information.The ones related to water velocity were derived from the limited literature available [16], while for the depth categorization the default values were maintained because of a lack of specific information.The river bed was considered uniform and constituted by sand.For the sake of simplicity, no vegetation was considered in this preliminary application.On the basis of these simplifications, to simulate the total time t = 1 × 10 5 s, as derived from the hydrodynamic simulations, a normal desktop PC requires around 12 h for each run.
Aiming to obtain a qualitative analysis of the effects of changing hydrology and because of a lack of specific field information, no sensitivity analysis on the EDD index was performed.

Results
The River2D solver was calibrated against measured data [44] and previous numerical simulations [35,45] in terms of (i) water elevation measured by means of a hydrometer at the gauging station of Revere (Figure 1); (ii) velocity field retrieved from an ADCP (Acoustic Current Doppler Profiler) survey performed along a 10-km reach close to Revere.Based on the outcomes of the hydrodynamics solver, the ecosystem dynamics can be computed by applying the DHABSIM solver.

River Hydrodynamics
For the sake of simplicity and to obtain a preliminary estimate of the flow field, the hydrodynamics of this reach of the Po River were simulated accounting for a fixed bed roughness in the main channel (0.035 s/m 1/3 ) and along the floodplains (0.055 s/m 1/3 ), without accounting for spatially variable patterns possibly related to vegetation or bedforms.Even if based on strong simplifications, the obtained results are in good accordance with previous simulations [36,46] and local field data recently acquired with an ADCP [44], and are therefore considered as representative of the overall hydrodynamics of the reach.
In particular, during mean flow conditions, the flow depth reaches a maximum of around 20 m along the thalweg, while the magnitude of the flow velocities ranges is below 1.5 m/s, comparable with field measurements [44].It is worth recalling here that such results represent an average over each cell, having a maximum area of around 250 m 2 , and therefore should be interpreted as a general trend of the overall river behaviour rather than very local estimates.
As visible in Figure 2, where the contours of the velocity fields are reported, in the main channel the higher the discharge, the higher the flow velocity, as one would expect.However, one can also see that, in the case of high water levels (Figure 2a), there are more active channels, especially in the lower part of the studied reach, suggesting that such conditions can be more suitable for the exchange of fluxes between the main channel and the adjacent floodplains.Moreover, as discussed in the following sections, these medium to high flow conditions are paramount for establishing good habitat conditions.In detail, a decrease of the flow velocity is observable going from upstream (Borgoforte) towards downstream (Sermide), because of a reduction of the bed slope and an increase in the river dimensions.Indeed, the first part of the river is quite narrow with a limited number of secondary channels and bars, while the middle part is characterized by several curves and islands, potentially creating a better habitat for the local fish community.The lower portion of the studied reach has a relatively large floodplain, which is also covered by water during low flow conditions (Figure 2c).
The levels in the main channel remain roughly constant along the whole reach, given that channel depth and width are compensating for each other, and therefore the active cross section varies only slightly.
Aiming to evaluate the effects of extreme flood and drought events, two major and one very minor discharges were simulated (Figure 3).As one can see, very high flow discharges involve very high velocity, over 1.5-2.0m/s along the entire reach (Figure 3a).On the other hand, very low discharges (Figure 3c) involve velocities of maximum 0.4 m/s, which can be detrimental to the fluvial habitat, because prolonged droughts can interrupt the water, sediment and biota continuity, ultimately affecting the overall ecological status of the watercourse.In detail, a decrease of the flow velocity is observable going from upstream (Borgoforte) towards downstream (Sermide), because of a reduction of the bed slope and an increase in the river dimensions.Indeed, the first part of the river is quite narrow with a limited number of secondary channels and bars, while the middle part is characterized by several curves and islands, potentially creating a better habitat for the local fish community.The lower portion of the studied reach has a relatively large floodplain, which is also covered by water during low flow conditions (Figure 2c).
The levels in the main channel remain roughly constant along the whole reach, given that channel depth and width are compensating for each other, and therefore the active cross section varies only slightly.
Aiming to evaluate the effects of extreme flood and drought events, two major and one very minor discharges were simulated (Figure 3).As one can see, very high flow discharges involve very high velocity, over 1.5-2.0m/s along the entire reach (Figure 3a).On the other hand, very low discharges (Figure 3c) involve velocities of maximum 0.4 m/s, which can be detrimental to the fluvial habitat, because prolonged droughts can interrupt the water, sediment and biota continuity, ultimately affecting the overall ecological status of the watercourse.Because the model is a fixed-bed one, only clear water is here considered, without describing the bed morphodynamics.Therefore, one cannot adequately capture the effect of high velocities on the sediment transport and bed morphology, but previous evaluations performed on the middle part of the reach [33,43] have suggested that the river is subject to general erosion and an overall redistribution of sand sediments during floods, with possible consequences for the local habitat.On the other hand, considering the effect of the normal hydrological flow regime on the sediment dynamics, [9] pointed out that, nowadays, the river results in a morphological quasi-equilibrium, which can be altered only in the case of flooding.

Eco-Environmental Diversity Index
Despite the non-detailed results obtained with the River2D solver, the simulations performed with the DHABSIM tool provided some qualitative indications of the impact of changing hydrology on the habitat diversity of this river reach.
Figure 4 reports the isolines of the Eco-Environmental Diversity (EDD) index along the studied reach, with a particular focus on the middle part, where several secondary channels are present.Because the model is a fixed-bed one, clear water is here considered, without describing the bed morphodynamics.Therefore, one cannot adequately capture the effect of high velocities on the sediment transport and bed morphology, but previous evaluations performed on the middle part of the reach [33,43] have suggested that the river is subject to general erosion and an overall redistribution of sand sediments during floods, with possible consequences for the local habitat.On the other hand, considering the effect of the normal hydrological flow regime on the sediment dynamics, [9] pointed out that, nowadays, the river results in a morphological quasi-equilibrium, which can be altered only in the case of flooding.

Eco-Environmental Diversity Index
Despite the non-detailed results obtained with the River2D solver, the simulations performed with the DHABSIM tool provided some qualitative indications of the impact of changing hydrology on the habitat diversity of this river reach.
Figure 4 reports the isolines of the Eco-Environmental Diversity (EDD) index along the studied reach, with a particular focus on the middle part, where several secondary channels are present.These modelling outcomes point out that the typical hydrological flow regime of the Po River is not highly influencing the fish habitat, at least at the studied scale.Comparing the three situations, indeed, one can see that the higher the discharge, the (slightly) higher the EDD index, which reaches a maximum of 0.5 points along the main channel.However, given that, as pointed out by the DHABSIM manual [22], if the EED index increases by 0.1 point, a species of fish may increase of about one species, and the optimum is reached for values higher than 0.8 points, it is possible to conclude that the present discharge variability does not influence the habitat in a relevant manner.indeed, one can see that the higher the discharge, the (slightly) higher the EDD index, which reaches a maximum of 0.5 points along the main channel.However, given that, as pointed out by the DHABSIM manual [22], if the EED index increases by 0.1 point, a species of fish may increase of about one species, and the optimum is reached for values higher than 0.8 points, it is possible to conclude that the present discharge variability does not influence the habitat in a relevant manner.Given the relatively large spatial scale adopted to discretize the domain and the low variability of the EDD index along the reach, even if relatively extreme events are considered, its averaged value over the whole reach can be assumed as a good proxy for describing the overall river behaviour (Figure 5).Given the relatively large spatial scale adopted to discretize the domain and the low variability of the EDD index along the reach, even if relatively extreme events are considered, its averaged value over the whole reach can be assumed as a good proxy for describing the overall river behaviour (Figure 5).Looking at the very long term of flow discharges measured along the Po River [45], one can notice a slight decrease in recent years, which could be explained by the combined effect of increased human pressure and climate change.While the maximum values are fairly constant in the long term, a decrease is observable in the minima, with more prolonged droughts during the summer months.
As recognizable in Figure 5, this decline in the flow discharge involves a decrease in the habitat conditions, here represented by the averaged EDD index, potentially threatening the overall river biota.Indeed, for very low discharges that last for long periods, no increase of the EDD index is observable anymore, meaning that no fish species can increase in number.

Discussion
The opportunity to combine two different solvers permits us to provide some insights about the possible changes on habitat conditions because of an alteration of the present river regime, which is paramount to driving future management.Indeed, since such alterations are related to both natural (climate change) and anthropogenic (dams, diversions) stressors, with different consequences on the fish community, they have to be adequately quantified for designing adaptation strategies at the regional and local scales.In this sense, the actual knowledge does not permit a reliable quantification of the impact of these external drivers on ecosystem services and habitat conditions [47,48].In fact, there are many shortcomings related to data availability and numerical coding that prevent the application of cascading modelling tools to gain reliable estimates.On the other hand, a lack of knowledge about water-sediment-biota interactions, especially in high anthropized environments, is a major challenge for scientists interested in adequately schematizing and reproducing these processes [49].
Regarding the preliminary outcomes reported in this paper, the use of two combined modelling tools provides qualitative indications of the response of the river habitat to different hydrological drivers, but the limitations of the codes clearly affect the final results.Even if the hydrodynamics are reproduced in quite a reliable manner with respect to the measured ones [44], there are many shortcomings that should be highlighted.On the one hand, as previously pointed out, a lack of consideration of the sediment transport and the fluvial morphodynamics can involve an underestimation of the habitat changes, especially in sandy-bedded rivers like the Po.Indeed, with the bed geometry being constant, the water depth is directly correlated to the flow discharge, and therefore fixed in time.However, it is well known that, because of stream morphodynamics, flow Looking at the very long term of flow discharges measured along the Po River [45], one can notice a slight decrease in recent years, which be explained by the combined effect of increased human pressure and climate change.While the maximum values are fairly constant in the long term, a decrease is observable in the minima, with more prolonged droughts during the summer months.
As recognizable in Figure 5, this decline in the flow discharge involves a decrease in the habitat conditions, here represented by the averaged EDD index, potentially threatening the overall river biota.Indeed, for very low discharges that last for long periods, no increase of the EDD index is observable anymore, meaning that no fish species can increase in number.

Discussion
The opportunity to combine two different solvers permits us to provide some insights about the possible changes on habitat conditions because of an alteration of the present river regime, which is paramount to driving future management.Indeed, since such alterations are related to both natural (climate change) and anthropogenic (dams, diversions) stressors, with different consequences on the fish community, they have to be adequately quantified for designing adaptation strategies at the regional and local scales.In this sense, the actual knowledge does not permit a reliable quantification of the impact of these external drivers on ecosystem services and habitat conditions [47,48].In fact, there are many shortcomings related to data availability and numerical coding that prevent the application of cascading modelling tools to gain reliable estimates.On the other hand, a lack of knowledge about water-sediment-biota interactions, especially in high anthropized environments, is a major challenge for scientists interested in adequately schematizing and reproducing these processes [49].
Regarding the preliminary outcomes reported in this paper, the use of two combined modelling tools provides qualitative indications of the response of the river habitat to different hydrological drivers, but the limitations of the codes clearly affect the final results.Even if the hydrodynamics are reproduced in quite a reliable manner with respect to the measured ones [44], there are many shortcomings that should be highlighted.On the one hand, as previously pointed out, a lack of consideration of the sediment transport and the fluvial morphodynamics can involve an underestimation of the habitat changes, especially in sandy-bedded rivers like the Po.Indeed, with the bed geometry being constant, the water depth is directly correlated to the flow discharge, and therefore fixed in time.However, it is well known that, because of stream morphodynamics, flow hydraulics and aquatic habitat are closely linked and sensitive to changes in the governing conditions [50,51], and any alteration in the natural variability of the flow conditions can potentially result in changes in both the local hydro-morphology and the fluvial habitat.Generally, in fact, both high and low flow conditions affect the aquatic habitat, with a different magnitude depending on the geomorphic structure of the riverbed (i.e., gravel or sand) and its channel forms (i.e., braided, single-thread, meandering, etc.) [52,53], the river hydrology [54,55] and the habitat resilience.On the other hand, the outcomes reported here should be considered only as a qualitative indication of the possible trends of the EDD index, provided that the parameters used for computing it were not calibrated against field evidence because of their lack.
In detail, regarding the habitat parameters assumed for describing the fish preferences in terms of flow depth and velocity, these values derive from a single study performed on the same basin [15] and some default values derived from Japanese streams.Obviously, a better estimate can be performed only using, as a basis, specific information about the fish species present along the reach.An extensive campaign, on the other hand, requires a lot of effort both economically and in terms of workforce, as demonstrated by the lack of field data.Therefore, it is very difficult to get detailed information about the fish population, especially for performing local-scale evaluations in large rivers.The present research represents not only the first application of a modelling technique to qualitatively estimate the changes in the fluvial habitat, but also a tool for pinpoint the need for additional research on fish population, crucial to future hydro-ecological studies.
This first application to a reach of the Po River, however, even if based on several simplifications, pointed out the potential for using numerical modelling as a support tool for river habitat studies.As shown here, there are many shortcomings that should be addressed to build a more reliable model for helping water managers.The steps planned for the future include: (i) the description of the domain by means of a finer grid, aiming to evaluate the very local hydrodynamics instead of the ones at the reach scale, and thus the habitat conditions at a smaller scale; (ii) the use of numerical codes able to reproduce the real hydro-morphodynamics, aiming to include the effect of sediments and solid transport on the fluvial habitat and not only the water-related impact; (iii) the consideration of riparian vegetation along the banks and on bars and islands, with possible effects on the local bed roughness.As for this last parameter, it can be changed spatially, but this requires very detailed spatially distributed information about the water depth, which was not available in this study, to accurately calibrate the 2D code over the whole reach.Even if the present example shows the potential of the iRIC suite in modelling the fluvial eco-hydrology, several shortcomings are also involved in applying this software.Among others, since the code is still under development, it is not user-friendly like other, more adopted, modelling tools.In particular, the post-processing of the DHABSIM solver is quite simple, and does not provide the user with tabular outcomes, which could eventually be used for elaborating ad hoc metrics for comparing the simulated scenarios.In addition, only a few parameters can be inspected, preventing the opportunity to perform a thorough sensitivity analysis.
Although among the scientific community the topic is still under discussion [56], there is a general agreement that, in the future, the global hydrological cycle will intensify, and extremes of droughts and floods will become more common, fostered by the combined effect of climate change and human pressure.In this context, many Mediterranean and Alpine rivers that used to be permanent tend to become more intermittent, with prolonged droughts, especially during the summer period.As pointed out in similar studies performed on benthic invertebrates [57], although some researchers have claimed that periodic droughts may be important for maintaining the diversity of freshwater systems [58], such shifts highly disturb the environment, with detrimental effects at multiple scales.Even if only from a qualitative point of view, the present study confirms these outcomes, showing that an extremely low discharge can negatively affect the EDD index at the reach scale.
Being a simple schematization of the real world, the modelling cascade proposed here considers just a few variables that can affect the local river habitat.For example, land use practices and climate change are not accounted for in the present application, but they could be very important, especially in highly anthropized environments, even if acting at different spatiotemporal scales.In addition, coupling several models is data-and labour-intensive, and the transferability of the obtained results, especially biological ones, is limited to watercourses and catchments having comparable characteristics [53] and located in similar climatic areas.
Given the need for improving knowledge on hydroecological interactions, the coupling of hydrodynamics and ecological models can be a powerful tool for future research, aiming to evaluate the effect of climate change on the local environment.In the past fluvial hydro-morphology and ecology/biology were generally studied as uncoupled phenomena, with scientists having different backgrounds and using manifold metrics.However, as recently pointed out and also highlighted by this simple example, river environments should be considered as a unicum [59], and studied and modelled as a whole to adequately characterize the interactions between the hydro-morpho-biodynamic processes involved.
On the one hand, the use of numerical codes should always be coupled with a detailed field campaign, aiming to obtain reliable and up-to-date input parameters for calibrating and validating the adopted modelling tool.On the other hand, water managers, practitioners and scientists should work together, aiming to choose the correct numerical model, populate it with the most accurate and field-based data, and derive the main conclusions to better manage the fluvial environment.Only a constant and timely revision of all the adopted numerical tools and monitoring techniques can assure reliable estimates of the future trends of highly anthropized alluvial watercourses, and allow one to plan adequate responses.

Conclusions
Using a 70-km reach of the Po River, in Italy, two numerical solvers were combined to first resolve the 2D fluvial hydrodynamics with grids having an area of around 250 m 2 , and then, based on that, provide a qualitative estimate of the Eco-Environmental Diversity index along the stream.Aiming to characterize the river's response to different flow discharges, several runs were performed, assuming, as input, constant water flows spanning from flooding to drought conditions, following the long-term hydrology measured along the reach.
Even if based on several simplifications under a modelling point of view and because of the lack of information about the local biota, the preliminary results reported here show, from a qualitative point of view, the dependence of the river habitat on the local hydrology, pointing out the need to improve our knowledge of eco-hydrological processes.Excessive water flow velocity and depth, correlated to flooding events, involves unsustainable conditions for the local fauna, preventing fish spawning and growth because of the high flow velocity.Meanwhile, a prolonged drought state like that recently observed along the Po River affects the community in a similar manner, decreasing the EDD index with detrimental effects on fish.
Thanks to its simplicity, the proposed methodology can be applied for a first, qualitative estimate of the local habitat as a response of changing hydrology, assuming the literature parameters for the biota to overcome the lack of information.However, for providing reliable and quantifiable values, proper calibration and validation of the models should be performed, in particular regarding the biological aspects but also the 2D distribution of the water flow.Based on these promising outcomes and for overcoming the modelling limitations highlighted in the present case study, in the future several improvements are planned, both to the numerical schematization and to the field-based evidence.

Water 2019 , 16 Figure 1 .
Figure 1.(a) Po River basin; (b) DEM of the studied area, from the upstream end of Borgoforte to the downstream end of Sermide; (c) yearly time series of flow discharge measured between 1986 and 2017 at the gauging station of Borgoforte (Mn).

Figure 1 .
Figure 1.(a) Po River basin; (b) DEM of the studied area, from the upstream end of Borgoforte to the downstream end of Sermide; (c) yearly time series of flow discharge measured between 1986 and 2017 at the gauging station of Borgoforte (Mn).

Water 2019 ,
11, x FOR PEER REVIEW 10 of 16

Figure 4 .
Figure 4. Simulated Eco-Environmental Diversity index for the (a) maximum, (b) mean and (c) minimum discharge.

Figure 4 .
Figure 4. Simulated Eco-Environmental Diversity index for the (a) maximum, (b) mean and (c) minimum discharge.

Figure 5 .
Figure 5. Averaged Eco-Environmental Diversity index along the reach.