Modelling of Pool-Type Fishways Flows: Efﬁciency and Scale Effects Assessment

: In this study, the 3D numerical modelling of ﬂow in a pool-type ﬁshway with bottom oriﬁces was performed using computational ﬂuid dynamics (CFD) software (FLOW-3D ® ). Numerical results were compared with experimental data obtained from acoustic Doppler velocimetry (ADV) and particle image velocimetry (PIV) measurements. Several hydrodynamic variables that inﬂuence ﬁshways efﬁciencies, such as ﬂow depths, ﬂow patterns, water velocity, turbulent kinetic energy, Reynolds normal stresses, and Reynolds shear stress parallel to the bottom component, were qual-itatively and quantitatively compared. The numerical model accurately reproduced the complex ﬂow ﬁeld, showing an overall good agreement between the numerical model predictions and the experimental data for the analysed variables. The importance of performing a numerical model validation for all the parameters under analyses was highlighted. Additionally, scaling effects were analysed by running an upscaled numerical model of the prototype ﬁshway. The model performed with similar accuracy for both physical model and prototype dimensions with no evidence of scale effects. The present study concludes that CFD models (namely FLOW-3D ® ) may be used as an adequate and efﬁcient design and analysis tool for new pool-type ﬁshways geometries, reducing and complementing physical model testing.


Introduction
Restoring the longitudinal connectivity of rivers remains a key issue in the recovery of freshwater ecosystems [1,2]. If well designed and constructed, fishways provide a path that allows fish to continue migrating past dams and weirs. In a review of fish passage efficiency, Noonan et al. [3] found that the design features of many existing fishways did not meet the needs of fish species adequately, although pool-type fishways presented the highest efficiency for all fish groups.
The flow field in pool-type fishways with bottom orifices is highly three-dimensional and much more complex than the flow field in vertical slot fishways (VSF), which is the design more frequently considered in previous studies on fishways numerical model validation [26][27][28][29]35]. To the authors' best knowledge, this is the first CFD study on pooltype fishways with bottom orifices, which also includes the most extensive comparison between experimental velocity data and three-dimensional numerical modelling results on pool-type fishways published. Two different measurement techniques (PIV and ADV), were used, thus allowing a detailed comparison and providing confidence in the results on CFD simulations for this type of flow field.
The study presents not only a qualitative comparison but also a detailed quantitative comparison, with statistical testing on the agreement between the numerical model results and measurements, including turbulence parameters, which was not presented in previous numerical model studies of other fishway types. Scale effects are addressed as well. Hence, this study will smooth the validation of CFD models of pool-type fishways, the most used worldwide [10], and it will hopefully encourage its use by the designers.
Furthermore, the use of CFD models (namely FLOW 3D ® ) as a design tool for new pool-type fishways geometries is discussed.

Experimental Setup and Procedure
The facility where the experimental work was undertaken is located at the Laboratory of Hydraulics and Environment of Instituto Superior Técnico (IST, Lisboa, Portugal). It is a 1:2.5 scaled fishway of a prototype size indoor pool-type fishway located at the National Laboratory for Civil Engineering (LNEC, Lisboa, Portugal). It includes a 5.7 m long, 0.4 m wide, and 0.5 m deep flume, with an adjustable slope, equipped with a recirculation circuit and a pump frequency converter. An electromagnetic flow meter, installed in the recirculation circuit, measures the discharge that is controlled by the pump frequency converter. A sluice valve controls the downstream flow depth. The flume has glass side walls, enabling flow visualization and PIV measurements ( Figure 1). type fishways with bottom orifices flow, two different sized numerical models were developed: a large-sized model with the prototype dimensions and a scaled small-sized one with the physical model dimensions.
The flow field in pool-type fishways with bottom orifices is highly three-dimensional and much more complex than the flow field in vertical slot fishways (VSF), which is the design more frequently considered in previous studies on fishways numerical model validation [26][27][28][29]35]. To the authors' best knowledge, this is the first CFD study on pool-type fishways with bottom orifices, which also includes the most extensive comparison between experimental velocity data and three-dimensional numerical modelling results on pool-type fishways published. Two different measurement techniques (PIV and ADV), were used, thus allowing a detailed comparison and providing confidence in the results on CFD simulations for this type of flow field.
The study presents not only a qualitative comparison but also a detailed quantitative comparison, with statistical testing on the agreement between the numerical model results and measurements, including turbulence parameters, which was not presented in previous numerical model studies of other fishway types. Scale effects are addressed as well. Hence, this study will smooth the validation of CFD models of pool-type fishways, the most used worldwide [10], and it will hopefully encourage its use by the designers.
Furthermore, the use of CFD models (namely FLOW 3D ® ) as a design tool for new pool-type fishways geometries is discussed.

