Prioritizing Design Parameters for Stepped Chutes and Shear Stress Distribution

Stepped chutes offer high efficiency in decreasing flow velocity due to roughness; however, negative impacts may still be experienced by the receiving water body into downstream. These effects might be mitigated if geometric and hydraulic parameters governing the structure are well addressed. Herein, five influential parameters were developed, i.e., longitudinal slope S (S = tan θ), discharge (Q), pool height above steps (hp), chute width (W), and chute height (H), employing a three dimensional (3D) numerical model. Through 600 simulations, two regression models were developed for predicting depth-averaged velocity at the last step Vd (m/s) and critical length Lc (cm) at the downstream where the maximum velocity occurs, using response surface methodology (RSM) based on the mixed-level full factorial design. The prediction data obtained by developed regression models were agreeable with actual data with coefficient determination (R2) of about 0.95, highlighting the accuracy and ability of the models for the prediction of Vd and Lc. Additionally, the analysis of variance (ANOVA) was used to prioritize the impact of the studied parameters on Vd and Lc. Results highlighted that among geometric parameters, W and S had a significant influence on Vd and Lc; however, the impact of W was more pronounced. Using a regression model for Lc, a cross section was obtained, and the shear stress distribution of the downstream was compared with that of the last step and sidewalls. The shear stress patterns showed that the maximum value shifted from the side walls to the downstream between the lower and higher slopes. Further, the longitudinal distribution of shear stress at the downstream revealed that geometric and hydraulic characteristics played a negligible role in the changing pattern of shear stress. The results of this study reveal the dynamic behavior of the given structure where different geometric options are available for structure design.


