Modelling boundary shear stress distribution in open channels using a face recognition technique

Boundary shear stress is the result of the tangential component of the hydraulic forces that act in the direction parallel to the channel's boundaries and transfer momentum to its bed and walls (Chow 1959). Excessive shear stress can undermine channel stability by eroding bank sides and cause changes in the river morphology by affecting the transport and deposition of sediments (Julien 1995). Erosion often results in higher levels of turbidity and lower water quality levels. Furthermore, an increase in sediment movement and deposition can cause a decrease in channel capacity and, consequently, higher flood risk. Computation of flow resistance, side-wall correction, sediment discharge, channel erosion or deposition, cavitation problems, and design of stable channels are among the problems which require accurate estimates of the boundary shear stress distribution (Yang & Lim 1997; Guo & Julien 2005; Blanckaert et al. 2010).

The distribution of boundary shear stress over the wetted perimeter of a channel cross-section is non-uniform. This is true even for steady flows in straight prismatic channels with a simple cross-sectional geometry. This non-uniformity is mainly due to the anisotropy of the turbulence which produces transverse gradients of Reynolds stresses and secondary circulations (Gessner 1973). Tominaga et al. (1989) and Knight & Demetriou (1983) showed that the boundary shear stress increases where the secondary currents flow towards the wall, and decreases when they flow away from the wall. Other factors that govern the distribution of shear stress are the geometry of the cross-section, lateral and longitudinal boundary roughness distributions (Blanckaert et al. 2010) and sediment concentration (Khodashenas et al. 2008).

To date, numerous investigations have been conducted and various mechanistic and empirical methods have been developed for understanding and estimating the magnitude and distribution of boundary shear stress. However, due to the complexities involved, boundary shear stress has proven to be one of the most challenging parameters to quantify and measure, even for simple smooth prismatic channels with uniform flow.

Preston's (1954) method is the most widely practised technique for measuring boundary shear stress in smooth channels. In this method, a Preston tube is used to infer the velocity of the water flow by recording the difference between static and total pressures. A non-dimensional calibration function is then established based on the ‘law of the wall’, Equation (3), and used to determine the boundary shear stress from the differential pressures. The simplicity of the experimental setup and its operation are the main reasons behind the popularity of this method. However, for rough boundaries, application of the technique is substantially more complicated, due to the absence of a viscous sublayer. A number of studies (Hwang & Laursen 1963; Ghosh & Roy 1970; Hollick 1976; Hollingshead & Rajaratnam 1980) have attempted to extend the use of this technique to rough surfaces, and have calibrated curves for the Preston tube by using Nikuradse's (1933) model of velocity distribution over rough boundaries. Although promising, these methods can only be applied when the sand equivalent roughness height of the surface is known, which makes them unsuitable for application to a variety of open channels. Other methods based on fitting the log law of the wall such as Clauser's (1956) method and the boundary characteristics method (Hinze 1975; Papanicolaou et al. 2012) have been developed and applied to gradually (Afzalimehr & Anctil 2000) and rapidly varying flows over spatially varying boundaries (Papanicolaou et al. 2012).

Geometrical methods for estimating shear stress distribution (Leighly 1932; Einstein 1942; Lundgren & Jonsson 1964; Yang & Lim 1997; Khodashenas & Paquier 1999; 2005; Yu & Tan 2007; Abderrezzak et al. 2008) consist of splitting the channel cross-section into sub-regions where the shear force along each segment of the boundary is calculated by balancing the forces against the weight of fluid in the corresponding subregion. In these approximations, mapping and discretising the wetted perimeter is often a complicated and sensitive task, however, they have the advantage of requiring relatively low computational effort.

Where abundant experimental data existed, researchers (e.g. Knight 1981; Knight et al. 1984, 1994; Flintham & Carling 1988; Pizzuto 1991; Olivero et al. 1999) have used regression and correlation analysis to derive empirical and semi-empirical equations for boundary shear stress. These equations are capable of only calculating mean, maximum and percentages of shear stress carried on the channel's walls and beds with relatively good accuracy, but are unable to provide the distribution of shear stress along the entire wetted perimeter. Some other researchers (e.g. Zheng & Jin 1998; Jin et al. 2004; Bilgil 2005; Guo & Julien 2005) have solved the governing energy transport, continuity, and momentum equations to formulate analytical and semi-analytical solutions for calculating the boundary shear stress. These methods often rely on a number of subjective and controversial assumptions and require a large amount of computing resources which make them impractical. With the advent of more powerful computers, computational fluid dynamic techniques have also been used (e.g. Christensen & Fredsoe 1998; De Cacqueray et al. 2009) to solve the referred set of equations and calculate the boundary shear stress distribution. Nonetheless, these methods are computationally expensive and the model outputs are extremely sensitive to mesh size, the turbulence closure model, and other internal parameters defined by the user.

