Numerical modelling as a support tool for river habitat studies

Numerical modelling is becoming a major tool for supporting environmental studies at different scales, thanks to the capability of up-to-date codes in reproducing the natural behaviour in a quite reliable manner. In evaluating the habitat diversity of anthropized rivers, however, many issues are rising because of the intrinsic complexity of the processes involved. Using a reach of the Po River in Italy as a case study, the present works aims to provide an estimate of the changes of the Eco-Environmental Diversity as a response to different constant flow discharges. The goals are achieved by means of two solvers of the iRIC suite, applied in sequence to firstly simulate the fluvial hydrodynamics and subsequently provide an estimate of the habitat conditions. Despite the simplifications intrinsically present in the models and the ones introduced for practical purposes, the results pointed out that the reduction of the flow discharge recently observed can threat the overall biological status of the river. Because of the modelling uncertainties, on the other side, these preliminary outcomes show the need for more research, both in terms of data acquisition and numerical schematization, for adequately evaluate the effects of transient hydrology on the river ecosystems. Moreover, additional field surveys are necessary to calibrate and validate the used model for having sufficiently reliable estimates.


Introduction
The response of fluvial ecosystems to external and internal drivers is a matter of study since more than a century [1,2], and involved scientists having very different backgrounds and applying manifold techniques.Research on the responses of river systems to external forcing has considered the influence of single [3], sequential [4] and multiple [5] stressors, showing that the cumulative effect produced by the latter is far for being completely understood, especially in large rivers.
Given that large watercourses are among the most productive ecosystems worldwide [6], they are also highly affected by anthropogenic activities such as dredging, bank stabilization or dam construction [2,[7][8][9], involving significant morphological modifications or species invasions [10,11].Nowadays, the use of numerical models for evaluating the effects of multiple natural and anthropogenic stressors in impacted ecosystems is becoming paramount thanks to the improvements in computational capacity.In fact, several modelling tools are available to evaluate the dynamics of river ecosystems at appropriate spatial and temporal scales [2,3,12].
To provide additional insights on the use of numerical models in assessing river habitat, the freeware iRIC suite was chosen because of its flexibility and the availability of different solvers which can be combined to adequately simulate both the eco-hydraulics of streams and large rivers [13].Being one of the most anthropized watercourses in Italy, the Po River represents a perfect case study for testing the capability of numerically reproduce the effects of different hydrological conditions on the Eco-Environmental Diversity (EDD) index.This index, computed on the basis of the Simpson's diversity index [14], is usually adopted to evaluate the fish diversity along watercourses [15], being one of a series of metrics frequently used to assess the hydraulic characteristics of in-stream habitats that need to be considered for successful stream habitat restoration [16].However, 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 just a qualitative indication of the habitat conditions along one of the most impacted Italian rivers.
After the description of the study site, the schematization used in the two numerical models will be addressed, pointing out the simplifications introduced.The results highlight the capability of the iRIC suite in working in a combined way, simulating firstly the fluvial hydrodynamics and then the changes of the habitat in terms of EDD index.The relationship between this latter parameter and the flow discharge shows the possible consequences of the 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 (Fig. 1a).With a drainage area of about 74,000 km², the watershed can be divided into three sectors based on the lithology and on 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 including the Po plain and the Po River delta [17].The middle and lower portions of the Po River are subjected to high flood-hazard and, consequently, heavily embanked since decades [18].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 is around 1000 m 3 /s [19], while the total annual sediment and freshwater discharges to the Northern Adriatic Sea are about 13x10 9 kg and 40-50 km 3 , respectively [20,21].
The iRIC models were applied to a reach of about 70 km between Borgoforte (Mn) and Sermide (Mn).As visible from the 2m Digital Elevation Model (DEM) acquired by the AIPO in 2005 (Fig. 1b), this reach has a mean width of 600 m and a bottom slope in the order of 0.05-0.1%.The present river course has been deeply affected by human interventions during the last century: gravel and sand mining from the riverbed, particularly intense during the 1960s, and creation of river embankments to prevent flooding of the surrounding floodplains and construction of hydropower dams in the Alpine tributaries, promoted a strong decrease in the sediment load and, consequently, an important degradation of the overall river bed, with an oversimplification of its path [8,22,23].Nowadays, indeed, the watercourse is like a single sandy-bedded canal flowing smoothly between fixed levees and vegetated banks, with a recognizable decrease in the habitat conditions.
The hydrological conditions were derived from the analysis of the water discharges measured at the Borgoforte gauging station between 1986 and 2017 (Figure 3c), located at the upstream end of the studied reach.

Numerical modelling
Several hydraulic models have been used to simulate instream habitat suitability at different temporal and spatial scale [16].Among others, the PHABSIM is a well-established hydro-ecological model that provides a suite of tools for the 1-D numerical modelling of hydraulic habitat suitability for fish and invertebrate species [24,25].However, complex environments require a 2-D or even 3-D representation of the flow field [26,27].Therefore, the present application will involve a complete suite composed by a 2-D solver for accurately evaluate the flow velocities and depths along the Po River, coupled with a 2-D tool for estimating the instream habitat suitability.
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 how to configure itself for the solver to be used, and a binary data structure that exchanges information [13].Thus, iRIC is only an interface, which permits to investigate the same case study at different spatial and temporal scales by means of many tools [28].In particular, the modelling of the habitat diversity along the Po River was performed applying two different solvers of the iRIC suite 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
River2D is a 2-d (depth-averaged) finite-element solver developed at the University of Alberta, Canada [29][30][31].The model assumes the flow is hydrostatic, while the bed roughness is computed using the White-Colebrook formulation of the Chézy coefficient and turbulence is included via the specification of the eddy viscosity.The computational time step is reduced thanks to an adaptive time step, and the overall simulations result quickly with respect to other unstructured grid models [13].
The River2D model distributes the flow across the inflow boundary proportional to 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 domain was covered using around 8x10 5 nodes, having a similar size regardless of being within the main channel or along the floodplains.The default upwind coefficient (0.5) was used in River2D's Petrov-Galerkin finite element scheme [32,33], while the groundwater parameters (maximum depth=1 m, transmissivity=0.1 m 2 /s, storativity=1) were accounted for in the wet/dry algorithm [30].The bed roughness ranged between 0.035 s/m 1/3 in the main channel and 0.055 s/m 1/3 along the floodplains and was calibrated in previous applications [17,28].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 is due primarily to the bed material roughness.
Provided that River2D simulates only constant discharges as input boundary conditions, six sets of scenarios were simulated.Three are corresponding to the long-term statistics as measured at the gauging station of Borgoforte in the period 1986-2017 (Figure 1c): minimum, mean and maximum discharges of 760 m 3 /s, 1300 m 3 /s and 2200 m 3 /s, respectively.Other three scenarios are assumed to reproduce the whole hydrograph, accounting for extreme events: 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 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) elevations were extracted from information retrieved from the local water authority.
The simulated time was 1x10 5 s, necessary to stabilize the code avoiding instabilities due to the initial conditions.For a sake of sensitivity, longer and shorter simulations were performed, pointing out that the obtained results are indicative of the hydrodynamic behaviour of the reach.

