Evaluation of Flow Resistance in Gravel-Bed Channels with Bed-Load Transport by using Multi- Gene Genetic Programming


 Evaluation of flow resistance is necessary for the computation of conveyance capacity in open channels. For the case of channels with bedload conditions, the significance of the friction factor is of paramount importance. The gravel-bed channel response for a bedload transport condition in terms of flow resistance is distinct from that of a fixed bed and requires a different evaluation technique. The paper studies the different empirical approaches by various authors to determine the friction factor under bedload transport conditions and proposes an expression using genetic programming to derive the same. Flow resistance in the bedload transport condition is affected by various hydraulic as well as geometric parameters. In the present study, the factors proposed to be influencing friction factor for such a flow condition are undertaken as bed slope, relative submergence depth, aspect ratio, Reynolds number, and Froude number. A wide range of experimental datasets is used to determine the effect of these influencing parameters. The predictability of the proposed model is compared to the various empirical equations in literature. The proposed model is observed to effectively predict the friction factor for a wide range of datasets, unlike the existing models. The evaluated value of the friction factor from different models is used to validate the conveyance capacity of a river. The developed Multi-Gene Genetic Programming (MGGP) model is observed to reasonably predict discharge in the river, signifying that the model is competent to be applied to field study, within the specified range of parameters.


Introduction
The bedload in an open channel usually refers to the sediment particles entrained into the turbulent ow of water and deposited again at a short distance in the ow direction. The movement of bedload depends on the strength of ow and the characteristics of the sediment particles. During movement, sediment particles offer resistance to the ow, often estimated as different roughness coe cients. A reliable estimation of the roughness coe cient is necessary for hydraulic analysis of the streams and evaluating change in conveyance for a particular ow condition. The measurement of ow resistance for bedload transport conditions is di cult as it depends on the strength of the ow, complex secondary ow structures, bed shear stress distribution, bed particle sizes, etc. Extensive research has been carried out on gravel-bed streams, with movable channel bed and various expressions for Darcy-Weisbach friction factor f, but for speci c geometrical and hydraulic conditions, limited to the respective experimental conditions. Assessment of the effect of bedload transport and its impact on the ow characteristic has acquired signi cance in the past few decades. Researchers such as Smart and Jaeggi (1983); Rickenmann (1990); Song et al. (1998); Omid et al. (2003); Gao and Abrahams (2004); Mahdavi and Omid (2004), and Campbell et al. (2005) are a few that have analyzed the effect of moving bedload on the ow resistance. Einstein and Barbarossa (1951), Van Rijn (1982), and Wu and Wang (1999) have considered the effect of bedload in terms of an extra resistance caused due to bedforms in the case of the alluvial rivers. Meyer-Peter and Muller (1948) proposed the use of correction factors to take into account such effects. A direct effect of bedload on the increase in ow resistance with slope (varying from 3-9%) was reported by Bathurst et al. (1982a) by calculating resistance encountered by ow over a mobile bed; compared with a xed bed condition. Song et al. (1998) observed that the movable bed increases the ow resistance and is dependent on the concentration and the size of the moving particles. The bedload transport was observed to increase the friction factor by up to 50%, causing an increase in depth by 14%. The widely accepted view is that change in ow resistance response in mobile bedload is based on the dissipation of kinetic energy due to the collision of sediment particles and bedload, extracting momentum from the ow. Wiberg and Rubin (1989) observed that the ow resistance associated with bedload transport could reach much higher values in upper plane bed conditions than those observed with xed bed conditions. Carbonneau and Bergeron (2000) concluded that the introduction of suspended sediment into a clean water ow can alter the turbulence intensity depending on the relative magnitude of ow and sediment transport variables. Abrahams (2004, 2005) compared bedload friction factor in mobile and xed beds and found that resistance to ow in moving bedload is more than that in xed beds. They observed that the friction factor due to the bedload movement was 22.06% of the total resistance due to the extraction of momentum from the ow. They reported that resistance to ow is in uenced by factors such as sediment concentration, relative submergence depth, particle diameter, bed slope, and Froude number. Recking et al. (2008) reported that friction factor f increases with the sediment concentration in mobile bed conditions. Movable bedload conditions affect the variations in the in uencing parameters, thus complicating the determination of ow resistance by employing analytical methods. In such a situation, evolutionary computing based on arti cial intelligence (AI) is increasingly used in various water resources engineering branches for modeling non- Gene programming is a type of arti cial intelligence-based method mostly applied to solve technical problems in different engineering elds. Genetic programming (GP) is based on the Darwinian theory of natural selection used as an alternate soft computing technique. According to the classi cation of modeling techniques based on colors by Giustolisi et al. (2008), it is de ned as the next generation of soft computing techniques, whose meaning is related to the three levels of prior information required, white-, black-, and grey-box models.
GP and its variant multi-gene GP (MGGP) can be classi ed as grey box techniques. Higher the physical knowledge used during the model development, the better the physical interpretation of the phenomenon that the model allows the user.
The MGGP modeling is acquiring importance and being used as an advancing computational technique in water resources engineering, for hydrological and hydraulic investigations. Azamathulla  (2019) discussed implementing MGGP modeling to predict ow parameters. The MGGP approach at present has become a widely accepted computational tool in many disciplines, including civil engineering (Gandomi et al. 2011), neurocomputing (Peng et al. 2013), and electrical engineering (Mekhamer et al. 2012).
The main advantage of genetic programming and its variant multi-gene genetic programming (MGGP) over traditional statistical or regression methods and even other soft computing techniques is its ability to develop a compact and explicit prediction equation in different models variables without assuming a prior form of the existing relationship. The main advantage of genetic programming over regression and other soft computing techniques is its ability to generate a simpli ed prediction equation without considering a prior form of the existing relationship.
This paper proposes a Multi-Gene Genetic Programming (MGGP) technique to develop a reliable model for predicting ow resistance expressed as the Darcy-Weisbach friction factor for moving bedload conditions. The model is developed using the observed experimental data of various researchers and compared with traditional friction factor expressions. The developed model is veri ed with the authors' experimentations at the hydraulics laboratory of the National Institute of Technology Rourkela, India, and even for a natural river.