Experimental Setup and Procedure
The facility where the experimental work was undertaken is located at the Laboratory of Hydraulics and Environment of Instituto Superior Técnico (IST, Lisboa, Portugal). It is a 1:2.5 scaled fishway of a prototype size indoor pool-type fishway located at the National Laboratory for Civil Engineering (LNEC, Lisboa, Portugal). It includes a 5.7 m long, 0.4 m wide, and 0.5 m deep flume, with an adjustable slope, equipped with a recirculation circuit and a pump frequency converter. An electromagnetic flow meter, installed in the recirculation circuit, measures the discharge that is controlled by the pump frequency converter. A sluice valve controls the downstream flow depth. The flume has glass side walls, enabling flow visualization and PIV measurements ( Five PVC cross-walls equipped with bottom orifices of 0.08 × 0.08 m 2 were used to create four fishway pools, each 0.76 m long, 0.40 m wide, and 0.50 m deep. The orifices were created using a thin (2 mm) perspex plate, and the edges were sharp-crested. Consecutive orifices were positioned on opposite sides of the cross-walls (Figure 1b), creating a sinusoidal flow path.
The flume was set at a slope of 8.5%, which is within the range of slopes commonly used in prototype fishways [48] and used in previous experiments with fish in the near- Five PVC cross-walls equipped with bottom orifices of 0.08 × 0.08 m 2 were used to create four fishway pools, each 0.76 m long, 0.40 m wide, and 0.50 m deep. The orifices were created using a thin (2 mm) perspex plate, and the edges were sharp-crested. Consecutive orifices were positioned on opposite sides of the cross-walls (Figure 1b), creating a sinusoidal flow path.
The flume was set at a slope of 8.5%, which is within the range of slopes commonly used in prototype fishways [48] and used in previous experiments with fish in the nearprototype fishway previously mentioned [7]. The discharge was set to 4.4 ± 0.4% Ls −1 , the head drop between pools (∆h) was 0.064 ± 5% m, and the pool mean water depth (h m ) was 0.352 ± 0.5% m. These values were scaled from the experiments performed by Silva et al. [7] and were controlled and maintained constant during all the experiments. The PIV and ADV measurements were carried out in the third pool in six planes: three planes parallel to the flume bottom at 11.4% h m (0.04 m-half the height of the orifice), 25% h m , and 50% h m , and three vertical planes parallel to the sidewalls at the midwidth of the pool and of each orifice ( Figure 2). The planes of the instantaneous velocities' measurement were chosen to provide an adequate characterization of the flow field in the range 0-50% h m , which according to Silva et al. [7] is the region where fish were mostly found. In Quaresma et al. [49] a complete outline of the PIV and ADV measurements can be found. prototype fishway previously mentioned [7]. The discharge was set to 4.4 ± 0.4% Ls -1 , the head drop between pools (Δh) was 0.064 ± 5% m, and the pool mean water depth (hm) was 0.352 ± 0.5% m. These values were scaled from the experiments performed by Silva et al. [7] and were controlled and maintained constant during all the experiments. The PIV and ADV measurements were carried out in the third pool in six planes: three planes parallel to the flume bottom at 11.4% hm (0.04 m-half the height of the orifice), 25% hm, and 50% hm, and three vertical planes parallel to the sidewalls at the midwidth of the pool and of each orifice ( Figure 2). The planes of the instantaneous velocities' measurement were chosen to provide an adequate characterization of the flow field in the range 0-50% hm, which according to Silva et al. [7] is the region where fish were mostly found. In Quaresma et al. [49] a complete outline of the PIV and ADV measurements can be found. The coordinate system has its origin in the bottom right corner of the third pool upstream orifice, with the longitudinal x-axis parallel to the bottom and the flume axis, the y-axis transversal and with the z-axis perpendicular to the bottom. The instantaneous velocity components u, v, and w correspond to the x-, y-, and z-axis, respectively.

Numerical Model
The numerical simulations were performed using FLOW-3D ® . This software is a general-purpose CFD program based on the finite volumes method that solve the equations of motion for fluids to obtain transient, three-dimensional solutions to multiscale and multiphysics flow problems [50]. In FLOW-3D ® , the computational domain is subdivided using Cartesian coordinates into a structured mesh grid of variable-sized hexahedral cells. The grid is defined independently of the geometry and, subsequently, the geometry is embedded in the grid by the Fraction Area/Volume Obstacle Representation (FAVOR TM ) technique [51]. Flow Science [50] presents additional details regarding the theoretical and numerical fundamentals of FLOW-3D ® . FLOW-3D ® has been used to model flows in fishways, hydraulic structures, and river reaches and bends (e.g., [45,[52][53][54][55][56].
In this study, the simulations were performed using the implicit solver generalized minimum residual method, which is highly accurate and efficient for a wide range of problems [50]. Several options of the software were investigated, namely the choice of the The coordinate system has its origin in the bottom right corner of the third pool upstream orifice, with the longitudinal x-axis parallel to the bottom and the flume axis, the y-axis transversal and with the z-axis perpendicular to the bottom. The instantaneous velocity components u, v, and w correspond to the x-, y-, and z-axis, respectively.

Numerical Model
The numerical simulations were performed using FLOW-3D ® . This software is a general-purpose CFD program based on the finite volumes method that solve the equations of motion for fluids to obtain transient, three-dimensional solutions to multiscale and multiphysics flow problems [50]. In FLOW-3D ® , the computational domain is subdivided using Cartesian coordinates into a structured mesh grid of variable-sized hexahedral cells. The grid is defined independently of the geometry and, subsequently, the geometry is embedded in the grid by the Fraction Area/Volume Obstacle Representation (FAVOR TM ) technique [51]. Flow Science [50] presents additional details regarding the theoretical and numerical fundamentals of FLOW-3D ® . FLOW-3D ® has been used to model flows in fishways, hydraulic structures, and river reaches and bends (e.g., [45,[52][53][54][55][56]).
In this study, the simulations were performed using the implicit solver generalized minimum residual method, which is highly accurate and efficient for a wide range of problems [50]. Several options of the software were investigated, namely the choice of the volume-of-fluid method, momentum advection method, turbulence model, and the roughness effect. To track the free surface, the VOF TM (volume-of-fluid) method [57] was used. The selected momentum advection method was the second-order monotonicity preserving method, which is based on a second-order, monotonicity-preserving upwinddifference method [58]. All the turbulence models available in FLOW-3D ® , namely, the Re-Normalization Group (RNG), k-ε, k-ω turbulence models and Large-Eddy Simulation (LES) model were investigated. The results presented in this paper were obtained with the LES model, which presented the best agreement with the experimental results.
The roughness effect was also investigated by testing a range of Manning's surface roughness coefficients that cover the experimental model materials possible values. The wall roughness was set using the roughness height computed from the Strickler formula. The Manning's surface roughness that led to the best agreement between simulated and measured discharges and flow depths was n = 0.01 sm −1/3 . The acceleration of gravity was applied in the negative z-direction and in the positive x-direction so that the x-axis was parallel to the flume bottom. The upstream and downstream boundaries were specified as pressure boundary conditions based on the water depths observed experimentally (0.33 m and 0.40 m, respectively). At the side and bottom no-slip wall boundaries were used and at the top, a symmetry boundary condition (zero value for normal velocity, zero gradients for the other quantities) was specified. In free surface simulations, FLOW-3D ® neglects the inertia of the gas adjacent to the liquid, and the volume occupied by the gas is replaced with an empty space, void of mass, represented only by uniform pressure and temperature. Therefore, free surface becomes one of the flow external boundaries. With this approach the computational effort is reduced, since the details of the gas motion are not important for the motion of the much denser liquid. To accurately capture the free-surface dynamics, the VOF TM method is employed. This method consists of three main components: (i) the definition of the volume of fluid function, (ii) a method to solve the VOF transport equation, and (iii) setting the boundary conditions at the free surface. The simulations were run for a number of time steps enough to attain a statistically stationary solution and obtain converged time-averaged values.
The fishway geometry was generated using AutoCAD 2014 (AutoDesk, San Rafael, CA, USA) based on the dimensions of the physical model and was imported into the code as a stereolithography (stl) file. The computational domain ( Figure 3) was discretized using a multiblock grid. The grid was constituted of five blocks; one contained the entire geometry with a uniform grid spacing of 0.02 m. The other four blocks were nested blocks with a uniform grid spacing of 0.01 m; one contained the third pool, and the other three blocks contained the remaining cross-walls. This cell size is in line with previous works, e.g., Fuentes-Pérez et al. [35] for a vertical slot fishway with a slot width of 0.281 m used a cell size of 0.03 m, approximately 10.7% of the slot width (in the present study, a cell size of 0.01 m represents 12.5% of the orifice width). A total of 380,650 grid cells were used to represent the fishway geometry in the coarser mesh hereafter, named mesh A model . volume-of-fluid method, momentum advection method, turbulence model, and the roughness effect. To track the free surface, the VOF TM (volume-of-fluid) method [57] was used. The selected momentum advection method was the second-order monotonicity preserving method, which is based on a second-order, monotonicity-preserving upwind-difference method [58]. All the turbulence models available in FLOW-3D ® , namely, the Re-Normalization Group (RNG), k-ε, k-ω turbulence models and Large-Eddy Simulation (LES) model were investigated. The results presented in this paper were obtained with the LES model, which presented the best agreement with the experimental results.
The roughness effect was also investigated by testing a range of Manning's surface roughness coefficients that cover the experimental model materials possible values. The wall roughness was set using the roughness height computed from the Strickler formula. The Manning's surface roughness that led to the best agreement between simulated and measured discharges and flow depths was n = 0.01 sm −1/3 . The acceleration of gravity was applied in the negative z-direction and in the positive x-direction so that the x-axis was parallel to the flume bottom. The upstream and downstream boundaries were specified as pressure boundary conditions based on the water depths observed experimentally (0.33 m and 0.40 m, respectively). At the side and bottom no-slip wall boundaries were used and at the top, a symmetry boundary condition (zero value for normal velocity, zero gradients for the other quantities) was specified. In free surface simulations, FLOW-3D ® neglects the inertia of the gas adjacent to the liquid, and the volume occupied by the gas is replaced with an empty space, void of mass, represented only by uniform pressure and temperature. Therefore, free surface becomes one of the flow external boundaries. With this approach the computational effort is reduced, since the details of the gas motion are not important for the motion of the much denser liquid. To accurately capture the freesurface dynamics, the VOF TM method is employed. This method consists of three main components: (i) the definition of the volume of fluid function, (ii) a method to solve the VOF transport equation, and (iii) setting the boundary conditions at the free surface. The simulations were run for a number of time steps enough to attain a statistically stationary solution and obtain converged time-averaged values.
The fishway geometry was generated using AutoCAD 2014 (AutoDesk, San Rafael, California, USA) based on the dimensions of the physical model and was imported into the code as a stereolithography (stl) file. The computational domain (Figure 3) was discretized using a multiblock grid. The grid was constituted of five blocks; one contained the entire geometry with a uniform grid spacing of 0.02 m. The other four blocks were nested blocks with a uniform grid spacing of 0.01 m; one contained the third pool, and the other three blocks contained the remaining cross-walls. This cell size is in line with previous works, e.g., Fuentes-Pérez et al. [35] for a vertical slot fishway with a slot width of 0.281 m used a cell size of 0.03 m, approximately 10.7% of the slot width (in the present study, a cell size of 0.01 m represents 12.5% of the orifice width). A total of 380,650 grid cells were used to represent the fishway geometry in the coarser mesh hereafter, named mesh Amodel.  Since the choice of the mesh element size is highly case-specific [55], to test the influence of the grid resolution on the results, a mesh sensitivity analysis was performed, based on the ASME criteria [55,59] with three more grids. One, with 2,714,000 grid cells, was obtained by considering a grid spacing of 0.01 m in the block that contained the entire geometry and a grid spacing of 0.005 m in the block that contained the third pool. This mesh is named mesh B model in this paper. Two more grids were considered using the grid overlay boundary condition available in FLOW-3D ® on the third pool mesh block. These grids had uniform grid spacings of 0.0025 m and 0.0020 m for mesh C model and D model , respectively.
To analyse the scale effects, a numerical model of a prototype size pool-type fishway (LNEC's flume) was used. This model (the prototype model) was geometrically similar to the IST flume model and the scale ratio of length was 2.5. Two grids were used, A prototype and B prototype . The computational meshes of the prototype and the physical model were also adjusted according to the geometric similarity to exclude the generation of numerical errors due to different scaled meshes. Table 1 presents the characteristics of the numerical meshes used in the present study and the computation time per 10 s of simulation time for each of the meshes. Table 1. Characteristics of the numerical meshes.

Grid Resolution and Quality Verification
LES models use a spatial filter to separate the turbulent flow field into two components: the larger and more energetic turbulent structures that are resolved by the numerical method on a given mesh (resolved scales); and the smaller structures that are not captured by the mesh (sub-grid scales). The influence of subgrid scales on resolved scales must be modelled [24]. In FLOW-3D ® , the subgrid scales' effects of turbulence are represented by an eddy viscosity, which is proportional to a length scale times a measure of velocity fluctuations on that scale. FLOW-3D ® uses the Smagorinsky method [60] which considers for the length scale the geometric mean of the grid cell dimensions (L = (∆x ∆y ∆z) 1/3 ). When using LES and implicit filtering, as the Smagorinsky method, since the filter size changes with the grid size and both the numerical discretization error and the sub-grid scale contributions are proportional to grid size, there is no grid-independent LES [61][62][63][64]. Considering this, Celik et al. [63] proposed the LES_IQ index to assess the quality of LES models. The LES_IQ k index is a measure of the percentage of the resolved turbulent kinetic energy (TKE) to the total. This index gives an indication of the resolution quality in LES models.
To apply the LES_IQ index it must be checked if the grids are in the asymptotic range. Therefore, the observed apparent order of the method p for the hydraulic parameters under analysis in this study was computed following the procedure outlined in Celik et al. [60]. Table 2 shows the average apparent order of the method, p, and the percentage of points showing oscillatory convergence for the ADV measurement grid points locations. The average observed apparent order of the method for all parameter under analysis is close to 2, which shows a good agreement with the formal order of the scheme used and indicates that the grids are within the asymptotic range [59]. The high percentage of points showing oscillatory convergence might be due to the fact that this flow field presents several recirculation regions [49]. Celik et al. [65] showed that the oscillatory convergence behaviour in finite difference solutions may be caused by the oscillatory velocity field that occurs in recirculation regions. This topic will be further addressed in the Appendix A.
Given the high number of points showing oscillatory behaviour, the modified version of LES_IQ k index [62] was used to assess the resolution quality of the grids used and was computed in all ADV measurement grid points locations. Mesh A model has an average LES_IQ k index of 0.72 with 45% of the points with values higher than 0.80. When using the theoretical order of the scheme (p = 2) instead of the observed apparent order of the method Mesh A model , the average LES_IQ k index increases to 0.75, with 49% of the points with values higher than 0.8. The average LES_IQ k index of the Mesh B model is 0.90, with 89% of the points presenting values higher than 0.80. When using the theoretical order of the scheme (p = 2), instead of the observed apparent order of the method Mesh B model , the average LES_IQ k index increases to 0.93, with 98% of the points with values higher than 0.80.
According to Pope [66], a good LES should resolve at least 80% of the TKE. Celik et al. [62] consider that a LES_IQ k of 75 to 85% can already be considered adequate for most engineering applications that typically occur at high Reynolds numbers. This flow field is highly complex and turbulent with a rapid mixing. Considering the mean velocity in the pool and the pool mean water depth, the Reynolds number is higher than 5 × 10 5 ; thus, a LES_IQ k of 75 to 85% can already be considered adequate. The values obtained for Mesh A model and B model indicate that sufficient grid resolutions were used. Therefore, the simulations can be qualified as LES. Mesh B model , with an average LES_IQ k index of 0.90, can already be considered a good LES simulation. It should, however, be emphasized that this index is a verification index which only assesses grid resolution quality. To assess model accuracy, a comparison with experimental data is still necessary.

Results
The numerical model accurately reproduced the experimental data with an average difference of 3.5% (3.6% for Mesh A model and 3.4 % for Mesh B model ) for the flow depths and 5% (5.4% for Mesh A model and 5.2% for Mesh B model ) for the discharge.
Subsequently, the numerical results are presented for the third pool (the ADV and PIV measurement pool) for several relevant hydraulic parameters, such as velocities, TKE, Reynolds stresses, and streamlines and are compared with the ADV and PIV measurements. The numerical results obtained for the ADV measurement grid points were compared with the corresponding ADV measurements. Concerning the PIV measurements, since the spatial resolution was higher than the used numerical meshes, the PIV data were interpolated within the mesh cells to obtain the values to be compared with the numerical results. Table 3 presents the values of the mean absolute difference (MAD), the coefficient of determination (R 2 ), and refined index of agreement (d r ) for both meshes of several relevant hydraulic parameters, namely time-averaged velocity components (u, v, and w), time-averaged velocity magnitude (U), Reynolds normal stresses u 2 , v 2 , and w 2 , turbulent kinetic energy (κ), and Reynolds shear stress component parallel to the bottom (τ uv ). Table 3. Comparison between numerical model results and ADV and particle image velocimetry (PIV) measurements. Since R 2 , the most widely used adjustment assessment coefficient, is oversensitive to outliers and insensitive to additive and proportional differences between model predictions and observations [67,68], d r (refined index of agreement) and MAD (mean absolute difference) were included in the analysis.

ADV-A model ADV-B model PIV-A model PIV-B model
where a CFD , a ADV are the hydraulic parameter under analysis obtained in each acquisition point with CFD and the experimental measurement techniques (PIV and ADV) respectively, and n is the number of compared points.
where a Exp is the arithmetic mean of the hydraulic parameter under analysis obtained with the experimental measurement techniques (PIV and ADV). The refined index of agreement is a refined statistical index of model performance proposed by Willmott et al. [69], which ranges from −1 to 1. Higher values of R 2 and d r indicate better agreement, with 1 indicating a perfect agreement between experimental data and numerical model predictions.
Meshes A model and B model showed similar levels of agreement with the experimental data (Table 3) for velocity (u, v, w, and U), turbulent kinetic energy (κ), longitudinal Reynolds normal stress u 2 , and Reynolds shear stress component parallel to the bottom (τ uv ). However, for the transversal and vertical Reynolds normal stresses v 2 and w 2 , a significant improvement in the agreement was achieved when doubling the number of cells in each direction. Therefore, mesh B model (the finer mesh) should be used when studying these turbulence parameters. Refining the mesh furthermore did not improve these results as shown in Appendix A. Figure 4 shows the streamlines computed for the time-averaged velocities obtained with the PIV and numerical model mesh A model in planes 2 and 5. The results from the numerical model are presented only for mesh A model since no significant differences could be found between mesh A model and B model results.
The flow patterns were well predicted by the numerical model, as shown in Figure 4. In plane 2, the jet flow trajectory and width were well captured, and the recirculation zones' location and size were in good agreement in both planes. Figure 5 shows the longitudinal variation in the time-averaged velocity components u, v, and w at two locations. One in the centre of the downstream orifice at planes 1 and 6 intersection (y = 0.36 m and z = 0.04 m), and the other one at the flume centre at planes 2 and 5 intersection (y = 0.20 m and z = 0.088 m). These profiles were chosen because of the high probability of fish being present at these regions, with the first location going through the entrance jet of the pool (considering fish moving upstream) and the second one going through the recirculation zones.
A very good agreement between the experimental and numerical results can be observed in Figure 5 for all the velocity components with the largest differences, that are still acceptable, occurring for values close to 0 of the transversal (Figure 5c) and vertical velocity components (Figure 5e,f). Figure 6 shows the longitudinal variation in the Reynolds normal stress components u 2 , v 2 , and w 2 and of the Reynolds shear stress component parallel to the bottom (τ uv ) at the same locations of Figure 5.  The flow patterns were well predicted by the numerical model, as shown in Figure  4. In plane 2, the jet flow trajectory and width were well captured, and the recirculation zones' location and size were in good agreement in both planes. Figure 5 shows the longitudinal variation in the time-averaged velocity components u , v , and w at two locations. One in the centre of the downstream orifice at planes 1 and 6 intersection (y = 0.36 m and z = 0.04 m), and the other one at the flume centre at planes 2 and 5 intersection (y = 0.20 m and z = 0.088 m). These profiles were chosen because of the high probability of fish being present at these regions, with the first location going through the entrance jet of the pool (considering fish moving upstream) and the second one going through the recirculation zones.
A very good agreement between the experimental and numerical results can be observed in Figure 5 for all the velocity components with the largest differences, that are still acceptable, occurring for values close to 0 of the transversal (Figure 5c) and vertical velocity components (Figure 5e,f).  Figure 6 shows the longitudinal variation in the Reynolds normal stress components u ' 2 , v ' 2 , and w ' 2 and of the Reynolds shear stress component parallel to the bottom (τuv) at the same locations of Figure 5.
Additionally, a good agreement between the experimental and numerical results was

Scale Effects
To analyse the scale effects, the prototype scale numerical model results were analysed and compared to the experimental data, scaled according to the Froude similarity.
The flow depths and discharges differences between the numerical model and the experiments slightly decreased for the prototype scale with average values of 3.4 and 4.8% for mesh Aprototype and 2.3 and 4.1% for mesh Bprototype, for flow depths and discharges, respectively. Table 4 shows the mean absolute difference (MAD), the coefficient of determination (R 2 ), and refined index of agreement (dr) for both meshes.
Comparing Table 3 with Table 4, there is no evidence of scale effects; the model performed with similar accuracy for both physical model and prototype dimensions. Only for Reynolds shear stress, a slightly better agreement exists at prototype scale, for the other parameters the agreement is similar. Additionally, a good agreement between the experimental and numerical results was found for the turbulence parameters, although not as good as that of the velocity components ( Figure 6). The largest differences occurred for the Reynolds shear stress. Nevertheless, for the correlation with fish behaviour in other experimental studies, these differences are considered acceptable. Regarding the transversal and vertical Reynolds normal stress components, a better agreement between Mesh B model results and the experimental data was observed (Figure 6d-f).

Scale Effects
To analyse the scale effects, the prototype scale numerical model results were analysed and compared to the experimental data, scaled according to the Froude similarity.
The flow depths and discharges differences between the numerical model and the experiments slightly decreased for the prototype scale with average values of 3.4 and 4.8% for mesh A prototype and 2.3 and 4.1% for mesh B prototype , for flow depths and discharges, respectively. Table 4 shows the mean absolute difference (MAD), the coefficient of determination (R 2 ), and refined index of agreement (d r ) for both meshes. Comparing Table 3 with Table 4, there is no evidence of scale effects; the model performed with similar accuracy for both physical model and prototype dimensions. Only for Reynolds shear stress, a slightly better agreement exists at prototype scale, for the other parameters the agreement is similar.

Discussion
Most of the previous numerical modelling studies on fishways relied only upon qualitative comparison of model predictions with experimental or field measurements (e.g., [25][26][27][28][29][30]) and only a few assessed turbulence parameters. As stated by Fuentes-Pérez et al. [35], finding numerical validation data in fishway studies is rather difficult, especially for turbulence metrics. However, Lane and Richards [70] stated that a qualitative assessment of the model predictions is not enough, as velocity measurements represent a sample of a much richer flow field. Tables 3 and 4 presented a quantitative statistical comparison between the results obtained with the numerical model and the experimental data obtained with the ADV and the PIV in a laboratory physical model.
The values obtained for MAD, R 2 , and d r showed an overall good agreement between the numerical results and the experimental data. The best agreement was achieved for the longitudinal and transversal velocity components (u and v) and velocity magnitude (U), with similar or better agreements relatively to previous studies on open channel flows (e.g., [71][72][73][74][75][76]). Additionally, good agreements were found in the vertical velocity component (w), turbulent kinetic energy (κ), and longitudinal Reynolds normal stress u 2 . A lower correlation was reported by several authors for w (e.g., [71][72][73][74][75][76]), since vertical velocities are typically very small. In addition, weaker agreements were found for the turbulent kinetic energy, since it involves second-order moment turbulence statistics [76].
For the transversal and vertical Reynolds normal stresses v 2 and w 2 , the agreement obtained for meshes A model and A prototype was relatively weak. For meshes B model and B prototype , the agreement was much higher and similar to the ones obtained in the vertical velocity, turbulent kinetic energy, and longitudinal Reynolds normal stress. This highlights the importance of validating the numerical model for all the analysed parameters and not only for the most common ones (e.g., discharge, flow depths, and mean velocity).
The weakest agreement, although still acceptable, was obtained for the Reynolds shear stress (τ uv ). This is not surprising since this is a third-order statistical moment, where measurements' uncertainties become more significant. At prototype scale, the agreement slightly improved. Reynolds normal stress and Reynolds shear stress comparison between experimental data and numerical results were seldom performed. It is worth mentioning that the agreement and differences observed for meshes A prototype and B prototype are in accordance with the values observed by Fuentes-Pérez et al. [35] for a vertical slot fishway flow.
The accuracy shown by the numerical model was similar for the prototype and for the physical model, showing no evidence of scale effects. Dargahi [45] used FLOW 3D ® to simulate the flow through a bottom outlet with a moving gate in two different scales, the prototype scale and a 1:22 physical model scale, and it also observed that the model performed with similar accuracy for both scales.
The results attest the ability of the numerical model to adequately characterize the complex flow field of a pool-type fishway with bottom orifices. Although no perfect match between the numerical model predictions and the experimental velocities and turbulence data was found, the agreement was generally better or comparable with what has been regarded as acceptable in previous studies on open channel flows (e.g., [71][72][73][74][75][76]) and fishways [35]. As the experimental data used in the present study was much more extensive (with two experimental measurement techniques used and a much higher number of points compared than in previous studies), with the flow field being more complex, the confidence to use such numerical models for modelling fish passes hydrodynamics is strengthened.
It is also important to mention that the numerical meshes used have a good compromise between accuracy and computational cost. The simulations were performed in a quite accessible computer (Table 1), and the run times for meshes A and B needed to obtain converged time-averaged values is of the order of some days, which is compatible with practical applications.
It is important to emphasize that to accurately reproduce this complex flow field, especially the turbulence parameters, LES model should be used. The lower performance of RNG, k-ε, and k-ω turbulence models in adequately reproducing this flow field, namely, the two last ones, to predict the sizes and shapes of the recirculation zones, may result from the isotropic eddy viscosity assumption and from improperly accounting for the interaction between swirl and turbulence. In this flow field the energy dissipation rate is higher adjacent to the jet. Energy is rapidly dissipated as the jet progresses along the pool. This rapid decay is caused by the entrainment of the recirculating flow on both sides of the jet.
It should also be mentioned that, although in this flow field wall friction does not play an important role, since turbulence production due to wall friction is low when compared to horizontal shear production, as found in other studies [26,27,29], the roughness effect of the walls and bed may be changed by considering different roughness heights. If substrate or bottom macro-roughness's are added, these can also be represented by using a smaller sized cell in those locations.
In designing an efficient pool-type fishway, one should consider several relevant hydrodynamic variables, such as velocity, flow depth, turbulence, and flow patterns in the pool [2,[4][5][6][7][8][9][10][11][12]. Results show that the numerical model accurately predicted these hydrodynamic variables. The flow field in this study mimics real conditions: the prototype numerical model is a real-scaled model based on an indoor full-scale prototype fishway configuration in which tests with fishes (the Iberian barbel) were performed [7,38]. In this study, we chose to analyse one of the configurations used in Silva et al. [38] that proved to be effective for the Iberian barbel, a large-bodied potamodromous benthic cyprinid. Pool-type fishway flows are subcritical flows controlled by the downstream water depth; thus, for a given configuration and slope, when fixing the pool mean water depth only one flow discharge assures uniform flow in the fishway guaranteeing the same hydraulic conditions in every pool, which are the ideal conditions [35] to maximize fishway efficiency. Since Silva et al. [7,39] found that the Reynolds shear stress was the turbulent parameter that most strongly influences the movements of Iberian barbel, the ability of the numerical model to accurately predict this parameter was also assessed in this study, which was not considered in previous fishway numerical modelling studies.
Considering that the flow parameters will be ultimately used for a correlation with fish behaviour, the differences are considered quite acceptable since fish behaviour uncertainties are normally much larger. Thus, it may be concluded that CFD models, namely FLOW 3D ® , can be an adequate and efficient tool to investigate and design pool-type fishways, even with complex flow patterns, such as the ones found in pool-type fishways with bottom orifices. Further research may include the 3D numerical modelling of different basin dimensions and cross-walls configurations.

Conclusions
In this study, two distinct measurement techniques, ADV and PIV, were used to measure instantaneous velocities in a laboratory pool-type fishway facility, which is a 1:2.5 physical model of an indoor prototype size pool-type fishway. The experimental data was compared with a 3D numerical model.
Several hydrodynamic variables relevant for fishway analyses and design, such as flow depths, time-averaged velocities, and turbulence parameters (namely, turbulent kinetic energy, Reynolds normal stresses, and Reynolds shear stress component parallel to the bottom), were determined and compared. The results highlighted the importance of validating the numerical model for all the parameters under analyses.
A good agreement between the numerical model and the experimental data was achieved, showing the ability of the numerical model to adequately reproduce the complex flow field of a pool-type fishway with bottom orifices. In addition, no scale effects were observed in the numerical modelling, which presented similar accuracy for near-prototype and physical model scales. Therefore, since the used numerical meshes have a good compromise between accuracy and computational cost, we consider that 3D CFD models, namely FLOW 3D ® , may be successfully used to simulate the pool-type fishway hydrodynamics, constituting a valuable tool to analyse and design fishways, allowing to eliminate and/or complement physical model testing. This could greatly improve fishways design process, as new configurations or different solutions for retrofitting existing nonoperational fishways could be tested and optimized in shorter periods, and in a less expensive way prior to their construction.

Data Availability Statement:
The data presented in this study is contained and it is available within the article and upon request to the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following symbols are used in this paper:

Appendix A
Given the significant improvement in the agreement achieved when doubling the number of cells in each direction, for the transversal and vertical Reynolds normal stresses, the mesh was refined using the grid overlay boundary condition and halving the grid spacing in each direction in pool 3. These mesh (C model ) characteristics were presented in Table 1. In addition, given the significant number of points that showed oscillatory convergence behaviour, an additional grid refinement (mesh D model ) was performed as recommended by Celik et al. [59]. It should be mentioned that the computation cost needed to obtain these meshes time-averaged results is much higher and not compatible with practical engineering applications. These results (C model and D model ) were compared with the ADV and PIV measurements (Table A1). Comparing Table 3 to Table A1, some reduction of the agreement with the experimental data can be observed for mesh C model , for all the analysed parameters, when compared to coarser meshes (A model and B model ). The best agreement with the measurements is found with the coarser meshes (A model and B model ), which is rather counterintuitive since the smaller the grid size is, the larger the range of eddy scales which are resolved by the LES model will be. Similar results can also be found in [61,64,77]. Meyers et al. [61] showed that counterintuitively to what is thought, an "improved simulation" at higher spatial resolution can yield results with a larger total error since the numerical and discretization errors in a coarse grid can have opposite signs counter-acting and partially cancelling each other. So, although the two error components (numerical and discretization errors) are higher in magnitude on a coarse grid, the sum (total error) can be higher on a finer grid. A discussion on this issue can also be found in Celik et al. [62].
Comparing mesh C model to mesh D model , an oscillatory behaviour is observed with mesh D model presenting a better agreement with the experimental data, although not as good as the coarser meshes agreement. Celik at al. [65] showed that this behaviour may be caused by the oscillatory velocity that occurs in recirculation regions.