Hydrodynamic instability of free river bars

Cite as: Phys. Fluids 33, 045105 (2021); https://doi.org/10.1063/5.0045530

Submitted: 27 January 2021 . Accepted: 10 March 2021 . Published Online: 02 April 2021

Rajesh Kumar Mahato, Sk Zeeshan Ali, and Subhasish Dey

Phys. Fluids 33, 045105 (2021); https://doi.org/10.1063/5.0045530 33, 045105

© 2021 Author(s).
 Physics of Fluids ARTICLE scitation.org/journal/phf

Hydrodynamic instability of free river bars

Cite as: Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530
Submitted: 27 January 2021 . Accepted: 10 March 2021 .
Published Online: 2 April 2021

Rajesh Kumar Mahato,1,a) Sk Zeeshan Ali,2,b) and Subhasish Dey1,3,c)

AFFILIATIONS
1
Department of Civil Engineering, Indian Institute of Technology Kharagpur, West Bengal 721302, India
2
Department of Civil Engineering, Indian Institute of Technology Hyderabad, Telangana 502285, India
3
Department of Hydraulic Engineering, State Key Laboratory of Hydro-Science and Engineering, Tsinghua University,
Beijing 100084, China

a)
[email protected]
b)
[email protected]
c)
Author to whom correspondence should be addressed: [email protected]

ABSTRACT
In this paper, we explore the hydrodynamic instability of free river bars driven by a weakly varying turbulent flow in a straight alluvial
channel with erodible bed and non-erodible banks. We employ linear stability analysis in the framework of depth-averaged formulations for
the hydrodynamics and the sediment transport. A significant fraction of the sediment flux is considered to be in suspension. The analysis is
performed for the alternate pattern of river bars at the leading order followed by the next order, covering the effects of flow regime. We find
that the unstable region bounded by a marginal stability curve depends significantly on the shear Reynolds number, which demarcates differ-
ent flow regimes, and the Shields number and the relative roughness (particle size to flow depth ratio). The results at the next order stabilize
the bars with longer wavenumbers. The variations of threshold aspect ratio with Shields number and relative roughness are studied for differ-
ent flow regimes. In addition, for a given Shields number and relative roughness, the diagram of threshold aspect ratio vs shear Reynolds
number is explained. Unlike the conventional theories of bar instability, the analysis reveals limiting values of Shields number and relative
roughness beyond which the theoretical results at the next order produce infeasible regions of instability. The limiting values of Shields num-
ber and relative roughness appear to reduce, as the shear Reynolds number increases.
Published under license by AIP Publishing. https://doi.org/10.1063/5.0045530

I. INTRODUCTION Formation of river bars has been studied extensively by means of

Riverbed erosion and deposition associated with the interaction theoretical analyses,2,4–12 laboratory experiments,3,13–17 field observa-
between a flowing fluid and an erodible bed form striking sedimentary tions,18,19 and numerical simulations.20,21 Recently, Crosato and
patterns over widely varying length scales.1 River bars are frequently Mosselman22 reported an in-depth review on this topic. Instability of
observed in both gravel and sand-bed rivers. They are large-scale sedi- river bars is conventionally analyzed using linear and nonlinear
mentary patterns with a wavelength of the order of several channel approaches.1,9,12 A linear stability analysis predicts the unstable region
widths.2,3 Bars are broadly grouped into two classes: free and forced on a plane formed by the channel aspect ratio (channel width to flow
bars.1,4 Instability of an erodible bed in a straight channel leads to the depth ratio) and the bar wavenumber. On the other hand, a nonlinear
formation of free bars.5 On the other hand, forcing effects (channel stability analysis quantifies the bar amplitude.
curvature and width variation) trigger the formation of forced bars.1 Sediment transport plays a crucial role toward the instability of
This study mainly puts into focus the hydrodynamic instability of free river bars. Some studies based on two-dimensional (2D) assumptions
bars, as they are the most common depositional features in a straight considered only the bedload flux in the linear stability analysis.2,9,10 It
fluvial system. The formation of river bars limits navigation, increases was found that the channel bed becomes unstable beyond a threshold
flood risk, and is linked with scour at banks and bridge piers.3,6 aspect ratio.9 Tubino et al.5 studied the instability of free river bars by
Moreover, bars play a key role in governing the dynamics of meander- means of a 3D framework, considering the effects of sediment suspen-
ing and braided rivers.1 Therefore, from an engineering perspective, sion. They concluded that the inclusion of sediment suspension leads
instability of river bars is a key to understanding their precise to a decrease in the threshold aspect ratio and an increase in the bar
formation. wavelength. Similar findings were reported by Federici and Seminara11

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-1

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

and Bertagni and Camporeale.12 In a slowly varying flow, Federici and respectively. The origin of the coordinate system lies on the channel
Seminara11 studied the formation of river bars applying the depth- centerline. Figure 1(b) displays an exemplary cross-sectional view of
averaged formulations for the hydrodynamics and the sediment the channel at a cross section S1–S2. In Fig. 1(b), D is the local flow
transport. Considering only the bedload flux, Federici and Seminara4 depth and H is the elevation of the free surface of flow from a horizon-
studied the convective nature of bar instability by means of branch point tal reference level. The broken line in Fig. 1(b) characterizes the undis-
singularities of the dispersion relation. Solving the nonlinear problem turbed channel cross section.
numerically, they concluded that groups of bars, emanating from either The depth-averaged momentum and continuity equations can be
an arbitrary distributed or a local bed topography perturbation, grow obtained by performing depth-averaging of the time-averaged
and propagate downstream, keeping the source area intact. Afterward, momentum and continuity equations, in addition to the appropriate
Federici and Seminara11 recognized that the inclusion of sediment sus- boundary conditions at the free surface and the bed. The time-
pension does not change the nature of bar instability. averaged pressure intensity is assumed to follow the hydrostatic law.
The formation of alternate bars was found to be influenced by As the wavelength of river bars is of the order of a few channel widths,
the flow unsteadiness. Tubino6 predicted the bar amplitude in an the use of depth-averaged hydrodynamic equations is justified. In
unsteady flow, performing a weakly nonlinear analysis close to the addition, the morphological timescale is much larger than the flow
neighborhood of the threshold condition for the bar formation. In timescale. Therefore, we consider the flow field to be a quasi-steady.
addition, Hall10 presented a linear stability analysis of alternate bars in Let the components of depth-averaged velocity, bed shear stress, bed-
an unsteady flow. It was found that the flow unsteadiness leads to a load flux, and suspended load flux in (x, y) be (U, V), (Tx, Ty), (Qbx,
reduction in the threshold aspect ratio and the bar wavelength. Qby), and (Qsx, Qsy), respectively.
Lanzoni and Tubino2 studied the effects of sediment heterogeneity on The depth-averaged momentum equations are expressed as23
the formation of alternate bars. They found that the growth rate, wave- ^ ^ ^ ^
length, and propagation speed of bars dampen owing to the presence ^ @U þ V
U ^ @ U þ @ H þ b T x ¼ 0; (1)
of sediment heterogeneity. @^x @^y @^x D^
The state-of-the-art of bar instability is quite mature as far as the ^ ^ ^ ^
U ^ @ V þ @ H þ b T y ¼ 0:
^ @V þ V (2)
flow is steady and the primary mode of sediment transport is the bed- @^x @^y @^y ^
D
load transport. It is well-known that the channel aspect ratio is the key
parameter controlling the formation of alternate bars. Earlier studies The depth-averaged continuity equation of the fluid phase is
reported that the sediment suspension plays a destabilizing role by expressed as23
reducing the threshold aspect ratio. This conclusion was drawn with
@ ^^ @ ^^
respect to the conventional treatment, where the sediment transport ðD U Þ þ ðD V Þ ¼ 0: (3)
occurs mainly as a bedload transport.5,11,12 However, this destabilizing @^x @^y
influence is yet to be explored in detail. Recently, Ali and Dey7,8 The continuity equation of the solid phase, given by Exner,24 is
reported a comprehensive analysis of the instability of large-scale river expressed as
bars considering the effects of sediment suspension and flow regimes.    
@ 2^ @Uby 1 @Usx @Usy
However, little is so far known about the effects of sediment suspen- ^ þ c @Ubx þ
ðF H  DÞ þ þ ¼ 0:
sion and flow regime on the bar instability at various orders of @^t @^x @^y 1  p @^x @^y
approximation. (4)
This study aims at exploring the effects of sediment suspension
and flow regime on the bar instability. We perform a linear stability In Eqs. (1)–(4), the following dimensionless variables are
analysis considering a weakly varying turbulent flow in a straight introduced:
channel with erodible bed and non-erodible banks. For the hydrody- ðx; yÞ ^ D
namics and the sediment transport, a depth-averaged formulation is ð^x ;^y Þ ¼ ; D¼ ; H ^ ¼ H ; F ¼ U0 ;
B D0 F 2 D0 ðgD0 Þ0:5
used. The analysis is capable to capture the bar instability over a wide
range of shear Reynolds numbers, Shields numbers and relative rough-
B
b ¼ ; ðU ^ ;V ^ Þ ¼ ðU; VÞ ; ^t ¼ tU0 ; ðT ^ y Þ ¼ ðTx ; Ty Þ ;
^ x; T
ness. It is revealed that the threshold aspect ratio, being the controlling D0 U0 B qf U02
ðQbx ; Qby Þ ðQsx ; Qsy Þ (5)
parameter for the bar formation, significantly varies with the shear ðUbx ; Uby Þ ¼ 0:5 ; ðUsx ; Usy Þ ¼ ;
Reynolds number, Shields number, and relative roughness. ½ðs  1Þgd 3 U0 D0
 0:5