Introduction
In the implementation of stepped chutes in natural waterways, there are usually different design options. Understanding the behavior of the design parameters and their influence on the structure performance can mitigate negative impacts such as erosion and scouring through bringing the flow condition close to the natural flow of the receiving water body [1]. Even with the efficiency of such structures in reducing flow velocity, the discharge into the downstream must be considered carefully owing to high uncertainty [2] in predicting the flow characteristics compared with smooth structures. This structure is widely used in water transmission systems ( Figure 1a) and in rivers for the purposes of fishways (Figure 1b). For these cases, contrary to the stepped chutes when they are used as spillways in embankment dams, the downstream is not equipped with the stilling basins to mitigate the vulnerability caused by the receiving water body. Studies have shown that geometric properties significantly influence chute p mance [3], and this has led to numerous laboratory works to better understand the e of gradual changes to geometric components. For example, in slope 26.6°, the st chute with flat steps, i.e., hp = 0, improved the performance of the chute in terms creasing flow velocity compared to the chute with hp = 3 cm [4,5]. For slope 8.9°, op results were reported by Felder et al. [6]. A comparison of the results obtained for 8.9° with data from Thorwarth [7] for slope 14.6° indicated that there was similar b ior in terms of how pools influence velocity. Kökpinar [8] and Andre [9] research stepped chutes with slopes of 30° and 18.6° on the basis of two-phase flows. T analyzed data of stepped chutes with slopes of 30° and 26.6° showed similar beh of the pools on the chute performance. Other publications that have presented de hydraulic analysis on stepped chutes include Gonzalez [10], Meireles and Mato Bung [12], and Zhang and Chanson [13,14]. The focus of these studies varies in ter the effects produced from changes to some geometric parameters or to changes in A comprehensive analysis of the effects produced from changes in multiple geom parameters or the influence that parameters have on each other as well as overall p mance has not been well addressed. Additionally, these studies did not focus on stress patterns in both steps and downstream; instead, some addressed this in the fo flow resistance [15,16]. In open channels such as stepped chutes in natural rivers, stress acts upon the sidewalls and bed surface, affecting the structural stability of the nel section [17].
To eliminate such limitations, a number of experimental studies have develope pirical formulas to estimate the effects of different characteristics and parameters i structure design [18]. The relationships and formulas determined in such works ar ited to the range or data set of a given study. Often, when these results are consi with data outside of the original study, the original relationship either no longer a or degrades, resulting in errors and potential problems in the design. Advances i chine learning and numerical modeling provide a base to these limitations of experim Studies have shown that geometric properties significantly influence chute performance [3], and this has led to numerous laboratory works to better understand the effects of gradual changes to geometric components. For example, in slope 26.6 • , the stepped chute with flat steps, i.e., h p = 0, improved the performance of the chute in terms of decreasing flow velocity compared to the chute with h p = 3 cm [4,5]. For slope 8.9 • , opposite results were reported by Felder et al. [6]. A comparison of the results obtained for slope 8.9 • with data from Thorwarth [7] for slope 14.6 • indicated that there was similar behavior in terms of how pools influence velocity. Kökpinar [8] and Andre [9] researched the stepped chutes with slopes of 30 • and 18.6 • on the basis of two-phase flows. The re-analyzed data of stepped chutes with slopes of 30 • and 26.6 • showed similar behavior of the pools on the chute performance. Other publications that have presented detailed hydraulic analysis on stepped chutes include Gonzalez [10], Meireles and Matos [11], Bung [12], and Zhang and Chanson [13,14]. The focus of these studies varies in terms of the effects produced from changes to some geometric parameters or to changes in slope. A comprehensive analysis of the effects produced from changes in multiple geometric parameters or the influence that parameters have on each other as well as overall performance has not been well addressed. Additionally, these studies did not focus on shear stress patterns in both steps and downstream; instead, some addressed this in the form of flow resistance [15,16]. In open channels such as stepped chutes in natural rivers, shear stress acts upon the sidewalls and bed surface, affecting the structural stability of the channel section [17].
To eliminate such limitations, a number of experimental studies have developed empirical formulas to estimate the effects of different characteristics and parameters in step structure design [18]. The relationships and formulas determined in such works are limited to the range or data set of a given study. Often, when these results are considered with data outside of the original study, the original relationship either no longer applies or degrades, resulting in errors and potential problems in the design. Advances in machine learning and numerical modeling provide a base to these limitations of experimental works. Jiang et al. [19], for example, were able to predict the energy dissipation of a stepped chute using support vector machine regression (SVR) with data available from different literature results. Khatibi et al. [20] employed artificial neural networks (ANNs) and genetic programming (GP) to predict energy dissipation over gabion stepped chutes. In addition, there are other methods that have been used in the field of hydraulic engineering to estimate discharge coefficient, scour, aeration efficiency, and energy dissipation (Baylar et al. [21]; Laucelli and Giustolisi [22]; Emiroglu et al. [23]; Azamathulla [24]; Dursun et al. [25]; Bagatur and Onen [26]; Roushangar et al. [27]). If properly applied, numerical models can also be employed to extend the results of different geometric changes to a given structure and their corresponding effects [28]. For example, Morovati et al. [28] and Li and Zhang [29] investigated further configurations for pools in addition to the configurations proposed by Felder et al. [15]. The same approach was adopted by Attarian et al. [30]. A comparison of the results reported by Morovati and Eghbalzadeh [31] with Li and Yang [32] indicated that pressure distribution on flat and inclined steps followed a virtually S-shaped distribution. Inception point location was also researched by Munta et al. [33], Bombardelli et al. [34], and Morovati and Eghbalzadeh [35]. Morovati et al. [3] investigated the effects of various configurations of openings on pooled stepped chutes, and they concluded that a staggered configuration of openings improved the chute's efficiency concerning the residual head. The proposed models rely on the laboratory models, which have limited conditions in terms of geometric changes investigated on the model. Hence, it seems that more research needs to be conducted to investigate the effects of gradual changes in geometric properties. Moreover, to the best of the authors' knowledge, no other research has addressed the prioritization of the influence of different parameters thus far.
The objective of this work is, therefore, to develop the geometric and hydraulic parameters to prioritize and predict their effects on chute performance. Herein, a 3D numerical model developed by Morovati et al. [3] is employed. The results of the 3D simulations were then used to develop two regression models for predicting depth-averaged velocity (V d ) at the last step and critical length (L c ) at the downstream where maximum velocity occurs (see Figure 2). Additionally, the influence of each parameter on the chute performance is prioritized based on the data yielded by the developed numerical and regression models. Then, by utilizing these results, shear stress pattern in the last step, sidewalls, and downstream are produced and analyzed. works. Jiang et al. [19], for example, were able to predict the energy dissipation of a stepped chute using support vector machine regression (SVR) with data available from different literature results. Khatibi et al. [20] employed artificial neural networks (ANNs) and genetic programming (GP) to predict energy dissipation over gabion stepped chutes.
In addition, there are other methods that have been used in the field of hydraulic engineering to estimate discharge coefficient, scour, aeration efficiency, and energy dissipation (Baylar et al. [21]; Laucelli and Giustolisi [22]; Emiroglu et al. [23]; Azamathulla [24]; Dursun et al. [25]; Bagatur and Onen [26]; Roushangar et al. [27]). If properly applied, numerical models can also be employed to extend the results of different geometric changes to a given structure and their corresponding effects [28]. For example, Morovati et al. [28] and Li and Zhang [29] investigated further configurations for pools in addition to the configurations proposed by Felder et al. [15]. The same approach was adopted by Attarian et al. [30]. A comparison of the results reported by Morovati and Eghbalzadeh [31] with Li and Yang [32] indicated that pressure distribution on flat and inclined steps followed a virtually S-shaped distribution. Inception point location was also researched by Munta et al. [33], Bombardelli et al. [34], and Morovati and Eghbalzadeh [35]. Morovati et al. [3] investigated the effects of various configurations of openings on pooled stepped chutes, and they concluded that a staggered configuration of openings improved the chute's efficiency concerning the residual head. The proposed models rely on the laboratory models, which have limited conditions in terms of geometric changes investigated on the model. Hence, it seems that more research needs to be conducted to investigate the effects of gradual changes in geometric properties. Moreover, to the best of the authors' knowledge, no other research has addressed the prioritization of the influence of different parameters thus far. The objective of this work is, therefore, to develop the geometric and hydraulic parameters to prioritize and predict their effects on chute performance. Herein, a 3D numerical model developed by Morovati et al. [3] is employed. The results of the 3D simulations were then used to develop two regression models for predicting depth-averaged velocity (Vd) at the last step and critical length (Lc) at the downstream where maximum velocity occurs (see Figure 2). Additionally, the influence of each parameter on the chute performance is prioritized based on the data yielded by the developed numerical and regression models. Then, by utilizing these results, shear stress pattern in the last step, sidewalls, and downstream are produced and analyzed.

