Modeling shear stress distribution in natural small streams by soft computing methods

1. Introduction

Stream flow is still a study area today for researchers who want to explore the complex characters of stream flows. The shear stress distribution is the one of these properties to be investigated. The particles of the fluid have different ve- locities. The velocity of the particle will vary from layer to layer as its distance varies from the boundary in open channel flow. The shear stress is the force per unit area in the flow direction. Shear stress in streams cannot be determined di- rectly but is predicted using measurements of flow velocity or flow properties and their relation with the shear stress. The transverse distribution of shear stress in streams is known to be affected by the geometry of the cross section, the bound- ary roughness distribution and the structure of secondary flows (Ghosh and Roy, 1970;Knight and Patel, 1985;Yang and Lim, 2005).

Determination of shear stress distribution in streams is important not only for the design of water resources but also for the definition of structure of the turbulent flow in streams. Although many researchers have already investigated for solution of the problem, the shear stress distribution in streams cannot be still fully understood. The importance of determination of transverse shear stress distribution is stressed by the use of local or average shear stress in many hy- draulic equations concerning the computation of flow resistance, sidewall correc- tion, sediment transport rate, channel erosion or deposition, and design of chan- nels. Currently, the mean bed and side wall shear stresses are used to designing and modeling river due to the shortage of knowledge of the transverse distribu- tion of shear stress. These approaches are not an economical and reliable (Hicks et al., 1990).

After the measurements on secondary flows have been started, investiga- tions on the shear stress distribution accelerated. An analytical method was im- proved to determine the bottom shear stress in a shallow, symmetrical channel byLundgren and Johnson (1960). Shear stress distributions in open channel flows have been investigated by a number of authors (Tracy, 1965;Knight et al., 1994;Zheng and Jin, 1998;Ardiclioglu et al., 2006).Knight et al. (1984) improved an empirical equation which gives the percentage of the total shear force as given in the following form

where B is the breadth of the channel and H is the depth of flow.

Then,Knight and Patel (1985) performed a series of experiments for rectan- gular duct flume and presented Eq. (2)

where b is the width of duct and h is the height of duct.

Determination of the local or average shear stress as sensitive is a quite dif- ficult study even using advanced turbulence models. Empirical, analytical or simplified approaches were presented to calculate shear stress distribution by many investigators such asChristensen and Fredsoe (1998),Berlamont et al. (2003). The distribution of shear stress in a V-shaped channel was investigated byMohammadi and Knight (2004). The distribution of boundary shear stress along the wetted perimeter in two dimensional flows was obtained using flow-net byYu and Tan (2007). Distribution of boundary shear stress and the average bed and wall shear stresses in prismatic open channel flows were calculated using six different methods and the results were compared byKhodashenas et al. (2008).Cacqueray et al. (2009) investigated the shear stress in smooth rectangular chan- nels using detailed Computational Fluid Dynamics (CFD) simulation results. They showed that their CFD simulations are consistent. All of these investiga- tions on secondary flows and shear stress distributions performed by researchers have showed that it is necessary and significant a simple, practical model to un- derstand the shear stress distribution in streams.

Ardiclioglu et al. (2012) stressed that one dimensional hydraulic equations was used to explain the flow in open channels. Because of river hydrodynamics is quite complicated, these hydraulic equations are insufficient to determine the flow properties. Recently, artificial neural networks (ANNs) and adaptive neuro- fuzzy inference system (ANFIS) techniques have been successfully used to solve problems of the water resources and hydraulic engineering.Yang and Chang (2005) used ANN for simulation of velocity profiles, velocity contours and dis- charges.Kocabas and Ülker (2006) modeled critical submergence for an intake in a stratified fluid media by using ANFIS.Dogan et al. (2007) estimated sediment concentration obtained by experimental study by ANN method.Cobaner et al. (2008) used ANN for modeling bridge backwater. Mamak et al. (2009) used both ANN and ANFIS models for analyzing bridge afflux through arched bridge con- strictions.Riahi-Madvar et al. (2009) predicted longitudinal dispersion coefficient in natural streams by using ANFIS.Kocabas et al. (2009) predicted critical sub- mergence by using ANN.