The paper is arranged as follows: The mathematical formulation ðs  1Þgd3
is presented in Sec. II. The linear stability analysis is performed in Sec. and c ¼ ;
ð1  pÞU0 D0
III. The results and discussion are furnished in Sec. IV. Finally, conclu-
sions are drawn in Sec. V. where F is the undisturbed flow Froude number, U0 is the undisturbed
average flow velocity, g is the gravitational acceleration, b is the chan-
II. MATHEMATICAL FORMULATION nel aspect ratio (channel half-width B to undisturbed flow depth D0
The schematic of free river bars is shown in Fig. 1. The bar wave- ratio), t is the time, s is the relative density (¼ qs/qf), qs is the mass
length is represented by k [Fig. 1(a)]. We consider a 2D turbulent flow density of sediment particles (¼ 2650 kg m3 for natural quartz sand),
of an incompressible fluid in a straight alluvial channel of constant qf is the mass density of fluid (¼ 1000 kg m3 for water), p is the sedi-
width 2B. The channel has an erodible bed and non-erodible banks. In ment porosity ( 0.3), d is the median sediment size, and c is the ratio
Fig. 1(a), x and y represent the streamwise and spanwise directions, of the scale of sediment flux to that of flow flux.

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-2

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

FIG. 1. Definition sketch of the physical system: (a) plan view (flow directed from left to right) and (b) cross-sectional view.

To solve the system of equations, closure relationships for the Rijn28 (see Appendix A). The flow Froude number F and flow
bed shear stress and the sediment flux vectors are required. The com- Reynolds number R can be expressed as follows:
ponents of the bed shear stress vector, as a function of dynamic pres-
sure, are expressed as " #0:5
8Hðs  1Þd^ ^  80:5
4DR
ðT ^ y Þ ¼ f ðU
^ x; T ^;V ^2 þV
^ ÞðU ^ 2 Þ0:5 ; (6) F¼ and R ¼ ; (11)
8 f ad^ f
where f is the Darcy–Weisbach friction factor, which describes the
where d^ is the relative roughness (¼ d/D0).
resistance due to frictional effects. The f can be determined from the
The bedload transport is defined as the transport of particles in
Colebrook–White equation. It is expressed as25
2 rolling, sliding, and saltating modes.27 When the fluid-induced bed
!1:1 3
^ shear stress just exceeds its threshold value, the particles are trans-
1 2:51 k
¼ 0:86 ln 4 0:5 þ 5;
s
(7) ported in sliding and/or rolling mode. With an increase in bed shear
f 0:5 Rf 14:8D ^
stress, the bedload transport occurs in a saltating mode.27,29 As the bed
where R is the flow Reynolds number (¼ 4UD/t), t is the coefficient shear stress increases further, finer particles are lifted up into suspen-
of kinematic viscosity of fluid, ^k s ¼ ks/D0, and ks is the bed roughness sion owing to the upward diffusion of near-bed turbulence.27 In this
height. The ks is expressed as a linear function of particle size, that is, situation, both the bedload and suspended load transport exist. The
ks ¼ ad, where a is a factor. Following the experimental observations direction of bedload transport deviates from that of bed shear
of Engelund and Hansen,26 we set a ¼ 2.5. It is worth mentioning that
stress.30–32 The components of bedload flux vector are expressed as
the Colebrook–White equation facilitates the estimation of friction fac-
tor over a wide range of flow regimes. The flow regimes are distin- ðUbx ; Uby Þ ¼ Uðcos v; sin vÞ; (12)
guished by the values of shear Reynolds number R as hydraulically
smooth (R  3), transitional (3 < R < 70), and rough (R  70) where U is the bedload flux intensity and v is the angle between the
flow regimes. The R is expressed as27 resultant bedload flux and the longitudinal direction. The v can be
obtained by imposing a dynamic equilibrium on a spherical particle
u ks moving along a plane tangential to the bed. The v is expressed as30
R ¼ ; (8)
t
^
V r @
where u is the shear velocity. The Shields number H, represen- sin v ¼ 2  ^  DÞ;
 ðF 2 H ^ (13)
tative of the dimensionless fluid-induced bed shear stress, is ^ 2 Þ0:5
^ þV
ðU bH0:5 @^y
expressed as27 where r is an empirical constant. The r varies over a certain range, as
u2 reported in the literature.33 However, following the study of Talmon
H¼ : (9) et al.,34 we consider r ¼ 0.56 in this study. For the bedload flux inten-
ðs  1Þgd
sity U, various power laws are available in the literature. Here, we
The coupling of R with H produces employ the empirical formula of Meyer-Peter and M€ uller.35 It is

R ¼ aðHD3 Þ0:5 ; (10) U ¼ 8ðH  Hc Þ1:5 ; (14)

2 1/3
where D is the particle parameter {¼ d[(s – 1) g/t ] }. As the sus- where Hc is the threshold Shields number. Note that Ali and Dey36
pended load transport is considered to be dominant in the analysis, obtained the above empirical exponent from the phenomenological
the Shields number H always exceeds the threshold of sediment sus- theory of turbulence. For the estimation of Hc, the force system acting
pension Hs. It implies that for H < Hs, the sediment suspension on a particle needs to be analyzed. Interested readers may refer to the
ceases. For the estimation of Hs, we use the empirical formula of van work of Ali and Dey37,38 and Dey and Ali39–41 in this regard. However,

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-3

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

to simplify the analysis in the present context, we use the following nonequilibrium effects originating from the spatial variation of the
empirical relationship of Hc(D), given by Wu and Wang:42 flow field. Therefore, we write
Hc ¼ KDn : (15) ðU ^ 1; H
^ 1; V ^ 1; D
^ 1 Þ ¼ ðU
^ 10 ; V
^ 10 ; H
^ 10 ; D
^ 10 Þ
In Eq. (15), the complete set (K, n) ¼ (0.126, 0.44) for D < 1.5, þ dðU ^ 11 ; V^ 11 ; H ^ 11 Þ þ Oðd2 Þ:
^ 11 ; D (23)
(0.131, 0.55) for 1.5  D < 10, (0.0685, 0.27) for 10  D < 20, The components of bed shear stress and bedload flux vectors are
(0.0173, 0.19) for 20  D < 40, (0.0115, 0.3) for 40  D < 150, and expressed as
(0.052, 0) for D  150.  
T ^ x11 ¼ f0 s1 ðU
^ x10 þ dT ^ 10 þ dU^ 11 Þ þ s2 ðD^ 10 þ dD ^ 11 Þ ; (24)
Following the work of Federici and Seminara11 and Bertagni and 8
Camporeale,12 we consider a 2D extension of the depth-averaged
model of Bolla Pittaluga and Seminara43 for the estimation of sus- T^ y10 þ dT^ y11 ¼ f0 ðV^ 10 þ dV ^ 11 Þ; (25)
8
pended load flux vector. Considering a slowly varying flow, Bolla  
Ubx10 þ dUbx11 ¼ U0 s3 ðU ^ 10 þ dU ^ 11 Þ þ s4 ðD^ 10 þ dD ^ 11 Þ ; (26)
Pittaluga and Seminara43 performed an asymptotic expansion of the (
exact solution of the advection-diffusion equation of sediment concen-
Uby10 þ dUby11 ¼ U0 ðV ^ 10 þ dV ^ 11 Þ  r  @
tration. In the perturbation approach, the effects of advective-cum- bH0:5 @^y
0
unsteadiness are considered to be smaller than those of gravitational )
settling and turbulent diffusion. Mathematically, it is expressed as  

