Comparison of West Balkan adult trout habitat predictions using a Pseudo-2D and a 2D hydrodynamic model

Flow regime is a key component for the ecological integrity of all biotic interactions among riverine ecosystems (Lai et al. 2016). Changes in river flow quantity and timing are more likely to cause perturbations to the instream biota and consequences to biodiversity that are more destructive than consequences caused by other choices (Richter & Thomas 2007). Nevertheless, the growing demands for water have increased the need to build dams, which in many cases alter the seasonal and interannual streamflow variability of rivers (Poff et al. 2007), while they also reduce the available amount of water downstream (Li et al. 2012; Mao et al. 2016). In order to mitigate the downstream impact of dams, several countries worldwide have established environmental flow rules in order to quantify the water required for ecosystem conservation and resource protection (Tharme 2003).

Today, there are several environmental flow estimation methods and they can be grouped into four categories: hydrological rules, hydraulic rating methods, habitat simulation methods and holistic methodologies (Dyson et al. 2003; Tharme 2003). The hydrological methods depend on statistical analysis of hydrological data, while the hydraulic methods rely on the estimation of a specific wetted perimeter or maximum depth, based on geomorphology and hydrodynamic behaviour of a river reach. Physical habitat methods quantify trade-offs between water use and environmental benefits of flowing water for particular species or groups of species (Dunbar et al. 2012). Finally, the incorporation of environmental, social and economic features, into a modified or regulated flow regime which would maintain the functional integrity of the riverine ecosystem, comprises the holistic methodologies (Arthington et al. 1992).

Among the previously described categories, physical habitat methods are commonly used to assess habitat quality and quantity of riverine ecosystems. This approach integrates a hydraulic model and a habitat suitability model, usually developed at a microhabitat scale for specific species. The basic idea of this approach is that the hydro-physical stream conditions formulate the abiotic background within which riverine biota have to adapt. In this concept habitat, suitability modelling provides information about the spatial requirements of specific organisms, usually by relating hydraulic parameters (e.g., depth, velocity) to a habitat suitability index (HSI) (Olsen et al. 2009). The weighted area by the composite HSI (a combined index including all the examined hydraulic parameters) within the entire domain of the hydrodynamic model corresponds to the weighted usable area (WUA) (Bovee et al. 1998), which is a well-known general indicator of habitat quality and quantity.

Although in situ studies of habitat selection can be made, to some extent, in a river reach, still it is impossible to cover all likely spatial and temporal heterogeneities (Heggenes et al. 1996). To overcome these limitations, prediction of ecological response to flow changes are made by hydraulic modelling integrated with habitat suitability modelling.

The hydraulic characteristics of the river are simulated by either one-dimensional (1D) (García et al. 2011) or two-dimensional (2D) hydrodynamic models (Leclerc et al. 1995) (in Benjankar et al. 2015). Flow properties in 1D models are estimated based on the physical characteristic of the cross sections (e.g., mean depth, average velocity). In some cases, it is possible to divide each cross section into sub-areas simulating Pseudo-2D situations. Regarding 2D models, water depths and velocities are calculated across a grid or mesh which defines the topographic information.

The selection of a 1D or 2D hydrodynamic model depends on the complexity of the river reach. 1D models are usually applied for long river lengths and 2D models are used over shorter representative river reaches (Katopodis 2012). Previous studies have shown that 1D and 2D models can provide comparable cross-sectional-averaged flow properties in simple uniform channels; nevertheless, different flow proprieties may be predicted in morphologically complex channels (Brown & Pasternack 2014). More specifically, complex flow patterns cannot be easily represented from cross-sectional-averaged properties predicted by 1D models (Mason et al. 2003), while 2D models predict depth, magnitude and direction (X,Y) of mean vertical velocity at points providing better habitat metrics (Bovee 1996).

Even though there are several sources of uncertainty incorporated into the modelling process including the input data, the required parameters, the structure of the models and the propagation of the uncertainty between several sub-models (Deletic et al. 2012), it is noted that the model structure is the most crucial to the model results (Refsgaard et al. 2006). In this study, a comparative assessment of the habitat quality spatial distribution, in terms of WUA, as simulated by a Pseudo-2D and a 2D model, has been attempted, focusing on the differences among model structures.