Cobaner et al. (2010) developed an ANN model to estimate the boundary shear stress distributions in smooth rectangular channels and ducts. They used 94 experimental data obtained from the studies of well-known researchers in this field. According to result of this study, the ANN model performed slightly better than the classical methods. Bilhan et al. (2010) used two different ANN techniques for modeling lateral outflow over rectangular side weirs.Guven (2011) modeled scour geometry downstream from hydraulic structures by using ANN.Dursun et al. (2012) estimated discharge coefficient of semi-elliptical side weir using ANFIS.Can et al. (2012) modelled the daily streamflow using autore- gressive moving average (ARMA) and ANNs in Çoruh Basin in Turkey. They selected nine gauging stations that are operated by Turkish General Directorate of Electrical Power Resources and Development Administration (EIE). Their study revealed that the historical time series have similar statistical parameters to those of the generated time series at 95% confidence level.Kisi et al. (2013) used ANFIS approach and successfully modeled the discharge capacity of rectan- gular side weirs. Marsili-Libelli et al. (2013) proposed a new stream flow assess- ment method based on fuzzy habitat suitability and large scale river modeling.Wolfs and Willems (2013) used ANFIS for computationally efficient lumped floodplain modeling.Rezaeianzadeh et al. (2013) investigated flood flow forecast- ing using ANN, ANFIS and regression models. They showed that nonlinear re- gression can be applied as a quick method for predicting the maximum daily flow.Tayfur et al. (2014) estimated hydraulic conductivity in heterogeneous aquifers by using ANN and ANFIS methods. Baghari et al. (2014) used ANN method in order to identify the most important parameters affecting the discharge coefficient of rectangular sharp-crested side weirs.

In this paper, the applicability of the ANN and ANFIS approaches for model- ing shear stress distributions on natural streams is investigated. The results are compared with the multiple-linear regression model.

2. Field measurements

Recent studies of flow measurements on Kızılırmak and Seyhan Basins in Central Turkey have broaden understanding of flow properties in the Kızılırmak and Zamantı River. Twentytwo (22) measurement studies were performed at four different stations on Kızılırmak and Seyhan Basins by a team that consists of first and third authors. The names of measured stations are Barsama, Bünyan, Şahsenem and Sosun. Sosun Station where is in Seyhan basin and tributary of Zamantı River was visited between 2009 and 2010. The Zamantı River which is one of the branches of the Seyhan River originates in Pınarbaşı, Kayseri and drains into the Mediterranean Sea in southern Turkey. Measurement studies on Barsama, Şahsenem and Bünyan Stations where are in Kızılırmak Basin and at tributary of Kızılırmak River was carried out between 2005 and 2010. Kızılırmak basin is the second biggest basin in Turkey. Kızılırmak is the longest river that originates in Sivas and drains into the Black Sea in northern Turkey (Fig. 1). InFig. 2, measurements at Barsama and Şahsenem stations are given.

Figure 1. Location of the study area and measurement stations at Bünyan, Şahsenem, Barsama and Sosun(Ardiclioglu et al., 2012).

Figure 2. Velocity and topographical measurements at (a) Barsama and (b) Şahsenem Stations(Ardiclioglu et al. 2012).

The technology of flow measurement devices has evolved rapidly in recent de- cades. Many flow measuring devices were discovered to determine the flow proper- ties. The Acoustic Doppler velocimeter, ADV which is one of these devices was uti- lized to measure three dimensional velocity data at the measured stations. The ADV is used to measure water velocity in a wide range of environments including laboratories, rivers, estuaries, and the ocean. An ADV operates by the principle of Doppler shift. All ADV measures the velocity samples by sending out short acous- tic waves from the transmitter probe. These waves travel to moving particulates in the liquid and three receiving probes record for the change in frequency of the re- turned waves. The velocity vector is transmitted to recorder of ADV. Three- dimensional flow velocities (u, v, w) for x, y, z directions can be measured by ADV in a sampling volume. Furthermore, ADV can be used to record one second velocity data for the specified averaging time, location, and water depth parameters (to document the data set), and a variety of statistical and quality control data. The ADV can measure sampling volume from 10 cm front of the probe head. Thus, the probe head itself does not much impact on the flow field surrounding the measure- ment volume. Velocity range is between ±0.001 m/s and 4.5 m/s, resolution 0.0001 m/s, accuracy ±1% of measured velocity (Sontek, 2002).