DHABSIM solver for the ecosystem dynamics
In Japan, where it was developed, DHABSIM (Diversity Habitat Simulation) is frequently adopted to compute the EDD index, in an attempt to quantify the diversity of habitat within a homerange area [34].Recent studies [15] founded that this index was positively correlated with fish species richness, and can, therefore, provide a good estimate of the river status.
This 2D solver imports geometry, velocity, and depth from other flow calculation solvers implemented in iRIC, like River2D, and then uses that information to calculate habitat conditions for small to medium-sized rivers without using Habitat Suitability Indices for specific fish species.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 [35], while the EED index of each point results as a combination of geometry, fish preferences on flow field (depth and velocity), sediment composition and vegetation cover.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>0.8 points, the habitat has good conditions for fish habitat.
For the present simulation, the grid was directly exported from River2D and has a maximum area of 250 m 2 , equal to the home-range.As for the values representing the fish preferences, the ones related to water velocity were derived from literature [11], while for the depth categorization the default values were maintained.For a sake of simplicity, no vegetation was considered in this preliminary application.
Based on these simplifications, to simulate 1x10 5 s, as derived from the hydrodynamic runs, a normal desktop pc needs around 12 hours for each run, pointing out the complexity of the code adopted.

Results
After calibrating the River2D solver against measured data and previous numerical simulations, the ecosystem dynamics can be computed by means of the DHABSIM solver.