Recently, information theory and machine-learning techniques have been used to tackle this problem. For instance, the principle of maximum entropy has been used (e.g. Sterling & Knight 2002; Li & Zhang 2008; Bonakdari & Moazaminia 2015) to establish relationships for the boundary shear stress. A comparison with experimental data has shown that these approximations provide relatively flat shear stress distributions which make them unreliable. The divergence between the numerical and experimental results increases at the regions around the corners of the sections where secondary flow structures are more pronounced. Cobaner et al. (2010) used a neural network with four hidden layers to predict the percentage of the shear force acting on the walls of smooth rectangular channels and ducts. The study concluded that the artificial neural network (ANN) predictions were less biased and slightly more accurate than the classic empirical models suggested by Knight et al. (1984) and Knight & Patel (1985).

Measuring the actual local shear stress along the channel's boundaries is difficult and costly owing to the complexity of the turbulent velocity field, presence of flow structures, and the small magnitude of the stress. Shear stress also represents a difficult parameter to calculate due to the variability of channel slope, geometry and flow structures, which are the main influencing factors in the complex flow process. To date, all the developed methods are inherently based on some sort of simplifying assumption, and therefore, the problem of accurately estimating these stresses has only been partially resolved (Zheng & Jin 1998).

The recent relative abundance in available experimental data has offered an opportunity to test and validate novel techniques for describing the boundary shear stress distribution. In this paper, an advanced pattern recognition technique is employed to predict the distribution of boundary shear stress in open channels. This technique, which results from merging two existing algorithms (recurrence plots (RPs) and Eigenfaces for Recognition), is combined with a standard regression model for the prediction of data series representing shear stress distribution of flows with known attributes (i.e. Froude numbers, flow depths and channel slopes).

In the following sections the RP analysis and its adaptation to the Eigenfaces for Recognition is explained. This is followed by a description of the experimental data used in the study and details of the proposed methodology. Next, the prediction of boundary shear stress distributions in trapezoidal and circular channels are presented and critically discussed. The paper concludes with a discussion on the advantages of the method and suggestions for improvement.

The proposed approach for predicting boundary shear stress distribution combines RP analysis (Eckmann et al. 1987) and Eigenfaces for Recognition (Turk & Pentland 1991). RP is used to transform one-dimensional data series into two-dimensional arrays which can be graphically represented. Eigenfaces for Recognition is then used as a means of identifying patterns in the arrays and to transform these into short-dimension vectors which can then be used to predict boundary shear stress distributions. It is to note that despite using a method that was originally developed for the recognition of human faces using 2D still images; no ‘recognition’ is involved in the proposed methodology. Instead, the technique is used to filter the original data and reduce their dimensionality whilst preserving intrinsic qualities. These reduced databases are then used to produce shear stress distribution for unseen cases.

By projecting the RP matrix on a Cartesian space, a map of the data can be generated. In this case each pixel on the map has the coordinates {j,k} with j, k = 1, 2, … , N, as well as a numerical value that is proportional to its associated distance, e_jk. In an 8-bit grayscale image representation, the values of e_jk can take values within the range 0 to 255, where the brightness intensity of each pixel indicates a larger e_jk. Through the calculation of the RP matrix, the correlation between all data points within the data series is established whilst preserving the basic structure of the database. Such transformation and visualization helps to make explicit features of data which otherwise would be difficult to observe in the original series.

As will be further explained in the following sections, a reliable estimation of the weight factors, i.e. vector Ω′, would be based on the known vectors obtained through Equation (15) together with the parameters that characterize the original data sets. These vectors and parameters can be typically related with the aid of a simple regression model or more complicated method such as an artificial neural network, and the output Ω′ can be considered to be reliable if the modelling error is less than 5%. Note that the validity of Equation (16) is provided in Appendix A (available with the online version of this paper).

In this study, laboratory measurements of flow velocity and boundary shear stress in trapezoidal and circular open channels were taken directly from the University of Birmingham's Flow Database (www.flowdata.bham.ac.uk).

To obtain homogenous subsets suitable for pattern recognition, and to test the sensitivity of the approach to the size of the training set, k-nearest neighbors (k-nn) clustering analysis (Fix & Hodges 1951) was first performed. The fundamental idea of the k-nn algorithm is to simply separate the data based on the assumed similarities between various clusters. Here, the Euclidean distance metric was used to measure the similarity between clusters, and shear stress data were non-dimensionalized by the average shear stress to eliminate the scale effects. K-nn was run with different k values, and consequently, three clusters (subsets) were identified by investigating the resultant dendrograms, i.e. graphical tree-structures that show the hierarchical relationships among clusters, ensuring highest similarity within each cluster (homogeneity) and lowest similarity between clusters were achieved.