According to the water surface width (T), the stream is divided into segments for each field measurements. As seen inTab. 1, wide channels (T / R ≥ 10) have 7–10 slices and narrow channels (T / R < 10) have 4–7 slices. Point velocity is mea- sured by ADV and cross-section area is measured using depth data measure- ments of distances from a fixed reference point on the riverbank. Point velocities were measured in the vertical direction starting 4 cm from the streambed for each vertical. Point velocity values were measured in the vertical direction start- ing 4 cm from the streambed for each vertical. Measurements were repeated ev- ery 2 cm from this point to water surface. The velocities of free water surface in all verticals were estimated using extrapolating the last two measurements of verticals. And also mean water surface velocities (u_ws) were measured at each visited stations. Water surface velocities were measured by how many seconds a tree branch passed distance of 10 meters using a chronometer. These measure- ments were repeated 10–15 times. In this way, we determined the average u_ws for all flow conditions and stations and given in the column 4 ofTab. 1.

The flow characteristics at each site are summarized inTab. 1. As shown inTab. 1, six measurements have been made at Barsama, Şahsenem and Bünyan stations in 2005–2006. Also, four measurement studies were undertaken at Sosun Station in 2009–2010.

The flow characteristics at each site are given inTab. 1. In this table, first and second columns show visit numbers and dates of stations, U_m (= Q / A) is the mean velocity, u_ws is the measured water surface velocity, with A being the area of the cross section, H_max is the maximum flow depth, and T / R is the aspect ratio, with T being the surface water width, R (= A / P) is the hydraulic radius, P is wetted perimeter. Re (= 4 U_m R / ν) is the Reynolds number, and ν is the kinematic 1/2 viscosity, Fr (= U_m / (g H_max) ) is the Froude number, where g is the gravitational acceleration, and S is water surface slope. According to the Froude and Reynolds numbers all the flow measurements are under subcritical and turbulent flow conditions according to Froude and Reynolds numbers. Froude numbers are be- tween 0.084–0.578, aspect ratios (T / R) vary between values 6.53–45.40. Reynolds numbers (Re × 10⁶) vary between 0.32–1.47.

Figure 3. The frequency histograms of the target variable, shear stress τ0.

Shear stress, τ (N/m²)

Figure 3 presents the scatter plots for frequency histogram of 145 shear stress, τ₀ values in measured vertical. As seen from the figure, the most of shear stress, τ₀ values change between the values 0–20 N/m². The minimum shear stress is 0.081 N/m² and the maximum shear stress is 78.4 N/m².

3. Shear stress distributions

Shear stress distribution basically depends upon the geometry of the cross section and the mechanism of the secondary flow cells. Shear force in uniform flow can be defined as in below

where γ is the specific weight or unit weight of water, A is the channel cross sectional area, L is the length of the control volume and S is the longitudinal slope of the channel. Unit shear force is calculated using Eq. (3) as following form.

where, τ₀ is the average value of the shear force per unit of the wetted area, R is the hydraulic radius (= A / P in which A is the wetted area and P is wetted perimeter).

It is clear that, shear force is not always uniformly distributed over the pe- rimeter. French (1985) stressed that for a sensitive calculation methodology to be determined, the distribution of the shear force on the perimeter of the channel should be predicted.

Schlicting (1987) informed that logarithmic relation between the shear velocity and the variation of velocity with height is used to determine the local bed shear stress.

where u is stream wise velocity at z, u_* = (τ₀/ρ)^1/2 is the shear velocity, ρ is water density, k_s is the Nikuradse’s original uniform sand grain roughness, χ is the von Kármán constant and z is distance from the bottom of the roughness elements.

4. Method

4.1. Artificial neural networks

ANNs are developed by getting inspired from real nervous system by neglect- ing most of the biological details. They are massively parallel systems consisting many processing elements, called neurons. Each layer in ANNs is fully connected to the proceeding layer by interconnection weights.Figure 4 represents a three layered ANN comprising layers i, j, and k, with the weights W_ij and W_jk between layers. Initial weight values are randomly assigned and then corrected during a training process. This process compares model outputs with measured outputs and back propagates any errors (from right to left inFig. 4). Thus, the final weights are obtained by minimizing the errors(Kisi, 2005).