We used a Mediterranean mountain reach, with limited anthropogenic disturbance and we examined several discharges in both Pseudo-2D and actual 2D models, with respect to the summer conditions of the study area. To analyse the differences in spatial distribution we used the error matrix which quantifies models’ differences (Benjankar et al. 2015). This analysis was performed on a cell-by-cell basis throughout the Pseudo-2D model inundated domain.

Habitat duration curves were constructed following the procedures described within the Instream Flow Incremental Methodology for environmental flow studies (Bovee et al. 1998) indicating the exceedance probability for the habitat area in the corresponding models (Pseudo-2D and 2D), with combined HSI higher than 0.5. The habitat duration curves were estimated based on the hydrological simulated summer flows (June to October) for the study reach, according to previous work (Papadaki et al. 2016).

The study was carried out in the upper part of Acheloos River located in the central western mountainous region of Greece. The mean summer (June–October) discharge of the upper Acheloos basin is 3.11 m³ s⁻¹. A typical native species of the basin is the Balkan trout (Salmo farioides) (hereafter W.B. trout). However, this species is now scarce, especially in the lower part of the Acheloos basin as a result of severe overfishing, even involving illegal spear fishing, electrofishing and habitat alteration.

Habitat mapping of a 1.5 km river stretch of the upper course of the Acheloos River (at 670 m a.s.l., 39.479443°, 21.326510°, WGS 84) was carried out in order to select a representative study site of 390 m (Figure 1). More specifically, several types of hydromorphological units (i.e., pools, runs, riffles), were identified, and their extent and physical attributes measured. Measurements of water depth and velocity at seven cross sections (Figure 1) perpendicular to the flow provided data for the validation of both models. Depth (m) was measured with a wading rod to the nearest cm and the mean flow velocity of the water column (hereafter velocity (m s⁻¹)) was measured with a propeller current meter (OTT^®).

The assessment of the flow requirements for the adult W.B. trout (>20 cm) was made by combining hydraulic simulation and habitat suitability modelling.

Riverbed topography was surveyed encompassing the main channel and banks with a GPS/GNSS Geomax Zenith 20 using geodesic references (i.e., GGRS ‘87 Greek Geodetic Reference System) to generate a digital elevation model as the base for the models. The fish microhabitat-use survey, as part of the habitat simulation method was conducted during summer 2014 in the Voidomatis River; a reference river with near-natural conditions within Greece's Northern Pindos National Park (Papadaki et al. 2016).

HEC-RAS (Version 4.1) was used to perform a Pseudo-2D hydrodynamic simulation for 14 flow scenarios with respect to monthly summer flows in the study reach. The model solves the 1D Saint-Venant equations, for steady state and gradually varied flow. In a HEC-RAS steady state simulation, Water Surface Elevation (WSE) profiles were computed from one cross section to the next using the standard step iterative procedure to solve the energy equation.

To account for the Pseudo-2D hydraulic simulation every cross section was subdivided into 12 cells both in the main channel and the overbank area. The number of cells was primarily a function of the number of water velocity measurements and substrate variation along the cross sections. For every cross section a single water stage was simulated by the standard solving procedure, while velocities were calculated separately for each cell for the simulated water stage (procedure described in details in HEC-RAS v4.1; Hydraulic Reference Manual 2010).

Pre-process of the geometry of the reach was carried out using Autodesk Civil 3D software. Simulations were performed at 27 cross sections along the river reach. Friction losses were calculated applying Manning's roughness coefficient to every cross section. Furthermore, roughness coefficient was horizontally varied to account for substrate variation. Manning's n initial values varied for the main channel between 0.023 and 0.044 while for the overbank areas between 0.055 and 0.08.

The FLOW-R2D model is an in-house numerical model (Tsakiris & Bellos 2014), which solves the fully dynamic form of the 2D shallow water equations (2D-SWE) through the finite difference method and the McCormack numerical scheme (McCormack 1969), in a cell-centred, non-staggered computational grid. The McCormack numerical scheme is explicit and therefore is stable under the Courant–Friedrichs–Lewy (CFL) condition.

