GEP- and MLR-based equations for stable channel analysis

A stable channel can be described by a balance of erosion and sedimentation processes that produce average dimensions, which do not change significantly over the years (Bonakdari et al. 2020). It is essential to estimate the dimension correctly to minimise the erosion and sedimentation effect (Garcia 2008). The prediction of a stable channel depends on various parameters and is described by the characteristics of the top width of the channel, the mean depth of the flow and the longitudinal slope (Bonakdari et al. 2020). For decades, researchers tried to correlate channel geometry to the geomorphological variables. According to Chang (1988) and Garcia (2008), a stable channel was determined by empirical and analytical methods. Blench (1957), Lacey (1930), Lindley (1919), Kennedy (1895) and Simons & Albertson (1960) investigated the regime concept by examining the channel stability parameters by using statistical methods. Leopold & Maddock (1953) introduced the concept of power function theory, and this theory has been the method of choice for many researchers. Channel geometry is described by the corresponding water and sediment discharge. The researchers divided this method into two categories, namely, downstream hydraulic geometry and at-a-station hydraulic geometry. By applying the hydraulic geometry method, Julien & Wargadalam (1995) extended this concept by introducing the relationship between the channel geometry and the full bank discharge (Q), the median size of bed material (d₅₀), the channel bed slope (S) and the Shields parameter (τ*). The width, mean depth and velocity depend on Q, d₅₀ and S, while the channel slope is dependent on Q, d₅₀ and τ*. Lee & Julien (2007) extended this research by using statistical analysis, and the new equations were found to be more accurate than the previous equations as a broader range of data was used in the study. Afzalimehr et al. (2010) investigated the relationship of channel geometry by using the at-a-station method and discovered that the median size of bed material and the Shields parameter did not contribute much to the width and depth. Other researchers, such as Kaless et al. (2014), Parker (1978), Parker (1979) and Pitlick & Cress (2002), use dimensionless width (B*), depth (y*) and dimensionless discharge (Q*) to describe the relationship of stable channel geometry. The variable B* is defined by B/d₅₀ and y* is presented by y/d₅₀, while Q* is defined by ⁠, where R is the submerged specific gravity of the sediment, assumed 1.65, and g is the standard gravity acceleration. The previous equations for stable channel design are shown in Table 1.

In the analytical approach, the design of the stable channel is controlled by the dominant process within the channel, such as flow resistance or sediment transport (Garcia 2008). The rational method proposed by Chang (1988), for instance, adopted three equations that are unique to the three dependent variables: width, depth and slope. Other methods derived from the analytical approach include those of Abou-Seida & Saleh (1987), Copeland (1994) and Yalin & da Silva (2000). The analytical method is very flexible, since the program enables the user to predict the stable channel geometry with the restriction in the depth or width of the river (Jang et al. 2016; Harun et al. 2020). The drawback of this method is its dependence on the accuracy of the sediment transport equation. Singh (2003) and Deshpande & Kumar (2012), however, use the concept of maximum entropy/minimum energy principle to examine the regime relationship for hydraulic geometry prediction. Gholami et al. (2018, 2019a, 2019b) also applied these concepts in predicting bank profile shape, and the results show that the method is capable of predicting geometry with better results. Bonakdari et al. (2020) asserted that rational approaches seem to be complicated and require a deep understanding of fluid flow, which is often time-consuming. Bonakdari et al. (2020) suggested including different key geomorphological variables to increase the prediction accuracy.

Currently, the use of machine learning in water resources engineering has been the method of choice for many researchers (Muhammad et al. 2018; Kargar et al. 2019; Montes et al. 2020; Roushangar & Shahnazi 2020; Sharafati et al. 2020). This method was extensively investigated by Gholami et al. (2019a, 2019b) and Shaghaghi et al. (2018) for use in a stable channel. Machine learning can be applied by using various methods such as artificial neural networks, adaptive network-based fuzzy inference systems, support vector machines, gene expression programming (GEP) and evolutionary polynomial regression (EPR) (Giustolisi & Savic 2009; Ebtehaj et al. 2019; Khosravi & Javan 2019; Roushangar & Ghasempour 2019; Yahaya 2019; Asheghi et al. 2020; Najafzadeh & Oliveto 2020). The research done by Bonakdari et al. (2020) showed that GEP and EPR had improved the accuracy of the model prediction; however, the improvement is not remarkable compared to the original equation.