In layers j and k, each neuron receives the x input which is the weighted sum of outputs from the previous layer. As an example, in layer j, y can be given as

in which θ_j is a bias for neuron j, O_pi is ith output of the previous layer, W_ij is the weights between the layers i and j. An output f(y) is obtained from each neuron in layers j and k by passing y value through a non-linear activation function. Logistic function is commonly used as an activation function

Figure 4. A three-layered ANN structure (Kisi, 2005).

The detailed theoretical information about ANNs can be obtained from the related references(Haykin, 2009).

4.2. Adaptive neuro-fuzzy inference system

ANFIS is first pioneered byJang (1993). As a universal approximator, it is capable of approximating any real continuous function. The structure of ANFIS is illustrated inFig. 5. As also seen from the figure that its structure is composed of a number of nodes connected through directional links and each node has a function consists of fixed or adjustable parameters(Jang, 1997).

Assume a fuzzy inference system having three x, y and z inputs and one f output and its rule base consist of two Takagi and Sugeno fuzzy IF-THEN rules

Rule 1: IF x is A₁, y is B₁ and z is C₁ THEN

Rule 2: IF x is A₂, y is B₂ and z is C₂ THEN

in which f₁ and f₂ respectively refer to the output function of rule 1 and rule 2.

Figure 5. The structure of an ANFIS.

Every node i in layer 1 includes an adaptive node function

where x is the ith node’s input and A_i is a linguistic label such as “small” or “big” associated with this node function. O_l,i refers the membership function of a fuzzy set A (in this example; A₁, A₂, B₁, B₂, C₁, or C₂). It denotes the degree to which the given input x satisfies the A_i. Gaussian function is usually chosen for the φA_i(x)

where a_i and b_i are the function’s parameters. When the values of these parameters change, the Gaussian function also varies accordingly, thus various forms of membership functions are exhibited on linguistic label A_i(Jang, 1993). The parameters in this layer are named as premise parameters.

In layer 2, the incoming signals are multiplied in each node. For instance,

The output of each node indicates the firing strength of a rule.

In layer 3, the ratio of the ith rule’s firing strength to the sum of all rules’ firing strengths is obtained by ith node as

In layer 4, each node has a function given in Eq. 14

in which $\stackrel{\_}{w_i}$ is the output of layer 3, and p_i, q_i, r_i and t_i are the consequent parameters.

In layer 5, the final output is computed as the summation of all incoming signals

Thus, the ANFIS network which is functionally equivalent to a first-order Sugeno fuzzy inference system has been constructed. Detailed information for ANFIS, the readers are referred to the related references (Jang, 1993).

5. Results and discussion

In the current study, two different program codes including MATLAB neural network and fuzzy toolboxes were prepared for simulation of ANN and ANFIS models. Water surface velocity u_ws, water surface slope S and y / T were used as inputs to the models to estimate shear stress distribution. One hundred forty five (145) field measured data were used in development of the ANN, ANFIS and MLR models. After randomly permutation, whole data were divided into training and testing parts. The first data set (116 values, 80% of the whole data) was used for training and the obtained models were tested using the second data set (29 values, 20% of the whole data). Before applying ANN to the data set, the training input and output values were normalized between 0.2 and 0.8 by using the following equation

in which x_max and x_min refer to maximum and minimum of the training/test data. Here, the c₁ and c₂ values were respectively assigned as 0.6 and 0.2. The optimal ANN models were obtained after trying various model structures. The evaluating criteria used in the study are root mean square errors (RMSE), mean abso- lute errors (MAE) and determination coefficient (R²). The equations of the RMSE, MAE and R² are:

where N and y_i refer to number of data set and entropy parameter, respectively. MSE is mean square errors.