The non-rectangular computational domains were approached with a pseudo-computational rectangular domain which encloses the real one. For the solid boundaries, a modified reflection boundary was used (Bellos & Tsakiris 2015a).

The wet/dry bed modelling was achieved using a water depth threshold (h_dry) which distinguishes wet and dry cells. As well, artificial viscosity was added in the context of the 2D-SWE discretization, through a diffusion factor (ω). For the friction modelling, the Manning empirical model was used, while several applications of the model for both experimental benchmark tests and real world applications are also available (Bellos & Tsakiris 2015a, 2015b).

In previous work, Pseudo-2D model roughness calibration was made in the upper section of the Acheloos River (Papadaki et al. 2016). Normal depth was used as downstream boundary condition, whereas the flow was considered as subcritical in the entire computational domain and therefore no specific upstream boundary condition was set. The energy slope required for the downstream boundary condition was assumed to be equal with the surface slope. Calibration of the model was conducted by adjusting Manning's coefficient (n) at seven critical cross sections for the discharge of 4 m³ s⁻¹, and by comparing graphically both the simulated water stage and velocities with the observed ones from field measurements. The adjusted Manning's n values through calibration were close enough to the starting, bibliographical values (Cowan 1956; Chow 1959; Barnes 1967) with few exceptions (e.g., n = 0.012) which were necessary in order to achieve velocity convergence with the measurements in some cross sections. This low roughness value was assigned in particular in two adjacent cross sections and found to be necessary due to their main channel complex geometry.

For the numerical simulation using the 2D model, a grid of 46,957 square cells with 1 × 1 m size was generated, in order to represent the computational domain. The lateral boundaries of the river were considered solid, in order to preserve the water mass balance. For the downstream boundaries, the open boundary condition was used. The wet/dry threshold was determined as h_dry = 1 cm, and the Courant number as CFL= 0.1, according to the authors’ experience in similar real-world case studies (Tsakiris & Bellos 2014; Bellos & Tsakiris 2015b, 2016). For the friction modelling, Manning's equation was used.

A slightly different procedure was used to calibrate the required parameters of the 2D model, in relation to the Pseudo-2D. However, conventional or evolution optimization methods cannot be implemented for the calibration phase, due to the computational burden of the 2D model. The parameters which were calibrated were the diffusion factor ω and the Manning's coefficients n, for the friction modelling. For the calibration phase, water depth and flow velocity measurements were used, derived from 100 points located across seven cross sections (Figure 1) of the study reach.

The trial and error method was used in order to calibrate the diffusion factor and the Manning's coefficient in the four friction zones of the computational domain (sand, gravel, cobble and boulder zone). Specifically, 32 scenarios were implemented, in which the diffusion factor ranged from 0.90 to 0.95, Manning's coefficient at the sand zone ranged from 0.03 to 0.04 sm^−1/3, at the gravel zone ranged from 0.04 to 0.06, at the cobble zone ranged from 0.05 to 0.07 sm^−1/3 and at the boulder zone ranged from 0.06 to 0.09 sm^−1/3. In a rough approximation, the diffusion factor simulates the eddy viscosities of the flow in all the scales (reach and sub-grid scale) (Tsakiris & Bellos 2014). Based on this approximation, it is generally recommended from the developers of the model, that in real-world conditions the diffusion factor takes the value ω = 0.90 (Bellos & Tsakiris 2016). However, when turbulence phenomena in reach or sub-grid scale are expected to be relatively small, this factor could be increased. The grid size of the present case study, which is relatively small, combining with the relatively low flow velocity values indicate that the interval of the diffusion factor could range from 0.90 to 0.95, according to the authors’ experience. For the intervals of Manning coefficients, the ranges were also predefined according to the authors’ experience (Bellos & Tsakiris 2015b, 2016).

