On the Difference of River Resistance Computation between the k − ε Model and the Mixing Length Model

: River resistance characteristics, which can be reﬂected by the resistance factor, have an impact on ﬂow and sediment transport. In the classical theory, Prandtl proposed the mixing length model for the simulation of the turbulence, and von K á rm á n established the logarithmic formula of the ﬂow velocity distribution. Based on that, the expression of the resistance factor can be derived. With the development of the numerical technology, the k − ε model has been widely applied in the channels computation. However, for the different closure ways of the Reynolds stress in turbulence equations, the outcomes of the k − ε model and the Prandtl mixing length model are not exactly identical. In this paper, both qualitative and quantitative studies are carried out on the difference between these two models, with respect to the resistance factor. This difference is evaluated by the ratio of the resistance factor computed with the two models. The result shows that with the increment of the relative ﬂow depth, the ratio ﬁrst increases and then decreases. Moreover, it is also affected by the bed slope. Therefore, the difference should be taken into account when a comparison is made between the simulation results of the k − ε model and the classical theory of river mechanics.


Introduction
In the (vertical) two-dimensional (2D) steady uniform open channel flow, as shown in Figure 1, the vertical distribution of velocity u can be expressed by the Prandtl-von Kármán logarithmic law, as follows: where z is the vertical coordinate; u is the time-averaged velocity in the x direction; u * = τ b /ρ is the shear velocity, τ b is the bed sheer stress; κ is the von Kármán constant, κ = 0.4; k s is the equivalent bed-roughness and is discussed in Section 3; B s = φ B s (Re * ) is the roughness function, Re * = u * k s /ν is the roughness Reynolds number; ν is the kinematic viscosity coefficient. B s can be determined by experimental curve [1,2] or calculated approximately by the following formula [1]: Remarkably, when ≥ * 70