Numerical Model
A 3D numerical model developed by Morovati et al. [3] was adopted for the present work purposes. Detailed information on the employed numerical model in terms of meshing, mesh convergence analysis, governing equations, and boundary conditions considered for turbulence model and walls is found in Morovati et al. [3], and thus only a summary is provided here. The conducted simulations enjoy RANS equations, which were solved using the finite-volume method, and all equations were formulated with area and volume porosity functions using the fractional area/volume obstacle representation (FAVOR) approach to model complex geometric regions [36]. The renormalization-group (RNG) method that was employed in this study has broad applicability and is known to describe more accurately the flows characterized with areas of strong shear forces [36]. Advection in the conducted simulations was completed using the second-order upwind method for the discretization scheme. The Courant-Friedrichs-Lewy (CFL) value and the Krylov subspace dimension were set at <0.6 and 15, respectively, and the CFL < 0.6 was maintained by automatically adjusting the lengths of the time steps [3]. Results are presented at the condition concerning the final steady-state throughout the simulations implemented in this paper. The simulation time for the implemented runs employing a workstation with 24 cores and 32G of RAM was one to three days.
The free surface boundary conditions must be clearly defined for an accurate assessment of the free-surface dynamics. Subsequently, considering the free surface or sharp interface option, the volume of fluid (VOF) method was used in this study to evaluate these conditions with a single-fluid approach [36]. In addition, the actual physical situation should be modeled as closely to the real-world conditions as possible. This is important to make sure the results predicted are accurate to what is observed. Subsequently, specific discharge with the upstream total head for inlet, the outflow boundary condition for outlet, specified pressure for top, and wall boundary for the bottom and side walls were applied.
The computational domain was defined using hexahedral cell and multi-block technique. Six blocks were defined in which a uniform distribution of cells was considered in the transverse direction, i.e., Z direction. For the chute blocks, a uniform distribution of cells was considered in the X and Y direction, while for the upstream and downstream blocks, variable-sized cells were applied in X and Y directions so that, at the shared planes with the chute blocks, the cell sizes remained constant and as the distance from the last step and chute crest increased, the cell size in the X direction increased. For more information, refer to Morovati et al. [3].
In employing k − ε(RNG) turbulence closure, the values for shear stress, turbulence kinetic energy k, and dissipation rate ε, were calculated using wall functions. Therefore, y + values remained at the logarithmic region for the cells at the boundary layer, i.e., y + ∈ [35,300]. Table 1 provides information about the adopted stepped chutes in the present work. For each step chute, the uniform discharges were considered at the beginning of the computational domain and reached the chutes through a broad-crested weir followed by identical steps. Table 1. Geometry properties of the stepped chutes adopted in the present study. NS is the number of steps.