F 2 ðH
^ 10 þ dH ^ 11 Þ  ðD^ 10 þ dD ^ 11 Þ ; (27)
U0 D0
dk ¼ 1; (16) where f0, U0, and H0 are the Darcy–Weisbach friction factor, bedload
Ws k
function, and Shields number, respectively, corresponding to the undis-
where Ws is the settling velocity of particles (see Appendix A) and k is turbed state. The coefficients si (i ¼ 1 to 4) are given in Appendix C. To
the scale of longitudinal variation of the flow field ( bar wavelength). perturb the suspended load flux, the w0 and w1 are expanded as follows:
As the dk depends on the bar wavelength, another small parameter d  
can be defined as11,12 w0 ¼ w00 f1 þ e w010 þ dw011 þ Oðd2 Þ g; (28)
k U0 w1 ¼ ew110 þ OðedÞ; (29)
d ¼ dk ¼ : (17) where w00 is the function w0 at the undisturbed state. The quantities
B bWs
w010, w011, and w110 are given in Appendix C.
The suspended load flux vector can be expressed as43 Substituting Eqs. (22)–(29) into Eqs. (1)–(4), we obtain the fol-
lowing linear differential problem at the leading order O(e):
^ U
ðUsx ; Usy Þ ¼ Dð ^;V
^ Þw; (18)
@U^ 10 @ H ^ 10 bf0  
where w is a function that can be expended in powers of d as follows: þ þ ^ 10 þ ðs2  1ÞD
s1 U ^ 10 ¼ 0; (30)
@^x @^x 8
w ¼ w0 þ dw1 þ Oðd2 Þ: (19) @V ^ 10 @ H ^ 10 bf0
þ þ V^ 10 ¼ 0; (31)
@^x @^y 8
In the above, w0 is the function for the uniform flow condition and w1
describes the O(d) correction to the function w. The correction at @U^ 10 @ V ^ 10 @ D ^ 10
þ þ ¼ 0; (32)
O(d) considers weak nonequilibrium effects owing to the spatial varia- @^x @^y @^x
tion of the flow field. The functions w0 and w1 are expressed as43 "
^ ^ ^
  @ 2^
ðF H 10  D^ 10 Þ þ cU0 s3 @ U 10 þ s4 @ D 10 þ @ V 10
w0 ¼ C0 K0 ; w1 ¼ K1 D ^ U ^ @C0 þ V
^ @C0 ; (20) @^t @^x @^x @^y
@^x @^y 2

r @
  2 ðF 2 H ^ 10  D ^ 10 Þ
where C0 is the depth-averaged concentration at the leading order, and bH0:50 @^y
K0 and K1 are the functions of relevant physical parameters (see Appendix B). !
w00 @ U ^ 10 @ D ^ 10 @ V ^ 10 @w010
Finally, the channel banks are considered to be impervious to the þ þ þ þ ¼ 0: (33)
flow and sediment fluxes. Hence, we write 1  p @^x @^x @^y @^x
^ j^y ¼61 ¼ Uby j^y ¼61 ¼ Usy j^y ¼61 ¼ 0:
V (21) Similarly, the linear differential problem at the next order O(ed)
is expressed as

III. LINEAR STABILITY ANALYSIS ^ 11 bf0 

^ 11 @ H
@U 
þ þ ^ 11 þ ðs2  1ÞD
s1 U ^ 11 ¼ 0; (34)
In the linear stability analysis, the primitive variables are @^x @^x 8
expanded as follows: @V^ 11 @ H ^ 11 bf0
þ þ V^ 11 ¼ 0; (35)
^;V
^;H
^ ; DÞ
^ ¼ ð1; 0; H
^ 0 ; 1Þ þ eðU
^ 1; V
^ 1; H
^ 1; D
^ 1 Þ; @^x @^y 8
ðU (22)
@U^ 11 @ V ^ 11 @ D ^ 11
where e is a small parameter. The perturbations need to be further þ þ ¼ 0; (36)
@^x @^y @^x
expanded in powers of d in order to account for the weak

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-4

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

"
^ ^ ^ In the above equation, t1 and t2 are given in Appendix C.
@ 2^
ðF H 11  D ^ 11 Þ þ cU0 s3 @ U 11 þ s4 @ D 11 þ @ V 11 At O(ed), a linear algebraic nonhomogeneous system is obtained as
@^t @^x @^x @^y
  
r @ 2
1 @U^ 11 AV ¼ C; (45)
2^ ^
 0:5  2 ðF H 11  D 11 Þ þ w00
bH0 @^y 1p @^x where V ¼ (u11, v11, h11, d11) , C ¼ (0, 0, 0, c41) , c41 ¼ ix1(F2h10
T T
  2
 d10) þ ^k (g1u10 þ g2d10), and g1 and g2 are given in Appendix C.
^ ^
@ D 11 @ V 11 @w011 @w
þ þ þ þ 110 ¼ 0: (37) At O(e), a dispersion relationship emerging from the solvability
@^x @^y @^x @^x
of Eq. (43) takes the form of
The solution of the problem is sought in the form of normal   
a11 a23 a32
modes. Therefore, we write ix0 F 2 ða11 a34  a14 a31 Þ  a13 a31 þ
  a22
^ 1j ; V
ðU ^ 1j ; H
^ 1j ; D
^ 1j Þjm¼odd ¼ u1j Sm ð^y Þ; v1j Cm ð^y Þ; h1j Sm ð^y Þ; d1j Sm ð^y Þ   
a41 a11 a23 a32
h i ¼ a13 ða11 a34  a14 a31 Þ þ a14 a13 a31 þ

exp ið^k^x  x^t Þ þ c:c:; ð38Þ a11 a22
!
 2
a23 a42 F crU0 p m 2 2
^ 1j ; V
ðU ^ 1j ; H
^ 1j ; D
^ 1j Þjm¼even ¼ u1j Cm ð^y Þ;v1j Sm ð^y Þ;h1j Cm ð^y Þ; þ ða11 a34  a14 a31 Þ 
h i a22 4bH0:5
0
d1j Cm ð^y Þexp ið^k^x x^t Þ þc:c:; (39)  "
a11 a23 a32 crU0 p2 m2
where m is the transverse Fourier mode, j is either 0 or 1 depending on  a13 a31 þ i^kcU0 s4 
a22 4bH0:50
the order of approximation, i is the imaginary unit [¼ (1)0.5], ^k is the 
w
dimensionless wavenumber (¼ 2pB/k), x is a complex quantity, whose þi^kð1 þ t2 Þ 00
: (46)
real and imaginary parts represent the growth rate and the dimensionless 1p
frequency of perturbations, respectively, when multiplied with –i, and We use Eq. (45) to evaluate x1. As the homogeneous part of Eq.
c.c. stands for the complex conjugate. The Sm and Cm are expressed as (45) admits a nontrivial solution, a solvability condition is required.
  The solvability condition suggests that det(A) ¼ 0 after replacing the
p last column of A by the column matrix C.9 Imposition of the solvabil-
Sm ð^y Þ ¼ sin m^y ; (40)
2 ity condition on Eq. (45) yields
 