Table 1 lists the geometric and hydraulic parameters of all the experiments. In each subset, one experiment (highlighted in Table 1) was randomly selected and excluded to be used for validation whilst the remaining were considered for training. Since the method requires all data series in the set to have the same number of measurements taken at relatively even distances, for each experiment, the horizontal coordinates of all data points were normalized using a perimetric distance defined as Pd=s/p, where s is the distance along the wetted perimeter starting at the left bank at the free surface, moving around the wetted perimeter, and p is the total length of the wetted perimeter. The shear stress measurements in each experiment were also non-dimensionalized by the average shear stress. Where required, linear interpolation was used to obtain shear stress values from adjacent neighboring points. It is notable that at regions where shear stress varied at higher rates, experimental measurements were taken at smaller increments, thus resulting in the standardization of the degree of accuracy across the wetted perimeter.

Table 1

Geometric and hydraulic parameters of trapezoidal experiments

.	.	1 .	2 .	3 .	4 .	5 .	6 .	7 .	8 .	9 .	10 .	11 .
.	ID .	2b (m) .	h (m) .	A (m2) .	S0 .	P (m) .	R .	V (ms−1) .	Q (m3 s−1) .	(Nm−2) .	Fr .	Re .
Set 1	#1	0.15	0.030	0.005	0.0040	0.235	0.023	0.565	0.003	0.894	1.12	11,392
	#2	0.15	0.058	0.012	0.0087	0.313	0.038	1.308	0.016	3.256	1.97	43,656
	#3	0.15	0.037	0.007	0.0234	0.255	0.027	1.843	0.013	6.223	3.35	43,800
(TVc-1)	#4	0.15	0.042	0.008	0.0234	0.269	0.030	1.901	0.015	6.871	3.27	50,251
	#5	0.15	0.050	0.010	0.0234	0.291	0.034	2.080	0.021	7.859	3.32	62,609
	#6	0.15	0.057	0.012	0.0234	0.310	0.038	2.190	0.026	8.625	3.32	72,546
Set 2	#7	0.45	0.044	0.022	0.0040	0.574	0.038	0.893	0.019	1.472	1.42	29,574
	#8	0.45	0.050	0.025	0.0010	0.591	0.042	0.398	0.010	0.414	0.60	14,743
	#9	0.45	0.056	0.029	0.0010	0.609	0.047	0.428	0.012	0.459	0.61	17,563
	#10	0.45	0.060	0.031	0.0010	0.620	0.049	0.439	0.013	0.484	0.60	18,998
	#11	0.15	0.029	0.005	0.0087	0.231	0.022	0.924	0.005	1.882	1.88	17,923
	#12	0.15	0.036	0.007	0.0087	0.250	0.026	1.010	0.007	2.244	1.87	23,347
	#13	0.15	0.041	0.008	0.0088	0.266	0.029	1.113	0.009	2.521	1.94	28,663
	#14	0.15	0.048	0.009	0.0087	0.284	0.033	1.231	0.012	2.815	2.01	35,702
(TVc-2)	#15	0.45	0.044	0.022	0.0087	0.574	0.038	1.381	0.030	3.228	2.19	45,751
	#16	0.45	0.059	0.030	0.0087	0.615	0.048	1.553	0.046	4.124	2.17	65,743
	#17	0.45	0.044	0.022	0.0145	0.574	0.038	1.822	0.040	5.384	2.89	60,371
	#18	0.15	0.029	0.005	0.0234	0.231	0.022	1.592	0.008	5.052	3.24	30,888
	#19	0.45	0.045	0.022	0.0234	0.576	0.038	2.272	0.050	8.752	3.59	76,145
	#20	0.45	0.059	0.030	0.0234	0.615	0.048	2.427	0.072	11.067	3.38	102,739
Set 3	#21	0.15	0.075	0.017	0.0040	0.362	0.047	0.960	0.016	1.812	1.29	39,299
	#22	0.15	0.107	0.028	0.0010	0.453	0.061	0.497	0.014	0.596	0.58	26,583