Re
, the flow is the rough turbulence, and thus ≡ 8.  where h is the depth, the Chézy resistance factor c is given by the following equation: which is also the classical formula for c . For the two-dimensional open channel uniform flow, there exists the following equation: (4) where S is bed slope along the channel. By combining Equations (3) with (4), u can be expressed as follows: = = * u cu c ghS (5) Equation (5) is in the form of the Chézy formula, in which the dimensionless resistance factor c has a relationship with the usual Chézy coefficient given by the following: The deduction of Equation (1) relies on the closure of the Reynolds stress in turbulence equations using the Prandtl mixing length model [2] (the details on how Equation (1) is derived from the Prandtl mixing length model can be seen in [3,4]). Hence, Equation (3) is also based on the Prandtl mixing length model. However, with the development of numerical technology, the ε − k turbulence model is more widely applied than the Prandtl mixing length model. In recent years, by employing the ε − k turbulence model, Wu et al. [5] simulated the flow and sediment transport in a 180° channel bend with movable bed in the laboratory; Zhang et al. [6] developed the two-dimensional ε − k turbulence model in curvilinear coordinates to simulate the behavior of turbulent flow in an open channel partially covered with vegetation; Thi et al. [7] investigated the turbulent flow around the local scour occurring around hydraulic structures with the ε − k turbulence model. On the basis of the three-dimensional geometry model with a geological model, Wang et al. [8] presented a threedimensional ε − k turbulence model considering Bingham liquid character. Besides these, the ε − k turbulence model is also applied in solving problems of other fields, such as air cooling technologies Remarkably, when Re * ≥ 70, the flow is the rough turbulence, and thus B s ≡ 8.5; the value is independent of Re * . As a semi-empirical solution, Equation (1) is widely applied in the calculation of open channel flows.
The dimensionless Chézy resistance factor c is defined as c = u/u * [1], where u is the depth-averaged velocity. Considering u being equal to u at the relative height z/h = e −1 = 0.368 [1], where h is the depth, the Chézy resistance factor c is given by the following equation: which is also the classical formula for c. For the two-dimensional open channel uniform flow, there exists the following equation: where S is bed slope along the channel. By combining Equations (3) with (4), u can be expressed as follows: u = cu * = c ghS (5) Equation (5) is in the form of the Chézy formula, in which the dimensionless resistance factor c has a relationship with the usual Chézy coefficient given by the following: The deduction of Equation (1) relies on the closure of the Reynolds stress in turbulence equations using the Prandtl mixing length model [2] (the details on how Equation (1) is derived from the Prandtl mixing length model can be seen in [3,4]). Hence, Equation (3) is also based on the Prandtl mixing length model. However, with the development of numerical technology, the k − ε turbulence model is more widely applied than the Prandtl mixing length model. In recent years, by employing the k − ε turbulence model, Wu et al. [5] simulated the flow and sediment transport in a 180 • channel bend with movable bed in the laboratory; Zhang et al. [6] developed the two-dimensional k − ε turbulence model in curvilinear coordinates to simulate the behavior of turbulent flow in an open channel partially covered with vegetation; Thi et al. [7] investigated the turbulent flow around the local scour occurring around hydraulic structures with the k − ε turbulence model. On the basis of the three-dimensional geometry model with a geological model, Wang et al. [8] presented a three-dimensional k − ε turbulence model considering Bingham liquid character. Besides these, the k − ε turbulence model is also applied in solving problems of other fields, such as air cooling technologies [9], automotive applications of on-board hydrogen storage [10], internal combustion engines scale-resolving simulation [11], pollutant removal [12], etc.
In general, the treatment of each layer in the k − ε turbulence model is shown in Figure 2, in which u 1 is the velocity near the bed; z 1 is a selected near-bed height. This treatment assumes that u conforms to the Prandtl-von Kármán logarithmic law only within the range of z 1 . Above this range, however, the velocity is computed by solving the turbulence equation closed by the k − ε model. On account of that, the resistance factor c corresponding to the k − ε model is not exactly identical to the one calculated by Equation (3). Note that Equation (3) is indeed the outcomes corresponding to the Prandtl mixing length model applied to the whole range of the depth. Thus, the fundamental reason causing the divergence is that there are two turbulence models using different ways in the closure of the Reynolds stress in turbulence equations. In this paper, c represents the outcome of the k − ε model and c represents that of the Prandtl mixing length model to distinguish them. In addition, the ratio c/c can be regarded as the quantitative parameter measuring the difference between the two models.
Water 2018, 10, x FOR PEER REVIEW 3 of 18 [9], automotive applications of on-board hydrogen storage [10], internal combustion engines scaleresolving simulation [11], pollutant removal [12], etc. In general, the treatment of each layer in the ε − k turbulence model is shown in Figure 2, in which 1 u is the velocity near the bed; 1 z is a selected near-bed height. This treatment assumes that u conforms to the Prandtl-von Kármán logarithmic law only within the range of 1 z . Above this range, however, the velocity is computed by solving the turbulence equation closed by the ε − k model. On account of that, the resistance factor c corresponding to the ε − k model is not exactly identical to the one calculated by Equation (3). Note that Equation (3) is indeed the outcomes corresponding to the Prandtl mixing length model applied to the whole range of the depth. Thus, the fundamental reason causing the divergence is that there are two turbulence models using different ways in the closure of the Reynolds stress in turbulence equations. In this paper, c represents the outcome of the ε − k model and ′ c represents that of the Prandtl mixing length model to distinguish them. In addition, the ratio ′ c c can be regarded as the quantitative parameter measuring the difference between the two models. In this study, the hydraulic parameters of natural rivers are collected, from which c and ′ c can be computed by using the above two models. On this basis, the variation pattern of ′ c c is studied.

Computational Methods of Resistance Factors by Using Different Models
Either way, the basic steps of computing the resistance factor are as follows:


Solve the turbulence equations to obtain the vertical distribution of u;  Average u over the depth to obtain u ; and  Calculate the resistance factor c or c′ by the ratio of u to * u .
Therefore, the core issue is how to solve the turbulence equations. The solving methods using the Prandtl mixing length model and the ε − k model are given separately in the below.

Prandtl Mixing Length Model
As stated in the introduction, solving the turbulence equations closed by the Prandtl mixing length model actually gives the Prandtl-von Kármán logarithmic distribution of u (i.e., Equation (1) in Section 1), which is in the continuous form. After integrating this equation, it can be identified that  In this study, the hydraulic parameters of natural rivers are collected, from which c and c can be computed by using the above two models. On this basis, the variation pattern of c/c is studied.

Computational Methods of Resistance Factors by Using Different Models
Either way, the basic steps of computing the resistance factor are as follows: • Solve the turbulence equations to obtain the vertical distribution of u; • Average u over the depth to obtain u; and • Calculate the resistance factor c or c by the ratio of u to u * .
Therefore, the core issue is how to solve the turbulence equations. The solving methods using the Prandtl mixing length model and the k − ε model are given separately in the below.

Prandtl Mixing Length Model
As stated in the introduction, solving the turbulence equations closed by the Prandtl mixing length model actually gives the Prandtl-von Kármán logarithmic distribution of u (i.e., Equation (1) in Section 1), which is in the continuous form. After integrating this equation, it can be identified that u is at the relative height z/h = e −1 = 0.368. For that, Equation (3) is used to compute c in this paper.

k − ε Model
Generally, a continuous form of velocity distribution cannot be obtained by closing and solving the turbulence equation with the k − ε model; thus, the discrete method is used. On the conditions of  [2] can be simplified into the following: where, u and w are the fluctuant velocity components in x and z directions, respectively; denotes the time-averaged value of the content within the angle brackets, and thus − u w actually denotes the Reynolds stress per unit mass.
According to the Boussinesq hypothesis, for the turbulent flow, ν t is defined as the eddy viscosity. Reynolds stress has a relationship with the time-averaged velocity gradient given by the following: Then, Equation (7) can be transferred to the following: Equation (9) is just the basic equation for the computation of c by using the k − ε model. Then, the finite volume method is adopted to discretize this equation.
The integral of Equation (10) on the grid of the i-th layer, as shown in Figure 3, gives the following: Water 2018, 10, x FOR PEER REVIEW 4 of 18 Generally, a continuous form of velocity distribution cannot be obtained by closing and solving the turbulence equation with the k ε − model; thus, the discrete method is used. On the conditions of the two-dimensional uniform flow in open channels, the RANS (Reynolds-averaged Navier-Stokes) equation in the x direction [2] can be simplified into the following: where, ′ u and ′ w are the fluctuant velocity components in x and z directions, respectively; denotes the time-averaged value of the content within the angle brackets, and thus − ′ ′ u w actually denotes the Reynolds stress per unit mass.
According to the Boussinesq hypothesis, for the turbulent flow, ν t is defined as the eddy viscosity. Reynolds stress has a relationship with the time-averaged velocity gradient given by the following: Then, Equation (7) can be transferred to the following: Equation (9) is just the basic equation for the computation of c by using the k ε − model. Then, the finite volume method is adopted to discretize this equation.
The integral of Equation (10) on the grid of the i-th layer, as shown in Figure 3, gives the following: In which the differential terms can be further discretized with the central differencing scheme, as follows: The substitution of Equations (11) and (12) into Equation (10) gives the following: In which the differential terms can be further discretized with the central differencing scheme, as follows: The substitution of Equations (11) and (12) into Equation (10) gives the following: which is the general discrete form of Equation (9).