2.1.1Gravel-Bed Laboratory Datasets
In total, 774 experimental ume datasets observed under equilibrium bedload transport conditions (i.e., no erosion or deposition of sediment particles on the bed occurs) have been undertaken in this study using the selection criteria of particle diameter higher than 2 mm and less or equal to 64 mm (i.e., gravel), sediment geometric standard deviation less than 1.6, and no formation of bedforms as per Liu's bedform classi cation (Julien, 2010). Gilbert (1914) performed experimental studies and provided 104 datasets by varying the particle size, ume dimensions, and bed slope. Cassey (1935) conducted 42 experimental observations by changing the particle size and channel bed slope. 174 combinations of experimental observations were performed by Mavis (1937) for ve different particle sizes on a xed channel width but varying longitudinal slope of the channel. Bogardi (1939) carried out 44 laboratory studies by varying width, the channel's slope, and sediment particles. Ho Pang Yung (1939) and Meyer-Peter and Muller (1948) provided 24 and 59 datasets of bedload transport conditions, respectively. Paintal (1971) conducted 37 sets of laboratory studies by changing the ume's bed slope for the different particle sizes of bed materials. Smart and Jaeggi (1983) carried out 24 experimental studies of bedload transport conditions by varying the bed-slope for three sizes of bed material. Rickenman (1990) performed 46 experiments with similar parameters as Smart and Jaeggi (1983), but for different slopes of the channel. Cao (1985) and Graf and Suszuka (1987) conducted 35 and 106 experiments respectively, on a 0.6 m wide ume by varying the bed-slope for different gravel sizes. Recking (2006) carried out ume experimentations on 79 combinations of discharge,bed-slope, and particle size of the bed and analysed the observational datasets of bedload conditions compiled from various laboratory studies. The different ow conditions for the above datasets have been tabulated in Table 2. These experiments have been performed in umes of different widths, and where umes are not su ciently wide, drag due to sidewall has a signi cant effect on ow resistance. Therefore it becomes necessary to delineate sidewall drag from overall ow resistance from each observation from different experiments to ensure that ow resistance is calculated from gravel bed only. Accordingly, each observational dataset was delineated for the sidewall correction on consideration of the distinct exchange of momentum by gravel bed ansidewallll following Vanoni and Brooks (1975) procedure wherein the friction factor for sidewall was calculated iteratively using the Prandtl's friction law (Schlitchting 1979).Thus,the ction factor corresponding to the ow cross-sectional area obtained after delineation is only due to the response from the gravel bed. Thus screened datasets from different ume experiments have been upscaled to represent the uniform condition and combined to make it a representative dataset for the purpose of analysis. width of 0.9 m. A moving bridge arrangement is provided across the ume's width to facilitate taking measurements throughout the cross-section and along the length of the ume. Water from the overhead tank at the upstream end ows into the ume, which consequently discharges into a volumetric tank, nally owing into an underground sump from where it is recirculated. The volumetric tank measures actual discharge owing into the ume for every experimental condition. The schematic diagram of the water recirculation system is illustrated in Figure 1, with the channel photograph shown in Figure 2.

