Formulae of Sediment Transport in Steady Flows (Part 1)

Abstract

This paper makes an attempt to answer why the observed critical shear stress for incipient sediment motion sometimes deviates from the Shields curve largely, and the influence of vertical velocity is analyzed as one of the reasons. The data with d50 = 0.016 ∼ 29.1 mm from natural streams and laboratory channels were analyzed. These measured data do not always agree with the Shields diagram’s prediction. The reasons responsible for the deviation have been re-examined and it is found that, among many factors, the vertical motion of sediment particles plays a leading role for the invalidity of Shield’s prediction. The positive/negative deviations are associated with the up/downward vertical velocity in decelerating/accelerating flows, and the Shields diagram is valid only when flow is uniform. A new theory for critical shear stress has been developed, a unified critical Shields stress for sediment transport has been established, which is valid to predict the critical shear stress of sediment with/without vertical motion.

1. Introduction

The incipient motion of sediment is one of the most important topics in sediment transport. Generally, two methods are available in the literature to express quantitatively it, the shear stress approach and velocity approach [1]. The latter assumes that if the mean velocity excesses its critical velocity, then the sediment motion can be observed. The former used by researchers represents the force acting on a particle. Shields [2] is the earliest one who used the shear stress approach, or Shields number τ/(ρ_s-ρ)gd₅₀ versus the Reynolds number, and he obtained a famous Shields curve to express sediment initiation. Francalanci et al. [3] interpreted the Shields number as the ratio of streamwise/vertical forces using the following form:

$${\tau }_{\ast c}=\frac{{\tau }_c}{\left({\rho }_s-\rho \right){\mathit{gd}}_{50}}=\frac{4}{3}\frac{{\tau }_c\frac{\pi }{2}{\left(\frac{d_{50}}{2}\right)}^2}{\left({\rho }_s-\rho \right)g\frac{4\pi }{3}{\left(\frac{d_{50}}{2}\right)}^3}$$E1

where ${\tau }_c$ (=ρu_*c²) is the critical shear stress for the median grain size of sediment, $d_{50}$; g is the gravitational acceleration; $u_{\ast c}$ is critical shear velocity; ${\rho }_s$ and $\rho$ are the sediment and fluid densities, respectively. The shear stress exerted by the fluid must be higher than the critical shear stress ${\tau }_c$ to initiate sediment motion at the bed. Based on available experimental data, Shields in 1936 found that the Shields number ${\tau }_{\ast }$ depends on the particle Reynolds number ($R_{\ast }$), i.e.,

$$R_{\ast }=\frac{u_{\ast c}{d}_{50}}{\nu }$$E2

where ν is the kinematic viscosity of the fluid.

The original Shields diagram has been reproduced and modified by many researchers. A comprehensive review has been done by many researchers [4, 5], in which some significant deviations of the observed critical shear stress from the standard Shields curve were observed. This has attracted extensive research by notable investigators, and some factors leading to the data scatter have been identified and discussed.

Some researchers believe that the definition of the incipient motion may cause the invalidity of Shields diagram, as the incipient motion depends more or less on the experimental observers’ subjective judgment. To address this, criteria like “individual initial motion”, “several grains moving” and “weak movement” has been introduced to express the incipient motion [6]. Subsequently, an error band has been included in the modified Shields diagram [7].

Other researchers attribute the large discrepancy to the stochastic nature of turbulence and sediment shape, its orientation, or exposure, protrusion [8, 9, 10, 11]. It is natural to expect that when sediment is non-uniform, the critical condition is very difficult to determine, as the larger particles could move relatively easily than the finer one that is sheltered [12].

Over the past eight decades, the incipient motion has been extensively studied again and again [4], because the Shields diagram has been found invalid to predict the critical shear stress of sediment transport in some circumstances. The invalidity is not fully explained, some researchers ascribe it to sediment’s characteristics, the other believe these deviations are caused by the flow conditions i.e., non-uniformity of flow [13].

Iwagaki [14] firstly linked the wide scatter in Shields diagram with flow’s non-uniformity based on his observation: for the same sediment by the same experimenter, the observed critical shear stress in non-uniform flows largely deviates from that in uniform flows. Afzalimhr et al. [13] confirmed Iwagaki’s results, they found experimentally that in decelerating flows, the critical shear stress is considerably below the Shields’ prediction, and their experimental data are in complete disagreement with the Shields diagram. Other experimental researchers [15, 16] obtain similar results as Afzalimhr et al.’s [13] who claimed “…there is no universal value for ${\tau }_{\ast }$”. Likewise, Buffington and Montgomery [4] also agreed “less emphasis should be given on choosing a universal ${\tau }_{\ast }$”.

Some researchers try to explain the large discrepancy between predicted and measured critical shear stress by considering channel’s characteristics, such as the channel shape and channel slope [16, 17, 18, 19, 20, 21, 22]. “the well-known Shields criterion is insufficient for large slope” was observed by Graf and Suszka [23], while Lamb et al. [24] comprehensively re-visited and examined almost all published datasets, and concluded that the critical shear increases with channel slope, this is totally different from the common sense that predicts increased mobility with increasing channel slope due to the added gravitational force in the downstream direction. But Chiew and Parker’s experiments [17] in very steep channels show that the critical shear stress is decreased, contrary to Lamb et al.’s [24] conclusion.

Therefore, the brief literature review shows that Shields diagram cannot predict the critical shear stress well and there are many different potential causes for the deviation, among them, it is necessary to clarify how the channel-bed slope and non-uniformity of flow affect the critical shear stress for sediment motion. The primary objectives of the present study are to

