Due to resistance of flow, a portion of the energy is lost to overcome it. Other factors, such as a change in velocity and flow depth, also lead to energy dissipation. The resistance models proposed by Chezy’s, Manning, and DarcyWeisbach uses this concept of energy loss due to resistance encountered in the direction of flow. The flow resistance is generally represented by the resistance coefficient inlieu of the bed condition, such as Manning’s roughness coefficient n, Chezy’s coefficient C, and the Darcy–Weisbach friction factor f.
The DarcyWeisbac relation is usually preferred in research experiments as it is nondimensional, contrary to C \(\left[{\text{L}}^{1/2}{\text{T}}^{1}\right]\) and n \(\left[ {\text{L}}^{1/3}\text{T}\right]\) and the results obtained for a given experimental set can be directly extended to other flow conditions.
The flow resistance equation for the DarcyWeisbach friction factor f in the condition of uniform flow in open channels is represented as,
\(\sqrt{\frac{8}{f}}\) = \(\frac{U}{\sqrt{gR{S}_{0}}}=\frac{U}{{u}_{*}}\) (1)
where \({u}_{*}\) is the shear velocity expressed as \({u}_{*}=\sqrt{gR{S}_{0}} =\sqrt{\frac{{\tau }_{o}}{\rho }}\), \({\tau }_{o}\) is the bed shear stress, U is the mean velocity, R is the hydraulic radius of the gravel bed, g is the acceleration due to gravity, S0 is the longitudinal slope of the bed, and \(\rho\) is the density of water.
Velocity distribution along the flow depth in a twodimensional channel is based on the consideration of boundary shear, also known as Prandtl–von Karman universal logarithmic velocity distribution, which is defined as
$$\frac{u}{{u}_{*}}=\frac{1}{{\kappa }} ln\frac{z}{{z}_{o}}$$
2
where u is the point velocity at a height z above the bed, \({\kappa }\) is von Karman’s universal constant, z0 is the height above the bed where velocity is considered as zero according to the law of the wall.
For rough turbulent flow in the twodimensional open channel when grain Reynold number, (Re*= \(\frac{\rho {u}_{*}D}{\mu })\) expressed in terms of the mean diameter of particle D exceeds 70, Keulegan (1938) expressed depthaveraged velocity by integrating Prandtl–von Karman universal velocity distribution along fulldepth and approximating full flow depth as hydraulic radius R,
$$U=\frac{{\int }_{Z0}^{R}udz}{{\int }_{Z0}^{R}dz}=\frac{{\int }_{Z0}^{R}\frac{{u}_{*}}{k} ln\frac{z}{{z}_{o}}dz}{R{z}_{o}}=\frac{{u}_{*}}{{\kappa }}\left[\left(\frac{R}{R{z}_{o}}\right)ln\left(\frac{R}{{z}_{o}}\right)1\right]$$
3
In the case of twodimensional open channel flow, hydraulic radius R is assumed to be full flow depth and \(R\gg {z}_{o}\) (height where velocity tends to zero) and converting base of natural logarithm to logarithm with the base as 10 (ln = 2.303*log), Equation (3) can be simplified as
$$\frac{U}{{u}_{*}}=5.75\text{log}\left(\frac{R}{{z}_{o}}\right)\frac{1}{{\kappa }}$$
4
Nikuradse (1933), verifying from experiments, has termed that surface absolute roughness is that height\({z}_{o}\) where velocity tends to zero and expressed \({z}_{o}\) in terms of the mean diameter of particle D, as \({z}_{o}=\frac{D}{30}\). Thus using, \({z}_{o}=\frac{D}{30}\), \({\kappa }=0.4\) and \(\sqrt{\frac{8}{f}}=\frac{U}{{u}_{*}}\)from Equation (1), the expression for friction factor expressed as
$$\sqrt{\frac{8}{f}}=5.75log\left(\frac{R}{D}\right)+5.94$$
5
Equation (5) can be written as
$$\sqrt{\frac{8}{f}}=5.75log\left(\frac{R}{D}\right)+E$$
6
where E is constant depending upon the crosssection of the channel.
3.1 Dimensional Analysis for Identifying Parameters Influencing Flow Resistance in GravelBed Channels
According to Yen (2002), for fluidcarrying cohesionless sediment flowing in a straight channel, the hydraulic resistance is a function of flow characteristics, the channel crosssection's geometry, and the properties of the fluid and sediment. For steady flow, the Darcy– Weisbach friction factor f for a mobile bed channel can be written as:
f = Function of (\(\nu\), \(\rho , {\rho }_{s }, D, \sigma , g , So ,\)U, H \(, \varsigma\)) (7)
where, \(\nu\)= kinematic viscosity, \(\rho\) = density of water, \({\rho }_{s}\) =density of sediment particle,, D = mean diameter of the particle, \(\sigma\) = gradation coefficient, g = acceleration due to gravity, S0 = longitudinal bed slope, H = Depth of flow, U = Mean velocity of flow and \(\varsigma\) = cross sectional shape factor .
An accurate prediction of the friction factor would require an explicit expression containing all the parameters in Equation (7). However, it can be simplified to make modeling more accessible. The ratio of specific density of sediment \({\rho }_{s }\)and water \(\rho\) is represented by the specific gravity, G. In this study, uniformly sized gravel particles are considered; hence, the gradation coefficient, σ, is almost constant and can be taken out of the above functional relation. The crosssectional shape factor of a rectangular channel can be represented by width, W, and hydraulic radius R. Thus, the functional relationship for the friction factor of steady uniform flow with uniformly sized spherical sediment in a rectangular channel can be expressed as
f = Function of (\(\nu , D ,G, g , {S}_{0} ,\)U, H\(, W, R\)) (8)
From the dimensional analysis, these parameters can be clubbed to form six important dimensionless terms, as (1) \(\frac{R}{D}\)= Relative submergence depth, (2) \(\frac{UR}{\nu }\)= Re (Reynolds number), (3) \(\frac{U}{\sqrt{gR}}\) = Fr (Froude number), (4) \(\frac{W}{H}\)=\(\delta\)(Aspect ratio ) (5) S0 = longitudinal slope, and (6) \(\frac{{S}_{0}R}{\left(G1\right)D}=\theta\) (Shields number\()\)
Froude, Reynolds, and Shields number are taken into consideration because flow in open channel with bedload transport condition is with a free surface, generally turbulent and moving bed material experiencing shear stress. However, it has been noticed that in the expression of Shields number, two already considered dimensionless terms \(\left(\frac{R}{D}\right)\) and S0 are appearing. Therefore, the Shields number, \(\theta\) has not been considered separately and thus, the above functional relation can be expressed as,
f = Function of (\(\frac{R}{D} ,{ S}_{0} ,\) \(\frac{W}{H} ,Fr , Re\)) (7)
In the present study, above selected five nondimensional terms i.e. longitudinal bed slope, S0; relative submergence height \(\left(\frac{R}{D}\right);\) aspect ratio\(\left(\delta =\frac{W}{H}\right)\); Reynolds number (Re); and Froude number (Fr) are considered to be influencing parameters affecting the variation of friction factor, f for the case of mobile bedload conditions. Table 3 presents a wide range of datasets facilitating the analysis to find the effect of various parameters on the friction factor for mobile bedload conditions.
Table 3
Influencing Parameters for the experimental datasets undertaken in the study
Datasets  Slope So (%)  Relative submergence height\(\left(\frac{R}{D}\right)\)  Aspect Ratio δ (W/H)  Froude no. Fr  Reynold’s No. Re  Discharge, Q (m3/s) 