In the next step of the calibration phase, the sum of the root mean square error (RMSE) was calculated, for both water depths and flow velocities at the 100 points where measurements exist, between the numerical results and the observed data. The above values were normalized through a utility function, in which the maximum sum of the RMSE took the value 0 and the minimum took the value 1. Finally, the normalized values of the water depths and the flow velocities were summed in order to determine the performance of the scenarios. The maximum score of the summed normalized values defined the best scenario. According to this procedure, the calibrated parameters took the following values: diffusion factor ω = 0.95, Manning's coefficient for the sand zone n = 0.03 sm^−1/3, for the gravel zone n = 0.06 sm^−1/3, for the cobble zone n = 0.07 sm^−1/3 and for the boulder zone n = 0.06 sm^−1/3. It is noted that the corresponding summed normalized value was also calculated for the results derived from the calibration phase of the Pseudo-2D model and it was found that the score obtained using the 2D model was better. For the upstream boundaries, the steady flow condition was used, determined by Manning's equation. The required parameters for these boundaries were the water elevation stage and the effective slope of the upstream cross section. The above parameters were optimized based on the flow inlet at the computational field which was defined at 4 m³ s⁻¹ (equal to the corresponding measured discharge), and the minimum RMSE between the observed data and the numerical results across the cross section for both water depth and flow velocity values. As previously mentioned, in the context of the calibration phase, two values for the Manning's coefficient were considered for each friction zone. The upstream cross section was exclusively located in the boulder zone and therefore two values for the effective slope and the water elevation stage were derived. Specifically, the effective slope values were determined as 0.00541 and 0.01224 for the Manning's coefficient 0.06 and 0.09 sm^−1/3, respectively, (boulder zone), whereas the water elevation stage was derived as 668.01 m a.s.l., in both of these cases.

Physical habitat was quantified using depth and velocity univariate HSI curves (Figure 2), developed according to Bovee (1986), representing generalized suitability for the W.B. trout; 103 fish adult (size >20 cm). The HSI curves relate the hydraulic variables (i.e., depth or velocity) with a suitability index (SI), ranging from 0 (unsuitable for the aquatic species) to 1 (excellent). The intermediate values represent the suitable range based on a specified hydraulic variable (i.e., depth or velocity) (Bovee et al. 1998). These two individuals’ suitability indices were then combined to form a composite SI for every cell of the hydraulic models, using the product mathematical operation (Vadas & Orth 2001). At the cell scale, Cell Suitability Index (CSI) indicates whether the physical parameters (depth and velocity in this case) are within those required by individual species and particular life stages (Bovee 1978).

WUA was calculated by multiplying the CSI for a given cell by the area assigned to that specific cell, which in this case is 1 m for both models. In order to study only the suitable conditions for the target species, WUA was estimated considering the cells with CSI higher than 0.5 only (hereafter WUA0.5). The whole procedure was carried out in R software (R Development Core Team 2014 R: A language and environment for statistical computing).

A comparison between the spatially distributed CSI, calculated from the Pseudo-2D and the 2D model using the error matrix method (Congalton & Green 2008) was made as a standard technique for quantifying the accuracy among maps, specifically designed for raster comparisons. The error matrix compares maps by calculating overall accuracy (OA) and the agreement index between the maps using the kappa statistic (K). A K value of 1 indicates perfect agreement, whereas a K value of 0 indicates agreement equivalent to chance. The OA is a ratio between the numbers of correctly predicted cells to total number of cells considered in the analysis. CSIs were separated into classes of 0 = 0, 0–0.19 = 1, 0.20–0.39 = 2, 0.4–0.59 = 3, 0.60–0.79 = 4 and 0.80–1.00 = 5 in order to estimate the K statistic using the error matrix.

The models’ performance indices indicated a relatively better performance of the FLOW-R2D model (RMSE: 0.11 and r: 0.49) in relation to HEC-RAS (RMSE: 0.12 and r: 0.36) for the depth simulations (Tables 1 and 2, Figure 1). The velocity model outputs were simulated with less accuracy by both models but HEC-RAS (RMSE: 33 and r: 0.39) illustrated a slightly better performance than FLOW-R2D (RMSE: 33 and r: 0.37).

The statistical comparisons between observations and simulated values indicated that both models underestimate the actual depth and velocity fluctuations while HEC-RAS presents a quite higher standard deviation in relation to the FLOW-R2D (Table 1, Figure 3). Moreover, both models have a significant degree of agreement in the depth simulated values (r: 0.95) while in the respective velocity values the agreement is quite lower (r: 0.65) with HEC-RAS presenting a very wide value distribution (Table 2, Figure 3).