of 18
As shown in Figure 4, for the situation that i = 1, Equation (10) is transferred to the following: Water 2018, 10, x FOR PEER REVIEW 5 of 18 which is the general discrete form of Equation (9). As shown in Figure 4, for the situation that = 1 i , Equation (10) is transferred to the following: The derivative term in the first term on the right-hand side can be discretized with the central differencing scheme as follows: According to the bed boundary condition, the second term on the right-hand side can be calculated as follows: Substituting Equations (15) and (16) into Equation (14), one obtains the following: As shown in Figure 5, for the situation that = i N , Equation (10) is transferred to the following:  The derivative term in the first term on the right-hand side can be discretized with the central differencing scheme as follows: According to the bed boundary condition, the second term on the right-hand side can be calculated as follows: Substituting Equations (15) and (16) into Equation (14), one obtains the following: As shown in Figure 5, for the situation that i = N, Equation (10) is transferred to the following: Water 2018, 10, x FOR PEER REVIEW 5 of 18 which is the general discrete form of Equation (9). As shown in Figure 4, for the situation that = 1 i , Equation (10) is transferred to the following: The derivative term in the first term on the right-hand side can be discretized with the central differencing scheme as follows: According to the bed boundary condition, the second term on the right-hand side can be calculated as follows: Substituting Equations (15) and (16) into Equation (14), one obtains the following: As shown in Figure 5, for the situation that = i N , Equation (10) is transferred to the following:  The derivative term in the second term on the right-hand side can be discretized with the center difference scheme as follows: According to the free surface boundary condition, the first term on the right-hand side can be calculated as follows: Substituting Equations (19) and (20) into Equation (18), one obtains the following: For the value of ν t , at the interface of the grid, the linear interpolation is resorted to the following: whereas the value of ν t at the center of the grid is calculated by the k − ε model, as follows: where k is the turbulent kinetic energy; ε is the dissipation rate; P z = (ν + ν t )( ∂u/∂z) 2 , is the production term; the standard values of the model coefficients are employed [5,13]: (25) and (26) are both differential equations, and they are discretized with the same scheme as Equation (9). In addition, for Equations (25) and (26), the bed boundary conditions [5] are as follows: and the free surface boundary conditions [5] are as follows: Combing Equation (13) with Equations (17) and (21), the following tridiagonal linear system is finally derived, as follows: where, However, the determinant of the coefficient matrix of this system turns out to be zero and the rank is N − 1 (< N, non-full rank), which means that there exists infinitely many solutions to this system. Furthermore, this system is derived as follows.
By means of accumulative summation from the top to the bottom, Equations (29) can be transferred to the following: The last equation of Equations (30) is as follows: Note that N∆z = h. Then, after transposition, Equation (31) can be transferred to which is the same equation as Equation (4). From Equations (30), it follows that only the difference between the velocities of any two adjacent points can be solved. Thus, in order to get the absolute value of the velocity, another supplementary equation is needed. As stated in the introduction, it is generally assumed that the near-bed velocity within the range of z 1 (z 1 ≤ 0.2h) [4,14] conforms to the Prandtl-von Kármán logarithmic law, which means the following: Then, the velocity u i (i = 1, 2, . . . , N) of any point can be obtained by successively solving Equation (33), and then Equations (30) from the top to the bottom. The depth-averaged velocity u can be evaluated as follows: Thus, the resistance factor c can be evaluated as c = u/u * . The above calculation still needs iteration for γ in the coefficient matrix of Equations (30), involving ν t = c ν k 2 /ε, which is determined by the differential Equations (25) and (26). The computation methods of c using the Prandtl mixing length model and the k − ε model are given, as above. Figure 6 shows the typical computations of velocity distribution by using these two models. The parameters of the Rio Grande River [15] are used in both computations: h = 0.332 m, k s = 0.028 m, and S = 0.830 × 10 −3 . The result corresponding to Prandtl mixing length model, which is calculated by Equation (1), is shown in the continuous form; in contrast, the result corresponding to the k − ε model, which is calculated by dividing the whole depth into ten layers, is shown in the discrete form. From Figure 6, an approximate coincidence can be detected between the two velocity distribution pattern computed by different models within the range of z ≤ 0.2h; however, the difference displays, to some extent, within the range of z > 0.2h. Indeed, the Prandtl mixing length model and the k − ε model treat the shear force between the flow layers in different ways. This is the main reason the resistance factors computed with these two models are not identical.
Then, the velocity i u (i = 1, 2, …, N) of any point can be obtained by successively solving Equation (33), and then Equations (30) from the top to the bottom. The depth-averaged velocity u can be evaluated as follows: Thus, the resistance factor c can be evaluated as = * c u u .
The above calculation still needs iteration for γ in the coefficient matrix of Equations (30), , which is determined by the differential Equations (25) and (26).
The computation methods of c using the Prandtl mixing length model and the k ε − model are given, as above. Figure 6 shows the typical computations of velocity distribution by using these two models. The parameters of the Rio Grande River [15] are used in both computations: h = 0.332 m, s k = 0.028 m, and S = 0.830 × 10 −3 . The result corresponding to Prandtl mixing length model, which is calculated by Equation (1), is shown in the continuous form; in contrast, the result corresponding to the k ε − model, which is calculated by dividing the whole depth into ten layers, is shown in the discrete form. From Figure 6, an approximate coincidence can be detected between the two velocity distribution pattern computed by different models within the range of 0.2 z h ≤ ; however, the difference displays, to some extent, within the range of 0.2 z h > . Indeed, the Prandtl mixing length model and the k ε − model treat the shear force between the flow layers in different ways. This is the main reason the resistance factors computed with these two models are not identical.

