Modeling open channel flow resistance with dune bedform via heuristic and nonlinear approaches

Determination of the flow resistance (i.e. roughness coefficient) in open channel hydraulics is a substantial and intransitive issue in the design and operation of hydraulic structures, the calculation of water depth, estimation of flow velocity and an accurate characterization of energy losses. Assessment of the roughness coefficient is not a trivial matter given the multitude of factors influencing roughness (e.g. bed material, bedforms, cross-sectional and plan form variability, vegetation, etc.). The complexities and uncertainties of bedform configurations in alluvial channels continue to pose challenges for engineers, which stem from a variety of bedform shapes that arise under different flow conditions. Starting from a plane bed without sediment transport, ripples, dunes and washed-out dunes develop in large experimental flumes as the flow intensity increases in magnitude over a bed of loose sand particles (Figure 1). Meyer-Peter & Mueller (1948), Einstein & Barbarossa (1952), Taylor & Brooks (1962), Raudkivi (1967), Richardson & Simons (1967), Smith (1968), Van Rijn (1984), Karim (1995), Yang et al. (2005), and Van der Mark et al. (2008) have all presented expressions for the total roughness coefficient due to bedform roughness. Therefore, a variety of analytical and semi-empirical approaches with both simple and complex structures have been developed to predict the roughness coefficient. Engel & Lau (1980) found that dune length to depth ratio has a significant effect on the roughness coefficient when the dunes are steep. Karim (1999) proposed a new method for predicting relative dune height in a sand-bed stream based on the concept of relating energy loss due to form drag to the head loss across a sudden expansion in open channel flows. Heydari et al. (2014) proved that by increasing the Shields number, the ratio of the Manning's coefficient related to dune bedforms and the total Manning's coefficient increased with a logarithmic trend. However, the existing equations rely on a limited database, untested model assumptions, and a general lack of field data, and do not show the same results under variable flow conditions. These issues cause uncertainty in the prediction of roughness coefficient; therefore, it is extremely critical to utilize methods which are capable of predicting roughness coefficient for the channels with bedform under varied hydraulic conditions.

On the other hand, artificial intelligence tools (e.g. artificial neural networks (ANNs), neuro-fuzzy models (NF), genetic expression programming (GEP) and support vector machine (SVM)) deal with the representation and generalization of datasets using data learning techniques. These methods are remarkable forecasting tools which in the past decade have been implemented in various fields of civil engineering, hydraulic and hydrological studies (e.g. Samui & Dixon 2011; Roushangar et al. 2014a, 2016a, 2016b; Sun et al. 2014; Roushangar & Koosheh 2015; Nourani et al. 2016). Such studies include the prediction of total bed load (Chang et al. 2012), prediction of suspended sediment concentration (Kisi et al. 2013), prediction of sediment transport in circular channels (Roushangar & Ghasempour 2017), modeling sediment transport (Bhattacharya et al. 2007), modeling the rainfall-runoff (Nourani et al. 2016), and prediction of total bed material load (Roushangar et al. 2014a).

The GEP as an artificial intelligence approach uses fundamental principles of the genetic algorithms (GA) and genetic programming (GP) and mimics the biological evolution to make a calculation algorithm for predicting a specified phenomenon. The GEP approach has been applied in modeling various components of water resource systems including: developing stage-discharge curves, predicting energy dissipation over spillways (Roushangar et al. 2014b), determining the rainfall-runoff process (Kisi et al. 2013), modeling bridge pier scour (Azamathulla et al. 2010), and examining suspended sediment modeling (Aytek & Kisi 2008).

Computation of the river stage and flow velocity relies on the determination of bedform roughness; therefore, knowledge of the bedform is very important. Predicting the roughness coefficient is a complex phenomenon due to the nonlinearity and uncertainties of the process. Considering the reviewed literature, there is a lack of research on the comprehensive study of predicting the roughness coefficient in alluvial channels with dune bedforms using artificial intelligence. Therefore, the present research used GEP as a heuristic approach to assess the total Manning's coefficient (which includes both grain friction and form resistance) in channels with dune bedform, and also to determine the effective parameters in the estimation of the roughness coefficient in open channels with dune bedform. This method can be successfully applied to predict any variable of interest where the interrelationships among the relevant variables are poorly understood, finding the size and shape of the ultimate solution is difficult, and conventional mathematical analysis methods do not (or cannot) provide analytical solutions (according to Banzhaf et al. 1998). The simulations were performed using four different data sets obtained from experimental studies in different laboratories (one original experiment was performed by the author and three other experiments by other researchers). Modeling for prediction of the Manning's coefficient was carried out under three scenarios and the impacts of flow, bedform, and sediment characteristics were assessed. Then, the outcome of GEP was compared with the results obtained from nonlinear equations. In addition to this analysis, since the knowledge about bedform geometry can provide useful information for hydraulic engineering research and also for investigating the capability of the GEP method in predicting a new output variable, the relative dune height (H/y) and dune steepness (H/L) were simulated.