Experimental Procedure
Three experiments are carried out in the ume by fabricating straight rectangular channels of width 0.9 m.
Sediment feeding is done using a feeding tank tted with a small conveyor belt equipped with an RPM controller for controlled feeding of sediment into the ume. The transported sediments are collected at the downstream end, dried, and recirculated back into the ume. The condition of plane bedload transport is achieved by varying the ume discharge, and bed slopes for a particular gravel size and hydraulic parameters are noted. Two uniformly graded bed materials are considered with mean diameters of 2.93 and 7.72 mm. With a grain diameter of 2.93 mm, plane bedload transport was observed with bed slope as 0.0015 and 0.002, represented as Series I and Series II, in this paper, respectively. The Series III represents the values of parameters for bed material of 7.72mm mean diameter for the condition of moving bedload observed at bed slope of 0.005. The hydraulic parameters are illustrated in Table 1.

Flow Resistance Equations For Gravel-bed Channels
Due to resistance of ow, a portion of the energy is lost to overcome it. Other factors, such as a change in velocity and ow depth, also lead to energy dissipation. The resistance models proposed by Chezy's, Manning, and Darcy-Weisbach uses this concept of energy loss due to resistance encountered in the direction of ow. The ow resistance is generally represented by the resistance coe cient in-lieu of the bed condition, such as Manning's roughness coe cient n, Chezy's coe cient C, and the Darcy-Weisbach friction factor f.
The Darcy-Weisbac relation is usually preferred in research experiments as it is non-dimensional, contrary to C L 1 / 2 T − 1 and n L 1 / 3 T and the results obtained for a given experimental set can be directly extended to other ow conditions.
The ow resistance equation for the Darcy-Weisbach friction factor f in the condition of uniform ow in open channels is represented as, where u * is the shear velocity expressed as u * = gRS 0 = τ o ρ , τ 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, S 0 is the longitudinal slope of the bed, and ρ is the density of water.
Velocity distribution along the ow depth in a two-dimensional channel is based on the consideration of boundary shear, also known as Prandtl-von Karman universal logarithmic velocity distribution, which is de ned as where u is the point velocity at a height z above the bed, κ is von Karman's universal constant, z 0 is the height above the bed where velocity is considered as zero according to the law of the wall.
For rough turbulent ow in the two-dimensional open channel when grain Reynold number, (Re*= ρu * D μ ) 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 full-depth and approximating full ow depth as hydraulic radius R, In the case of two-dimensional open channel ow, hydraulic radius R is assumed to be full ow depth and R ≫ 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) (1), the expression for friction factor expressed as Equation (5) can be written as where E is constant depending upon the cross-section of the channel.