Gilbert (1914)  0.62  6.38  1.18  0.86  14325  0.0051 
2.25  27.32  12.5  1.67  73814  0.0316 
Cassey (1935)  0.119  10.31  1.39  0.49  10473  0.0047 
0.509  75.2  14.29  0.93  101812  0.099 
Mavis (1937)  0.155  6.61  6.25  0.58  7036  0.006 
1.01  50.57  50  1.2  71782  0.078 
Bogardi (1939)  1.04  4.24  2.17  0.87  23681  0.016 
2.45  10.87  25.00  1.61  89381  0.064 
Ho Pang Yung (1939)  0.128  13.81  1.61  0.53  29593  0.016 
0.504  37.36  5.56  0.93  89833  0.069 
MeyerPeter and Muller (1948)  0.128  6.06  1.43  0.42  11211  0.005 
2.27  33.53  11.11  1.47  1281043  4.613 
Paintal (1971)  0.13  5.94  4.35  0.54  24853  0.02605 
1.00  61.76  14.29  1.04  193151  0.25484 
Smart and Jaeggi (1983)  3.00  3.67  2.38  1.42  19919  0.005 
20.00  14.14  8.33  4.57  93467  0.03 
Cao (1985)  0.50  1.30  2.78  0.85  55815  0.04 
9.00  15.26  10.00  1.52  201308  0.19 
Graf (1987)  0.50  3.84  2.32  0.84  53764  0.040 
2.50  17.43  8.33  1.30  185294  0.205 
Rickenmann (1990)  7.00  3.04  2.33  1.50  34472  0.01 
20.00  7.85  6.25  3.93  97059  0.03 
Recking (2006)  1  1.97  2.86  1  4588  0.006 
9  28.19  8.33  1.8  37500  0.015 
Figure 3 illustrates the effect of these parameters on friction factor, f. The four subfigures show the direct variation of friction factor for four influencing parameters with different bed slope values.
Figure 3(a) shows that for lower bed slopes, there is not much change in the friction factor for mobile beds with the change in relative submergence height\(\left(\frac{R}{D}\right)\); but for higher slopes, even a slight variation of\(\left(\frac{R}{D}\right)\), causes significant variations in the friction factor. On the other hand, the friction factor appears to show little variation with the aspect ratio, δ for lower slope, as seen in Figure 3(b). Similar observations are made on the effect of Reynolds number on friction factor at low bed slope, but the variations of f at lower slopes due to change in Re is not much as in the case of relative submergence height as seen in Figure 3(c); where the friction factor increases with the decrease in Re and the change in variation is observed to be similar for different conditions of bed slope. It is observed that there is a powerlaw variation of friction factor for the Froude number, as seen in Figure 3(c). Thus, from Figure 3, it can be inferred that friction factor f is influenced by selected geometric and hydraulic parameters, i.e., S0; R/D; W/H; Re; and Fr.
3.2 Empirical Flow Resistance Equations.
In the case for open channel flow where turbulence is fully developed and without sediment transport, Keulgan (1938) proposed the value of constant E in equation(6) as 6.25 for rough channel flow and the expression of friction factor formulated as,
$$\sqrt{\frac{8}{f}}=6.25+5.75log\left(\frac{R}{D}\right)$$
8
Manning and MeyerPeter Muller (1948) examined the flow resistance characteristics in case moving bedload and the friction factor expressed as,
$$\sqrt{\frac{8}{f}}=\frac{U}{{u}_{*}}=6.75{\left(\frac{R}{D}\right)}^{1/6}$$
9
Engelund and Hansen (1966) proposed the following formulation for finding out friction factor with \(\left(\frac{R}{{d}_{65}}\right)\) as a factor,
\(\sqrt{\frac{8}{f}}=4.27+5.75log\left(\frac{R}{{d}_{65}}\right)\), they defined f as \(\frac{U}{{U}_{*}}=\sqrt{\frac{2}{f}}\) (10)
Cao (1985) carried out experiments on different combinations of bedslope and gravel size and proposed logarithmic formulation for friction factor as
$$\sqrt{\frac{8}{f}}=3.75+5.91log\left(\frac{R}{D}\right)$$
11
Smart and Jaeggi (1993) proposed the following formulation wherein slope was also incorporated for finding out friction factor as:
$$\sqrt{\frac{8}{f}}=2.5\left[1\text{e}\text{x}\text{p}(0.05{Z}_{90}/{So)}^{0.5}\right]\text{ln}\left(8.2{Z}_{90}\right) with {Z}_{90}=R/{d}_{90}$$
12
Julien (2006) proposed the following formulation for friction factor as:
$$\sqrt{\frac{8}{f}}=5.75 log\left(\frac{2R}{D}\right)$$
13
Recking (2006) carried out laboratory flume studies by considering various combinations of bedslope varying from 1–9% and uniformly sized particle size of 2.3, 4.9, 9, and 12.5 mm to evaluate friction factor, f under equilibrium condition of moving bedload. In no bedload transport condition ,it was observed that\(\sqrt{\frac{8}{f}}\) increases with the relative submergence depth and becomes constant when bedload appears. In the case of high bedload transport, both the bedload and \(\sqrt{\frac{8}{f}}\) increases with the relative submergence depth. Different flow resistance equations were studied for a wide range of flow conditions, and new equations for two different range of relative submergence depth were formulated as:
$$\sqrt{\frac{8}{f}}=1+9.5\text{log}\left(\frac{R}{D}\right) for \frac{R}{D}cript>$$
14
$$\sqrt{\frac{8}{f}}=3.6+5.75\text{log}\left(\frac{R}{D}\right) for \frac{R}{D}>16.9$$
15
On further investigations, Recking et al. (2008) produced expressions for flow without and with bedload transport as:
\(\sqrt{\frac{8}{f}}=3.6+3.2\text{ln}\left(\frac{R}{D}\right)\) (without bedload) (16)
\(\sqrt{\frac{8}{f}}=0.67+3.2\text{ln}\left(\frac{R}{D}\right)\) (with bedload) (17)
The various expressions for friction factor have undertaken relative depth, shear velocity, hydraulic radius, bed slope, grain characteristics, etc. established according to the respective datasets. This paper attempts to customize as many parameters as possible into a single expression of friction factor by undertaking 774 laboratory flume datasets, as presented in Table 3.
3.3 MultiGene Genetic Programming Based Model for Predicting of Flow Resistance Equations
Gene Expression Programming (GEP), suggested by Azamathulla, Ahmad, and Ghani (2013), is a search technique that involves computer programs such as mathematical expressions, decision trees, polynomial constructs, logical expressions, etc. GEP was developed by Ferreira (2001) based on generation and evaluation of its suitability and has offered a wide possibility to solve complex scientific problems.
In addition to the standard GP evolution mechanisms discussed, some unique MGGP crossover mechanisms allow the exchange of genes between individuals, called Twopoint Highlevel crossover and Lowlevel crossover. MGGP provides six mutation methods for the genes, such as subtree mutation; mutation of constants; the substitution of a randomly selected input node; substitute a randomly selected constant; setting randomly selected constant to zero or as one.
A functional set may include arithmetic operators (+;; /;x; etc), mathematical functions (sin(); cos(); tan(); ln()), Boolean operators (AND, OR, NOT, etc.), logical expressions (IF or THEN) or any other suitable functions defined by the user. Since the trigonometric, Boolean, and logical operators would not help in defining a proper mathematical explicit expression in a generalized sense, only the arithmetic, logarithmic and exponential functions are used with constants and the input parameters as variables. GeneXproTools 5.0 is used for modeling the gene expression programming in this study.
The model developed in this study designates friction factor f as the output and the five influencing dimensionless terms in Eq. (7) as input. Four basic arithmetic operators (+, , *, /) and some basic mathematical functions (sqrt, ex, ln) were used as the function set in the model development. MultiGenic programming, i.e., with three genes and addition as the linking function, is used. A large number of generations were tested. The functional set and operational parameters used in MGGP modeling during this study are listed in Table 5. Bedload experimental datasets used in this study were distributed into both training and validation purposes. Out of the total 774 series of datasets, 545 series of data sets were used as training data, and 229 series of data sets were used for validation purposes. To decide about the most appropriate or best MGGP model based on the possible combination of an individual parameter, predicted f derived from five different models or equations are computed using MGGP. In the first model, all the five parameters are considered input to predict the friction factor, f. Similarly, four different models are developed by neglecting one input parameter in each case. The performance of each of the models has been evaluated by conducting statistical error analysis for different variations of five models. For this purpose, standard statistical indicators Mean Absolute Percentage Error (MAPE), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Percentage Error (MPE), and Coefficient of Determination (R²) for different models have been determined.
Table 5
Functional set and operational parameters used in MGGP Model
Description of Parameter