S l (m) h (m) H (m) W (m) NS
Felder et al. [

Design of Experiment
Design of experiment (DOE) is an effective method to understand a process as opposed to observing a process. Studying the influence of input parameter(s) on the response(s) is highly essential in a process, especially when different options are available for design [37][38][39][40]. Hence, DOE was applied by selecting the response surface methodology (RSM) based on the mixed-level full factorial design. This efficient method of DOE can also be used for modeling and optimizing processes [40,41]. Therefore, some important geometric and hydraulic parameters dominating over stepped chutes were developed in relation to experimental models employing high accurate 3D numerical model [3]; then, results were used as input parameters for the DOE model ( Table 2). The input parameters are chute slope S, discharge (Q), pool height (h p ), chute height (H), and chute width (W). The responses are depth-averaged velocity V d (m/s) at the last step and critical length L c (cm). Table 2 lists these input parameters along with their respective levels. The full factorial design based on the user-defined option resulted in 600 CFD runs 5 2 × 4 1 × 3 1 × 2 1 . Table 2. The input parameters used for DOE together with their levels.

Input Parameters Level 1 Level 2 Level 3 Level 4 Level 5
Chute slope S (S = tanθ) The levels given in Table 2 were selected based on real-world constructed examples. The range considered for H, Q, S, and h p is representative of some stepped chutes used as waterways in rivers, drop structures, and transmission channels (see Guenther et al. [42]). For a chute width W that is greater than the range considered, there are practical examples that exist; however, such cases are more applicable for embankment dams and spillways, which are not the focus of the present study. Additionally, larger pool heights, i.e., h p > 0.1 m, can cause problems associated with sediment accumulation and stability in the long term. Regarding discharge Q, since the flow regime passing on stepped chutes may change over time, the discharge ranges were chosen in a way that three types of flow regimes governing the chute, including nappe, transition, and skimming, are formed. It should be noted that the developed cases based on chute height H were performed in such a way that only the step height h was changed, i.e., l remained unchanged. In addition, in all cases, the number of steps in the 3D numerical model implemented based on the levels presented in Table 2 was similar to the experimental ones given in Table 1. For example, for H = 1.5 in the stepped chute with a slope of 26.6 • , the step height was changed from 0.1 m to 0.15 m.
This design allows researchers to analyze responses at whole combinations of the input parameters levels. Here, all significant parameters and their interactions have a 95% confidence interval level or a probability value less than 0.01 (p-value < 0.05) towards the responses. The analysis of variance (ANOVA) was also used to investigate the statistical significance of input parameters on the responses. Design Expert 10.0.4.0 and Minitab ® 18.1 were employed to analyze the results.

Model Development
Two-factor interaction (2FI) models suggested by the software were developed to predict V d and L c as responses. All main parameters and their interactions with the probability value less than 0.05 (p-value < 0.05) were considered as significant parameters. According to ANOVA results presented in Table 3 for the V d response, all main parameters are highly significant. The ANOVA results in Table 4 for L c show that all main parameters are significant towards the response except for the H parameter with a p-value of 0.0992.
Since the interaction of this main parameter is significant, it is not possible to remove the H parameter from the initial model [43]. As can be seen in Table 4, the interaction of S-H has a p-value of less than 0.0001; then, the H parameter was kept in the final model for the prediction of L c . For both models, insignificant interactions were removed from the initial models. Therefore, the final improved 2FI models were obtained according to the following equations: The Equations (1) and (2) developed are based on real-world factors and conditions. they can be used to predict responses for given levels of each factor but cannot determine the relative significance of each factor. To find the relative importance of the design factors, all analysis and model fittings need to be performed using coded design variables [43]. In the coded design variable, i.e., −1 ≤ x i ≤ +1, the magnitudes of the model coefficients are directly comparable, that is, they are all dimensionless, and they measure the effects of changes in a design factor over a one-unit interval. In this study, each independent parameter was coded as x i according to Equation (3) [40]: where X 0 is the value of X i (selected parameters) at the center point, and ∆X denotes the step change. In addition, depth-averaged velocity V d and critical length L c as response data were normalized by dividing their actual value into a mean value of 2.3 (m/s) and 58 (cm), respectively. The final developed models obtained in terms of coded factors and normalized response are as follows. It is noteworthy that all parameters are unitless.
The coded equation helps to detect the relative impact of the factors by comparing the absolute value of factor coefficients [37]. In this regard, the descending order from highest to lowest significant factors on the response is as follows: Prioritizing the effects of parameters is of great importance in recognizing the structure's behavior and the design aspects. Tables 3 and 4 give information regarding the impact of each parameter employing ANOVA results. Regarding V d response, the highest important parameter is Q with an F-value of 4416.80, followed by W, S, H, and h p , as presented in Table 3. In other words, the Q and h p parameters have the highest and lowest impacts on the V d response, respectively. This order also remained constant for L c response except for H and h p , in which H has the lowest impact. For both responses, W as a geometric parameter plays a significant role in V d and L c . For both responses, H and h p parameters have an insignificant impact compared to the other three parameters; however, this effect is more pronounced for the L c response. According to ANOVA results presented in Tables 3 and 4, the model's F-values of 1157.31 and 1012.23 (p-value < 0.0001) for V d (m/s) and L c (cm), respectively, show that the final 2FI models are highly significant and fairly appropriate for the present results.

Statistical Validation of Regression Models
After identifying the significance of parameters and developing the regression models, it is important to validate the models and evaluate their accuracy. Hence, four main approaches were selected: (1) plot of the predicted data versus simulated data, (2) normal probability plot, (3) plot of residuals vs. run order, and (4) plot of residuals versus predicted data. Figure 3a,b illustrates the predicted data by 2FI models against the simulated data together with the absolute relative error (ARE) and root mean square error (RMSE) to verify the ability of the models for the prediction of the response. The formulas for the statistical measures, coefficient of determination (R 2 ), as well as adjusted and predicted R 2 are listed in Table 5.  As an additional step for model validation, it is necessary to compare the probability to residuals and residuals to the run order. The normal probability plots are used to identify whether the data follows a normal distribution. As shown in Figure 4, the points fall    Figure 3a,b shows the close proximity of the model prediction with the simulated data demonstrating the validity of 2FI models. From Figure 3 and Tables 3 and 4, the R 2 was found to be almost 0.95, which is close to 1, indicating that the models can explain 95% of the variability in response, and the predicted data from the 2FI models agree well with the simulated data. Sometimes a model includes several independent parameters and the interaction terms; in such cases, in order to measure the quality of the model, it is fair to consider the adjusted R 2 [43]. In addition, the difference between the predicted R 2 and the adjusted R 2 values should be less than 0.2. As seen in Tables 3 and 4, both V d and L c results show reasonable agreement between the predicted R 2 and the adjusted R 2 with a difference of less than 0.001.
As an additional step for model validation, it is necessary to compare the probability to residuals and residuals to the run order. The normal probability plots are used to identify whether the data follows a normal distribution. As shown in Figure 4, the points fall along the reference line, and there is no kurtosis or significant variation in the sample distributions, indicating the normal distribution of data.
The plot of externally studentized residuals is an effective method to detect the outlier data since externally studentized residuals are more sensitive to detect outliers as compared to internally studentized residuals. According to Figure 5a, the discrepant data and potential outliers are not detected in the actual data for V d . Figure 5b also displays the externally studentized residuals versus predicted values for L c . The plot ideally exhibits a constant variation from the left (low level of response) to the right (highest predicted level), and no outlier data are detected. Similar results were obtained for the externally studentized residuals versus predicted values for V d , and externally studentized residuals versus run for L c . along the reference line, and there is no kurtosis or significant variation in the sample distributions, indicating the normal distribution of data. The plot of externally studentized residuals is an effective method to detect the outlier data since externally studentized residuals are more sensitive to detect outliers as compared to internally studentized residuals. According to Figure 5a, the discrepant data and potential outliers are not detected in the actual data for Vd. Figure 5b also displays the externally studentized residuals versus predicted values for Lc. The plot ideally exhibits a constant variation from the left (low level of response) to the right (highest predicted level), and no outlier data are detected. Similar results were obtained for the externally studentized residuals versus predicted values for Vd, and externally studentized residuals versus run for Lc. igure 5. The plots of (a) residuals vs. ran for depth-averaged velocity Vd (m/s), and (b) residuals vs. predicted data for ritical length Lc (cm).
In order to evaluate the performance of the model further, the data available from the literature that have different geometry and hydraulic characteristics in relation to what were considered as input parameters (see Table 2) were compared with the data predicted by the 2FI model. Figure 6 shows the scatter plot of the results obtained from the 2FI model  The plot of externally studentized residuals is an effective method to detect the outlier data since externally studentized residuals are more sensitive to detect outliers as compared to internally studentized residuals. According to Figure 5a, the discrepant data and potential outliers are not detected in the actual data for Vd. Figure 5b also displays the externally studentized residuals versus predicted values for Lc. The plot ideally exhibits a constant variation from the left (low level of response) to the right (highest predicted level), and no outlier data are detected. Similar results were obtained for the externally studentized residuals versus predicted values for Vd, and externally studentized residuals versus run for Lc. In order to evaluate the performance of the model further, the data available from the literature that have different geometry and hydraulic characteristics in relation to what were considered as input parameters (see Table 2) were compared with the data predicted by the 2FI model. Figure 6 shows the scatter plot of the results obtained from the 2FI model In order to evaluate the performance of the model further, the data available from the literature that have different geometry and hydraulic characteristics in relation to what were considered as input parameters (see Table 2) were compared with the data predicted by the 2FI model. Figure 6 shows the scatter plot of the results obtained from the 2FI model and the actual data reported by Felder et al. [6,15]. The values of the statistical measures presented in Figure 6 indicate that the given 2FI model enjoys a high capability of predicting V d , even with the geometry and hydraulic changes; however, these changes should be in the range presented in Table 2. and the actual data reported by Felder et al. [6,15]. The values of the statistical measures presented in Figure 6 indicate that the given 2FI model enjoys a high capability of predicting Vd, even with the geometry and hydraulic changes; however, these changes should be in the range presented in Table 2.

Main Impact of Factors and Their Interaction Effects
In Tables 3 and 4, the parameters' impacts were prioritized, although the effects of the parameters on the responses are difficult to understand. Therefore, the 3D surface plots are provided to better illustrate the parameter's impact on the Vd and Lc responses. Here, these plots are presented through data obtained by the 2FI models with two parameters (the three other parameters were fixed at their specific values, i.e., = 20°, Q = 0.1342 m 3 /s, hp = 0.05 m, H = 1.5 m, and W = 0.75 m). Figure 7a,b demonstrates the effects of Q, , and H parameters on the Vd. In both figures, the simultaneous increase of the parameters increases the responses, with the difference that the slope of is higher than H, resulting in a more significant increase. There is a different behavior pattern for the impact of and W on the Vd response; for example, for the minimum value of Vd, W should increase while the needs to be decreased (Figure 7c). Such a similar pattern also exists in the case where Lc forms near the last step; however, this seems not to be a cautious approach since high velocity exiting the area after the last step can cause scouring and breach formation, undermining the stability of the structure. Moreover, a comparison of W and parameters in Figure 7f indicates that W plays a more significant role in controlling the location of maximum velocity at the downstream. Therefore, according to the results presented in Tables 3 and 4 and Figure 7, W can be regarded as the most important geometric parameter influencing Vd as well as Lc.

Main Impact of Factors and Their Interaction Effects
In Tables 3 and 4, the parameters' impacts were prioritized, although the effects of the parameters on the responses are difficult to understand. Therefore, the 3D surface plots are provided to better illustrate the parameter's impact on the V d and L c responses. Here, these plots are presented through data obtained by the 2FI models with two parameters (the three other parameters were fixed at their specific values, i.e., S = 20 • , Q = 0.1342m 3 /s, h p = 0.05 m, H = 1.5 m, and W = 0.75 m). Figure 7a,b demonstrates the effects of Q, S, and H parameters on the V d . In both figures, the simultaneous increase of the parameters increases the responses, with the difference that the slope of S is higher than H, resulting in a more significant increase. There is a different behavior pattern for the impact of S and W on the V d response; for example, for the minimum value of V d , W should increase while the S needs to be decreased (Figure 7c). Such a similar pattern also exists in the case where L c forms near the last step; however, this seems not to be a cautious approach since high velocity exiting the area after the last step can cause scouring and breach formation, undermining the stability of the structure. Moreover, a comparison of W and S parameters in Figure 7f indicates that W plays a more significant role in controlling the location of maximum velocity at the downstream. Therefore, according to the results presented in Tables 3 and 4 and Figure 7, W can be regarded as the most important geometric parameter influencing V d as well as L c .
When there is more than one variable in a process, the interaction between the factors may occur. There is interaction when one variable's impact on the response relies on the level of another variable. The DOE method's strength is the management of interactions amongst variables [37].
As shown in Equations (1) and (2) and as evidenced by Tables 3 and 4, there are interactions between the factors. The significant interaction terms are S × Q, h p , Q × h p , and Q × H in Equation (1), while in Equation (2), Q, S × h p , S × H, S × W, Q × h p , and Q × W are significant interaction terms. The interactions between the factors can be established using an interaction graph to display any two-factor interactions. These figures illustrate the response as a function of one variable with the two different levels of the other variable while the remaining independent variables are kept constant. In such plots, the nonparallel lines demonstrate that an interaction exists between the two factors (the details on the interaction evaluation are provided in Mathews [37]). Two of such interaction terms are highlighted and discussed in the subsequent paragraphs. When there is more than one variable in a process, the interaction between the factors may occur. There is interaction when one variable's impact on the response relies on the level of another variable. The DOE method's strength is the management of interactions amongst variables [37].
As shown in Equations (1) and (2) and as evidenced by Tables 3 and 4, there are interactions between the factors. The significant interaction terms are × , ℎ , × ℎ , and × in Equation (1), while in Equation (2), , × ℎ , × , × , × ℎ , and × are significant interaction terms. The interactions between the factors can be established using an interaction graph to display any two-factor interactions. These figures illustrate the response as a function of one variable with the two different levels of the other variable while the remaining independent variables are kept constant. In such plots, the nonparallel lines demonstrate that an interaction exists between the two factors (the details on the interaction evaluation are provided in Mathews [37]). Two of such interaction terms are highlighted and discussed in the subsequent paragraphs.
The interaction between the discharge (Q) and slope ( ) and its impact on the depthaveraged velocity is shown in Figure 8. Depending on the discharge value Q, the stepped chutes experience fully developed or partially developed hydraulic jumps in lower discharges, coherent stream in higher discharges, and cavities above the steps and the corner created between the last step's height and downstream bottom [44,45]. These features play a significant role in passing flow in terms of decreasing flow velocity and can be changed if the geometry and hydraulic characteristics are altered. Figure 8a shows that Vd varies as the Q changes on each slope, with more pronounced variations attributed to slope 8.9°. On slope 8.9° and for the discharge range given in Table 1, i.e., 0.049 < Q (m 3 /s) < 0.2195, the chute experiences nappe, transition, and skimming flow regimes while for stepped slope, i.e., 30°, a skimming flow regime forms (see Chanson [44]). This implies that the more significant variation of Vd is attributed to the changes between various flow regimes, i.e., nappe to skimming, than the changes in Q within the limit of a flow regime, i.e., skimming. The interaction between the chute width (W) and discharge (Q) and its impact on critical length is illustrated in Figure 8b. As seen, changing the chute width W as a factor changing the flow regime causes significant differences for Lc, especially in greater discharges. According to Equation (6), as the chute width W increases at a constant discharge Q, the flow depth passing over the chute decreases. As a result, cavities and hydraulic The interaction between the discharge (Q) and slope (S) and its impact on the depthaveraged velocity is shown in Figure 8. Depending on the discharge value Q, the stepped chutes experience fully developed or partially developed hydraulic jumps in lower discharges, coherent stream in higher discharges, and cavities above the steps and the corner created between the last step's height and downstream bottom [44,45]. These features play a significant role in passing flow in terms of decreasing flow velocity and can be changed if the geometry and hydraulic characteristics are altered. Figure 8a shows that V d varies as the Q changes on each slope, with more pronounced variations attributed to slope 8.9 • . On slope 8.9 • and for the discharge range given in Table 1, i.e., 0.049 < Q (m 3 /s) < 0.2195, the chute experiences nappe, transition, and skimming flow regimes while for stepped slope, i.e., 30 • , a skimming flow regime forms (see Chanson [44]). This implies that the more significant variation of V d is attributed to the changes between various flow regimes, i.e., nappe to skimming, than the changes in Q within the limit of a flow regime, i.e., skimming. The interaction between the chute width (W) and discharge (Q) and its impact on critical length is illustrated in Figure 8b. As seen, changing the chute width W as a factor changing the flow regime causes significant differences for L c , especially in greater discharges. According to Equation (6), as the chute width W increases at a constant discharge Q, the flow depth passing over the chute decreases. As a result, cavities and hydraulic jumps' properties are likely influenced because their formation directly depends on hydraulic and geometric properties [6,15,44,45].
where H 1 and L creast represent the flow depth above chute crest and crest length, respectively. The results presented in Figure 8b, where the chute width increases from 0.5 to 1 m, indicate that the critical length L c is almost halved for all discharge. Such a decrease appears to be due to (1) the flow depth decreasing as discussed earlier, and thus the critical length occurs at a shorter distance from the last step, and also (2) the decreased flow depth, which likely results in changes in the flow regime as well as the size of the vortices above the steps and at the downstream.
The perturbation plot compares the effects of all parameters at a specific value as well as at the optimum conditions in the design space. Design-Expert sets the reference point default at the middle of the design space (the coded zero level of each parameter). The response is plotted by changing only one parameter over its range while holding all the other parameters constant. The perturbation plots for the V d and L c are displayed in Figure 9a,b, respectively. All parameters experience a linear effect on the V d and L c responses. When the slope of the line is sharper, it means that the parameter has a more significant impact on the response. Both the V d and L c responses are controlled mainly by Q and W parameters, which significantly impact these responses. Two slopes are observed for lines, positive and negative, which means that if the slope of the line is positive, increasing the corresponding parameter to that line would increase responses, but for negative slope, as the parameter increases, the values for responses decreases (for example, see W and Q for L c ). According to the parameters' values shown, both h p and H have an insignificant impact in determining the critical length at the downstream. Additionally, the discharge Q of natural waterways varies constantly as it is affected by many environmental factors. If we assume the geometric parameters are constant, the length of the critical area of the downstream can be calculated based on a given range of discharge values, for example, the minimum and maximum values. jumps' properties are likely influenced because their formation directly depends on hydraulic and geometric properties [6,15,44,45].
where H1 and Lcreast represent the flow depth above chute crest and crest length, respectively. The results presented in Figure 8b, where the chute width increases from 0.5 to 1 m, indicate that the critical length Lc is almost halved for all discharge. Such a decrease appears to be due to (1) the flow depth decreasing as discussed earlier, and thus the critical length occurs at a shorter distance from the last step, and also (2) the decreased flow depth, which likely results in changes in the flow regime as well as the size of the vortices above the steps and at the downstream. The perturbation plot compares the effects of all parameters at a specific value as well as at the optimum conditions in the design space. Design-Expert sets the reference point default at the middle of the design space (the coded zero level of each parameter). The response is plotted by changing only one parameter over its range while holding all the other parameters constant. The perturbation plots for the Vd and Lc are displayed in Figure  9a,b, respectively. All parameters experience a linear effect on the Vd and Lc responses. When the slope of the line is sharper, it means that the parameter has a more significant impact on the response. Both the Vd and Lc responses are controlled mainly by Q and W parameters, which significantly impact these responses. Two slopes are observed for lines, positive and negative, which means that if the slope of the line is positive, increasing the corresponding parameter to that line would increase responses, but for negative slope, as the parameter increases, the values for responses decreases (for example, see W and Q for Lc). According to the parameters' values shown, both hp and H have an insignificant impact in determining the critical length at the downstream. Additionally, the discharge Q of natural waterways varies constantly as it is affected by many environmental factors. If we assume the geometric parameters are constant, the length of the critical area of the downstream can be calculated based on a given range of discharge values, for example, the minimum and maximum values.

Shear Stress
With the help of Equation (2), the critical length at the downstream was measured. The shear stress pattern at the edge of the last step, the sidewalls, and downstream where the maximum velocity occurs (Lc) is displayed in Figure 10 for different slopes and dis-

Shear Stress
With the help of Equation (2), the critical length at the downstream was measured. The shear stress pattern at the edge of the last step, the sidewalls, and downstream where the maximum velocity occurs (L c ) is displayed in Figure 10 for different slopes and discharges. The shear stress pattern is almost symmetrical in all cases, although there are some insignificant differences. For lower slopes, i.e., 8.9 • and 14.6 • , the maximum shear stress occurs on the sidewalls, while in the higher slopes, sidewalls and step surfaces experience approximately identical shear stress. The shear stress pattern on the sidewalls differs between steeper and lower slopes. For lower slopes, the pattern on the sidewalls forms a parabolic distribution with the maximum stress occurring near the middle of the flow depth. For steeper slopes, the pattern on the sidewalls forms an approximately linear distribution with the maximum stress near the base where the sidewalls and step surface intersect. The shear stress distribution at the downstream showed some similarities and small changes in each case. The most significant difference observed was the location of the maximum shear stress value. With higher slopes, the maximum shear stress occurred at the downstream, while for lower slopes, the maximum value was recorded on the sidewalls. Therefore, with increasing slope, the location of the maximum shear stress changes from the sidewall of the last step to the downstream. It should be noted that the results shown in Figure 10 correspond to the flat stepped chutes, i.e., h p = 0. Figure 11 illustrates the shear stress distribution along downstream immediately after the last step at Z/W = 0.5 for different slopes and discharges. The significant variation of shear stress occurs at the area close to the last step, although the values of shear stress in this area are much smaller than the rest of the downstream. After this area, the shear stress sharply increases, followed by an insignificant variation for the rest of the downstream. A comparison of the results highlights the point that the changing pattern of the shear stress seems to be independent of geometric and hydraulic characteristics since almost a similar pattern is seen for all cases; however, these characteristics may impact shear stress values.

Conclusions
This work was conducted to investigate and prioritize the influence of geometric and hydraulic parameters on chute behavior. Five important parameters were selected: slope, discharge, pool height, chute height, and chute width. These parameters were then utilized to develop 600 cases, employing a highly accurate 3D numerical model. Results were used for RSM design, which resulted in two developed regression models for predicting Vd and Lc with coefficient determination (R 2 ) of about 0.95, highlighting the high accuracy of the models. Employing the analysis of variance (ANOVA), prioritization of the influence of the studied parameters on the dynamic behavior of the investigated chutes was

Conclusions
This work was conducted to investigate and prioritize the influence of geometric and hydraulic parameters on chute behavior. Five important parameters were selected: slope, discharge, pool height, chute height, and chute width. These parameters were then utilized to develop 600 cases, employing a highly accurate 3D numerical model. Results were used for RSM design, which resulted in two developed regression models for predicting V d and L c with coefficient determination (R 2 ) of about 0.95, highlighting the high accuracy of the models. Employing the analysis of variance (ANOVA), prioritization of the influence of the studied parameters on the dynamic behavior of the investigated chutes was accomplished in which Q, W, θ, H, and h p had the most significant influence on V d , respectively. The same order was also reported for L c ; however, h p had greater influence than H. Results revealed that the H parameter played an insignificant role in detecting L c and can be neglected. By employing the regression model for L c , the shear stress distribution at the critical section of the downstream in terms of the maximum velocity was investigated and compared to that of the last step and sidewalls. For θ ≤ 14.6 • , the greatest value of shear stress was formed on the sidewalls compared to the last step and downstream; however, the difference between downstream and side walls was not significant as compared to the sidewalls and the last step's edge. For θ ≥ 21.8 • , the maximum shear stress was formed at the downstream, and the difference between the shear stress on the last step's edge and side walls was not significant. Additionally, the shear stress in the longitudinal direction of the downstream showed that the behavior pattern of shear stress did not depend on the geometric and hydraulic parameters, and so an approximately similar pattern was observed in all cases. The findings in this study prove that these methods could effectively investigate and prioritize the parameters of a structure to improve the design and performance. Data Availability Statement: All data, models, or codes that support the findings of this study are available from the primary author of this paper upon reasonable request.