Determination of Equivalent Bed-Roughness
In natural rivers, bed forms are induced by sediment transport. This means that the bed resistance to flow includes not only the grain resistance, but also the bed-form resistance-For some special cases, the resistance of other factors (e.g., the emergent vegetation) should also be taken into account [16], but they are not considered in this study. In this study, the van Rijn method [17] is employed to determine s k , which takes the influence of these two kinds of resistance into consideration.
The expression of s k proposed by van Rijn [17] is as follows:

Determination of Equivalent Bed-Roughness
In natural rivers, bed forms are induced by sediment transport. This means that the bed resistance to flow includes not only the grain resistance, but also the bed-form resistance-For some special cases, the resistance of other factors (e.g., the emergent vegetation) should also be taken into account [16], but they are not considered in this study. In this study, the van Rijn method [17] is employed to determine k s , which takes the influence of these two kinds of resistance into consideration.
The expression of k s proposed by van Rijn [17] is as follows: where D is the characteristic bed material diameter; ∆ and δ are the bed-form height and the bed-form steepness, respectively, which are determined as follows: where T is a transport stage parameter, where (τ b ) f is the bed shear stress related to grains; (τ b ) cr is the critical bed-shear stress for initiation of motion of grains, which can be determined by experimental curve [1,18] or calculated approximately by the following formula: where Ξ is the dimensionless material number, which is defined as the following: where ρ s is the density of the bed material, ρ s = 2650 kg/m 3 ; ρ is the density of water, ρ = 1000 kg/m 3 .
In the work of van Rijn [17], D was identified with D 90 ; in this study, however, as a result of the lack of information of D 90 reported in the references, D is thus identified from all the available data within a range of D 50 ∼D 90 . Practices have indicated that although the exact value of c/c varies with the D selected, the pattern of c/c versus D does not change in a noticeable manner. This is further discussed and evidenced in Section 5.1.