GEP as a nonlinear approach was developed by Ferreria (2001), using fundamental principles of the GA and GP. GEP is a heuristic search and optimization technique that mimics the biological evolution to make a computer program to predict a specified phenomenon. The problems are encoded in linear chromosomes with fixed-length as a computer program (Ferreria 2001) which are then expressed or translated into expression trees (ETs). One strength of the GEP approach is its unique, multigenic nature, which allows the evolution of more complex programs composed of several subprograms. A GEP algorithm begins by selecting the five elements: the function set, terminal set, fitness function, control parameters, and stopping condition. Figure 2 shows the overall process of GEP algorithm. According to Figure 2, the process begins with the random generation of the chromosomes of the initial population. Then the chromosomes are expressed and the fitness of each individual is evaluated. The individuals are then selected according to fitness to reproduce with modification, leaving progeny with new traits. The individuals of this new generation are, in turn, subjected to the same developmental process: expression of the genomes, confrontation of the selection environment, and reproduction with modification. The process is repeated for a certain number of generations or until a solution has been found (Ferreria 2001). In fact, in GEP the first major step is to identify the set of terminals to be used in the individual computer programs. The major types of terminal sets contain the independent variables of the problem, the state variables of the system and the functions with no arguments. The second step is to determine the set of functions (+, –, ×, /, x^a, log(x), ln(x), sqrt(x), etc.). The third step is fitness measure, which identifies the method of evaluating how well a given program solves a particular problem. In the current study, root mean square error (RMSE) is taken as fitness function. The fourth step is the selection of certain parameters to control the runs. The control parameters contain the size of the population, the rate of crossover, etc. The last step is the determination of the criteria to terminate the run. In a subsequent step, a comparison between the predicted and actual values is completed. When the desired results are in accordance with the error criteria initially selected, the GEP process is terminated. If the desired error criteria could not be achieved, some chromosomes are chosen by a method called ‘roulette wheel sampling’, and the selected chromosomes are then mutated to obtain new chromosomes. After the desired fitness score is found, this process terminates and the chromosomes are then decoded for the best solution to the problem (Teodorescu & Sherwood 2008).

GEP models utilized. The quality of the data driven model does not depend only on the selected input features. It can also be related to the parameters of the model (GEP). In this study, GEP has been trained for the Manning's coefficient prediction in alluvial channel with dune bedforms. One benefit of genetic programming is the possibility of generating explicit formulae. In this study, four basic arithmetic operators (+, –, ×, /) and some basic mathematical functions (⁠⁠, ⁠, ⁠, ⁠, e^x) were utilized as the GEP function in order to develop simple equations. The architecture of the chromosomes, including number of chromosomes (25-30-35), head size (7–8) and number of genes (3–4), were selected and different combinations of the mentioned parameters were tested. The model was run for a number of generations and was stopped when there was no significant change in the fitness function value and coefficient of correlation. It is observed that the model with number of chromosomes of 30, head size of seven, and number of genes of three yielded better results. Addition and multiplication were also tested as linking functions and it was found that linking the sub-ETs by addition represented better fitness values. One of the important steps in preparing the GEP model is to choose the set of genetic operators. In the current study a combination of all genetic operators (recombination, mutation, transposition, and crossover) was used for this aim. Mutations can occur anywhere in the chromosome. However, the structural organization of chromosomes must remain intact. In GEP there are three kinds of recombination: one-point, two-point, and gene recombination. In all cases, two parent chromosomes are randomly chosen and paired to exchange some material between them. The transposable elements of GEP are fragments of the genome that can be activated and jump to another place in the chromosome. There are three kinds of transposable elements in the GEP model: short fragments with a function or terminal in the first position that transpose to the head of genes, except to the root (IS elements); short fragments with a function in the first position that transpose to the root of genes (root IS elements or RIS elements); and entire genes that transpose to the beginning of chromosomes (Ferreria 2001). Also, the inversion operator is restricted to the heads of genes. Here any sequence might be randomly selected and inverted. The certain rates of genetic operators which define a certain probability of a chromosome were determined based on a trial and error process. Each GEP model was evolved (trained model) until there was no significant change in the fitness function value, then the program was stopped (fixed model). Therefore, the GEP model optimized parameters were determined. The order of the data sets was selected in such a way that the training data set contains a representative sample of all the behavior in the data in order to obtain the model with higher accuracy. To find a good training set which can give a good accuracy both in training and testing sets, one method is the instance exchange (Bolat & Yildirim 2004). The process starts with a random selected training set. Parameters of the optimized GEP model are shown in Table 1. It should be noted that in this study, for all models, 75% of the entire dataset was used for training goals and the remaining 25% of the data were used for testing goals. The same training set and the same testing set were used for all models.