1. investigate the mechanism that causes the invalidity of the Shields curve for the incipient motion of sediment transport;

2. examine why the Shields number depends on the water depth’s variation or channel slope;

3. establish a universal Shields diagram that is valid for all data available in the literature; and

4. verify the newly established equations using data from the literature.

2. Theoretical considerations of influence of vertical velocity on the critical shear stress

The author has been systematically investigating the role of vertical velocity on the mass and momentum transfer and has obtained a series of important and interesting conclusions [25, 26, 27]. It is found that omission of vertical velocity in our existing theorem of sediment transport makes many phenomena unexplainable. For example the presence of vertical velocity in non-uniform flows leads to the deviation of measured Reynolds shear stress from the linear distribution from the free surface to the bottom, consequently the upward velocity causes the positive deviation of velocity from the log-law or the wake-law is needed to express the velocity distribution, and the downward velocity results in the dip-phenomenon, or the maximum velocity is submerged and does not occur at the free surface as the log-law predicts. As the momentum and mass transfers are closely related to each other, it is interesting to investigate how the vertical motion affects sediment transport.

As a continuous effort, this study investigates the influence of upward/downward velocity on sediment incipient motion and the validity of Shields’ diagram. Figure 1 shows how a river flow interchanges with groundwater and the Darcy law tells that vertical velocity is proportional to the hydraulic gradient, i.e., the suctions and injections inside groundwater can be expected in flood/dry seasons alternatively. The upward flow or injection flow may increase the sediment particles’ mobility, or the required critical shear stress is reduced due to the “buoyant effect”, which reduces the net settling velocity, mathematically

Figure 1.

The upward and downward vertical velocity generating from seepage face injection seepage.