Data Sources and Range
The data of natural rivers were collected from previous publications [15,[19][20][21][22][23], including the discharge Q, the water depth h, the bed slope S, and the characteristic bed material diameter D.
There are 213 groups data in total. The detailed information is shown in Table 1.

Analysis of Influential Factors
According to Sections 2 and 3, c/c is influenced by the following seven dimensional parameters: ρ, g and D are selected as the basic dimensionally independent variables that represent the dynamic, kinematic, and geometric characteristics, respectively. Then, the Buckingham π theorem is employed for the dimensional analysis. One determines for c/c that is influenced by the following four dimensionless parameters: within these parameters, ρ s /ρ is considered as constant (=2.65) in this paper. Thus, the influence of ρ s /ρ is not considered in the following analyses. The product of D gD/ν, √ h/D and √ S gives the following: where Re D is the grain size Reynolds number. According to Equation (35), the roughness Reynolds number Re * can be expressed as follows: Here, Re ∆ = u * ∆/ν, which is defined as the bed-form Reynolds number. According to Equation (44), Re D is included in Re * . Also, note that Re * only appears in the roughness function B s in Equations (3) and (33). For the data of natural rivers collected in this study, it can be proved that Re * ≥ 70, meaning that the flow is an invariably rough turbulent, for which B s ≡ 8.5; as a consequence, in this study, B s is independent of Re * and thus Re D . Then, the influence of D gD/ν is neglected.
In the final analysis, c/c is a function solely depending on h/D and S, for which one can write the following: Note that h/D is referred to as the relative water depth [1]. The function form of f will be determined in the following subsection.

Computational Results and Analyses
First, by clustering the data of S and h/D corresponding to all the collected natural rivers, the relation of S versus h/D is plotted in Figure 7 with log−log axes. From Figure 7, it can be seen that all the (h/D, S) points are approximately distributed around the following five horizontal lines: S = 10.050 × 10 −3 , 1.540 × 10 −3 , 0.806 × 10 −3 , 0.190 × 10 −3 , and 0.072 × 10 −3 . Note that despite a considerable scatter, the (h/D, S) points around S = 0.190 × 10 −3 still follow this pattern. Thus, a classification will be made in the following contents for all the collected natural rivers on the basis of these five values of S. Water 2018, 10, x FOR PEER REVIEW 11 of 18 It is assumed that curves 1, 2, 4, and 5 can be obtained by translating curve 3 on the (h/D, c/c ) plane. Thus, two translation parameters ∆x and ∆y (see Figure 8) Apparently, both of ∆x and ∆y should be functions of S. Apparently, both of Δx and Δy should be functions of S . ( )  The values of ∆x and ∆y for curves 1, 2, 4, and 5 can be obtained by fitting Equation (47) to the corresponding (h/D, c/c ) points. Taking also into account ∆x = 0 and ∆y = 0 for curve 3, the relations of ∆x versus S, and ∆y versus S are plotted in Figure 9a,b, respectively, with lin−log axes. It is indicated that both ∆x and ∆y are linear with the logarithm of S (ln S). With the premise that the fitting lines pass through the points corresponding to curve 3 in Figure 8 This is the functional form sought for f . Note that there exists a vertical asymptote on the right side of each c/c ∼h/D curve, as follows: Thus, Equation (51) is only applicable to the following: According to Equation (50), together with all the curves in Figure 8, the following aspects can be summarized: (a) For a specific S, c/c first increases and then decreases with the increment of h/D; the variation pattern of c/c versus h/D is shown as a convex upward curve; (b) The peak value of the convex upward c/c ∼ h/D curve increases with the increment of S; Figure 9b, along with Equation (49), indicates that the peak value is approximately linear with ln S; and (c) The abscissa of the peak value decreases with the increment of S; Figure 9a, along with Equation (48), indicates that the abscissa of the peak value is approximately linear with ln S as well.
Moreover, from the (h/D, c/c ) points in Figure 8, it follows that the values of c/c are generally greater than unity. This means that for the collected natural rivers, the resistance factors computed with the k − ε model are generally greater than those computed with the Prandtl mixing length model. However, it cannot be ignored that the circumstances where the values of c/c are smaller than unity can also be detected, as long as the corresponding values of h/D are sufficiently large (see Figure 8, the points on the lower right sides of curves 3 and 4). According to Equation (50), together with all the curves in Figure 8, the following aspects can be summarized: Moreover, from the ( ) , h D c c′ points in Figure 8, it follows that the values of c c′ are generally greater than unity. This means that for the collected natural rivers, the resistance factors computed with the k ε − model are generally greater than those computed with the Prandtl mixing length model. However, it cannot be ignored that the circumstances where the values of c c′ are smaller than unity can also be detected, as long as the corresponding values of h D are sufficiently large (see Figure 8, the points on the lower right sides of curves 3 and 4). In practice, if the data collected in a natural river include the water depth h, the bed slope S, and the characteristic bed material diameter D, all of which satisfy Equation (52), then the flow discharge Q corresponding to the k − ε turbulence model can be simply calculated with the following steps, instead of by running a complex numerical model:

•
Calculate the resistance factor c corresponding to the Prandtl mixing length model by employing Equation (3)

Discharge Verification
The computed flow discharge Q and Q corresponding to the k − ε model and the Prandtl mixing length model, respectively, are verified by comparing them with the field-measured flow discharge Q m . The results are shown in Figure 10. It can be seen that both the (Q m , Q) and (Q m , Q ) sets of points are evenly distributed around the line y = x.
The computed flow discharge Q and ′ Q corresponding to the k ε − model and the Prandtl mixing length model, respectively, are verified by comparing them with the field-measured flow discharge m Q . The results are shown in Figure 10. It can be seen that both the m ( , ) Q Q and  Therefore, each of two turbulence models has its own advantages and disadvantages, if considered from different aspects.
Considering the function f proposed in Section 4.3, which quantitatively describes the difference between two turbulence models in the computation of river resistance, Q and ′ Q corresponding to the ε − k model and the Prandtl mixing length model, respectively, seem interconvertible, as follows: Furthermore, the frequency distributions of the computational errors of Q and Q m relative to Q m are shown in Figure 11. Both approximately conform to the normal distribution as follows: Q relative error ∼ N µ = 21.65%, σ 2 = 28.79 2 (53) Q relative error ∼ N µ = 7.78%, σ 2 = 31.11 2 (54) where N is the normal distribution; µ is the mathematical expectation of normal distribution; and σ 2 is the variance.
As above, it is further indicated that the simulation results corresponding to the ε − k model and the Prandtl mixing length model are not identical. When comparing them, the difference between the river resistance characteristics reflected by these two models should be taken into account. Figure 11. Frequency distribution of computed discharge relative error.

Influence of the Selection of the Characteristic Bed Material Diameter
As stated in Section 3, the calculation of the equivalent bed-roughness s k involves a certain characteristic bed material diameter D. In this study, the available data of D are within a range of D50~D90. Indeed, with regards to Figure 8, even though the a different type of D is selected for each point, the same trend of the computed ′ c c h D  curves can still be obtained. The points are still From the values of µ in Equations (53) and (54), the mean relative errors of Q and Q are 21.65% and 7.78%, respectively. This means that the (Q m , Q ) points in Figure 10b are closer to the line y = x than the (Q m , Q) points in Figure 10a. From the values σ 2 of in Equations (53) and (54), the variance of the relative errors of Q and Q are 28.79 2 and 31.11 2 , respectively. This means that the (Q m , Q ) points in Figure 10b are more scattered than the (Q m , Q) points in Figure 10a. Therefore, each of two turbulence models has its own advantages and disadvantages, if considered from different aspects.
Considering the function f proposed in Section 4.3, which quantitatively describes the difference between two turbulence models in the computation of river resistance, Q and Q corresponding to the k − ε model and the Prandtl mixing length model, respectively, seem interconvertible, as follows: As above, it is further indicated that the simulation results corresponding to the k − ε model and the Prandtl mixing length model are not identical. When comparing them, the difference between the river resistance characteristics reflected by these two models should be taken into account.

Influence of the Selection of the Characteristic Bed Material Diameter
As stated in Section 3, the calculation of the equivalent bed-roughness k s involves a certain characteristic bed material diameter D. In this study, the available data of D are within a range of D 50 ∼D 90 . Indeed, with regards to Figure 8, even though the a different type of D is selected for each point, the same trend of the computed c/c ∼ h/D curves can still be obtained. The points are still distributed around the same one out of the five curves, although positioned at different positions.
Only a few references report various types of D for a certain river; (e.g., from Rio Grande River [15]) four groups of data corresponding to different river reaches of both D 50 and D 84 can be extracted. In Figure 8, the (h/D, c/c ) points corresponding to these river reaches are computed with D 50 ; these points, together with curve 3, which they are fitted to, are replotted in Figure 12.
To show the influence of the selection of the type of D, the values of c/c corresponding to the same river reaches are further recomputed with D 84 , and the corresponding (h/D, c/c ) points are plotted in Figure 12, too. It can be seen from Figure 12  However, according to Rodi [13], in a case where the production of the turbulence kinetic energy is much larger than its dissipation, the value of ν c can be as small as about 0.05. Figure 13

Influence of the Coefficients in the k − ε Model
As stated in Section 2.2, the k − ε turbulence model involves five coefficients and their standard values are used, namely, c v = 0.09, σ k = 1.0, σ ε = 1.3, c ε1 = 1.44, and c ε2 = 1.92. However, according to Rodi [13], in a case where the production of the turbulence kinetic energy is much larger than its dissipation, the value of c v can be as small as about 0.05. Figure 13 displays a partial result of c/c versus h/D, with different values of c v : c v = 0.09 and c v = 0.05. From Figure 13, it follows that the computed values of c/c with c v = 0.05 are somewhat larger than those with c v = 0.09. Consequently, the result presented in this paper still depends on the values of c v in the k − ε model. Nevertheless, it should also be noted that with either value of c v , c/c appears invariably as a function of h/D; that is, it can never be constant.
For lack of knowledge on how the other coefficients change according to different cases, their influences are not included in this paper. How these coefficients, including c v , affect the resistance factor ratio c/c needs further studies. However, according to Rodi [13], in a case where the production of the turbulence kinetic energy is much larger than its dissipation, the value of ν c can be as small as about 0.05. Figure 13  model. Nevertheless, it should also be noted that with either value of v c , c c′ appears invariably as a function of h D; that is, it can never be constant.
For lack of knowledge on how the other coefficients change according to different cases, their influences are not included in this paper. How these coefficients, including v c , affect the resistance factor ratio c c′ needs further studies.

Conclusions
With the data of natural rivers, the computational results corresponding to the k − ε model and the Prandtl mixing length model are compared and analyzed. The following has been suggested: (a) The resistance factors corresponding to these two models are not identical. A new expression Equation (50), has been derived to quantify it. The difference of these two models, evaluated by the ratio c/c , first increases and then decreases with the increment of the relative water depth h/D. This variation pattern can be represented by a convex upward curve, of which the peak value increases and the abscissa of the peak value decreases with the increment of the bed slope S. (b) Both turbulence models can generate reasonable results with acceptable errors with respect to the flow discharges in natural rivers. Each model has its own advantages and disadvantages. Within the range of the collected natural rivers, the mean relative error of the computed flow discharges corresponding to the Prandtl mixing length model is less than that corresponding to the k − ε model; in contrast, the error distribution corresponding to the former is more scattered than that corresponding to the latter. Furthermore, the flow discharges computed with these two models can be converted into each other by taking into account Equation (55). (c) The velocity distribution patterns corresponding to these two models within the range of z ≤ 0.2h are approximately consistent, whereas the difference displays to some extent within the range of z > 0.2h.
(d) For river computation issues, the simulation results corresponding to the k − ε model and the Prandtl mixing length model are not exactly identical, although both can well simulate the flow. Accounting for that, the difference between the river resistance characteristics reflected by these two models should be considered.