There are different concepts and approaches that are used in the derivation and extraction process of roughness coefficient formulae, where some are simple and some are complex, such as linear and nonlinear approaches. Based on the linear superposition approach, the roughness coefficient is separated into the friction resistance and form resistance components. The components of the roughness coefficient vary according to the bed conditions (Yen 2002). In linear modeling of Manning's coefficient in movable beds, it is necessary to calculate both plane bed resistance and form resistance. This is a complex and time-consuming process. In the nonlinear approach, the resistance factor is not separated into grain roughness (friction resistance) and bedform roughness (form resistance), as in the linear superposition approach. Instead, it is kept as a single factor. Most of the available nonlinear approaches are developed based on dimensional analysis and statistical fitting of data to the parameters considered in the functional relationships. According to the reviewed literature, there are no comprehensive roughness coefficient formula specified only for alluvial channels with dune bedform. Due to this fact, the nonlinear approach is selected based on the classification of Yen (2004) which can be divided into four groups: (a) those that consider the roughness coefficient as a dependent variable; (b) methods based on shear stresses; (c) equations based on the mean velocity; and (d) the energy approach derived from Bagnold's (1966) stream power concept for expressions of velocity that can be used to calculate resistance. Most of these methods implicitly assume the flow to be steady, reachwise uniform with equilibrium sediment transport. Furthermore, among the four groups, only group (a) considers the roughness coefficient explicitly, whereas the other three groups consider the roughness indirectly and implicitly.

Table 2 represents the used nonlinear equations in this study. The approach developed by Strickler (1923) considers only the median grain diameter (D₅₀). The formula developed by Bruschin (1985) was designed by considering the grain and total Darcy–Weisbach roughness coefficient (f₀ and f), median grain diameter, the energy slope gradient (S_f), and hydraulic radius (R). Both approaches consider the D₅₀ as the main parameter. Similar to Bruschin's formula, Karim (1995) developed a formula without using S_f valid for Froude numbers below 0.4. For values higher than 0.4, different formulae were used. On the other hand, the approach developed by Camacho & Yen (1989) took advantage of the Reynolds (Re) and Froude numbers together with R and D₅₀. The formula varies for Froude numbers between 0.4 and 0.7 and between 0.7 and 1. Griffiths (1981) studied the hydraulic resistance in coarse gravel-bed rivers. Based on dimensional analysis and statistical justifications, he proposed relations for gravel rivers with rigid and movable beds using 136 field data. Yalin (1977) proposed that for subcritical flows, the roughness coefficient is a function of the grain Reynolds number, mobility number and grain size to depth ratio (see Table 2). Garde & Ranga Raju (1966) developed a roughness coefficient equation based on velocity. In this approach, the mean velocity was viewed as an implicit representation of the resistance. White et al. (1987) investigated the roughness coefficient in alluvial channels based on mixed momentum through an energy approach, resulting in a proposed equation for calculating the Darcy–Weisbach roughness coefficient. In this study, the output variable is parameter n, therefore, the equation was used to convert the dependent variables.

From extensive laboratory experiments at Colorado State University by Simons & Richardson (1963, 1966), a classification of several bedforms types was performed in terms of governing non-dimensional parameters. Flat bed, or plane bed, refers to a bed surface without bedforms. Ripples are small bedforms with wave heights less than a few cm (∼3.048 cm) and only seen in cases of hydraulically smooth bed conditions. Ripple shapes vary from nearly triangular to almost sinusoidal. Dunes are much larger than ripples and are out of phase with the water surface waves. Dunes are often triangular with relatively gentle upstream slopes and downstream slopes approaching the angle of repose of the bed material (Julien 2010).