$${\omega }^{\text{'}}=\omega -V_b$$E3

where ω = particle’s settling velocity in still water and ${\omega }^{\text{'}}$ = the net settling velocity subject to the vertical velocity of groundwater, $V_b$. The submerged weight in Eq. 1 can be represented by a drag force with the falling velocity ω in still water ($V_b=0$) as:

$$C_d\pi \frac{d^2}{4}\frac{{\mathit{\rho \omega }}^2}{2}=\pi \frac{4}{3}{\left(\frac{d}{2}\right)}^{3}g\left({\rho }_s-\rho \right)$$E4

where d is the particle diameter, C_d is the drag coefficient.

If the upward velocity $V_b$ of groundwater is so high and $V_b$ = ω, the net settling velocity of the particle becomes zero, thus the particle is neutrally suspended, i.e., liquefaction state. This is often observed during earthquakes. In such case, saturated soil loses its strength and stiffness, it is natural that the Shields diagram cannot predict the particle’s critical shear stress. Similarly, if the groundwater in Figure 2 is downward, then the net falling velocity ${\omega }^{\text{'}}$ should be higher than ω, the threshold critical shear stress should be unpredictable using the existing Shields diagram.

Figure 2.

Schematic diagrams showing interaction of streamwise and vertical motions after Ladson (2008), p99.

The above discussion clearly demonstrates that velocity $V_b$ in a sediment layer may cause the invalidity of Shields diagram, which is supported experimentally by many researchers [28] who conducted experiments by observing the critical shear stress subject to injection and suction flows. Lu et al. [29] has reviewed these experimental results comprehensively. The influence of vertical motion on the critical shear stress has been discussed by many researchers, the parameters used to express the vertical motion include (i) the hydraulic gradient of seepage, e.g., Cheng and Chiew [30]; (ii) the pressure variation in flows [3]. But, there is no research available to investigate the role of time-averaged vertical velocity on the incipient motion of sediment transport.

The introduction of apparent sediment density is similar to Francalanci et al.’s treatment [3]. Instead of modifying the sediment density, they modified the water’s density to eliminate the effect of pressure variation over time and space (like pressure induced by waves or bridge piers) on sediment’s critical shear stress. Their results show that higher pressure yields higher “apparent water density”, and lower pressure corresponds to lower “apparent water density”. They found that the Shields number shown in Eq. 1 is actually the ratio of friction force in the streamwise direction (i.e., τ_cπd²/2) to the net force in the vertical direction, i.e., downward gravitational force (=ρ_sgπd³/6) minus the upward buoyant force (=ρgπd³/6), so that the effective density of the particle is reduced from ρ_s to ρ_s-ρ. When a particle is experienced in the environment with upward velocity V_b, additional upward force will be generated, and the problem is how to determine the additional upward force F_vb induced by the upward velocity.

Eq. 4 shows that the denominator of Shields number, i.e., (ρ_s-ρ)gd can be replaced by 3C_dρω²/4, this is why the parameter of settling velocity ω can be used in the study of sediment incipient motion. Some researchers believe that ω is a parameter for suspended load and it should not be used to express the sediment initial motion. The first one who uses ω to discuss the critical velocity is Yang [1], and this treatment significantly simplifies the problem. But this treatment has been blamed by many researchers who believe that the initiation problem does not involve any settling process. Now, Eq. 4 clearly shows that it is logical to express the incipient motion with sediment settling velocity.

For the case shown in Figure 2, if the upward velocity is zero, this is a static problem and the net force balance in vertical direction is expressed in Eq. 4. When the upward velocity is non-zero, this becomes a dynamic problem where the lift force F_vb should be included, i.e., submerged weight minus F_vb must be balanced by the drag force with settling velocity ω’.

$$C_d^{\text{'}}\pi \frac{d^2}{4}\frac{{\mathit{\rho \omega }}^{\text{'}2}}{2}=\pi \frac{d^3}{6}g\left({\rho }_s-\rho \right)-F_{\mathit{vb}}$$E5

where F_vb = C_d,liftπd²ρV_b²/8. Let ${\rho }_s^{\text{'}}={\rho }_s-F_{\mathit{vb}}/\left(\pi d^{3}g/6\right)$ and inserting the apparent density of sediment into Eq. 5, then the net force in the vertical direction can be alternatively expressed as:

$$C_d^{\text{'}}\pi \frac{d^2}{4}\frac{{\mathit{\rho \omega }}^{\text{'}2}}{2}=\pi \frac{d^3}{6}g\left({\rho }_s^{\text{'}}-\rho \right)$$E6

In this study, apparent sediment density is introduced, and it depends on the vertical velocity of groundwater. Therefore, it is expected to have a relationship between the apparent sediment density and the settling velocity, similar to Einstein’s relativity theory that the length/time depends on velocity if the light’s speed is assumed to be constant. Therefore, the effect of vertical motion caused by pressure variation or seepage on sediment transport is eliminated after the introduction of apparent density. In other words, real lightweight particles motion must be the same as those with reduced settling velocity ${\omega }^{\text{'}}$ in terms of mobility when both have the same settling velocity. Hence, the apparent density can greatly simplify the mathematical treatment for the complex F_vb induced by vertical motions. From Eqs. (4) and (5), the relationship between the modified settling velocity and the apparent density can be expressed by:

$$\frac{{\rho }_s^{\text{'}}-\rho }{{\rho }_s-\rho }=\alpha {\left(\frac{\omega -V_b}{\omega }\right)}^2$$E7

where α is a coefficient (= C_d′/C_d) and α = 1 is assumed in this study to simplify the mathematical treatment. Because the drag coefficient depends on Re, this assumption is approximately correct if the particle size is very coarse or the Re does not change significantly with/without V_b. Eq. 7 tells that if $V_b$ is equal to zero, then ${\rho }_s^{\text{'}}$ is the same as the natural sediment; if $V_b$ is positive or upwards then ${\rho }_s^{\text{'}}$ is less than the density of natural sediment ${\rho }_s$, and the particles behave like “lightweight sand” or plastic sands; if $V_b$ = ω, then ${\rho }_s^{\text{'}}$ is same as the density of water or similar to neutrally buoyant milk; if $V_b$ is negative or downward, the higher apparent density of sediment behaves like heavy metals.

The vertical velocity $V_b$ in Figure 2 has the similar effect for the particles’ stability as the buoyancy effect, i.e., the submerged weight of the particles is no longer ${\rho }_s-\rho$, but ${{\rho }_s}_{\text{'}}^{}-\rho$, one may give the general expression of Shields number

$${\tau }_{\ast }^{\text{'}}=\frac{{\tau }_c^{\text{'}}}{\left({\rho }_s^{\text{'}}-\rho \right){\mathit{gd}}_{50}}$$E8

Substituting Eq. 7 into Eq. 8, one obtains:

$${\tau }_{\ast }^{\text{'}}=\frac{{\tau }_c^{\text{'}}}{\left({\rho }_s-\rho \right){\mathit{gd}}_{50}}{\left(\frac{\omega }{\omega -V_b}\right)}^2$$E9

or

$${\tau }_{\ast }^{\text{'}}=\frac{{\tau }_c^{\text{'}}}{\left({\rho }_s-\rho \right){\mathit{gd}}_{50}}{\left(\frac{1}{1-Y}\right)}^2$$E10

where Y = V_b/ω.

Eq. 9 or 10 generally expresses the influence of vertical velocity $V_b$ on the critical shear stress. It is clear that V_b can be induced by seepage in the sediment layer, it can be also inferred that V_b can be estimated by the Darcy Law using the hydraulic conductivity and hydraulic gradient. Obviously, Y = 1 means that the particles can be suspended in water, i.e., liquefaction. If Y > 1, it means that particles flow in the upward direction with a net velocity of V_b –ω, this may have a devastating impact on dikes in flood defense as it may cause piping failure. In the following section, the analysis shows that the vertical velocity, V_b is ubiquitous in open channel flows, which is induced by non-uniform flows.

3. Influence of non-uniform flow on the critical shields stress

Ideal uniform flow is very rare in natural conditions, flow rate and water depth/channel width keep always changing, i.e., non-uniform, as shown in Figure 3. It is interesting to discuss how accelerating or decelerating flows generate the vertical velocity. To simplify the discussion, it is assumed that the flow rate is constant, i.e., dh/dx (≠ 0). The 2-D continuity equation is:

Figure 3.

Non-uniform flows in open channel and the variation of water depth, in which u and v are mean velocities in x and y direction, respectively.

$$\frac{\partial \overline{u}}{\partial x}+\frac{\partial \overline{v}}{\partial y}=0$$E11

where $\overline{u}$ and $\overline{v}$ are the mean local velocity at any point in x and y directions, respectively. By integrating Eq. 11, one has

$$\overline{v}=-{\int }_0^y\frac{\partial \overline{u}}{\partial x}\mathit{dy}$$E12

∂$\overline{u}$/∂x > 0 means accelerating, thus Eq. 12 tells that accelerating flows yield a negative or downward vertical velocity; but decelerating flows generates the positive or upward vertical velocity, i.e., ∂$\overline{u}$/∂x < 0. Hence, the vertical velocity can be generated by in non-uniform flows.

For a channel with a constant width, its discharge per unit width can be expressed by:

$$Q/b=\mathit{Uh}$$E13

where $Q$ = discharge; $U$ is the depth-averaged velocity, $b$ is the channel width and h is the water depth. If Q/b could be constant in x direction, one has:

$$\frac{d\left(\mathit{Uh}\right)}{\mathit{dx}}=\frac{d}{\mathit{dx}}{\int }_0^h\overline{u}\mathit{dy}=0$$E14

The vertical velocity ${\overline{v}}_h$ at the free surface can be obtained from Eq. 12 using Leibniz’s rule, i.e.

$${\overline{v}}_h=-{\int }_0^h\frac{\partial \overline{u}}{\partial x}\mathit{dy}=-\frac{d}{\mathit{dx}}{\int }_0^h\overline{u}\mathit{dy}+{\overline{u}}_h\frac{\mathit{dh}}{\mathit{dx}}$$E15

where ${\overline{u}}_h$ is the horizontal velocity at the surface in the x direction. By inserting Eq. 14 into 15, one obtains:

$${\overline{v}}_h={\overline{u}}_h\frac{\mathit{dh}}{\mathit{dx}}$$E16

Eq. 16 shows that $\mathit{dh}/\mathit{dx}$ > 0, i.e., a decelerating flow generates ${\overline{v}}_h$ > 0, but $\mathit{dh}/\mathit{dx}$ < 0 or an accelerating flow yields the negative ${\overline{v}}_h$. Therefore, the vertical velocity in the main flows can also generate vertical velocity. Its interaction with groundwater can be obtained as V_b = V_g + V_s, or the groundwater velocity at the bed is jointly caused by Darcy velocity V_g and the vertical velocity caused by the main flow, V_s. Even the velocity on the solid–liquid interface may be very small, its importance for sediment transport should not be underestimated [3], and Eq. 16 shows that the vertical velocity has similar amplitude like the secondary current, i.e., about 1% of mean velocity.

Julien [5] replaced the Reynolds number in Shields’ diagram by dimensionless particle diameter:

$$d_{\ast }={\left[\frac{{\rho }_s-\rho }{\rho }\frac{{\mathit{gd}}_{50}^3}{{\nu }^2}\right]}^{1/3}$$E17

Similarly, d_* needs modification by introducing the apparent density with the following form:

$$d_{\ast }^{\prime }={\left[\frac{{\rho }_s^{\prime }-\rho }{\rho }\frac{{\mathit{gd}}_{50}^3}{{\nu }^2}\right]}^{1/3}$$E18

Inserting Eq. 7 into Eq. 18, one has

$$d_{\ast }^{\text{'}}={\left[\frac{{\rho }_s-\rho }{\rho }{\left(1-Y\right)}^2\frac{{\mathit{gd}}_{50}^3}{{\nu }^2}\right]}^{1/3}$$E19

or

$$\frac{d_{\ast }^{\prime }}{d_{\ast }}={\left(1-Y\right)}^{2/3}$$E20

Therefore, the empirical equation of Shields curve by Yalin and Silva [31] can be modified with the following form:

$$\tau {\text{'}}_{\ast }=0.13d_{\ast }^{\text{'}-0.392}exp\left(-0.015d_{\ast }^{\text{'}2}\right)+0.045\left[1-exp\left(-0.068d_{\ast }^{\text{'}}\right)\right]$$E21

For the fall velocity, many empirical equations are available in the literature. Julien [5] related c_d in Eq. 4 with the particle diameter d_* and obtained the following empirical equation:

$$\frac{\omega d_{50}}{\nu }=8\left[\sqrt{1+0.0139d_{\ast }^3}-1\right]$$E22

The incipient motion in uniform flows has been extensively investigated, but no one investigates the influence of vertical velocity on incipient motion, probably because this vertical flow may not large enough to induce discernible seepage, thus it is useful to estimate V_b using some measured parameters. The depth-average vertical velocity can be determined by,

$$V=U\frac{\mathit{dh}}{\mathit{dx}}$$E23

where U, the average streamwise velocity, and both U and dh/dx are measurable parameters, thus Eq. 23 is convenient to use.

The vertical velocity is jointly induced by either the groundwater or the surface variation, the joint effect can be assumed as the proportional V and the nominal seepage velocity V_s, i.e., λV + λ_sV_s, or:

$$V_b=\frac{\mathit{\lambda V}+{\lambda }_s{V}_s}{1-{\epsilon }_0}$$E24

where λ and λ_s are the coefficients to relate V_b with the mean vertical velocity V and nominal seepage velocity (V_s) defined by Darcy (V_s = ki, k = hydraulic conductivity, i = hydraulic gradient), ε₀ = porosity of granular materials.

Generally in laboratory flumes, the second term of Eq. 24 is negligible (i.e., V_s = 0), but in natural streams both the river flow and underground water flow can generate the velocity at the river bed, thus two terms co-exist in Eq. 24.

4. Re-analysis of data on the original shields diagram

To verify whether Eq. 21 is applicable to non-uniform flows, 329 data points are comprehensively compiled [13, 17, 23, 29, 32, 33, 34, 35, 36, 37, 38, 39, 40]. The hydraulic conditions of the used data are summarized in Table 1, and the experimental conditions are briefly outlined as follows:

Researchers	$d_{50}$(mm)	$S$	h (m)	U (m/s)	$u_{\ast }$(m)	${{\tau }_c}^{\text{'}}/{\tau }_c$	$\omega$ (m/s)	No. of data points	remark
Neil (1967)	5–29.1	0.01	0.03–0.192	0.28–0.35	0.029–0.165	0.685–1.18	0.13–0.62	59	uniform
Gaucher et al. (2010)	0.91–4.36	0.01	0.125–0.14	0.29–0.56	0.021–0.038	0.405–0.907	0.105–0.245	6	Uniform
Carling (1983)	62; 77	0; 0	0.213; 0.226	0.163; 0.124	0.141; 0.315	—	0.94; 1.05	2	uniform
Shvichenko & Pender (2000)	1.5–12	0.0019–0.0287	0.002–0.65	0.1–1.07	0.026–0.1157	—	0.14–0.41	21	uniform
White (1970)	0.016–2.2	0.02	0.02–0.07	0.0018–0.232	0.0062–0.045	—	0.00023–0.174	26	Non-uniform
Sarker & Hossain (2006)	0.64–1.02	0.00026–0.00063	0.089–0.214	0.3–0.59	0.019–0.024	—	0.084–0.1135	19	Non-uniform
Afzalimhr et al. (2007)	8	0.0075, 0.015	0.13–0.21	0.726–0.86	0.05–0.061	0.362–0.535	0.338	9	Non-uniform
Graf & Suszka (1987)	12.2,23.5	0.0075, 0.025	0.102–0.2	0.23–1.6	0.087–0.155	—	0.41–0.58	9	Non-uniform
Emadzadeh et al. (2010)	0.8, 1.3, 1.8	$\pm$0.7, $\pm$ 0.9, $\pm$ 1.25, $\pm$ 1.5	0.146–0.25	0.15–0.44	0.007–0.021	0.078–2.9	0.097–0.156	72	Non-uniform
Everts (1973)	0.09–1.79	0.005	0.0094–0.09	0.1312–0.38	0.018–0.043	0.39–1.79	0.007–0.156	35	Non-uniform
Liu & Chiew (2012)	0.9	0.01	0.12–0.14	0.28–0.35	0.0215	0.98–1.71	0.105	5	seepage
Cheng & Chiew (1999)	0.63–1.95	0.01	0.027–0.076	0.09–0.399	0.017–0.032	0.02–1.048	0.08–0.163	50	Seepage
Kavcar & Wright (2009)	0.16, 0.5, 1.2	0.01	0.23–0.29	0.23–0.412	0.013–0.022	0.67–1.84	0.019–0.124	16	seepage

Table 1.

Summary of experimental conditions by previous researchers.

Neil conducted his experiments in a flume 0.9 m wide and 5 m long by using sands with different particle sizes and densities [32]. Among the data sets, 11 data points are obviously above the Shields curve. White collected his data from a recirculating flume 6 m long and 0.3 m wide, uniform sediment was used with diameter between (0.016–2.2) mm [33]. The experimental datasets by Everts included 35 runs with size $d_{50}$ from 0.127 to 1.79 mm and specific gravity of 2.65, and 11 runs having $d_{50}$ from 0.09 to 0.18 and specific gravity of 4.7 [34]. Figure 4 shows that almost all his data points are located below the Shields’ prediction. Carling’s data [35] were collected from a narrow natural stream and in a broad stream. Graf and Suszka measured the critical shear stress in a flume 16.8 m long, 0.6 m wide and 0.8 m high, gravel sediment with uniform size was used [23]. Shvidchenko and Pender [36] used a flume to study the effect of relative depth on the incipient motion of coarse uniform sediments. Gaucher et al.’s [40] experiments were conducted in a horizontal, rectangular glass walled flume with dimensions of 6 m long, 0.5 m wide and 0.7 m deep, different types of non-cohesive materials were used ranged from $d_{50}$ = 0.91 to 4.36 mm. Cheng and Chiew [30] investigated the influence of upward seepage on the critical conditions of incipient motion, the experiments were conducted in a horizontal flume 7.6 m long, 0.21 m wide and 0.4 m deep, with particle sizes of $d_{50}=$ 0.63, 1.02 and 1.95 mm, and the seepage velocity (injection) was measured with a range between (0–0.0138) m/s. They found that the upward seepage reduces significantly the critical shear stress required by Shields curve. Kavcar and Wright [38] conducted experiments in a 7.5 m long, 0.6 m wide flume with both injection and suction seepage using sediment particle of $d_{50}$ =0.16, 0.5 and 1.2 mm and the observed value of seepage velocity, i.e. $V_s$ is range between (−0.0026–0.00223) m/s. Liu and Chiew [29] examined the critical shear stress for sediment with $d_{50}$ = 0.9 mm subject to downward seepage with velocity between (−0.00314–0) m/s. Their glass-sided flume was 30 m long, 0.7 m wide and 0.6 m deep, they observed that the upward seepage (injection) decreases the critical shear velocity while the downward seepage (suction) increases it.

Figure 4.

Measured critical shear stress versus d_* and its comparison with shields curve or Eq. 21 at Y = 0 (i.e., uniform flow) and 100% error band.

Nineteen flume experiments from Sarker and Hossain [37] are also included in Figure 4. They investigated the initiation of sediment motion under non-uniform sediment mixtures. Afzalimhr et al. [13] conducted experiments to investigate the effect of non-uniformity of flow on the critical shear stress in a channel (14 m long, 0.6 m width and 0.5 m depth), the sediment size of $d_{50}=8$ mm was used for their observation. Different from Lamb et al’s [24] prediction, their experimental data reveal that the value of critical shear stress is smaller than Shields’ prediction by at least 50%. Similarly, Emadzadeh et al. [39] conducted experiments in accelerating and decelerating flow conditions, his flume was 14 m long, 0.6 m wide and 0.6 m deep. The sediment size used were $d_{50}$ = 0.8 and 1.3, 1.8 mm for a total of 72 data sets. The decelerating/accelerating flows were obtained by adjusting negative and positive bed slope (±0.7%, ±0.9%, ±1.25% and ± 1.5%). It is found that the critical shear stress and Shields parameter for incipient motion in accelerating flow are higher than those predicted by Shields in uniform flow while their values in decelerating flow are considerably lower than that in accelerating flow.

These data mentioned are plotted in Figure 4, where the observed critical shear stress highly deviates from the standard Shields curve. All has been noticed and commented by many researchers [4, 24]. The consensus is that this discrepancy cannot be simply attributed to measurement errors or methodological bias. In Figure 4, the three lines are the Eq. 21 ($Y=0$) $\pm 100\%$ error band.

5. Dependence of critical shields stress on channel slope

Many researchers have noticed that high channel’s slope can cause the deviation of data from the Shields curve. For example Chiew and Parker [17] proposed that

$$\frac{{\tau }_{\ast }^{\text{'}}}{{\tau }_{\ast }}=cos\phi \left(1-\frac{tan\phi }{cos\theta }\right)$$E25

where ϕ = angle of streamwise bed slope, θ = angle of repose. Eq. 25 shows that the Shields number decreases with the increase of channel slope.

However, the formula given by Lamb et al. (2008) shows that the steep channel has a higher Shields number with the following form:

$${\tau }_{\ast }^{\text{'}}=exp\left[0.0249X^4+0.107X^3+0.199X^2+0.476X-3.57\right]$$E26

where X = 0.407ln(142 $S$), and the slope S is in the regime 10⁻⁴< $S$ <0.5.

Figure 5 demonstrates the comparison of the measured data from Table 1 and Eqs. 25 and 26. Obviously these equations do not agree the data points well. The measured τ_* could be largely different even the same type of sediment and channel slope are used. Therefore, the invalidity of Shields prediction cannot be simply explained by the dependence of channel slope, and there are some physics inside for the discrepancy.

Figure 5.

Dependence of critical shear stress on the channel slope.

6. Seepage on critical shields stress

Figure 5 demonstrates that for the same particle size in the same channel slope, the data points behave largely different, which cannot be explained by any existing theory. Beyond other factors, Eq. 24 shows that the scatter could be induced by either groundwater or the main flow’s non-uniformity, or both of them. The effect of seepage on the critical shear stress is discussed first, the experimental data [29, 30, 38] are showed in Figure 6.

Figure 6.

Comparison of measured and predicted critical shear stress subject to seepage.

The modified Shields number in Eq. 8 (i.e., that with seepage) will be the same as that obtained from the Shields curve if one uses both the apparent sediment density and the apparent critical shear stress (that with seepage), i.e.

$${\tau }_{\ast }^{\text{'}}=\frac{{\tau }_c^{\text{'}}}{\left({\rho }_s^{\text{'}}-\rho \right){\mathit{gd}}_{50}}=\frac{{\tau }_c}{\left({\rho }_s-\rho \right){\mathit{gd}}_{50}}$$E27

Using Eq. 7, one obtained the ratio of critical shear stresses with/without V_b in the following form:

$$\frac{{\tau }_c^{\text{'}}}{{\tau }_c}={\left(1-Y\right)}^2$$E28

Figure 6 shows the critical shear stress predicted by Eq. 28 and the empirical factor λ_s is found to be 8.5, the experimenters determined the critical shear stress without seepage using the Shields curve. The good agreement between the measured and predicted critical shear stress indicates that the introduction of apparent sediment density is acceptable.

Similarly, local scour by large vortices (e.g., scour holes around bridge piers) is not caused by higher velocity or higher boundary shear stress, but the upward velocity Y. The mechanism is similar to the helicopter whose rotor blades generate the upward velocity by “vortices”. Consequently, low water pressure induces the seepage or upward velocity, large particles like stones/helicopter can be lift. One can easily infer the relationship between the upward velocity and “vortices” in front of an electricity fan. Likewise, by observing how tornadoes damages large particles like cars, houses on surface, one can easily concluded that the upward velocity or lift force is the cause, by no means the shear force.

7. Effect of non-uniformity of flow on the critical shear stress

Figure 4 shows that the Shields’ curve could be totally invalid sometimes, these noticeable deviations imply that the non-uniformity of flow could affect the predictability of Shields curve, for example, Afzalimhr et al.’s [13] data points locate below the curve when the flow was decelerating, Emadzadeh et al.’s data points [39] were far from the Shields’ prediction, and his data points were obtained from both decelerating and accelerating flows. Hence, the large deviations from Shields curve shown in Figure 4 can be used to verify Eq. 24, i.e., the vertical velocity induced by flow’s non-uniformity is responsible for the invalidity of Shields curve.

To confirm whether the invalidity of Shields curve is caused by the non-uniformity of flow, the data without seepage in Table 1 are used, and the water depth variation $\mathit{dh}/\mathit{dx}$ is calculated using the following formula:

$$\frac{\mathit{dh}}{\mathit{dx}}=\frac{S-S_f}{1-U^2/\mathit{gh}}$$E29

where $\mathit{dh}/\mathit{dx}$ is the water depth’s variation, $S$ and $S_f$ are the bed and energy slopes, respectively. Manning coefficient ($n$) can be assessed using the Strickler’s formula:

$$n=\frac{d_{50}^{1/6}}{21.1}$$E30

The energy slope $S_f$ in Eq. 29 can be determined from the Manning equation using the hydraulic radius R, i.e.,

$${\mathrm{S}}_f=\frac{n^2{U}^2}{R^{4/3}}$$E31

In Table 1, the calculated $\mathit{dh}/\mathit{dx}$ could be either negative or positive and the data points in Figure 4 are included and replotted in Figure 7, where the data point is represented by the sign “+” if the obtained dh/dx is positive, otherwise the data point is marked by “-” for all negative $\mathit{dh}/\mathit{dx}$ cases. Figure 7 clearly shows that nearly all data points above the Shields curve have “-” signs, indicating the flows were accelerating, whilst almost all data points below the Shields curve have the sign of “+”, or decelerating. Therefore, the non-uniformity of flow can play an important role for the deviation of measured critical shear stress from the Shields curve. Figure 7 reveals that the presence of vertical velocity is one of the main causes responsible for the deviation of observed critical shear stress from the Shields curves for these data, the accelerating flow enhances particles’ stability, and decelerating flow enables sediment’s mobility. In Figure 7, the calculated positive $\mathit{dh}/\mathit{dx}$ ranges from 0.000237 to 0.0526 and − 0.024 to −0.00073. It remains necessary to investigate whether the higher $\mathit{dh}/\mathit{dx}$ has the higher deviation, and its analysis is shown below.

Figure 7.

The variation of water depth $\mathit{dh}/\mathit{dx}$ has different values based on the influence of vertical velocity on the initial motion, where (−0.024< $\mathbf{dh}/\mathbf{dx}$ <0.0526) for all data sets from Figure 3.

8. Modification of shields diagram

To examine whether data points without seepage shown in Figure 4 can be expressed by Eq. 21, we can analyze the datasets without artificial seepage or with negligible groundwater effects, only those data are analyzed in which V_b is caused by the non-uniformity of flow in the main flow. Therefore, Eqs. 23 and 24 can be simplified as follows:

$$Y=\frac{V_b}{\omega }=\frac{\mathit{\lambda U}}{\left(1-{\epsilon }_0\right)\omega }\frac{\mathit{dh}}{\mathit{dx}}$$E32

Experiments [13, 34, 39, 40] are analyzed first. They reported that their measured critical shear stress is lower than Shields’ prediction. Besides, the datasets [32, 39] are examined; they claimed that higher values of critical shear stress were observed.

In these studies, the experimental data sets from non-uniform flows are plotted in Figure 8 where the empirical factor λ is found to be 8.5 for both decelerating and accelerating flows. The comparison of the predicted and measured critical shear stress in Figure 8 shows that the agreement is reasonably good. Better agreement can be obtained if λ is calibrated as a function of sediment gradation and shapes, turbulence. Here, the assumption is that sediment particle size is uniform and can be represented by d₅₀.

Figure 8.

Comparison of experimental results on threshold condition without seepage with Eq. 32.

Figure 9 shows the comparison of measured and the predicted critical shear stress for the datasets [23, 33, 35, 36, 37]. Obviously, the observed critical shear stress largely deviates from the solid line, i.e., Shields curve (Y = V_s/ω = 0), all data points can be covered by Eq. 9 or 10 when the parameter Y is used. In other words, Figure 9 suggests that the scatter might be explained by variation of Y.

Figure 9.

Influence of vertical velocity on critical shear stress, the solid line is the original shields curve (or Y = 0) and other lines are calculated from Eq. 21 with different Y.

9. Discussion on slope’s influence

As mentioned, some researchers have found the dependence of the critical shear stress on the channel slope, but it is still an open question about the validity of Shields curve, especially when the bed slope is large, thus it is worthwhile to discuss this dependence.

This study reveals that the deviation from the Shields curve could be caused by the vertical velocity, the Shields curve is approximately valid only when the flow is uniform, when the vertical velocity is almost zero. As the true uniform is very rare in laboratory or nature, thus it is understandable why Shields curve is invalid to express most of observed critical shear stress. Hence, one needs to answer whether the dependence of τ_* on the channel slope is also caused by the flow’s acceleration.

Eqs. 23 and 24 show that in almost all cases, there always exists the vertical velocity caused groundwater and flow’s non-uniformity. Therefore, the widely observed dependence by Lamb et al. [24] may be also caused by the parameter Y (≠ 0). Obviously, they assumed that the data used in their analysis were collected from uniform flows, thus these data can be used to compare the data with the Shields curve, and conclusion of the slope-dependence can be drawn. It is useful to examine this assumption by checking whether Lamb et al’s data [24] are observed from uniform flows. Their data are listed in Table 2, in which only the laboratory data are included as their field data were certainly collected from non-uniform conditions. The last column of Table 2 shows the length of flumes, and from it one can see that almost half of the flumes were less than 10 m. Kirkgöz and Ardiçlioğlu [41] measured the minimum length to form a uniform flow in a flume and found that a channel should be longer than 10 m as there is a transition zone from non-uniform flow to uniform flow. Even for those data from flumes longer than 10 m, the flow still could be non-uniform also when the parameter $\mathit{dh}/\mathit{dx}$ is calculated using Eq. 6. Paola and Mohrig [42] suggest that uniform flow can only be assumed when the channel length is longer than h/S, if the water depth is 0.1 m, and slope is 1%₀, this means that the channel length should be longer than 100 m. Therefore, one can conclude that it is likely that the data were generated in non-uniform conditions, this may lead to the different interpretation of the dependence of critical shear stress on the channel slope. Figure 3 shows the accelerating flows in steep channels, thus it is likely that the downward velocity increases particles’ stability.

Researchers	$d_{50}$ (mm)	${\tau }_{\ast }$	$R_{\ast }$	Flume length (m)
Neil (1967)	6.2, 8.5, 10.6, 20, 23.8, 29.1, 5, 16, 6.4	0.04–0.06	184.3–4800	5
Paintal (1971)	7.95, 2.5	0.05, 0.05	638, 112	15
Everts (1973)	3.57, 1.79, 0.895, 0.508, 0.395, 0.254, 0.127, 0.18, 0.09	0.018–0.07	1.3–162	16.8
Ashida & Bayazit (1973)	22.5, 12, 6.4	0.0386–0.1178		20
Fernandez Luque & Van Beek (1976)	0.9, 1.5, 1.8, 3.3	0.021–0.047	12–127	8
Ikeda (1982)	0.42; 1.3	0.02; 0.047	8.7; 72	4
Graf & Suszka (1987)	12.2, 23.5	0.05–0.07	800–5000	16.8
Wilcock (1987)	1.83; 1.83; 0.67; 5.28	0.03; 0.036; 0.023; 0.037	61; 12; 332; 115	23
Wilcock & Mcardell (1993)	5.3	0.02	219	7.9
Shvidchenko & Pender (2000)	1.5, 2.4, 3.4, 4.5, 5.65, 7.15, 9, 12	0.025–0.065	40–2000	6.5

Table 2.

Previously reported data selected from lamb et al. (2008).

Chiew and Parker’s data [17] is used as an example, their observation is opposite to Lamb et al’s prediction [24], they also found the dependence of critical shear stress on the channel slope based on their own data. In their experiments, the channel slope was specially adjusted from -10^o to 31^o, their channel lengths used were 4 m and 2 m only. Obviously, their experiments were conducted in the non-uniform flow conditions as the 2 ∼ 3 m length is too short to form a uniform flow. In other words, both conclusions drawn by Lamb et al. and Chiew and Parker [17, 24] are not very convincing as they did not check the parameter of $\mathit{dh}/\mathit{dx}$, and the data they used may be generated from non-uniform conditions.

10. Conclusions

This paper investigates why the observed critical shear stress widely deviates from the Shields curve, its discrepancy or validity could be caused by many factors like sediment shapes, gradation, measurement errors, turbulence and channel-bed slopes. However, this study reveals that the vertical motion also plays an important role, and the vertical velocity could be induced by non-uniformity of flow and seepage turbulence alike. After re-examining 329 data points from the literature, the following conclusions can be drawn:

1. The upward velocity increases sediment mobility and downward velocity increases sediment stability. The mobility or stability can be equivalently expressed by its apparent sediment density which is able to eliminate the effect of vertical velocity as shown in Eq. 7. This shifts a dynamic problem into a simplified static problem.

2. There exists vertical velocity on the channel bed and this vertical velocity could be induced by seepage or non-uniformity of flow, similar to the secondary currents, the small vertical velocity’s influence on sediment incipient should not be underestimated. The joint effect is expressed by Eq. 24. For non-uniform flow, the sediment tends to move in decelerating flows, but it becomes more difficult to move in accelerating flows.

3. The Shields curve is valid only when the flow is nearly uniform, but a general Shields curve can be obtained by introducing the apparent sediment density, thus the modified Shields curve could be extended to express complex flows, this modified relationship for critical shear stress has been established.

4. A new parameter Y can be used to express the influence of non-uniformity of flow or seepage, this parameter should be included in the models of sediment transport. According to available experimental data in the incipient motion in non-uniform flows or in the seepage cases, good agreements between the measured and predicted values can be achieved if Y is included in the existing model, but more research is needed to determine the coefficients λ_s and λ in Eq. 24, they could be a function of sediment gradation and shapes, and turbulence.

All in all, high horizontal motion can make a plane (a big particle) to fly, high vertical velocity can also make the same particle called helicopter to fly. Two mechanisms are totally different. It is wrong to ascribe all sediment transport phenomena to the horizontal motion only, without considering the vertical motion.

Notations

b = channel widthC_d = drag coefficient;d₅₀ = median size of sediment particles;F_vb = force induced by the vertical velocity;g = gravitational acceleration;h = water depth;i = hydraulic gradient;k = hydraulic conductivity;n = Manning coefficient;Q = discharge;R_* = Reynolds number;S_f = energy slope.U = mean velocity;u_* = shear velocity;u_*c = critical shear velocity (τ_c = ρu_*c²);$\overline{u}$ and $\overline{v}$ = time-averaged velocity in the streamwise and vertical directions;${\overline{u}}_h$, ${\overline{v}}_h$ = horizontal and vertical velocities at the surface;V = vertical velocity;V_b = vertical velocity at the bed;V_s = nominal seepage velocity at the bed;X = 0.407ln(142S);y = distance normal to the wall;Y = V_b/ωε₀ = porosity of granular materialsθ = angle of repose.λ and λ_s = coefficients;ν = kinematic viscosity;ρ = fluid density;ρ_s = sediment density;ρ_s’ = apparent density of sediment;τ = boundary shear stress;τ_c = critical boundary shear stress;τ_* = Shields number;τ_*’ = modified Shields number subject to vertical velocity;ϕ = angle of streamwise bed slope,ω = particle fall velocity;ω’ = net falling velocity subject to vertical velocity.