River hydrodynamics
The hydrodynamics of this reach of the Po River were simulated accounting for a simplified description of the bed roughness, but the obtained results are in good accordance with previous simulations [28] and field data recently acquired with an Acoustic Doppler Current Profiler [36].
In particular, during mean flow conditions, the flow depths reaches a maximum of 20 m along the thalweg, while the magnitude of the flow velocities ranges is below 1.6 m/s.It is worth to recall, here, that such results represent an average over each cell, having a maximum area of around 250 m 2 , therefore should be interpreted as a general trend of the overall river behaviour rather than very local estimates.As visible from Figure 2 and as one can expect, along the main channel the higher the discharge the higher the flow velocity.But one can also notice that, in the case of high water levels, 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 section, these high flow conditions are paramount in 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 of the river dimensions.Indeed, the first part results quite narrow with a very limited number of secondary channels and bars, while the middle part is characterized by several curve and islands, creating a better habitat for the local fish community.The lower portion of the studied reach has a relatively large floodplain, which is covered by water also during low flow conditions.
The levels in the main channel remain quite constant along the whole reach, given that channel depth and width are compensating each other, and therefore the active cross section results similar.Aiming to evaluate the effects of extreme flood and drought events, two major and one very minor discharges were simulated (Figure 3).As observable (Fig. 3a,b), very high discharges involve very high velocity, over 1.5 m/s along the entire reach.Being the model a fixed-bed one (i.e., only clear water is considered, without describing the morphodynamics) one cannot catch the effect of these high velocities on the sediment transport and the bed morphology, but previous evaluations [28,37] pointed out that the river is subjected to a general erosion and an overall redistribution of sand sediments during floods.Probably such morphological alterations affect the habitat, but, because of the modelling approach adopted, are not adequately represented here.On the other part, very low discharges (Fig. 3c) are worst for the fluvial habitat in a similar manner, because prolonged droughts can interrupt the water, sediment and biota continuity, ultimately affecting the overall ecological status of the watercourse Peer-reviewed version available at Water 2019, 11, 482; doi:10.3390/w11030482

Eco-Environmental Diversity
Despite the non-detailed results obtained with the River2D solver, the simulations performed with the DHABSIM tool provided some indications about the impact of changing hydrology on the habitat diversity.Figure 4 reports the contours of the Eco-Environmental Diversity (EDD) index in the studied reach, pointing out that the typical hydrological flow regime of the Po River is not highly influencing the fish habitat.Comparing the three situations, one can notice that the higher the discharge, the 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 [15], if the EED index increases by a 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, the present discharge variability does not influence the habitat conditions 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 extreme events are considered, its averaged value can be assumed as a good proxy of describing the overall river behaviour (Figure 5).Looking at the very long-term of flow discharges measured along the Po River [37], one can notice a slight decrease in the recent years, which could be addressed by the combined effect of increased human pressure and the climate change.While the maximum values result rather constant at the longterm, a statistically significant decrease is observable in the minima, with prolonged droughts during the summer months.As recognizable in Figure 5, this decline of the flow discharge involves a decrease in the habitat conditions, potentially threatening the overall biota.