The relationship between observed and simulated water velocities and depths derived from the Pseudo-2D and the 2D models (FLOW-R2D) indicated that both models do not capture satisfactorily the variance of the observations since all linear R² values are relatively low (Figure 3). Nevertheless, the FLOW-R2D illustrates a slightly better performance for depth with a R² value of 0.24.

The comparison of the observed and estimated depths and velocities in four characteristic cross sections for the two models endorsed the aforementioned respective results. In most cases the depth is relatively accurately estimated by both models while the velocity illustrates significant discrepancies (Figure 4). In the particular four cross-sections, HEC-RAS illustrated a slightly better performance regarding velocity in comparison to the 2D model since its estimated values follow more closely the observations’ fluctuations (Figure 4).

From the spatial distribution of the model outputs, it was observed that relatively small differences in the simulated water depths derived from the Pseudo-2D and the 2D model exist in all discharge scenarios (Figure 5). The most distinctive differences between the two models were in the average depth values (0.6–1 m) that were more spatially extended in the FLOW-R2D outputs than in HEC-RAS. Moreover, these differences became more pronounced as the discharge increased from 2 to 6 m³ s⁻¹ (Figure 5).

However, the estimated flow velocities between the two models had greater differences (Figure 6). HEC-RAS indicated significantly higher velocities in the discharge scenarios of 2 and 4 m³ s⁻¹ in relation to the FLOW-R2D while in the discharge scenario of 6 m³ s⁻¹ the differences between the models were much lower (Figure 6). The study reach has an inherent behaviour of a run-type river with turbulence. This can be justified by the fact that flow velocity varies rapidly in magnitude and direction, along space and time, which underlines the capability of a 2D model to simulate complex flow distribution that a Pseudo-2D model cannot.

After transforming the models’ estimated velocities and depths for each discharge scenario into WUAs by using the West Balkan adult trout CSI, a comparison diagram was created (Figure 7). The results of the two models’ habitat analysis show that the Pseudo-2D (HEC-RAS) model underestimates WUA_0.5 for discharges over 3 m s⁻¹ in comparison with the 2D model, while the opposite occurs for lower discharges (Figure 7). Moreover, following the peak WUA values (Figure 7), the estimated best minimum ecological flow differs significantly between the two models with HEC-RAS indicating a value close to 2 m³ s⁻¹ and FLOW-R2D a value close to 6 m³ s⁻¹. This is a very big difference from a water management perspective and further in-depth investigation should be performed in order to identify the potential causes.

In Figure 8, habitat duration curves for summer conditions (June–October) for the years 1986 to 2004 are presented for both models. The discharge range was from 4 m³ s⁻¹ to 45 m³ s⁻¹. The results indicate that 30% of the time, habitat area according to the 2D model is equalled or exceeded in comparison with the Pseudo-2D model (HEC-RAS) results. In contrast, 70% of the time the Pseudo-2D model (HEC-RAS) indicated higher habitat area than the FLOW-R2D model.

Based on a graphical comparison (Figure 9) of the spatially distributed CSI, the largest differences were near the channel banks, where the 2D model generally predicted higher CSI classes than those predicted with the Pseudo-2D model. Again, the differences between the two models are higher in the low discharge scenario (2 m³ s⁻¹) and become smaller in the higher discharge scenario (6 m³ s⁻¹).

An evaluation of the spatial similarities between raster maps of the two models was based on pixel-by-pixel comparison technique. The kappa statistic defined a similarity measure between the Pseudo-2D model and the 2D model for 14 discharge scenarios in respect to summer conditions. Heat maps of four discharge scenarios (2, 4, 6, 8 m³ s⁻¹) are presented in Figure 10 showing the kappa statistics results. Generally, kappa coefficient showed low agreement between the two models for all the examined discharges except from discharges 2 and 6 m³ s⁻¹ where kappa coefficient indicated higher similarity.