(TVc-3)	#23	0.15	0.125	0.034	0.0010	0.504	0.068	0.544	0.019	0.669	0.59	32,599
	#24	0.15	0.150	0.045	0.0010	0.574	0.078	0.584	0.026	0.768	0.59	40,170
	#25	0.45	0.075	0.039	0.0010	0.662	0.060	0.514	0.020	0.583	0.64	26,845
	#26	0.15	0.073	0.016	0.0087	0.356	0.046	1.468	0.024	3.896	2.00	58,886
	#27	0.15	0.099	0.025	0.0087	0.430	0.057	1.667	0.041	4.891	2.00	84,008
	#28	0.15	0.030	0.005	0.0145	0.235	0.023	1.296	0.007	3.272	2.58	26,146
	#29	0.15	0.075	0.017	0.0145	0.361	0.046	1.943	0.033	6.598	2.62	78,915
	#30	0.15	0.099	0.025	0.0234	0.430	0.057	2.690	0.066	13.129	3.23	135,513
	#31	0.45	0.044	0.004	0.0234	0.574	0.007	1.679	0.007	5.264	3.13	9,997
	#32	0.45	0.058	0.006	0.0234	0.614	0.010	1.785	0.011	6.513	2.96	15,045

.	.	1 .	2 .	3 .	4 .	5 .	6 .	7 .	8 .	9 .	10 .	11 .
.	ID .	2b (m) .	h (m) .	A (m2) .	S0 .	P (m) .	R .	V (ms−1) .	Q (m3 s−1) .	(Nm−2) .	Fr .	Re .
Set 1	#1	0.15	0.030	0.005	0.0040	0.235	0.023	0.565	0.003	0.894	1.12	11,392
	#2	0.15	0.058	0.012	0.0087	0.313	0.038	1.308	0.016	3.256	1.97	43,656
	#3	0.15	0.037	0.007	0.0234	0.255	0.027	1.843	0.013	6.223	3.35	43,800
(TVc-1)	#4	0.15	0.042	0.008	0.0234	0.269	0.030	1.901	0.015	6.871	3.27	50,251
	#5	0.15	0.050	0.010	0.0234	0.291	0.034	2.080	0.021	7.859	3.32	62,609
	#6	0.15	0.057	0.012	0.0234	0.310	0.038	2.190	0.026	8.625	3.32	72,546
Set 2	#7	0.45	0.044	0.022	0.0040	0.574	0.038	0.893	0.019	1.472	1.42	29,574
	#8	0.45	0.050	0.025	0.0010	0.591	0.042	0.398	0.010	0.414	0.60	14,743
	#9	0.45	0.056	0.029	0.0010	0.609	0.047	0.428	0.012	0.459	0.61	17,563
	#10	0.45	0.060	0.031	0.0010	0.620	0.049	0.439	0.013	0.484	0.60	18,998
	#11	0.15	0.029	0.005	0.0087	0.231	0.022	0.924	0.005	1.882	1.88	17,923
	#12	0.15	0.036	0.007	0.0087	0.250	0.026	1.010	0.007	2.244	1.87	23,347
	#13	0.15	0.041	0.008	0.0088	0.266	0.029	1.113	0.009	2.521	1.94	28,663
	#14	0.15	0.048	0.009	0.0087	0.284	0.033	1.231	0.012	2.815	2.01	35,702
(TVc-2)	#15	0.45	0.044	0.022	0.0087	0.574	0.038	1.381	0.030	3.228	2.19	45,751
	#16	0.45	0.059	0.030	0.0087	0.615	0.048	1.553	0.046	4.124	2.17	65,743
	#17	0.45	0.044	0.022	0.0145	0.574	0.038	1.822	0.040	5.384	2.89	60,371
	#18	0.15	0.029	0.005	0.0234	0.231	0.022	1.592	0.008	5.052	3.24	30,888
	#19	0.45	0.045	0.022	0.0234	0.576	0.038	2.272	0.050	8.752	3.59	76,145
	#20	0.45	0.059	0.030	0.0234	0.615	0.048	2.427	0.072	11.067	3.38	102,739
Set 3	#21	0.15	0.075	0.017	0.0040	0.362	0.047	0.960	0.016	1.812	1.29	39,299
	#22	0.15	0.107	0.028	0.0010	0.453	0.061	0.497	0.014	0.596	0.58	26,583