Parameter Setting

Function Set

+, , *, /, square, exp, ln

Number of Chromosomes

30

Head Size

8

Number of Genes

3

Linking Function

Addition

Fitness Function

RMSE

Program Size

32

Number of Generations

1197

Constants per Gene

10

Data Type

Integer

Block Mutation

0.00138

Table 4
Model Sensitivity of different MGGP Models for f
MGGP Models for f

MAPE

RMSE

MAE

MPE

R²

f =F(R/D,So ,W/H,Fr,Re)

0.6047

0.0008

0.0005

0.3555

0.9999

f=F(R/D, So ,Fr,Re)

57.0932

0.0427

0.0363

51.9099

0.9045

f =F(R/D, So ,W/H,Fr)

42.2142

0.0446

0.0342

15.8183

0.8478

f =F(R/D,W/H,Fr,Re)

24.5616

0.0482

0.0269

3.3520

0.8124

f=F(R/D, So ,W/H,Re)

33.1842

0.0514

0.0329

17.6761

0.7967

Table 4 indicates the error analysis of the prospective MGGP models for overall data sets. It is observed that higher MAPE, MAE, MPE, RMSE, and lower R2 values are obtained on excluding any one of the dimensionless terms. Thus signifying that each of the five selected dimensionless terms has a significant influence on the friction factor. Hence, the first model in Table 4 is taken as the best model, and Figure 4 represents its expression tree.
The simplified functional relationship is obtained as,
f = \(\frac{So}{\left(R/D\right)*(3.67 {Fr}^{2}+\delta +Fr)}+\frac{49\left(\frac{R}{D} So \right)}{\delta \left\{{\left(\frac{R}{D}\right)*Re+\delta }^{2}\right\}}+ \frac{7.93*So}{{Fr}^{2}}\) (18)
The MGGP has been used to predict the friction factor f. The influential parameters derived from the data sets were used as inputs. Figure 5(a) shows the comparison of predicted f versus the target f at the end of the training and found that the coefficient of determination (\({R}^{2}\)) is 0.9999. This indicates that the training of the MGGP model was satisfactory. Figure 5(b) plot the performance of the MGGP model for f at the end of the validation stage; the coefficient of determination, (\({R}^{2}\)), is 0.9999.