p   
Cm ð^y Þ ¼ cos m^y : (41) a11 a23 a32
2 ix1 F 2 ða11 a34  a14 a31 Þ  a13 a31 þ
a22
The complex quantity x is expanded in powers of d as (   
2 a11 a23 a32 g1
¼ ^k g2 a13 a31 þ  a13 ða11 a34  a14 a31 Þ
x ¼ x0 þ dx1 þ Oðd2 Þ; (42) a22 a11
 #)
where x0 is the value of x at the leading order and x1 is the O(d) cor- a11 a23 a32
rection to the x. þa14 a13 a31 þ : (47)
a22
When Eqs. (38)–(42) are substituted into the linear differential
problem at O(e), a linear algebraic homogeneous system is obtained. It is worth mentioning that by substituting x1 into Eq. (42), the
In matrix form, this can be expressed as marginal stability (MS) curve correct to O(ed) can be obtained.
The MS curves at O(e) and O(ed) are defined as Re(–ix0) ¼ 0 and
AU ¼ 0; (43)
Re(–ix) ¼ 0, respectively, where Re represents the real part of the
4
4 T complex function.
where A ¼ aij 2 R and U ¼ (u10, v10, h10, d10) . The matrix coeffi-
cients aij are as follows:
IV. RESULTS AND DISCUSSION
bf0 In this section, computational results are presented for the alter-
a11 ¼ i^k þ s1 ; a12 ¼ a21 ¼ a24 ¼ a33 ¼ 0;
8 nate pattern of river bars (m ¼ 1). In order to cover various flow
bf0
a13 ¼ a31 ¼ a34 ¼ i^k; a14 ¼ ðs2  1Þ; regimes, for instance, smooth (R  3), transitional (3 < R < 70), and
8
bf0 pm pm rough (R  70) flow regimes, we consider R ¼ 1 in a smooth flow,
a22 ¼ i^k þ ; a23 ¼ ð1Þ1þm ; a32 ¼ ð1Þm ; R ¼ 3 at the extremity of a smooth flow, R ¼ 12 and 30 (two charac-
8  2  2
w00 teristic values) in transitional flows, and R ¼ 100 in a rough flow. It
a41 ¼ cU0 s3 þ ð1 þ t1 Þ i^k;
1 p  (44) might appear to the readers that the resolution of the shear Reynolds
pm w number selection is limited, as only five characteristic values are con-
a42 ¼ ð1Þm cU0 þ 00 ;
" 2 1p sidered, and therefore the corresponding simulation results may not
 2 #
crU 0 pm capture the entire hydrodynamic instability. However, this is not the
a43 ¼ F 2 ix0 þ ;
bH0:5 2 case. It may be emphasized here that the selected shear Reynolds num-
 0 2 bers correspond to a specific flow regime and therefore, they offer an
crU0 pm w
a44 ¼ ix0 þ i^kcU0 s4  þ i^kð1 þ t2 Þ 00 : idea of how the stability pattern evolves with an increase in shear
bH0:5
0 2 1 p
Reynolds number. In the subsequent discussion, it is shown that the

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-5

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

above-mentioned consideration of shear Reynolds numbers is justified

in order to identify the stability pattern for an intermediate shear
Reynolds number. The MS curves, for a given shear Reynolds number
R, are presented on the b(^k) plane for different values of Shields num-
ber H and relative roughness d^ (Figs. 2–5). A large value of H signifies
a considerable amount of sediment suspension, whereas a large value
of d^ for a given particle size suggests a shallow flow depth. A typical
MS curve demarcates the stable (exterior of the curve) and the unsta-
ble (interior of the curve) regions. In a stable (unstable) region, the
perturbations decay (amplify) with time.
The MS curves corresponding to a smooth flow regime (R ¼ 1)
are furnished in Fig. 2 on a semi-logarithmic scale. The solid and the
dashed lines represent the MS curves at O(e) and O(ed), respectively
(also in succeeding Figs. 3–5). At O(e), the MS curves (solid lines) on
the b(^k) plane appear to be U-shaped, embracing the region of insta-
^ there exists a single MS curve with
bility. For a given triplet (R, H, d),
local minima (^k c , bc), which correspond to db/d ^k ¼ 0. This suggests
that for b < bc, the formation of alternate bars is inhibited. Beyond
the threshold aspect ratio bc, the alternate bars are formed over a wide
range of wavenumbers.
Figure 2(a) shows the MS curves for a given relative roughness d^
and different Shields numbers H. As the shear Reynolds number is
kept constant, an increase in Shields number corresponds to a decrease
in particle size [Fig. 2(a)]. This is because the particle parameter fol-
lows an inverse scaling with the Shields number as D / H1/3, in
accord with Eq. (10).
At O(e), the unstable region bounded by an MS curve tends to
reduce, as the Shields number increases [Fig. 2(a)]. It is attributed to
the presence of considerable sediment suspension that imparts a stabi-
lizing effect. The stabilizing effect is more encouraging for a longer
wavenumber. This reveals that even for the same aspect ratio, a small
Shields number (large particle size) may favor the formation of alter- FIG. 2. The MS curves in a smooth flow regime, R ¼ 1 [solid lines: O(e) and
nate bars with longer wavenumbers, while a large Shields number dashed lines: O(ed)]: (a) ^
d ¼ 0.005 and different values of H ¼ 0.15, 0.25, and
(small particle size) may not favor. Similar conclusion can be drawn 0.62 and (b) H ¼ 0.25 and different values of ^d ¼ 0.001, 0.005, and 0.023. The
for the MS curves at O(ed) [see dashed lines in Fig. 2(a)]. However, it curves of bc vs H and bc vs ^
d are shown in the insets.
^ the MS curve at O(ed) holds
appears that for a given triplet (R, H, d),
a smaller unstable region than that at O(e). Therefore, the longer
wavenumbers at O(ed) become more stable than those at O(e). This Figure 2(b) shows the MS curves for a given Shields number H
leads to an enhanced stable region, an observation being in conformity and different values of relative roughness d.^ As the Shields number is
with Federici and Seminara.11 Note that for shorter wavenumbers (^k kept constant in Fig. 2(b), it suggests that the particle size remains a
< 0.04), the MS curves at O(e) almost coincide with those at O(ed). In constant. Therefore, an increase in relative roughness suggests a
Fig. 2(a), the maximum Shields number is taken as 0.62. It has been decrease in flow depth. At O(e), the unstable region appears to dimin-
found that for H > 0.62, the MS curve at O(ed) produces an infeasible ish strikingly, as the relative roughness increases [Fig. 2(b)]. This is
region of instability, which contradicts the classical results. Therefore, owing to a change in the friction factor, in accord with Eq. (7). The sta-
Hm ¼ 0.62 represents the limiting Shields number for the formation bilizing effect remains effective for both shorter and longer wavenum-
of alternate bars in a smooth flow regime. Notably, this subtle conclu- bers. It indicates that even for the same aspect ratio, a small relative
sion was not drawn in the previous studies. roughness is favorable for the formation of alternate bars with both
It is interesting to find how the threshold aspect ratio bc varies shorter and longer wavenumbers, while a large relative roughness may
with the Shields number H. To this end, the bc(H) curves at O(e) and favor their formation at intermediate wavenumbers. However, at
O(ed) are plotted in the inset of Fig. 2(a). The bc(H) curves at both O(ed) [see dashed lines in Fig. 2(b)], dissimilar results are found.
orders show similar trend. For a given order, the threshold aspect ratio Unlike the behavior of the MS curves at O(e), for a given triplet (R,
increases with an increase in Shields number attaining a peak and then ^ the unstable region bounded by the MS curve at O(ed) appears
H, d),
reduces, as the Shields number increases further. In addition, for a to enhance with an increase in relative roughness. Note that the MS
given Shields number, the threshold aspect ratio at O(ed) is larger than curves at O(ed) nearly overlap with those at O(e) for shorter wave-
that at O(e). Moreover, the peak of bc(H) curve at O(e) appears early numbers. In Fig. 2(b), the maximum value of relative roughness is set
than that at O(ed). as 0.023, beyond which the MS curve at O(ed) yields an unrealistic

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-6

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

FIG. 4. The MS curves in a transitional flow regime, R ¼ 12 [solid lines: O(e) and
FIG. 3. The MS curves at the extremity of a smooth flow regime, R ¼ 3 [solid dashed lines: O(ed)]: (a) ^
d ¼ 0.005 and different values of H ¼ 0.15, 0.25, and
lines: O(e) and dashed lines: O(ed)]: (a) ^
d ¼ 0.005 and different values of H 0.38 and (b) H ¼ 0.25 and different values of ^ d ¼ 0.001, 0.005, and 0.009. The
¼ 0.15, 0.25, and 0.47 and (b) H ¼ 0.25 and different values of ^
d ¼ 0.001, 0.005, curves of bc vs H and bc vs ^
d are shown in the insets.
and 0.013. The curves of bc vs H and bc vs ^d are shown in the insets.

region of instability. Therefore, d^ m ¼ 0.023 is the limiting relative ¼ 0.47 and d^ m ¼ 0.013. Unlike Fig. 2(a), the bc(H) curves at O(e) and
roughness of the bar formation in a smooth flow regime. O(ed) show that the threshold aspect ratio decreases with an increase
Akin to Fig. 2(a), the variations of threshold aspect ratio bc with in Shields number [see the inset of Fig. 3(a)]. On the other hand, the
relative roughness d^ at both orders [O(e) and O(ed)] are shown in the ^ curves at O(e) and O(ed) remain quite simi-
qualitative trends of bc(d)
inset of Fig. 2(b). The behavior of bc(d)^ curves at O(e) and O(ed) is lar [see the insets of Figs. 2(b) and 3(b)].
similar. For a given order, the threshold aspect ratio reduces almost Figure 4 presents the MS curves at O(e) and O(ed) in a transi-
linearly with an increase in relative roughness reaching a minimum tional flow regime (R ¼ 12). In Fig. 4(a), the MS curves are shown for
value and then, it increases in a linear fashion with a further increase a given relative roughness d^ and different Shields numbers H, whereas
in relative roughness. For a given relative roughness, the threshold in Fig. 4(b), those are shown for a given Shields number H and differ-
aspect ratio at O(ed) appears to be larger than that O(e). Note that the ent values of relative roughness d. ^ It appears that the MS curves
minima of bc(d)^ curves appear at d^  0.004. change strikingly, as the shear Reynolds number increases (compare
Figure 3 displays the MS curves at O(e) and O(ed) at the extrem- Figs. 2 and 4). Unlike the MS curves at O(e) in a smooth flow regime
ity of a smooth flow regime (R ¼ 3). In Fig. 3(a), the MS curves are [see Fig. 2(a)], Fig. 4(a) displays that the unstable region confined to
shown for a given relative roughness d^ and different Shields numbers an MS curve increases with an increase in Shields number. Therefore,
H, while in Fig. 3(b), they are shown for a given Shields number H in a transitional flow regime (R ¼ 12), the sediment suspension intro-
and different values of relative roughness d.^ Although the qualitative duces a destabilizing effect, which appears to be pronounced for
patterns of MS curves remain similar (see Figs. 2 and 3), it appears shorter wavenumbers. Thus, for the same aspect ratio, a large Shields
that for a given relative roughness at O(e), the unstable region number (small particle size) is favorable for the bar formation at
bounded by an MS curve changes considerably for intermediate wave- shorter wavenumbers. The MS curves at O(ed) also show that that the
numbers as the Shields number increases [Fig. 3(a)]. The limiting val- unstable region bounded by an MS curve increases toward shorter
ues of Shields number and relative roughness are found to be Hm wavenumbers, as the Shields number increases. Similar to a smooth

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-7

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

[Fig. 2(b)], the unstable region at O(ed) in a transitional flow regime

increases, as the relative roughness increases [Fig. 4(b)]. For a given
^ the MS curve at O(ed) produces smaller unstable
triplet (R, H, d),
region than that at O(e). The MS curves at O(e) perfectly coincide
with those at O(ed) for shorter wavenumbers. It also appears that for a
given relative roughness, the unstable region decreases substantially, as
the shear Reynolds number increases [compare Figs. 2(b) and 4(b)]. In
Fig. 4(b), d^ m ¼ 0.009 represents the limiting relative roughness of the
formation of bars in a transitional flow regime (R ¼ 12). In the inset
of Fig. 4(b), the bc(H) curves at both orders [O(e) and O(ed)] display
a near-linear decreasing trend of the threshold aspect ratio with rela-
tive roughness. Note that for a given relative roughness, the threshold
aspect ratio at a given order increases with an increase in shear
Reynolds number [compare Figs. 2(b) and 4(b)].
As the transitional flow regime spans over a wide range of shear
Reynolds numbers (3 < R < 70), it is interesting to explore the behav-
ior of the MS curves with a further increase in shear Reynolds number
in this flow regime. Figure 5 shows the MS curves at O(e) and O(ed)
in a transitional flow regime for R ¼ 30. The MS curves for a given
relative roughness d^ and different Shields numbers H are shown in
Fig. 5(a), while those for a given Shields number H and different val-
ues of relative roughness d^ are shown in Fig. 5(b). It is noticeable that
^ the unstable region for a given order
for a given triplet (R, H, d),
reduces with an increase in shear Reynolds number (compare Figs. 4
and 5). The trends of the MS curves at O(e) and O(ed) for R ¼ 30
remain similar to those for R ¼ 12. Unlike Fig. 4(a), it appears that
for a small Shields number (say, H ¼ 0.15), the MS curves at both
O(e) and O(ed) coincide, as the shear Reynolds number increases [Fig.
5(a)]. However, for H > 0.15, the MS curve at O(ed) holds slightly
smaller unstable region than that at O(e). As the shear Reynolds num-
ber increases, the MS curves at both O(e) and O(ed) match for shorter
and intermediate wavenumbers [Fig. 5(a)]. The limiting Shields num-
FIG. 5. The MS curves in a transitional flow regime, R ¼ 30 [solid lines: O(e) and ber is found to reduce from Hm ¼ 0.38 to Hm ¼ 0.36, as the shear
dashed lines: O(ed)]: (a) ^
d ¼ 0.005 and different values of H ¼ 0.15, 0.25, and Reynolds number increases from R ¼ 12 to 30. The bc(H) curves [see
0.36 and (b) H ¼ 0.25 and different values of ^d ¼ 0.001, 0.005, and 0.0083. The the inset of Fig. 5(a)] indicate that for a given Shields number, the
curves of bc vs H and bc vs ^
d are shown in the insets.
threshold aspect ratios at both orders remain the same. For a given
Shields number, the threshold aspect ratio at a given order for R ¼ 30
is larger than that for R ¼ 12 [compare insets of Figs. 4(a) and 5(a)].
flow regime, for a given triplet (R, H, d), ^ the MS curve at O(ed) Unlike Fig. 4(b), as the shear Reynolds number increases, the MS
embraces smaller unstable region than that at O(e). As expected, the curve at O(ed), for a small relative roughness (d^ ¼ 0.001), occupies an
MS curves at both O(e) and O(ed) match for shorter wavenumbers. It unstable region being of the order of the unstable region formed by
is worth noting that for a given Shields number, the unstable region the MS curve at O(e) [Fig. 5(b)]. For d^ > 0.001, the MS curves at O(e)
reduces significantly, as the shear Reynolds number increases [com- almost coincide with those at O(ed), as the shear Reynolds number
pare Figs. 2(a) and 4(a)]. In Fig. 4(a), the limiting Shields number of increases [compare Figs. 4(b) and 5(b)]. Therefore, for d^ > 0.001, the
the formation of bars in a transitional flow regime (R ¼ 12) corre- unstable regions created by the MS curves at O(e) and O(ed) become
sponds to Hm ¼ 0.38. The bc(H) curves at both orders [O(e) and nearly equal for R ¼ 30. Comparison of Figs. 4(b) and 5(b) reveals
O(ed)] show that for a given order, the threshold aspect ratio reduces that the limiting relative roughness reduces from d^ m ¼ 0.009 to
with an increase in Shields number [inset of Fig. 4(a)]. Note that for a 0.0083 with an increase in shear Reynolds number from R ¼ 12 to
given Shields number, the threshold aspect ratios at both orders are 30. The inset of Fig. 5(b) shows that for a given relative roughness, the
nearly equal. Further, it appears that for a given Shields number, the threshold aspect ratios at both orders are almost equal. Comparison of
threshold aspect ratio at a given order in a transitional flow regime the insets of Figs. 4(b) and 5(b) suggests that for a given relative rough-
is larger than that in a smooth flow regime [compare Figs. 2(a) ness, the threshold aspect ratio at a given order for R ¼ 30 is larger
and 4(a)]. than that for R ¼ 12.
Figure 4(b) also shows that at O(e), as the relative roughness Figure 6 displays the MS curves at O(e) and O(ed) in a rough
increases, the unstable region reduces for shorter and longer wave- flow regime (R ¼ 100). In Fig. 6(a), the MS curves are shown for a
numbers, while that increases for intermediate wavenumbers. Similar given relative roughness d^ and different Shields numbers H, whereas
to the features of the MS curves at O(ed) in a smooth flow regime in Fig. 6(b), those are shown for a given Shields number H and

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-8

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

different values of relative roughness d.^ The MS curves contract strik-

ingly with an increase in shear Reynolds number, producing smaller
unstable region (compare Figs. 5 and 6). Similar to a transitional flow
regime, the MS curves at O(e) in a rough flow regime shows that the
unstable region bounded by an MS curve increases with an increase in
Shields number, destabilizing shorter wavenumbers owing to the pres-
ence of sediment suspension. At O(ed), as the Shields number
increases, the MS curves reveal that that the unstable region increases
toward shorter wavenumbers. Unlike smooth and transitional flow
regimes, the MS curve at O(ed) in a rough flow regime embraces the
same unstable region as that at O(e) [Fig. 6(a)], as the MS curves at
both orders coincide. The limiting Shields number in Fig. 6(a) corre-
sponds to Hm ¼ 0.34. Similar to a transitional flow regime, the bc(H)
curves in the inset of Fig. 6(a) at both orders [O(e) and O(ed)] overlap
perfectly. However, for a given Shields number, the threshold aspect
ratio at a given order in a rough flow regime appears to be larger than
those in smooth and transitional flow regimes. Akin to a transitional
flow regime, the MS curves at O(e) for different values of relative
roughness in a rough flow regime show that the unstable region, for a
given relative roughness, increases toward smaller aspect ratios with
intermediate wavenumbers, as the relative roughness increases [Fig.
6(b)]. The MS curves at O(e) and O(ed) in a rough flow regime, for a
given relative roughness, produce an identical unstable region. The
limiting relative roughness in Fig. 6(b) is found to be d^ m ¼ 0.007 8. In
the inset of Fig. 6(b), the bc(H) curves at both orders [O(e) and O(ed)]
almost overlap. Note that for a given relative roughness, the threshold
aspect ratio at a given order in a rough flow regime is larger than those
in smooth and transitional flow regimes.
Figs. 2–6 do not show how an MS curve, for a given relative
roughness and Shields number, evolves with the shear Reynolds num-
ber. It is therefore interesting to explore the sensitivity of an MS curve
to various flow regimes. To this end, Fig. 7 presents the MS curves for
different shear Reynolds numbers (R ¼ 1, 3, 12, 30, and 100) belong- FIG. 6. The MS curves in a rough flow regime, R ¼ 100 [solid lines: O(e) and
ing to hydraulically smooth, extremity of smooth, transitional, and dashed lines: O(ed)]: (a) ^
d ¼ 0.005 and different values of H ¼ 0.15, 0.25, and
rough flow regimes. The relative roughness and the Shields number 0.34 and (b) H ¼ 0.25 and different values of ^
d ¼ 0.001, 0.005, and 0.007 8. The
are considered as d^ ¼ 0.005 and H ¼ 0.25, respectively. The unstable curves of bc vs H and bc vs ^
d are shown in the insets.
region appears to contract with an increase in shear Reynolds number.
At O(e), the MS curves stabilize shorter and longer wavenumbers, as As discussed previously, the maximum values of Shields number
the shear Reynolds number increases. By contrast, at O(ed), the MS and relative roughness for a given flow regime in Figs. 2–6 indicate
curves destabilize longer wavenumbers, as the shear Reynolds number their limiting values Hm and d^ m , respectively. The limiting values cor-
increases. In smooth (R ¼ 1), extremity of smooth (R ¼ 3), and tran- respond to a critical flow, for which the flow Froude number is unity.
sitional (R ¼ 12 and 30) flow regimes, the unstable region bounded Beyond the limiting values of Hm and d^ m , the theoretical results at
by the MS curve at O(e) appears to be larger than that at O(ed). O(ed) produce two unstable regions. Among them, the first unstable
However, in a rough flow regime (R ¼ 100), the MS curves at both region favors the formation of bars beyond a threshold aspect ratio, as
orders become identical. In essence, Fig. 7 offers an understanding of furnished in Figs. 2–6. The second unstable region (not shown here)
how the MS curves evolve with the shear Reynolds number. In addi- suggests the formation of bars at all aspect ratios, thus contradicting
tion, it provides a qualitative prediction of the MS curves for an inter- the first unstable region in addition to a large corpus of theoretical and
mediate shear Reynolds number between any two neighboring shear experimental observations. Therefore, the second unstable region
Reynolds numbers, as shown in Fig. 7. The inset of Fig. 7 shows the appearing at longer wavenumbers represents an infeasible solution on
bc(R) curves at both orders [O(e) and O(ed)]. The threshold aspect the b(^k) plane.
ratio for a given order increases with an increase in shear Reynolds It is further interesting to explore the variations of limiting values
number. Note that for a given order, a point of inflection is noticeable of Shields number and relative roughness, Hm and d^ m , respectively,
in the bc(R) curve. For the chosen values of relative roughness and with shear Reynolds number R. Figure 8 shows the curves of Hm(R)
Shields number, the threshold aspect ratio, for a given shear Reynolds (solid lines) and d^ m (R) (dashed lines) for different values of relative
number, at O(ed) appears to be slightly larger than that at O(e) up to roughness and Shields number, respectively. The Hm(R) curves indi-
R ¼ 20. However, in a rough flow regime, both the orders predict the cate that for a given relative roughness, the limiting Shields number
same threshold aspect ratio. reduces with an increase in shear Reynolds number. In addition, for a

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-9

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

unsteadiness.6,10 In this regard, Hall10 reported linear and weakly non-

linear stability analyses of alternate bars that develop in an unsteady
flow. However, the effects of sediment suspension were overlooked in
Hall’s10 mathematical formulation. Therefore, a direct comparison of
this study with the linear results of Hall10 is not feasible. Note that
some preliminary conclusions of this study are in agreement with the
previous work of Federici and Seminara,11 as discussed before. It is
however important to highlight that the key results of this study may
offer guidelines for new research pathways, including the sensitivity of
free bars to the combined effects of flow unsteadiness, flow regimes,
and bedload and suspended load transport.

V. CONCLUSIONS
FIG. 7. Sensitivity of the MS curves to flow regimes for different values of R ¼ 1, A linear stability analysis of free river bars in weakly varying
3, 12, 30, and 100 [solid lines: O(e) and dashed lines: O(ed)]. The ^d ¼ 0.005 and flows with dominant sediment suspension is presented. The analysis is
H ¼ 0.25 are considered. The bc vs R is shown in the insets. done at leading and next orders of approximation. The behavior of sta-
bility curves is explored in different flow regimes. The results of linear
given shear Reynolds number, the limiting Shields number reduces, as stability analysis are discussed for different values of Shields number
the relative roughness increases. The d^ m (R) curves show that for a and relative roughness.
given Shields number, the limiting relative roughness reduces with an In a smooth flow regime, the threshold aspect ratio for a given
increase in shear Reynolds number. On the other hand, for a given order increases with an increase in Shields number reaching a peak
shear Reynolds number, the limiting relative roughness reduces, as the value and thereafter reduces with a further increase in Shields number.
Shields number increases. In general, Fig. 8 displays that the variations By contrast, the threshold aspect ratio reduces with an increase in rela-
of Hm(R) and d^ m (R) are more prominent in a smooth flow regime tive roughness attaining a minimum value and then, it increases
as compared to those in transitional and rough flow regimes. monotonically, as the relative roughness increases further. On the
It is pertinent to mention that the major focus of this study is to other hand, in transitional and rough flow regimes, the threshold
explore the stability of free river bars from a linear stability perspective aspect ratio for a given order reduces, as the Shields number and rela-
in weakly varying turbulent flows carrying significant load of sediment tive roughness increase. For a given shear Reynolds number, the
suspension. The key finding of this study is that the behavior of stabil- Shields number and the relative roughness play a decisive role in gov-
ity curves at the leading and next orders of approximation largely erning the bar instability. The former reflects the role of significant
depends on the values of shear Reynolds number, Shields number, sediment suspension, while the latter takes into account the effects of
and relative roughness. Another important outcome of this study, frictional resistance.
being unidentified in the previous studies, is that the analysis limits the At the leading order, as the shear Reynolds number increases,
values of Shields number and relative roughness, beyond which results shorter and longer wavenumbers are stabilized. However, at the next
at the next order show that the alternate bars are unable to develop. order, the marginal stability curves destabilize longer wavenumbers
However, this study does not take into account the effects of flow with an increase in shear Reynolds number. The unstable region for a
given order contracts significantly with an increase in shear Reynolds
number. In smooth and transitional flow regimes, the unstable region
at the leading order is larger than that at the next order. By contrast, in
a rough flow regime, the leading and next orders produce the same
region of instability.
In essence, this study offers an insight into the instability of free
bars in weakly varying flows, highlighting the role of shear Reynolds
number at different orders of approximation. The analysis is capable
to address the effects of Shields number and relative roughness on the
stability curves in different flow regimes. However, the analysis is pri-
marily focused on the free bar instability from a linear perspective.
Therefore, predicting the bar amplitude in an unsteady flow covering
different flow regimes by using a nonlinear approach remains an inter-
esting aspect as a future scope of research. In addition, analytical pre-
dictions of bar wavelength and bar amplitude in gravel-bed rivers
remain a challenging task. This arises owing to the complex fluid-
sediment interaction in the near-bed flow layer44–47 together with the
sediment transport.48,49 To resolve this issue, researchers need to use a
FIG. 8. Limiting Shields number Hm vs R for ^ d ¼ 0.001 and 0.005 (solid lines) 3D framework in modeling the near-bed flow, turbulent characteris-
and limiting relative roughness ^
d m vs R for H ¼ 0.15 and 0.25 (dashed lines). tics, and sediment transport.

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-10

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

ACKNOWLEDGMENTS Ws Sc
f¼ ; (B5)
ju
S.D. acknowledges the J C Bose Fellowship Award [Funded by
DST j Science and Engineering Research Board (SERB), Grant No. where Sc is the turbulent Schmidt number ( 1) and j is the von
JCB/2018/000004] in pursuing this work. Karman coefficient (¼ 0.41).
The function K0 in Eq. (20) is expressed as11
 
Cf0:5 I2
APPENDIX A: ESTIMATION OF HS K0 ¼ ð1  ^a Þ þ K2 ; (B6)
j I1
According to van Rijn,28 the condition required for the initia-
where Cf is the conductance coefficient (¼ f/8), I2 is the integral
tion of sediment suspension is the saltating length of a particle function, and K2 is a function that depends on Cf. The I2 and K2 are
exceeding 100 times the particle size. He proposed following empir- expressed as
ical relationships:
ð1  f
us 4 ^a 1  ^z
ð1 < D  10Þ ¼ ; (A1) I2 ¼ ðln ^z þ 1:84^z 2  1:56^z 3 Þ  d^z ; (B7)
Ws D 1  ^a ^z
^a
us j
ðD > 10Þ ¼ 0:4: (A2) K2 ¼ 0:777 þ : (B8)
Ws Cf0:5
In the above, us is the shear velocity at the threshold of sedi-
In addition, the function K1 in Eq. (20) is expressed as11
ment suspension. With us, the Shields number Hs at the threshold
of sediment suspension can be obtained from Eq. (9). For the esti- Cf0:5
mation of settling velocity Ws, Wu and Wang50 proposed an empir- K1 ¼ K3 ; (B9)
j
ical formula as follows: where K3 is a function of ^a and f. The K3 has the following form:43
8" 9
<  1=q #0:5 =
P t 1 4Q 3 1 ð1 ð1
Ws ¼  : þ D  ;; (A3)
Q d 4 3P  2 K3 ¼ C12 ðln ^z þ 1:84^z 2  1:56^z 3 Þd^z  ln ^z 0 C12 d^z ; (B10)
where ^z 0 ^a¼z0/D, z0 is the zero-velocity level where
^a streamwise
where P and Q are the coefficients and q is an exponent. They are velocity (x component) vanishes according to the logarithmic law,27
expressed as P ¼ 53.5exp(0.65Sp), Q ¼ 5.65exp(2.5Sp), and and C12 is a function. The C12 is obtained from the solution of the
q ¼ 0.7 þ 0.9Sp, where Sp is the Corey shape factor ( 0.7 for natural following differential equation:
particles).   
@C12 1 @ @
þ ^z ð1  ^z Þ C12
@^z f @^z @^z
APPENDIX B: EXPRESSIONS FOR C0, K0, AND K1   f
ð1  ^a ÞCf0:5 ^z ^a 1  ^z
¼ ln þ 1:84^z 2  1:56^z 3  :
The depth-averaged concentration at the leading order C0 in jI1 ^z 0 1  ^a ^z
Eq. (20) is expressed as11 (B11)
1 The boundary conditions associated with Eq. (B11) are expressed as
C0 ¼ C a I1 ; (B1)
1  ^a
C12 ð^z ¼ 1Þ ¼ 0; (B12)
where ^a ¼ a/D, a is the reference level that acts as the interface sep-
arating the bedload and the suspended load particles, Ca is the refer- @C12
ð^z ¼ ^a Þ ¼ 0: (B13)
ence concentration, and I1 is the integral function. For a and Ca, @^z
several empirical formulas exist in the literature.28,51,52 Here, we use Equation (B11) is solved numerically using a shooting technique.53
van Rijn’s28 relationships as follows:
 
d H APPENDIX C: COEFFICIENTS sj AND
Ca ¼ 0:015  1 D0:3 ; (B2)
a Hc PERTURBATIONS w010, w011, AND w110
aðks < 0:01DÞ ¼ 0:01D; aðks  0:01DÞ ¼ ks : (B3) The coefficients sj (j ¼ 1 to 4) in Eqs. (24) and (26) are
The integral function I1 in Eq. (B1) is expressed as expressed as
   
ð1  f H0 @f 1 1 @f H0 @f 1
^a 1  ^z s1 ¼ 2 1   ; s2 ¼  1  ;
I1 ¼  d^z ; (B4) f0 @H ^
f0 @ D f0 @H
1  ^a ^z (C1)
^
a H0 @U H0 @U 1 @U
s3 ¼  s1 ; s4 ¼  s2 þ  :
where ^z ¼ z/D, z is the vertical distance from the channel bed, and U0 @H U0 @H ^
U0 @ D
f is the Rouse number that quantifies the settling velocity relative to We express the perturbations w010, w011, and w110 in terms of
the shear velocity. The f is expressed as the perturbations of the primitive variables as

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-11

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

9
^ 10 þ t2 D
w010 ¼ t1 U ^ 10 ; (C2) M. Colombini, G. Seminara, and M. Tubino, “Finite-amplitude alternate bars,”
J. Fluid Mech. 181, 213–232 (1987).
^ 11 þ t2 D
w011 ¼ t1 U ^ 11 ; (C3) 10
P. Hall, “Alternating bar instabilities in unsteady channel flows over erodible
^ 10
@U ^ 10
@D beds,” J. Fluid Mech. 499, 49–73 (2004).
11
w110 ¼ g1 þ g2 : (C4) B. Federici and G. Seminara, “Effect of suspended load on sandbar instability,”
@^x @^x Water Resour. Res. 42(7), W07407, https://doi.org/10.1029/2005WR004399 (2006).
12
In the above, t1, t2, g1, and g2 are expressed as M. B. Bertagni and C. Camporeale, “Finite amplitude of free alternate bars with
 suspended load,” Water Resour. Res. 54(12), 9759–9773, https://doi.org/
H0 @f 1:5H0 10.1029/2018WR022819 (2018).
t1 ¼ s 1  þ  0:5ðI20 þ K20 I10 Þ1 13
R. Kinoshita, “Investigation of channel deformation in Ishikari River,” Report
2f0 @H H0  Hc
" #9 of the Bureau of Resources, Department of Science and Technology, Tokyo,
jH0 I10 @f = 14
Japan, 1961.

f0 ðI21 þ K20 I11 Þ þ  ; (C5) Y. Fujita, “Fundamental study on channel changes in alluvial rivers,” D.Eng.
f0 Cf0:5
0 @H ; thesis (Kyoto University, Kyoto, Japan, 1980).
 15
M. Iguchi, “Tests for fine gravel transport in a large laboratory flume,” Report
1 1 1:5H0
t2 ¼ s 2  þ  0:5ðI20 þ K20 I10 Þ1 of the National Science Foundation, School of Earth Science, Tsukuba
2 s2 H 0  H c University, Ibaragi, Japan, 1980.
9 16
" # S. Ikeda, “Prediction of alternate bar wavelength and height,” J. Hydraul. Eng.
jI10 = 110(4), 371–386 (1984).

f0 ðI21 þ K20 I11 Þ þ 0:5 ; (C6) 17
A. Crosato, F. B. Desta, J. Cornelisse, F. Schuurman, and W. S. J. Uijttewaal,
Cf 0 ;
“Experimental and numerical findings on the long-term evolution of migrating
  alternate bars in alluvial channels,” Water Resour. Res. 48(6), W06524, https://
s1 Ca0 K10 1:5H0
g1 ¼ I10  0:5f0 I11 ; (C7) doi.org/10.1029/2011WR011320 (2012).
1  ^a H0  Hc 18
J. P. C. Eekhout, A. J. F. Hoitink, and E. Mosselman, “Field experiment on alter-
  
s2 Ca0 K10 1:5H0 1 nate bar development in a straight sand-bed stream,” Water Resour. Res.
g2 ¼  I10  0:5f0 I11 ; (C8) 49(12), 8357–8369, https://doi.org/10.1002/2013WR014259 (2013).
1  ^a H0  Hc s2 19
M. T. Ramirez and M. A. Allison, “Suspension of bed material over sand bars
where f0, Cf0, Ca0, I10, I20, and K20 are the values of Rouse num- in the Lower Mississippi River and its implications for Mississippi delta envi-
ronmental restoration,” J. Geophys. Res.: Earth Surf. 118(2), 1085–1104,
ber, conductance coefficient, reference concentration, I1, I2, and
https://doi.org/10.1002/jgrf.20075 (2013).
K2 for the undisturbed flow, respectively. The I11 and I21 are 20
A. Defina, “Numerical experiments on bar growth,” Water Resour. Res. 39(4),
expressed as 1092, https://doi.org/10.1029/2002WR001455 (2003).
21
A. Siviglia, G. Stecca, D. Vanzo, G. Zolezzi, E. F. Toro, and M. Tubino,
ð1  f0  
“Numerical modelling of two-dimensional morphodynamics with applications
^a 1  ^z ^a 1  ^z
I11 ¼  ln  d^z ; (C9) to river bars and bifurcations,” Adv. Water Resour. 52, 243–260 (2013).
1  ^a ^z 1  ^a ^z 22
A. Crosato and E. Mosselman, “An integrated review of river bars for engineer-
^a
ing, management and transdisciplinary research,” Water 12(2), 596 (2020).
ð1  f   23
B. de Saint-Venant, “Theorie du mouvement non permanent des eaux, avec
^a 1 ^z 0 ^a 1 ^z
I21 ¼ ðln^z þ 1:84^z 2  1:56^z 3 Þ  ln  d^z : application aux crues de rivieras et a l’introduction des marces dans leur lit,” C.
1  ^a ^z 1  ^a ^z R. Acad. Sci. 73, 147–154 (1871).
^
a 24 €
F. M. Exner, “Uber die wechselwirkung zwischen wasser und geschiebe inin
(C10) fl€ussen,” Sitzungsber Akad. Wiss. 134, 165–203 (1925).
25
C. F. Colebrook and C. M. White, “Experiments with fluid friction in rough-
DATA AVAILABILITY ened pipes,” Proc. R. Soc. A 161, 367–381 (1937).
26
F. Engelund and E. Hansen, A Monograph on Sediment Transport in Alluvial
The data that support the findings of this study are available Streams (Danish Technical Press, Copenhagen, 1967).
from the corresponding author upon reasonable request. 27
S. Dey, Fluvial Hydrodynamics: Hydrodynamic and Sediment Transport
Phenomena (Springer-Verlag, Berlin, Germany, 2014).
28
L. C. van Rijn, “Sediment transport, Part II: Suspended load transport,”
REFERENCES J. Hydraul. Eng. 110(11), 1613–1641 (1984).
29
1
S. Dey and S. Z. Ali, “Fluvial instabilities,” Phys. Fluids 32(6), 061301 (2020). H. A. Einstein, “The bed-load function for sediment transportation in open
2 channel flows,” in Technical Bulletin (United States Department of Agriculture,
S. Lanzoni and M. Tubino, “Grain sorting and bar instability,” J. Fluid Mech.
393, 149–174 (1999). Soil Conservation Service, Washington DC, USA, 1950), p. 1026.
3 30
M. Redolfi, M. Welber, M. Carlin, M. Tubino, and W. Bertoldi, S. Ikeda, “Lateral bed load transport on side slopes,” J. Hydraul. Div. 108(11),
“Morphometric properties of alternate bars and water discharge: A laboratory 1369–1373 (1982).
31
investigation,” Earth Surf. Dyn. 8(3), 789–808 (2020). G. Parker, “Discussion on lateral bed load transport on side slopes,” J. Hydraul.
4 Eng. 110(2), 197–199 (1984).
B. Federici and G. Seminara, “On the convective nature of bar instability,”
32
J. Fluid Mech. 487, 125–145 (2003). S. Z. Ali and S. Dey, “Hydrodynamic instability of meandering channels,”
5 Phys. Fluids 29(12), 125107 (2017).
M. Tubino, R. Repetto, and G. Zolezzi, “Free bars in rivers,” J. Hydraul. Res.
33
37(6), 759–775 (1999). A. W. Baar, J. de Smit, W. S. J. Uijttewaal, and M. G. Kleinhans, “Sediment
6
M. Tubino, “Growth of alternate bars in unsteady flow,” Water Resour. Res. transport of fine sand to fine gravel on transverse bed slopes in rotating annular
27(1), 37–52, https://doi.org/10.1029/90WR01699 (1991). flume experiments,” Water Resour. Res. 54(1), 19–45, https://doi.org/10.1002/
7
S. Z. Ali and S. Dey, “Interfacial instability of sand patterns induced by turbu- 2017WR020604 (2018).
34
lent shear flow,” Int. J. Sed. Res. (in press) (2021). A. M. Talmon, N. Struiksma, and M. C. L. M. Van Mierlo, “Laboratory mea-
8
S. Z. Ali and S. Dey, “Instability of large-scale riverbed patterns,” Phys. Fluids surements of the direction of sediment transport on transverse alluvial-bed
33(1), 015109 (2021). slopes,” J. Hydraul. Res. 33(4), 495–517 (1995).

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-12

Published under license by AIP Publishing
 Physics of Fluids ARTICLE scitation.org/journal/phf

35 45
E. Meyer-Peter and R. M€ uller, “Formulas for bed-load transport,” in E. Padhi, N. Penna, S. Dey, and R. Gaudio, “Spatially averaged dissipation rate
Proceedings of the 2nd Meeting of International Association for Hydraulic in flows over water-worked and screeded gravel beds,” Phys. Fluids 30(12),
Research, Stockholm, Sweden (1948), Vol. 3, pp. 39–64. 125106 (2018).
36 46
S. Z. Ali and S. Dey, “Origin of the scaling laws of sediment transport,” Proc. E. Padhi, N. Penna, S. Dey, and R. Gaudio, “Near-bed turbulence structures in
R. Soc. A 473, 20160785 (2017). water-worked and screeded gravel-bed flows,” Phys. Fluids 31(4), 045107
37
S. Z. Ali and S. Dey, “Hydrodynamics of sediment threshold,” Phys. Fluids (2019).
47
28(7), 075103 (2016). N. Penna, E. Padhi, S. Dey, and R. Gaudio, “Structure functions and invariants
38
S. Dey and S. Z. Ali, “Stochastic mechanics of loose boundary particle transport of the anisotropic Reynolds stress tensor in turbulent flows on water-worked
in turbulent flow,” Phys. Fluids 29(5), 055103 (2017). gravel beds,” Phys. Fluids 32(5), 055106 (2020).
39 48
S. Dey and S. Z. Ali, “Mechanics of sediment transport: Particle scale of entrainment R. Maurin, J. Chauchat, B. Chareyre, and P. Frey, “A minimal coupled fluid-
to continuum scale of bedload flux,” J. Eng. Mech. 143(11), 04017127 (2017). discrete element model for bedload transport,” Phys. Fluids 27(11), 113302
40
S. Dey and S. Z. Ali, “Review article: Advances in modeling of bed particle (2015).
49
entrainment sheared by turbulent flow,” Phys. Fluids 30(6), 061301 (2018). D. Berzi and L. Fraccarollo, “Intense sediment transport: Collisional to turbu-
41
S. Dey and S. Z. Ali, “Bed sediment entrainment by streamflow: State of the sci- lent suspension,” Phys. Fluids 28(2), 023302 (2016).
50
ence,” Sedimentology 66(5), 1449–1485 (2019). W. Wu and S. S. Y. Wang, “Formulas for sediment porosity and settling veloc-
42
W. Wu and S. S. Y. Wang, “Movable bed roughness in alluvial rivers,” ity,” J. Hydraul. Eng. 132(8), 858–862 (2006).
51
J. Hydraul. Eng. 125(12), 1309–1312 (1999). F. Engelund and J. Fredsøe, “A sediment transport model for straight alluvial
43
M. Bolla Pittaluga and G. Seminara, “Depth-integrated modeling of suspended channels,” Hydrol. Res. 7(5), 293–306 (1976).
52
sediment transport,” Water Resour. Res. 39(5), 1137, https://doi.org/10.1029/ M. Garcia and G. Parker, “Entrainment of bed sediment into suspension,”
2002WR001306 (2003). J. Hydraul. Eng. 117(4), 414–435 (1991).
44 53
E. Padhi, N. Penna, S. Dey, and R. Gaudio, “Hydrodynamics of water-worked and R. S. Esfandiari, Numerical Methods for Engineers and Scientists Using
screeded gravel beds: A comparative study,” Phys. Fluids 30(8), 085105 (2018). MATLABV R , 2nd ed. (CRC Press, Boca Raton, 2017).

Phys. Fluids 33, 045105 (2021); doi: 10.1063/5.0045530 33, 045105-13

Published under license by AIP Publishing