Discussion
The opportunity to combine two different solvers permits to provide some insights about the possible changes on habitat conditions because of an alteration of the present river regime.Such alteration could be related to both natural (climate change) and anthropogenic (dams, diversions) stressors, with different consequences on the fish community.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 [39,40].In fact, on the one side, there are many shortcomings related to data availability and numerical coding that prevent the application of the modelling tools.On the other side, the lack of knowledge about water-sediment-biota interactions, especially in high anthropized environments, prevents scientists to adequately schematize and reproduce these processes.
Regarding the preliminary outcomes reported in this paper, the use of two combined modelling tools results paramount in providing an estimate of the response of the river habitat to different hydrological drivers, but the limitations of the codes clearly affect the final results.Indeed, as previously pointed out, the lack of consideration about the sediment transport and the morphodynamics can involve an underestimation of the habitat changes, being the water depth just directly correlated to the flow discharge, and therefore fixed in time.
Because of stream morphodynamics, flow hydraulics and aquatic habitat are closely linked and sensitive to changes in the governing conditions [41,42], any alteration of the natural variability of the flow conditions results in a variation of both the local hydro-morphology and the fluvial habitat.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.), the river hydrology [43] and the habitat resilience.
This first application to a reach of the Po River, however, even based on several simplifications, pointed out the opportunity to use numerical modelling as a support tool for river habitat studies.In this case, to overcome some assumptions towards a more reliable model, the planned steps to be performed in the future comprise: i) the description of the domain by means of a finer grid, aiming to evaluate the very local hydrodynamics 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; iii) the consideration of riparian vegetation along the banks and on the island, with possible effects on the local bed roughness.As for this last parameter, it can be changed spatially, but this requires very detailed and spatially-distributed information about the water depth, which was not available in this study, to accurately calibrate the code.
Regarding the habitat parameters assumed for describing the fish preferences in terms of flow depth and velocity, these values derive from a study performed on the same basin and some default values derived from Japanese streams.Obviously, a better estimate can be derived only from specific information about the fish species present along the reach.An extensive campaign, on the other side, requires a lot of effort both economically and in terms of workforce, and is, therefore, difficult to have detailed data about the fish population, especially in the case of large rivers.
Being a simple schematization of the real world, the modelling cascade proposed here considers just a few variables that can affect the 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.In addition, coupling several models is data and labour intensive, and the transferability of the obtained results, especially the biological ones, is limited to watercourses and catchments having comparable characteristics [43].However, provided the need for improving the knowledge on hydro-ecological interactions, the coupling of hydrodynamics and ecological models can be a powerful tool for future research.

Conclusions
Using a 70-km reach of the Po River, in Italy, two numerical solvers were combined to resolve firstly the fluvial hydrodynamics and, based on that, to provide an index of the fish diversity along the stream.Several runs were performed aiming to characterize the river response to a variable flow discharge, spanning from flooding conditions towards very low water depths.
Even if base on several simplifications both under a modelling point of view and because of the lack of information about the local biota, the preliminary results reported here show the dependence of the river habitat from the local hydrology.On the one side, in fact, an excessive water flow velocity involves unsustainable conditions for the local fauna, preventing the fish spawning.On the other part, a prolonged drought state affects the community in a similar manner, decreasing the EDD in a relevant manner.
Thanks to its simplicity, the proposed methodology can be applied for a first estimate of the habitat conditions as a response of changing hydrology but, for providing reliable and quantifiable values, a proper calibration and validation of the models should be performed, in particular regarding the biological aspects.Based on these promising outcomes, for the future, several improvements are planned, both on the numerical schematization and on the field-based evidence.

Figure 1 .
Figure 1.(a) Po River basin; (b) DEM of the studied area; (c) hydrological time series between 1986 and 2017 measured at the gauging station of Borgoforte (Mn).

Figure 2 .
Figure 2. Simulated velocity profiles for the (a) maximum; (b) medium and (c) minimum discharge.

Figure 4 .
Figure 4. Simulated diversity index for the (a) maximum; (b) medium and (c) minimum discharge.