The geometry of dune bedforms, such as their height and resulting decrease of water depth for navigation, can be a concern. Also their roughness might increase water levels which increase flood risk. The geometry of dune bedform refers to the representative dune height (H) and length of dunes (L) as a function of the average flow depth (y), median bed particle diameter (D₅₀), and other flow parameters, such as the grain shear Reynolds number. Figure 3 shows the dune bedform geometry.

The laws of conservation of energy and momentum must be considered for hydraulic resistive forces in the calculation of open channel hydraulics. Steady and uniform flow conditions are needed for the driving and resisting forces to be balanced, where flow is not accelerating nor decelerating, so the average channel cross-section, slope, and velocity are assumed to be constant under constant discharge conditions. In natural streams, velocity or discharge must often be predicted or calculated using other flow parameters, most commonly hydraulic radius, energy slope (or some approximation), and an estimate of channel roughness. Even though natural streams do not strictly comply with uniform flow assumptions, uniform flow conditions are often assumed to simplify velocity and discharge computations. The roughness coefficient in river engineering applications is commonly calculated using one of the three equations: Manning, Chézy, or Darcy–Weisbach (Yen 2002).

In this paper, the following four kinds of data corresponding to the roughness coefficient in open channels with dune bedform selected from published literature are used:

• (I)

  Williams (1970) organized several experiments that were made in channels with different widths and water depths in laboratories in Washington, DC. Sediment transport rates, grain size, water depth, and channel width were measured. Furthermore, water discharge, mean velocity, slope (energy gradient), and bedform characteristics were considered as the dependent variables.

• (II)

  As a part of the research programme of the Water Resources Division of the US Geological Survey, a project was organized at Colorado State University between 1956 and 1961 to determine the effects of the size of bed material, temperature of flow, and the fine sediment in the flow of the hydraulic and transport variables. The investigations for each set cover flow phenomena ranging from a plane bed with no sediment movement to violent anti-dunes (Guy et al. 1966).

• (III)

  The US Army Corps of Engineers Waterways Experiment Station in 1935 organized several experiments to research sediment transport.

• (IV)

  Roushangar (2010) studied the roughness coefficient in alluvial channels. New original experiments were carried out in the hydraulic laboratory of Caen University using a 5 m long, 0.15 m wide, and 0.25 m high rectangular flume. Natural quartz sand with a specific gravity of 2.65 and uniform average diameters of 0.15 and 0.4 mm were used as sediment particles in the experiments. By changing the flow depth and discharge, different average velocities, Froude numbers, dune height, and other hydraulic parameters were measured or calculated. The ranges of selected parameters used in the tests are listed in Table 3, where b represents the channel width. It should be noted that the measured values of the Manning's coefficient are calculated from Manning's famous equation which is simple in form and the formula is well proven by much practical experience.