These earlier findings suggest that a new equation should be developed to increase the accuracy of stable channel prediction. To date, many researchers have focused on parameters Q, d₅₀S_o and Q* to determine stable channel geometry. This paper will focus on introducing a new equation for stable channel geometry by adopting sediment transport parameters that resemble the hydraulic characteristic of the river. The width/depth (B/y) parameter will be introduced to the equation for predicting stable channel geometry to solve the problem when there is a restriction in width or depth, as discussed earlier in the analytical approach. The proposed equation will be further enhanced by using the machine learning method to increase the model prediction accuracy.

This study uses ordinary multiple regression in linearised form after log transforming all the variables included in the equation. Since all the possible test cases are described in one power law equation, the fitting of all regressions can be determined using statistical analysis software, Statistical Product and Service Solutions (SPSS). The dependent variable used for this equation is B*, y*, S_o and B/y.

GEP has been developed as an extension to genetic programming (Koza 1994). This program involved several search techniques, such as decision trees and logical expression, that help to evolve computer programs. GEP is encoded in linear chromosomes that can be represented by an expression tree (ET). With the aim to analyse particular data through genetic modification, the population of an ET will adapt and discover the traits of the data. The decision tree and logical expression have to be initially carefully determined to ensure the program yields the best results (Ferreira 2001b). The chromosome of each individual is randomly selected and evaluated by the fitness function and selected by the genetic operator to reproduce with modification (mutation). The same process occurs in the new individual for whom the trends are repeated until the required accuracy is obtained (Ferreira 2001a, 2001b). Figure 1 summarises the steps taken to execute the GEP. The data are first randomised before continuing to the next step. Approximately, 80% of the data are used for training and the remaining data for testing purposes. The setting of the function and chromosome architecture includes chromosome length, gene number, head size, linking function and genetic operators. The root-mean-square error (RMSE) was used as a fitness function to fit the curve to the target value. The stopping criteria were set at 25,000 generations for all the runs of the various models. The summarised parameters used in the GEP model are shown in Table 2.

The river data used in this study were obtained from the Malaysian Department of Irrigation and Drainage. Three rivers were chosen for the study to represent different hydraulic characteristics of a large river (Muda River), a medium river (Langat River) and a small river (Kurau River). The locations of these rivers are shown in Figure 2.

The Muda River originates in the highland areas of Kedah, a northern state of Malaysia, which borders Thailand. It is the largest river in Kedah and is important in supplying water to the three states of Kedah, Perlis and Pulau Pinang (Sim et al. 2015). The Muda River Basin covers a drainage area of 4,210 km². In terms of length, the Muda River is approximately 180 km long, with a slope of 1/2,300 (or 0.00043). The channel width of the Muda River is typically around 10 m upstream, 100 m in midstream and widest at its estuary, averaging 300 m (DID 2009). In terms of depth, bathymetric surveys show that the shallowest point in the river is located 2.5 km upstream from the river mouth, resulting in difficulty for navigation during low tides (Julien et al. 2010).

The Langat River forms one of the four major river systems in the state of Selangor. The Langat River is a medium-sized river, approximately 180 km long with an average annual flow of 35 m³/s and a mean annual flood of 300 m³/s. The Langat River system flows from the north-eastern state of Selangor to Negeri Sembilan and the Federal Territory of Putrajaya, finally emptying into the Straits of Malacca. The Langat River Basin has a total catchment area of 2,396 km² (DID 2009).