For estimating shear stress distribution of the streams, three different input combinations were used. The Pearson correlations between the inputs u_ws, S, y / T and output are 0.804, 0.772 and 0.031, respectively. The lowest linear relation- ship exists between the y / T and output and the u_ws has the highest correlation. According to the correlation values, u_ws seems to be the most effective variable on shear stress distribution. The optimal hidden node numbers were obtained for each ANN model by using simple trial and error method. The hidden node num- ber tried for each ANN model ranges between 1 and 10. As an example, the vari- ation of hidden node number vs. mean square error (MSE) in test stage for the ANN model with three inputs is illustrated inFig. 6. The ANN models were trained by using Conjugate Gradient algorithm which is more powerful than the classical gradient descent technique. The sigmoid activation functions were used for the hidden and output nodes. The ANN training was stopped after 1000 iterations. Training and test results of the optimal ANN models are shown inTab. 2 in respect of RMSE, MAE and R2 statistics. The optimal hidden node numbers are also provided in the second column of this table. It is apparent from the table that the ANN model comprising three inputs corresponding to u_ws, S_ws and y / T, 6 hidden and 1 output nodes has the lowest RMSE, MAE and the highest R² both in training and test periods. FromTab. 2, it is clear that the accuracy of ANN grid partition method was used by applying genfis1 in MATLAB.Table 2 represents training and test results of the ANFIS model for each input combination. The second column of this table gives the optimal number of membership func- tions. It is clearly seen fromTab. 2 that the ANFIS model comprising 3 Gaussian membership functions for the inputs, u_ws, S and y/T outperforms the other ANFIS models for the both periods. Comparison of ANN and ANFIS models clearly indicates that the ANFIS models generally provide better shear stress estimates than the ANN models. The RMSE, MAE and R2 statistics of the MLR models are also shown inTab. 2. Second column of this table gives the regression coefficients of each model. The MLR model with the inputs of u_ws, S and y / T gives inferior results in relative to the ANN and ANFIS models. This implies the strong non- linear relationship between u_ws, S and y / T inputs and shear stress distribution. The optimal ANFIS model respectively reduced the root mean square errors and mean absolute errors by 47% and 50% with respect to the optimal MLR model. t-statistics of each coefficient of the MLR model comprising three inputs calcu- lated using Excel program are reported inTab. 3. According to the table, all three coefficients seem to be significant at the 5% significance level.

Figure 6. Variation of hidden node number vs MSE in test stage for the ANN model with three inputs.

The shear stress estimates of the ANN, ANFIS and MLR models in test pe- riod are illustrated inFig. 7 for the test period. It is clear from the figure that the estimates obtained using ANFIS and ANN models are much closer to the corre- sponding measured shear stress values than those of the MLR model. Significantly under/over-estimations are clearly seen for the MLR.

Figure 8 demonstrates the test estimates of each model in the scatterplot form. It is clear from the figure that the shear stress estimates of the ANN and ANFIS models are closer to the corresponding measured values than those of the MLR model. The MLR having highly scattered estimates seems to be insufficient in estimating shear stress distribution of the natural streams. The ANFIS model performs better than the ANN model especially for the low and high shear stress values.

Figure 7. Time variation of the measured and estimated shear stress values by ANN, ANFIS and MLR models.

Figure 8. The scatterplots of the measured and estimated shear stress values by ANN, ANFIS and MLR models.

The results are tested by using one way ANOVA for verifying the robustness (the significance degree of differences between the model estimates and mea- sured shear stress distribution) of the models. The test is set at a 95% significant level.Table 4 reports the test statistics. It is clear from the table that the ANFIS and ANN models yield small testing values with high significance levels. FromTable 4, the ANFIS and ANN models seem to be more robust (the similarity between the measured entropy parameters and model estimates are significantly high) in shear stress estimation than the MLR model.

6. Conclusion

This study investigated the accuracy of ANFIS and ANN methods for esti- mating shear stress distribution in streams. Overall 145 field data gauged from four different cross-sections at four sites on the Sarımsaklı and Sosun streams in Turkey were used. The water surface velocity u_ws, water surface slope S and y / T were used as inputs to the ANFIS and ANN models for estimating shear stress distribution. The estimates obtained from the ANFIS and ANN models were compared with multiple-linear regression model. The results revealed that the ANFIS and ANN models performed much better than the MLR model in esti- mating shear stress distribution. The ANFIS was found to be slightly better than the ANN. The best ANFIS model was respectively reduced the root mean square errors and mean absolute errors by 47% and 50% with respect to the best MLR model. The study recommends that the ANFIS and ANN techniques can be successfully used for estimating shear stress distribution of the natural streams.