(TVc-3)	#23	0.15	0.125	0.034	0.0010	0.504	0.068	0.544	0.019	0.669	0.59	32,599
	#24	0.15	0.150	0.045	0.0010	0.574	0.078	0.584	0.026	0.768	0.59	40,170
	#25	0.45	0.075	0.039	0.0010	0.662	0.060	0.514	0.020	0.583	0.64	26,845
	#26	0.15	0.073	0.016	0.0087	0.356	0.046	1.468	0.024	3.896	2.00	58,886
	#27	0.15	0.099	0.025	0.0087	0.430	0.057	1.667	0.041	4.891	2.00	84,008
	#28	0.15	0.030	0.005	0.0145	0.235	0.023	1.296	0.007	3.272	2.58	26,146
	#29	0.15	0.075	0.017	0.0145	0.361	0.046	1.943	0.033	6.598	2.62	78,915
	#30	0.15	0.099	0.025	0.0234	0.430	0.057	2.690	0.066	13.129	3.23	135,513
	#31	0.45	0.044	0.004	0.0234	0.574	0.007	1.679	0.007	5.264	3.13	9,997
	#32	0.45	0.058	0.006	0.0234	0.614	0.010	1.785	0.011	6.513	2.96	15,045

View Large

The initial step in using RPs and Eigenface Recognition for predicting boundary shear stress distribution is forming a training set from available experimental data. As mentioned in the previous section, the raw experimental data sets typically consist of a number of local boundary sheer stress measurements, taken along the wetted perimeter of a channel, and presented as a data series. If measurements are not taken at the same relative locations across the different channels, then, to make the raw data suitable for use in the data mining algorithm, interpolation is carried out to find shear stresses at the same relative distances along the wetted perimeter, for all data series.

As the Eigenfaces for Recognition algorithm accepts two-dimensional arrays of numbers, the original (one-dimensional) data series in the training set has to be pre-processed. This can be done through the RP algorithm. To this end, the dimension of the vectors identified in Equation (4) is set to be equal to 1, so that each vector would contain numerical differences of shear stress between consecutive points in the data series, effectively resulting in each row of the RP matrix becoming a normalized version of the original data series. In line with the image recognition algorithm, once the RP matrix is formed, Equation (6) is used to construct a unique Г_i vector representing the i-th experimental data set.

Figure 7 shows the steps involved in the process of encoding the original data series. This process can be reversed to obtain a new set of weighting values, ω_ɩ, e.g. for test cases not included in the training set. Equation (17) enables finding those weighting factors which can then be stored in the short-dimension vector (Ω′). Following, the vector Γ′ can be predicted by applying Equation (16).

To help the reader better understand the entire modelling process, a simple step-by-step guide to using the proposed approach is presented in Appendix B (available with the online version of this paper). Furthermore, a copy of the code written in C ++ is available at: http://shear-stress-using-face-recog.sourceforge.net.

Finally, the predicted series were rescaled by multiplying their ordinates by the estimated average shear stress obtained by the Slope method (Equation (1)). For practical purposes, this parameter could be estimated using any other prediction model, such as the ones suggested by Knight (1981), Knight et al. (1984, 1994) and Flintham & Carling (1988).

RP analysis and Eigenface for Recognition were used to predict the distribution of boundary shear stress in trapezoidal and circular channels. In this approach, first, the RPs of all training set members are constructed and the differences between them and the average RPs are computed. PCA is then performed and weight factors proportional to the eigenvectors are obtained. To obtain predictions of boundary shear stress, a simple regression equation is established to relate the weight factors to non-dimensional attributes of the training set's hydraulic and geometric characteristics. For each validation case, corresponding weights are obtained by the regression equations, and the reverse of the process is performed to obtain the distribution of the boundary shear stress.

The method was applied to two trapezoidal data sets and one circular data set. The results showed that: The technique is capable of capturing the intrinsic patterns of the RPs which makes it suitable for the prediction of shear stress distributions; the method is valid for both sub- and supercritical flow condition; the average error obtained across all predicted series is 2.09% and the cross correlation is within 92% of accuracy for all trapezoidal and circular verification cases.

The accuracy of the predictions was found to be somewhat higher for trapezoidal channels compared to circular channels. This can be a reflection of the consistency of the input information which in the case of circular channels is less, i.e. the distributions are less uniform. The variation of the shape of the wetted section with the increase of the water level appears to be the reason of such variability.

The present investigation was based on a database formed by a limited number of experimental test cases, particularly for the case of circular sections. Nevertheless, the prediction results were satisfactory. The robustness of the methodology should be further tested with a larger training database containing further combinations of hydraulic parameters and section dimensions, and additional validation cases. This would help to ensure the generality of the weighting factors, and therefore, the overall accuracy of the prediction models. Based on the analysis presented here it is clear that the method works for a relatively low number of input data series which in this research ranged between 5 and 13 data series in the clusters. Furthermore, the linear regression models were demonstrated to be adequate estimators for the relatively smooth bed shear stresses studied, as they were able to capture the rates of shear stress variation with accuracy. It is notable that such simple estimators might not be accurate when modelling more complex configurations and patterns, and a more robust estimator (e.g. artificial neural network) may be more suitable depending on the degree of non-linearity observed in the input data.