Selection of various parameters as input combinations can affect the accuracy of the results throughout the modeling process. Thus, selection of appropriate parameters and input combinations plays a key role during modeling. Therefore, in order to determine the output variables (Manning's coefficient) under three different scenarios, a set of dimensionless variables as input combinations were attempted to develop the process. The scenarios considered in this study are as follows.

• (I)

  Models based on flow characteristics; Scenario 1

• (II) Models based on flow and bedform characteristics; Scenario 2

• (III) Models based on flow and sediment characteristics; Scenario 3

In the first scenario, models were defined according to the flow characteristics. Table 5 sums up the performance evaluation indices of the models with flow characteristics. According to Table 5, the first three combinations were included with only one variable. Results revealed that the parameter Re led to the best results amongst the single parameter models for the Manning's coefficient. The next two combinations were defined with two variables. The combination including Re and y/b proved to be the best double parameter GEP model with flow characteristics for predicting the roughness coefficient. According to the results of Table 5, it could be seen that using the Reynolds number as the only input parameter yielded an acceptable accuracy, while using only the Froude number as input did not lead to a desired prediction performance. Adding parameter y/b to input parameters caused a fair increase in model efficiency. It could be deduced that Re and y/b were the most efficient parameters in modeling the Manning's coefficient based on flow characteristics. Figure 5 presents the scatter plot of the observed data versus model results that uses these two flow characteristics as input. According to Figure 5, the lower values are predicted more accurately compared to the higher values of the Manning's coefficient. Additionally, based on trial and error, data sets were divided into two groups in terms of the Reynolds number (Re < 80,000 and Re > 80,000). Then, the best input combination (Re, y/b) was rerun for both data categories. Results are represented in Table 5 (last two rows) and Figure 6. These findings demonstrate that when the Reynolds number is less than 80,000, the model tends to be more accurate in predicting the Manning's coefficient. Figure 6 highlights the better model performance for Reynolds numbers less than 80,000 (Re < 80,000 category), in comparison to the higher range of the Reynolds number (Re > 80,000).

In order to consider both flow and bedform features in modeling processes, the models were defined based on flow and bedform characteristics. Considering the Manning's coefficient simulations with bedform characteristics, the GEP was fed with single, double and triple input variables. The first input combination with only Re as an input variable led to the best results for the prediction of the roughness coefficients. Moreover, adding H/L, H/y and L/H (⁠ respectively are the dune's height and length) to the Reynolds number (Re) had less impact on the results of Manning's coefficient predictions. The next three input datasets consisted of a double input combination. Among the mentioned input combinations, the one including Re and L/y ended with a better outcome. The last four datasets with more than two input combinations did not perform as well as the single input (Re). Results of models in the verification step for the Manning's coefficient with bedform datasets are presented in Table 6. Based on the outcome, the developed models with flow and bedform characteristics performed relatively weakly in comparison to developed models with only flow characteristics. Figure 7 presents the scatterplot obtained from developed models using flow and bedform characteristics for the best performing model (the first model).

Developed models based on flow characteristics along with modeling via flow and bedform characteristics demonstrated that the flow characteristics were capable of more precise predictions so far. Moreover, to consider the flow and sediment features in roughness coefficient estimation, flow and sediment characteristics were employed in the establishment of input combination. To measure the significance of the sediment features in the modeling of flow resistance, sensitivity analysis was performed for predicting the variables based on the performance of models. Results in Table 7 indicate that the input combination including Re, R/D₅₀ and led to the best consequence for the prediction of the Manning's coefficient. In other words, it could be inferred that the ratio of hydraulic radius to the sediment diameter and Reynolds number as the most effective input combination led to the best results in predicting the Manning's coefficient. Therefore, both hydraulic and sediment characteristics have a significant effect on the prediction process. From the results it is clear that adding (relative discharge) to Re and R/D₅₀ caused an improvement in the results of the Manning's coefficient. Looking into the data presented in Table 7, using Re and R/D₅₀ as the individual input variables shows almost equal performance as the three combination. However, the use of R/D₅₀ seems to be a more practical option since velocity measurements in the field situations are typically very difficult. The scatter plot of the observed and predicted Manning's coefficient for the best input combination with sediment characteristics is shown in Figure 8.

So far, several different studies have been implemented in the field of predicting the roughness coefficient. In many practical applications, it was mentioned that it is desirable and sometimes necessary to predict the ratio of dune height to flow depth (H/y). For this purpose, in the present study different parameters were used as input combinations to investigate their role in the modeling of relative dune height. Results of models are presented in Table 8. According to the listed performances in Table 8, it was observed that the combination of Re, y/b, and led to the best outcome for prediction of H/y (with the highest R and DC and the lowest RMSE and MAPE). Furthermore, these three parameters were also dominant parameters for predicting the ratio of dune height to its length (H/L).

By comparing the various input combinations, it can be concluded that adding to other parameters does not cause an increase in the accuracy of the results. According to the results, it can be inferred from Table 8 that the use of as input parameter caused an increment in models efficiency. Finally, according to the results, it is almost impossible to draw a sensible relation between H/y and each one of the input parameters solely. Resultantly, it can be assumed that the relative dune height has a suitable correlation with a combination of ⁠, and ⁠. In Figure 9, the scatter plot of the simulated relative dune height vs. the observed values for this best performing model for dune height prediction is presented. It can be seen that the model incidentally underestimates the high values of relative dune height. It is noteworthy to mention that other models listed in Table 8 yielded similar results for the high H/y values.

The key point of gene expression programming is that it is able to give the explicit expression of the relationship between the variables. It should be noted that in each run of the GEP method a unique formula will be obtained. Therefore, every model was run several times using different settings of GEP parameters and the best models (with minimum error) were selected. The mathematical expressions of the GEP best models for all cases are presented in Table 9. In this table, the d parameters represent the input parameters and c parameters are constant values which are determined by GEP. In this study, two constants per gene were selected for each model. Looking into Table 9, the equation of scenario 3 model is not only the best performing equation, but also the simplest relation compared to the other equations. Furthermore, it does not include any constants.