In general, there is low agreement between the two models in most of the CSI values categories in scenarios 4 and 8 m³ s⁻¹ (Figure 10), while in scenarios 2 and 6 m³ s⁻¹ there is a significantly better agreement especially for the medium and low CSI values. FLOW-R2D indicated more pixels with low CSI values for all discharge scenarios in relation to the HEC-RAS, while in discharges 2 and 6 m³ s⁻¹ an extension of the area with the larger HEC-RAS number of pixels was illustrated for the medium CSI values (Figure 10).

In the area of high CSI values, both models indicated a similar number of pixels in the model domain. The percentage of agreement between the two models (pixels that indicated CSI values in the same category) followed a decreasing trend with increasing discharge, ranging from 73% in the discharge scenario of 2 m³ s⁻¹ down to 54% in the scenario of 8 m³ s⁻¹.

The goal of instream habitat models is to formulate relationships between the area of suitable habitat and discharge. Nevertheless, data requirements for these models are very demanding in terms of biological, geomorphological and hydrological data and accuracy. Several studies in the past concluded that the geomorphology of the river reach is one of the most important factors for the output of hydrodynamic models (Tarbet & Hardy 1996; Gard 2009; Boyraz & Kazezyılmaz-Alhan 2014; Benjankar et al. 2015; Greco 2015) which influence the output of the habitat models. The Mesochora reach resembles a relatively deep run with high velocities and medium depths creating complex flow patterns and for this reason the two models’ outputs seem to induce differences on the WUA results.

Even though, in the past, the computational power requirements were possible constraints for the application of 2D models, nowadays they can be applied on desktop PCs (Tonina & McKean 2010; Pasternack & Senter 2011; Tonina et al. 2011), and they are expected to become the most common tool in aquatic habitat modelling (Tonina & Jorde 2013). 1D models are still very useful as they are computationally efficient and they can simulate first-order conditions over much larger stream domains and over much longer periods than 2D models (Benjankar 2009; Burke et al. 2009).

Results deriving from our study highlight the fact that 1D models’ performance in a variety of flow and channel conditions requires establishment before using them broadly for purposes of habitat prediction. 2D models’ meshes must be capable of reproducing the spatial flow patterns created by the topographic features at the resolution important to the aquatic organisms under study according to Crowder & Diplas (2000).

For the specific case study, significant discrepancies were observed for the results between the Pseudo-2D model and the 2D hydrodynamic model, as far as the hydraulic variables are concerned (water depth and flow velocities). It seems that the structure of the models has a significant impact on the results derived from the combined use of hydraulic/habitat model. More specifically, differences between the Pseudo-2D and the 2D model have been examined and found to be larger for velocities and significantly lower for depth. The majority of velocity differences were observed along the water edges and channel banks. This probably occurred because the 2D model incorporated the transverse velocity component and it was able to estimate secondary flows created at the flow boundaries, however this statement should be further investigated in order to be verified. In general, the 2D hydrodynamic model was able to simulate the local complexities of the flow due to the microtopography effects and the spatial distribution of roughness, using a sophisticated, physical-based approach, whereas the Pseudo-2D hydrodynamic model approaches these complexities by implementing interpolations between the derived results to each cross section. However, the results of the two models indicate the usage of the Pseudo-2D approach as a reliable alternative for habitat predictions.

There is no specific pattern with which one model underestimates or overestimates the hydraulic variables in comparison with the other, and thus extensive future research should focus on the determination of a methodological framework to quantify the uncertainties created from models’ structure. It seems that as far as the water depths is concerned, the differences between the Pseudo-2D and 2D results were smaller (from ±0.30 m), with a decreasing trend as the discharge increases. Regarding flow velocities, the models’ outputs indicated higher differences (from ±0.60 m s⁻¹), but no significant trend was observed in relation to the discharge. It should be noted that the local extreme differences which were observed were excluded from the analysis since they were probably caused by various inherent errors and biases, such as the possible numerical oscillations in 2D modelling and the interpolations implemented for the Pseudo-2D modelling.

Further research should be conducted to understand the mechanisms associated with the biological responses from the hydrodynamic behaviour of rivers. More parameters, such as the substrate, the cover, the temperature and the food availability should be incorporated in future similar works in order to adapt a holistic approach in the ecological flow estimations.