The Kurau River Basin is a typically small river basin, draining an area of approximately 682 km². In terms of elevation, elevations at the river headwaters are moderately high, being 900–1,200 m. In terms of slopes, the upper 6.5 km of the river averaged 12.5%, while those lower down the valleys are significantly lower, in the order of 0.25–5%. The average velocity of the Kurau River ranges from 0.45 to 0.636 m/s with the highest sediment load being 0.878 kg/s (Saleh et al. 2018). A reservoir was constructed at the middle section of the river, approximately 65 km upstream. Two major systems, the Kurau River system and the Merah River system, are upstream of this reservoir, and both rivers drain into the reservoir (DID 2009).

River surveys, flow measurements and field data collection were conducted based on the Guidelines for Field Data Collection and Analysis of River Sediment by Ab. Ghani et al. (2003). The data collection includes flow discharge, bed and bank material, suspended load, bedload and water surface slope. In addition, bed elevation, water surface and thalweg measurements were also conducted at selected cross-sections. The range of the data for this study is shown in Table 3.

The new equations were formulated by using dimensionless hydraulic geometry relations as proposed by Kaless et al. (2014), Parker (1978), Parker (1979) and Pitlick & Cress (2002). Along with the dimensionless stable geometry equation in terms of B*, y* and S, a new dimensionless B/y equation was introduced to solve the restriction issue in either width or depth prediction. The equation also uses dimensionless full bank discharge Q* as one of the parameters in the development of the equation.

Researchers have divided the significant parameters that contribute to sediment transport into four-parameter classes, namely, mobility, transport, sediment and flow resistance (Ebtehaj & Bonakdari 2016; Harun et al. 2020). This concept has been widely used by researchers to define the sediment transport in the river (Sinnakaudan et al. 2006; Harun et al. 2020), in closed channels or in pipes (Ebtehaj et al. 2019; Danandeh Mehr & Safari 2020). These parameters were also found to be prevalent in the application of sewer sedimentation analysis (Ebtehaj & Bonakdari 2017; Kargar et al. 2019). The details of the parameters are shown in Table 4.

These parameters were combined to get the best relationship with regard to the dimensionless hydraulic geometry relations to width, depth, slope and width/depth. The parameters governing the stable channel geometry are presented in a dimensionless form, where ψ is the flow parameter, R is the hydraulic radius, U* is the shear velocity, S_o is the bed slope, S_s is the specific gravity of the sediment, is the fall velocity of the bed material, V is the cross-sectional velocity, g is the acceleration due to gravity, A is the area of cross-section, v is the kinematic viscosity of water, ϕ is the transport parameter, D_gr is the dimensionless grain size and is the friction factor.

The significance of each parameter in Table 4 was assessed to ascertain which parameter is influential to the stable channel geometry. The correlation analysis determines the significance of the parameter. The probability (p) value < α = 0.05 indicates convincing evidence that the parameter is influencing the stable channel geometry. The results from the correlation analysis show that the stable channel geometry (B/y, B*, y*, S) were significant to the parameters Q*, ψ, ⁠, ⁠, ⁠, C_v, ϕ, ⁠, ⁠, R/d₅₀, and ⁠. The selected best combination of characteristic parameters, together with the related equations from multiple linear analyses, were further assessed in the study of variance (ANOVA) approach. Figure 3 explains the different input combinations for the determination of dimensionless width/depth, width, depth and the slope of a stable channel.

The B/y parameter was evaluated with a different combination of inputs, as shown in Figure 3. This study shows that B/y was significant to the combination of parameters Q*, ψ, ⁠, ⁠, ⁠, D_gr, and R/d₅₀. In total, there were 37 combinations that were significant to the prediction of B/y. The best five combinations of the parameters are given in Table 5.

The R² value for model one is the highest among all the models (0.890) followed by model two (0.638), model three (0.598), model four (0.542) and model five (0.453). In contrast, MSE for model one yields the lowest MSE (45.549) followed by model two (157.900), model three (187.70), model four (208.298) and model five (253.270). As for the NSE, model one yields the highest value (0.889), followed by model two (0.614), model three (0.542), model four (0.492) and model five (0.382). A comparison between all models shows that model one yields the highest DR value which ranges between 0.8 and 1.2 (79.06%), and model five yields the lowest DR value (31.62%). The results of this study indicate that the B/y is best described by the combination parameters of Q*, and R/d_50. The above findings conclude that model one is superior in comparison to the other models.