To investigate the impacts of different parameters of the GEP-best performing models on the Manning's coefficient, sensitivity analysis was performed. The significance of each parameter was evaluated by eliminating the parameter. Table 10 reveals that in the prediction of the Manning's coefficient with flow characteristics, the Reynolds number is the most efficient parameter. This also applies to the prediction of bedform and sediment feature characteristics given that the best performing model has Re as the only input parameter. Finally, in predicting the relative dune height for both dimensionless forms H/y and H/L, is the dominant parameter.

The experimental data presented in Table 2 was used to evaluate the applicability of several existing nonlinear equations for the Manning's coefficient in the channels with dune bedform. The overall performance of each equation is presented graphically in a plot of observed Manning's coefficient against their computed values. Four evaluation criteria (R, DC, RMSE, and MAPE) were used as indications of the accuracy of the equations. Then a comparison was performed among the GEP-best performing model and the equations. The nonlinear equations were divided in four groups: those that considered the roughness coefficient as dependent variable, methods based on shear stresses, equations based on the mean velocity, and the energy approach. The results of the comparisons are plotted in Figure 10. From the obtained results among nonlinear equations, Karim's formula that considered the roughness coefficient as a dependent variable predicted the roughness coefficient better in comparison to the others. It should be noted that in the Karim equation D₅₀, f and f₀ parameters were used as input variables. The parameters f and f₀ were taken from the measurements (Table 3), therefore, this equation presented good performance in comparison to the other equations. However, in most cases, the outcome of the equations did not provide reasonable results. On the other hand, the obtained results by GEP-best performing model were close to the measured data and were in good agreement. It should be noted that the existing equations are developed based on special flow conditions, therefore, their application is limited to those special cases of their development and they did not show uniform results under different conditions. The above-mentioned issue can be seen in Figure 10, which shows that the obtained results from equations differ from each other, and have less accuracy than GEP model which, due to its multigenic nature, allows the evolution of more complex programs composed of several subprograms. Based on Figure 10, GEP had the highest R and DC and the lowest RMSE and MAPE. This confirms the workability of GEP as an efficient approach in the modeling of the Manning's coefficient in channels with dune bedform.

In the present study, the performance of nonlinear approaches as semi-empirical equations and GEP model were evaluated in order to predict the Manning's coefficient in open channels with dune bedform. In order to develop GEP models, three scenarios were considered and the impact of different modeling for the estimation of the roughness coefficient was evaluated. Then, the outputs from the proposed GEP methods were compared to the nonlinear approaches. Results verified the capability of GEP in predicting the Manning's coefficient and it was found that for prediction of the roughness coefficient in channels with dune bedform, the third scenario that took the advantages of flow and sediment characteristics performed more successfully than the others. The obtained results revealed the accurate prediction via the GEP approach in comparison with nonlinear equations. According to the outcome, when predicting the Manning's coefficient by using flow characteristics, including the ratio of flow depth to channel width (y/b) together with the Reynolds number (Re), there was an improvement in the results compared to the other parameters (Fr, Re and y/b) when they were used as the individual input variables. Also, it was deduced that the prediction of the Manning's coefficient for Reynolds values of less than 80,000 was more precise than in those with higher Reynolds values. It was found that in the prediction of the Manning's coefficient using flow and bedform characteristics together, the single input of the Reynolds number was the highest influential parameter and plays a key role in the process. Regarding the Manning's coefficient with flow and sediment characteristics, Re, R/D₅₀ and played the most influential parameter role. Also, the results showed that using Re and R/D₅₀ as the individual input variables shows almost equal performance as the three combination. Since velocity measurements in field situations are typically very difficult, the use of R/D₅₀ seems to be a more practical option. In the estimation of the relative dune height (H/y or H/L), using the Re, y/b and produced improved GEP results. According to sensitivity analysis, the Reynolds number (Re) played the most important role in predicting the Manning's coefficient in channels with dune bedform. It should, however, be noted that the GEP is a data-driven model and the GEP-based models are data sensitive, so further studies using data ranges out of this study and preferably field data should be used to find the merits of the models to predict roughness coefficient in the real condition of flow.