Dimensional Analysis for Identifying Parameters In uencing Flow Resistance in Gravel-Bed Channels
According to Yen (2002), for uid-carrying cohesionless sediment owing in a straight channel, the hydraulic resistance is a function of ow characteristics, the channel cross-section's geometry, and the properties of the uid and sediment. For steady ow, the Darcy-Weisbach friction factor f for a mobile bed channel can be written as: f = Function of (ν, ρ, ρ s , D, σ, g, So, U, H , ς) (7) where, ν= kinematic viscosity, ρ = density of water, ρ s =density of sediment particle,, D = mean diameter of the particle, σ = gradation coe cient, g = acceleration due to gravity, S 0 = longitudinal bed slope, H = Depth of ow, U = Mean velocity of ow and ς = 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 simpli ed to make modeling more accessible. The ratio of speci c density of sediment ρ s and water ρ is represented by the speci c gravity, G. In this study, uniformly sized gravel particles are considered; hence, the gradation coe cient, σ, is almost constant and can be taken out of the above functional relation. The cross-sectional 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 ow with uniformly sized spherical sediment in a rectangular channel can be expressed as f = Function of (ν, 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) In the present study, above selected ve non-dimensional terms i.e. longitudinal bed slope, S 0 ; relative submergence height R D ; aspect ratio δ = W H ; Reynolds number (Re); and Froude number (Fr) are considered to be in uencing 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 nd the effect of various parameters on the friction factor for mobile bedload conditions.  Figure 3 illustrates the effect of these parameters on friction factor, f. The four sub-gures show the direct variation of friction factor for four in uencing parameters with different bed slope values. 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 power-law 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 in uenced by selected geometric and hydraulic parameters, i.e., S 0 ; R/D; W/H; Re; and Fr. Julien (2006) proposed the following formulation for friction factor as: Recking (2006) carried out laboratory ume studies by considering various combinations of bed-slope 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 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 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 ume datasets, as presented in Table 3.

Multi-Gene 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 scienti c problems.
In addition to the standard GP evolution mechanisms discussed, some unique MGGP crossover mechanisms allow the exchange of genes between individuals, called Two-point High-level crossover and Low-level crossover. MGGP provides six mutation methods for the genes, such as sub-tree 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 de ned by the user. Since the trigonometric, Boolean, and logical operators would not help in de ning 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 ve in uencing 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. Multi-Genic 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   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 R 2 values are obtained on excluding any one of the dimensionless terms. Thus signifying that each of the ve selected dimensionless terms has a signi cant in uence on the friction factor. Hence, the rst model in Table 4 is taken as the best model, and Figure 4 represents its expression tree.
The MGGP has been used to predict the friction factor f. The in uential parameters derived from the data sets were used as inputs. Figure 5

Assessment of Equations for Calculating Flow Resistance in Gravel-Bed Channels
The coe cient of determination is possibly a suitable statistical parameter to evaluate the performance of different models. As can be observed in Figure 6, the R 2 value for almost all the researchers is less than 0.688, even as low as 0.396 for the case of Julien(2006). Whereas, the developed model has a high R 2 value of 0.999. The models by other researchers are established by taking into consideration their respective experimental observations. On the other hand, the proposed model in this paper considers the experimental observation of a vast range of datasets from the literature. Therefore, the proposed model has feasibility over a wider range of data sets.
To evaluate the overall performance of different models, four error analysis and one e ciency criteria were considered. The MAPE, MAE, MPE, RMSE, and R 2 , computed for each respective friction factor model, and presented in Table 6.
( ) ranged from 4.588x10 3 to 1.281x10 6, and Froude number ranged from 0.42 to 4.57. The MGGP model provided a better prediction of friction factor f compared to other models because the MGGP is capable of mapping the complex non-linear relationship between the various parameters involved in the phenomenon. Other models' performance could not found to be satisfactory for the global geometrical and hydraulic input datasets. However, compiling all the in uencing geometrical and hydraulic, the MGGP model proved to be the best for friction factor f for bedload conditions. Nevertheless, when all the data are analysed together and compared to the MGGP model, the model is observed to provide promising predictions.
The above analysis provided an overall understanding of the feasibility of different models. It was observed that the expressions for friction factor from the literature have a large error when used in the datasets undertaken. It can be postulated that these models would be valid for a speci c range of datasets. Therefore, to understand every model's validity to various datasets, error analysis is carried out for all the assumed datasets of this paper, using the different models. Such a study would provide a generalized understanding of the viability of every model.