Results from the correlation analysis revealed that B* was significant with the parameters Q*, ϕ, ψ,⁠, ⁠, D_gr, ⁠, and R/d₅₀. There were 51 combination models altogether for the B* prediction. Table 6 lists the best five combinations of the parameters that are significant to the B prediction.

Model one has the highest DR (93.16%), followed by model two (64.10%), model four (55.56%), model three (46.15%) and model five (37.61%). In terms of R², model one has the highest value (0.964), and model five has the lowest value (0.745). The trend was also the same for NSE, with model one showing the highest value (0.964), thereby indicating that the model predicts better compared to the other models. Moreover, model one tends to yield the lowest MAE (3,031.499); meanwhile, model five has the highest MAE among the models (10,167.036). Model one is considered superior compared to the other models because it has the highest R² and NSE and is able to give better DR that has 93.16% accuracy.

Significant parameters were combined to yield the best equation that can predict the average depth of the river. In the present study, it is observed that parameter Q*, C_v, ϕ, ψ,⁠, D_gr, and R/d₅₀ were significant to the development of depth prediction. The analysis showed that 51 combination parameters were significant to the y* prediction. Table 7 lists the best five combinations of parameters that are significant to the depth prediction.

Models one and two have the highest value of DR square (92.31%) followed by model three (91.88%), model four (91.45%) and model five (91.03%). Among all the models, Model two has the lowest in terms of MSE (58,451.679) and the highest in terms of R² (0.977) and NSE (0.977). Model two is selected to predict the depth instead of model one because model two gives a higher R² and NSE and the lowest MSE compared to the other models.

The prediction of the slope was calculated by combining the significant parameter, as discussed in Figure 3. In this study, parameters C_v, ϕ, ψ,⁠, D_gr and R/d₅₀ were observed significant to the slope prediction. The study revealed that 57 models were significant to the slope prediction. The best five combinations for the slope prediction are listed in Table 8.

Model one has the highest R² (1) followed by model three (0.807), model two (0.795), model four (0.789) and model five (0.763). The same trend also applies to the NSE, but for MSE, the trend is in the reverse order. Model one has the lowest MSE, and model five has the highest MSE. Model one is found to be superior compared to the other models, as it has the highest DR ratio (100%) followed by model two (55.56%), model three (54.70%), model four (53.85%) and model five (50.43%).

The summary of the stable channel geometry prediction model in dimensionless form by using the MLR method is shown in Table 9.

The best model for every stable channel parameter was further determined by using GEP. The function for the stable hydraulic geometry was analysed based on the significant parameter, as shown in Table 9. The formulations of GEP for stable hydraulic geometry are shown in Table 10.

The performance of the training and testing of the GEP model is shown Table 11. The GEP model for all stable channel geometries has a high correlation for training and testing data. The same trend was also observed for the errors in terms of MSE and MAE.

The performance of the GEP model in comparison to the MLR model is shown in Table 12. Both models were analysed based on R², NSE, MSE and MAE. In all parameter predictions, GEP models outperformed MLR models, in which the GEP model yields better R² and NSE and produces the least errors in terms of MSE and MAE. It can be concluded that the GEP model is a far superior model compared to the traditional MLR method.

The performance of the currently developed equation by using MLR and GEP was compared to the previous equation (Figures 4–6), and the results are shown in Tables 13 and 14. The channel width, depth and slope of the stable channel are estimated by using the existing equation. As the equations were presented in dimensionless form, calculations of the width and the depth were done by multiplying the equation with d_50. For the B/y equation, the width is calculated by multiplying the equation by the depth of the channel. However, the equation will be divided by the width to obtain the depth of the channel.

In the prediction of channel width (B), the majority of the previous equations had moderate R² values but with low or negative values in terms of NSE. This shows that several of the equations were underpredicted, which explains the higher observed values compared to the predicted data. Only Pitlick & Cress (2002), Lee & Julien (2007) and Haron et al. (2019) have a positive value for NSE. Lee & Julien (2007) has the highest NSE of all, which is 0.61. Bonakdari et al. (2020) have the highest values of R² in EPR and GEP, namely, 0.685 and 0.655; however, these lack accuracy in the prediction, as the NSE value transpired to be negative. The NSE for both equations is −0.06 and −0.54, respectively. In terms of the error, Julien & Wargadalam (1995) have the highest MSE and MAE values, which are 1,335.68 and 24.76, respectively. Lee & Julien (2007) have the lowest MSE and MAE of all the equations (MSE = 206.52, MAE = 9.64). With the introduction of the new parameter, all newly developed equations were determined to give better stable channel width prediction. By using the MLR method, the current study MLR (B/y) and current study MLR (B) equations give higher R² and NSE values, and the revised equation, using GEP for both equations, returns improved R² and NSE values compared to the MLR method. The current study GEP (B) equation has the highest R² and NSE values, which are 0.965 and 0.965, followed by the current study MLR (B) equation (R² = 0.953, NSE = 0.953), the current study GEP (B/y) equation (R² = 0.863, NSE = 0.862) and the current study MLR (B/y) equation (R² = 0.835, NSE = 0.832). The trend for error indices is also the equivalent but in the reverse order. The current study GEP (B) equation has the lowest MSE and MAE values, which are 18.472 and 2.734, followed by the current study MLR (B) equation (MSE = 24.773, MAE = 2.674), the current study GEP (B/y) equation (MSE = 72.651, MAE = 4.811) and the current study MLR (B/y) equation (MSE = 88.308, MAE = 5.085).

In the prediction of flow mean depth y, many of the existing equations have high R² and NSE values, indicating that most of the equations are able to predict the depth of the river with low error (MSE and MAE). Hey & Thorne (1986), Lee & Julien (2007) and Bonakdari et al. (2020) EPR were particularly excellent in predicting the depth of the river with the DR (0.8–1.2) 39.32, 37.6 and 52.99%, respectively. Hey & Thorne (1986) have the lowest MSE (0.122) and MAE (0.282) among all the existing equations. Conversely, Simons & Alberstone (1960) have the highest MSE (40.49) and MAE (6.311). The newly developed equation is able to predict with better accuracy, coupled with less error. The current study GEP (y) equation has the highest R² and NSE values (R² = 0.960, NSE = 0.960) followed by the current study MLR (y) equation (R² = 0.955, NSE = 0.952), the current study GEP (B/y) equation (R² = 0.839, NSE = 0.819) and the current study MLR (B/y) equation (R² = 0.839, NSE = 0.783). The MSE and MAE of the newly developed equation were reduced to the ranges of 0.027–0.148 (MSE) and 0.095–0.195 (MAE).

In the slope prediction, many existing equations are unable to predict the slope accurately, and many of the existing equations had low R² and NSE values. Julien & Wargadalam (1995) and Lee & Julien (2007) equations were found to be excellent in predicting the slope with an extremely high R² value (1, 0.982) and an high NSE value, which are 1.0 and 0.861, respectively. The MSE and MAE for both equations are also fairly low, being 3.131 × 10⁻¹¹ and 4.985 × 10⁻⁸ (MSE), and 3.8467 × 10⁻⁶ and 2.054 10⁻⁴ (MAE). In all the equations, the EPR and GEP equations of Bonakdari et al. (2020) were found to have the highest error. This may be due to the sensitivity of the equations used and is, therefore, not suitable for studying Malaysian rivers. The newly developed slope equation using GEP and MLR is able to predict with high R² and NSE values (R² = 1, NSE = 1). The MSE for both equations is fairly low, namely, 1.773 × 10⁻¹² and 6.6217 × 10⁻¹², respectively, and MAE is 1.092 × 10⁻⁶ and 5.766 × 10⁻⁶, respectively. Both Julien & Wargadalam (1995) and the current study equation (MLR and GEP) with fewer variables are determined to be suitable to predict the stable channel slope.