4.2Application of the proposed MGGP model in a River
An explicit mathematical formulation can be considered as enhanced if it could reasonably predict the parameters on the eld. The proposed MGGP model, developed for predicting friction factor for mobile bedload conditions, is applied to a eld study to observe its performance.
River Koina is a tributary of river South Koel in Odisha. This river originates at Bhangaon in the state of Odisha and ows through West Singhbhum district in the state of Jharkhand, India. The river covers a total distance of 83km before it gives off into the South Koel River. A stretch of the river Koina located approximately 0.5km west of Manoharpur Railway Station. is undertaken for the eld study. The aerial photograph of the river has been shown in Figure 8. The river, during its course, transports quite a lot of bedload of varying sizes. The mean diameter of the transported bedload is analysed and found to be around 40mm, which lies in the range of the present study. The channel geometry and other hydraulic features are analogous to the present study for which the stretch of this river is undertaken for the eld study.
The eld survey was conducted at a straight planform of the river site to determine the discharge at three sections approximately 50m apart. A wire rope was attached across each of the cross-sections (S1, S2, and S3) along which the experimental observations undertaken with the help of a boat (Figure 9). At 2m intervals, along the cross-section, the depth of ow and depth-averaged velocity was measured using a measuring staff and a Flow tracker. For the measurement of depth-averaged velocity, the average velocity readings at 0.2 and 0.8 times the depth of ow were considered. The discharge capacity of the channel was consequently determined by the Velocity-Area method. The longitudinal slope of the river was measured by taking the reduced levels of various points along the centreline and determined to be 8.046x10 4 . Sieve analysis was carried out to determine the mean diameter of the bed particles and found to be 40 mm. The measured geometric and hydraulic parameters are illustrated in Table 8. Field measurements in rivers require the evaluation of discharge capacity, and any proposed method must reasonably perform well in these conditions. Therefore it is necessary to check whether the proposed model can predict the discharge capacity of a stream. Discharge is calculated by employing the developed MGGP model (Eq. 24) for friction factor, f in the Darcy-Weisbach Equation.
It is assumed that the modeled friction factor is no longer dependent only on the roughness factor but the other geometric and hydraulic factors incorporated into the MGGP model. Different models analysed in the previous section is applied to the three surveyed sections of river Koina. Table 9 illustrates the error analysis of discharge capacity in terms of MAPE, RMSE, MAE, and MPE.  Cao(1985) has over predicted the discharge by 9.79%. The models suggested by Recking(2006), and Recking et al. (2008) have performed well, but the MGGP model is observed to successfully predict the discharge capacity for a natural section, much better than the other traditional models. It is relevant to state that the proposed model is developed by undertaking ve geometric and hydraulic parameters and appropriately applied to real cases within the prescribed ranges.

Conclusions
The present paper proposed a model based on MGGP for the prediction of friction factor in bedload transport conditions. Data from previous ume experimental studies covering a wide range of ow conditions were utilized to develop a model for predicting ow resistance in two-dimensional open-channel ows moving over a mobile bed. The model developed was veri ed for laboratory experiments and real eld situations. This study aims to assess the impact of in uential parameters and develop a convincing model for the prediction of friction factor f using MGGP modeling. From the statistical analysis, the effectiveness of MGGP modeling was demonstrated. The main reason for the high level of accuracy lies in genetic programming's capability, which uses a genetic algorithm to solve optimization problems by creating a group of possible solutions to the problem. The comparison of correlation plots for different models shows that predicted values from the MGGP model are tted very closely with the actual value and with better accuracy. The coe cients of determination, R 2 for prediction of friction factor f using equation suggested by Meyer-Peter and Muller (1948), Cao (1985), Recking (2006) can be concluded that the other models could not satisfactorily predict frictional ow resistance. Although the Equation suggested by Cao (1985) and Julien (2006) had a reasonable better RMSE and R 2 values, but for the entire global and individual data sets, the MGGP proved to be a more promising model. Therefore MGGP modeling can be used to investigate the framework of one of the important real-world sedimentation problems involved in the river.
The proposed model is valid for the range of datasets undertaken in the study. The application of the proposed model, beyond the range of the parameters undertaken, needs to be performed taking proper precaution. The present model can be further improved by incorporating sinuosity and eddy viscosity as geometric and hydraulic parameters to take a more generalised form for solving real eld problems.
In engineering applications, errors in determining governing parameters affect the accuracy of the prediction of ow resistance. Experimentation has inherent errors that cannot be neglected. Thus, to obtain an accurate prediction of the friction factor, precise determination of ow properties such as slope, velocity, and cross-section is required. The variability of the friction factor with different parameters is also essential to an understanding of the process.

Availability of Data and Materials
All data, models, and code generated or used during the study appear in the submitted article. In detail, data can be found from the relevant literature.        Coe cient of Determination for friction factor by different approaches  Sample photographs during velocity measurement and cross-section at Koina River, India Highlightsforreview.docx