The presently developed equations were validated by using the data from Saleh et al. (2017, 2018), Chang et al. (2008), Sinnakaudan et al. (2006) and Ariffin (2004).

The results from Table 15 show that both MLR and GEP equations can predict the width of the river precisely. The R² and NSE values are high for all the validation data except for the current study GEP (B/y) equation and for the Ariffin (2004) data, in which the values for R² and NSE are 0.880 and 0.543, respectively. The predicted data are highly correlated but lack accuracy. This could be due to the fact that the GEP equation developed is sensitive to the characteristics of the current river data. The data from the Ariffin (2004) study are generally based on small river data. Conversely, the present study MLR (B/y) equation is able to give far better results with better R² (0.954) and NSE (0.706) values. Overall, the new equation developed in this study is better in predicting the width compared to the existing equation.

Validation for the depth prediction in Table 16 shows that the developed MLR and GEP equations have an excellent degree of depth prediction accuracy based on the evidence of high R² and NSE values with low MSEs and MAEs. The current equation also outperformed the existing equations, as many of the existing equations yielded underpredicted values. However, the current study MLR (B/y) and current study GEP (B/y) equations yield low accuracy on the Saleh et al. (2017) data, whereby low R² and negative NSE values were observed. The model prediction is poorly correlated and inaccurate, as the predicted model performed worst. This could be due to the characteristics of the river that has a higher discharge compared to the current study data. This indicates that the equation is not suitable to apply to rivers with high discharge.

The slope validation prediction as shown in Table 17 revealed an equation from this study that is also able to predict the slope with high accuracy, particularly current study S (GEP), which yields better R² and NSE values and lower errors. This study also applies to the much-improved parameter that only consists of two parameters, namely, flow parameter and R/d₅₀.

As shown in Table 18, the new equation B/y also performed well in the validation data. The current study GEP (B/y) and and current study MLR equations have high R² and NSE values in all the validation data with low MSEs and MAEs. This indicates that parameter B/y is suitable to represent the geometry of the river. In the case of restriction either in the width or depth of the river, this equation is fundamental in helping decision-makers to determine the suitable channel dimension.

The current study focuses on developing a new equation for stable channel geometry that is suitable for the median size of bed material between 0.29 and 3.0 mm. Different key geomorphological variables such as bank profile and vegetation types need to be considered in future studies to enhance researchers' understanding of the effect of changes to the stable channel geometry. Applications of different machine learning techniques with higher accuracy can be conducted to increase the accuracy of the model prediction. For the B/y model prediction, a different variable that resembles the hydraulic characteristic of the river should be further explored to increase the efficiency of the equation, especially in rivers with a high volume of discharge.

The current study set out to develop a new equation that can improve the stable channel geometry prediction accuracy. This study also aims to establish an equation in the dimensionless form of B/y to provide an alternative solution for the width and depth prediction. Existing equations predict stable channel geometry with low accuracy, as many model predictions use the same parameter. This research has enhanced our knowledge in predicting stable channel geometry by associating the sediment transport parameters to predict stable channel geometry. The newly proposed equations by using MLR and GEP give a better prediction for stable channel hydraulic geometry with a far better R², NSE and DR. Moreover, it was established that B/y is a good parameter that can be used to predict the geometry of the river. This result also suggests that the proposed MLR and GEP models are robust in predicting stable channel geometry, with lower MSEs and MAEs compared to the existing equations.

The authors acknowledge the technical support provided by REDAC, USM, for the completion of this study. The first author also acknowledge the Public Service Department of Malaysia for the scholarship provided under the Hadiah Latihan Persekutuan (HLP) program for the duration of this study.

Data cannot be made publicly available; readers should contact the corresponding author for details.