Mega riverbed-patterns: linear and weakly nonlinear perspectives

In this paper, we explore the mega riverbed-patterns, whose longitudinal and vertical length dimensions scale with a few channel widths and the flow depth, respectively. We perform the stability analyses from both linear and weakly nonlinear perspectives by considering a steady-uniform flow in an erodible straight channel comprising a uniform sediment size. The mathematical framework stands on the dynamic coupling between the depth-averaged flow model and the particle transport model including both bedload and suspended load via the Exner equation, which drives the pattern formation. From the linear perspective, we employ the standard linearization technique by superimposing the periodic perturbations on the undisturbed system to find the dispersion relationship. From the weakly nonlinear perspective, we apply the centre–manifold-projection technique, where the fast dynamics of stable modes is projected on the slow dynamics of weakly unstable modes to obtain the Stuart–Landau equation for the amplitude dynamics. We examine the marginal stability, growth rate and amplitude of patterns for a given quintet formed by the channel aspect ratio, wavenumber of patterns, shear Reynolds number, Shields number and relative roughness number. This study highlights the sensitivity of pattern formation to the key parameters and shows how the classical results can be reconstructed on the parameter space.

SZA, 0000-0003-0763-7437; SD, 0000-0001-9764-1346 In this paper, we explore the mega riverbed-patterns, whose longitudinal and vertical length dimensions scale with a few channel widths and the flow depth, respectively. We perform the stability analyses from both linear and weakly nonlinear perspectives by considering a steady-uniform flow in an erodible straight channel comprising a uniform sediment size. The mathematical framework stands on the dynamic coupling between the depth-averaged flow model and the particle transport model including both bedload and suspended load via the Exner equation, which drives the pattern formation. From the linear perspective, we employ the standard linearization technique by superimposing the periodic perturbations on the undisturbed system to find the dispersion relationship. From the weakly nonlinear perspective, we apply the centre-manifold-projection technique, where the fast dynamics of stable modes is projected on the slow dynamics of weakly unstable modes to obtain the Stuart-Landau equation for the amplitude dynamics. We examine the marginal stability, growth rate and amplitude of patterns for a given quintet formed by the channel aspect ratio, wavenumber of patterns, shear Reynolds number, Shields number and relative roughness number. This study highlights the sensitivity of pattern formation to the key parameters and shows how the classical results can be reconstructed on the parameter space.  mega riverbed-patterns was found to be convective [29,31], as any tiny local disturbance is convected downstream without disturbing the flow domain over a long duration [30]. However, the convective nature of patterns does not appear to have changed even in the presence of sediment suspension [30]. The mega riverbed-patterns have been also analysed theoretically considering the effects of unsteadiness [28,31], vegetation [38,39] and sediment heterogeneity [27]. The sediment heterogeneity causing the selective transport of particles was found to dampen the wavelength, growth rate and migration speed of patterns [27].
The weakly nonlinear analysis for bedload-dominated channels was carried out by using the multiple-scale technique [3] and the centre-manifold-projection technique [39], while that for combined bedload and suspended load transport was performed by using the centre-manifoldprojection technique [32], extending the previous linear analysis [30]. In the weakly nonlinear analysis, the amplitudes of patterns for both plane and dune-covered beds were predicted [32]. It was found that the presence of dunes is to enhance the equilibrium amplitude of patterns.
Despite magnificent advances in studying the mega riverbed-patterns theoretically [4,8], little is known about how the patterns are sensitive to a class of parameters pertinent to flow and particle transport. Although a recent study has clearly shown that the inclusion of sediment suspension changes the stability and amplitude of patterns significantly [32], more insights are to be gained regarding how and why these changes evolve with the characteristic parameters. Moreover, the behaviours of the growth rate and amplitude of patterns in transitional and rough flow regimes need to be examined.
In this paper, we analyse the mega riverbed-patterns, being free from an external forcing, from both linear and weakly nonlinear perspectives. Employing the depth-averaged formulations for the flow and particle transport, we follow the standard mathematical techniques to analyse the pattern formation. From the linear perspective, we apply the standard linearization technique, whereas from the weakly nonlinear perspective, we apply the centre-manifold-projection technique [32]. The main objective of the linear stability analysis is to explore how the mega riverbed-patterns are sensitive to the characteristic parameters, e.g. channel aspect ratio, shear Reynolds number, Shields number and relative roughness number. In addition, the primary objective of the weakly nonlinear stability analysis is to explore the sensitivity of the finite amplitude of patterns to the key parameters.
The paper is organized as follows. In §2, the theoretical analysis of flow is presented. The theoretical analysis of particle transport is described in §3. The stability analyses from the linear and weakly nonlinear perspectives are furnished in §4. The computational results are discussed in §5. Finally, the conclusion is drawn in §6.

Theoretical analysis of flow
We consider mega riverbed-patterns, having a wavelength λ * and a finite amplitude A * , driven by a free-surface turbulent flow in an infinitely long straight channel of width 2B * ( figure 1a,b). Henceforth, a star superscript is used to represent a dimensional variable. We employ a Cartesian coordinate system (x * , y * , z * ), where x * , y * and z * are the longitudinal, transverse and vertical distances, respectively. The origin of the coordinate system lies on the centreline of the channel. Two well-recognized patterns are highlighted in figure 1. In figure 1a, the alternate bars are shown, where the channel embraces a single row of bars. On the other hand, in figure 1b, the central bars are shown, covering two parallel rows of alternate bars. The presence of bars causes channel-bed oscillations in both longitudinal and transverse directions. For alternate bars (figure 1a), a typical section S 1 -S 2 across the channel shows the channel-bed erosion and deposition near the right and left channel-banks, respectively. In addition, for central bars (figure 1b), the cross-sectional view shows the channel-bed erosion and deposition near the channel-banks and at the channelbed centreline, respectively. At a given time t * , D * (x * , y * , t * ) is the local flow depth, which is the vertical distance between two characteristic points P 1 and P 2 at a given transverse distance ( figure 1a,b). In addition, z * = R * (x * , y * , t * ) denotes the elevation of the channel-bed from a fixed   horizontal datum. Note that in figure 1, the dashed rectangular section represents the undisturbed cross section of flow.
In the theoretical analysis, we employ the quasi-steady and the depth-averaged assumptions. The former is appropriate here, as the time scale of the morphological patterns is quite larger than that of the flow. In addition, the latter assumption emerges from the fact that the wavelength of the mega riverbed-patterns scales with a few channel widths [3]. Now, we introduce the pertinent variables in dimensionless form as follows: where D * 0 and U * 0 are the undisturbed flow depth and flow velocity (i.e. the flow depth and flow velocity corresponding to the uniform flow), respectively, U * = (U * x , U * y ) is the depth-averaged velocity vector, T * = (T * x , T * y ) is the bed shear stress vector and ρ f is the mass density of fluid. Moreover, we introduce the channel aspect ratio β and the undisturbed flow Froude number Fr as where g is the acceleration due to gravity. The depth-averaged momentum and continuity equations of flow are expressed as In order to close the hydrodynamic equations, an estimation of the bed shear stress vector is essential. The bed shear stress vector is expressed herein as a function of dynamic pressure.
Hence, we write where C f is the conductance factor. It can be linked with the Darcy-Weisbach friction factor f as f = 8C f . The C f and f can be obtained from the Colebrook-White equation as follows [40]: 1 where k * s is the bed roughness height (= αd * ), α is a coefficient (= 2.5) [41], d * is the particle size, a r and a s are the constants (= 14.8 and 2.51, respectively), Re is the flow Reynolds number and υ is the coefficient of kinematic viscosity of fluid. The above equation provides an accurate estimation of the friction factor over a wide range of flow regimes. The shear Reynolds number R * is commonly used to demarcate the flow regimes [42]. The R * is defined as where u * f is the shear velocity.

Theoretical analysis of particle transport
The erodible channel-bed comprises uniform cohesionless sediment particles. The flow over the channel-bed exerts hydrodynamic force on the particles. The particles are entrained into the flow when the dimensionless fluid-induced bed shear stress crosses its threshold value Θ c [43]. Both the modes of particle transport as bedload and suspended load transport are considered herein. The bed shear stress is quantified by introducing the Shields number Θ as follows: where s is the relative density (= ρ p /ρ f ) and ρ p is the mass density of particles. It is pertinent to mention that the suspended load transport prevails when the Shields number Θ exceeds its threshold value Θ s to bring the particles into the suspension. Hence, for Θ c < Θ < Θ s , the bedload remains the dominant mode of particle transport, while for Θ > Θ s , both the bedload and suspended load transport take place. For the estimation of Θ s , we use van Rijn's [44] empirical formulae as In the above, D * = d[(s − 1)g/υ 2 ] 1/3 , called the particle parameter and w * s is the settling velocity of particles. where In the above, is the suspended load flux vector and ρ 0 is the porosity of particles. We consider the channel-banks to be impermeable to the flow and the sediment fluxes, so that the quantities U y , Φ by and Φ sy vanish at the channel-banks (y = ±1).
The bedload flux vector Φ b is expressed as where χ is the angle that the resultant bedload flux makes with the longitudinal direction and Φ is the bedload flux intensity. The χ is expressed as follows [46,47]: where ϑ is the transverse slope coefficient. For the bedload flux intensity Φ, we use Meyer-Peter & Müller's formula as follows [48]: For the estimation of Θ c , we use the Θ c (D * ) relationships of Cao et al. [49]. Bolla Pittaluga & Seminara [37] presented an analytical solution for the estimation of suspended load flux in slowly varying flows by means of asymptotic expansion of the exact solution of the advection-diffusion equation of sediment concentration. In this study, to obtain the suspended load flux vector, we employ the two-dimensional extension of the depth-averaged formulation of Bolla Pittaluga & Seminara [37], as was done by Federici & Seminara [30] and Bertagni & Camporeale [32]. The perturbation approach assumes that both the advective and unsteady effects are trivial compared with the gravitational settling and turbulent diffusion. This suggests δ p = [U * 0 D * 0 /(w * s λ * )] to be much smaller than unity. As the δ p depends on the priori unknown wavelength of riverbed-patterns, we define another small parameter δ as δ = δ p λ * /B * = U * 0 /(βw * s ) [30,32]. The suspended load flux vector is expressed as follows [37]: where the function ψ is expanded as ψ = ψ 0 + δψ 1 + O(δ 2 ). Here, ψ 0 corresponds to the uniform flow condition and ψ 1 is the correction to ψ at O(δ). This correction takes into account the weak non-equilibrium effects due to the spatial change in the flow field. The ψ 0 and ψ 1 are expressed as follows [37]: where C 0 is the depth-averaged suspended particle concentration at the leading order, and K 0 and K 1 are functional parameters.
The C 0 is expressed as follows [30]: In the above, a = a * /D * , z = z * /D * , a * is the reference level, C a is the reference concentration, Z is the Rouse number and κ is the von Kármán coefficient (= 0.41). In this study, we use van Rijn's [44] empirical formulae to obtain a and C a . The K 0 in equation (3.12) is expressed as follows [30]: (3.14) In the above, K 2 = 0.777 + κC −1/2 f . In addition, the K 1 in equation (3.12) is expressed as follows [30,37]: where z 0 = z * 0 /D * , z 0 is the zero-velocity level at which the time-averaged longitudinal velocity becomes zero (in accord with the classical logarithmic law) and C 12 is a function. The C 12 can be found by imposing the boundary conditions, namely ∂C 12 /∂z = 0 at z = a and C 12 = 0 at z = 1, to the following differential problem: (3.16)

Stability analysis
We expand the set of variables where R 0 refers to the value of R at the undisturbed state. We express the components of bed shear stress and bedload flux vectors as where C f 0 , Φ 0 and Θ 0 are the conductance factor, bedload flux intensity and Shields number, respectively, at the undisturbed state and the coefficients c 1−4 are given in appendix A. Further, in order to expand the suspended load flux vector, we express ψ 0 and ψ 1 as where ψ 00 refers to the value of ψ 0 at the undisturbed state. The coefficients c 5−8 are given in appendix A. Substituting equations (4.1)-(4.4) into equations (2.3)-(2.5) and (3.6), we obtain where L 0 is the 4 × 4 matrix whose elements are zero except the lower right entry (being equal to unity), L m is the 4 × 4 matrix differential operator and the term N(V) contains all the second-order nonlinearities. The perturbations can be expanded in terms of the eigenfunctions as follows [ where A m,p corresponds to the amplitude of the longitudinal mode p and the transverse mode m, k is the dimensionless longitudinal wavenumber and i is the unit imaginary number. In the above expansion, the eigenfunctions are the product of eigenvectorsv m and the exponential term.
We write the eigenvectorsv m as follows [32]: where c.c. stands for the complex conjugate. The eigenvectorsv m are the product of another vector u m depending on the wavenumber of perturbation and the corresponding Fourier harmonics satisfying the boundary conditions. Note that the subscript m denotes the variable associated with the transverse mode.

(a) Linear perspective
Linear analysis leads to the prediction of the favourable domain of the pattern formation.
In the linear analysis, we consider only the fundamental longitudinal mode, p = 1. Therefore, each transverse mode leads to an independent problem. The amplitude is considered to grow exponentially with time as where ω m is the complex number, whose real and imaginary parts represent the growth rate and the dimensionless frequency, respectively. With p = 1, equations (4.6)-(4.8) are substituted into equation (4.5) and the nonlinearities are ignored. The resultant linear system is where elements l ij of the matrix L m are expressed as The solvability condition of equation ( [32], this analysis aims at obtaining the finite amplitude solution for the fundamental mode (m = 1). Initially, we solve the following linear adjoint eigenvalue problem: where the superscript † stands for the adjoint operator and the overbar denotes the complex conjugate. In accordance with the boundary conditions, the inner product defining the adjoint operator is expressed as 1 −1 (La 1 ) ·ā 2 dy = 1 −1 a 1 · (L † a 2 ) dy, (4.13) where L is the linear operator, and a 1 and a 2 are the generic vectors. In this study, L 0 = L † 0 and L † m represents the conjugate transpose of L m . It is worth mentioning that for L = L 0 , the eigenvectorŝ v m and adjoint eigenvectorsv † m become orthonormal after the proper normalization. Hence, we write where δ ij is the Kronecker delta function, and i and j are the transverse modes. Unlike the linear analysis, the nonlinear interactions of different modes are to be accounted for in the weakly nonlinear analysis. In this study, we consider that the contributions from the higher-order modes are negligible. Therefore, the expansion of equation (4.6) is considered up to p = 2 and m = 2.
We substitute the expansion of equation (4.6) into equation (4.5), take the inner product of the resulting equation with the adjoint eigenfunctions and subsequently collect the terms of the same Fourier modes. This operation yields the amplitude equation of the fundamental and two other stable modes as follows:  (4.17) and The centre-manifold-projection technique facilitates to project the amplitudes of the stable modes on the unstable ones. Therefore, the stable amplitudes A m,2 can be written as a nonlinear combination of A 1,1 andĀ 1,1 as where p 1 , p 2 and p 3 are the projection coefficients. After neglecting the higher-order terms, the time derivative of equation (  To solve the above equation, we consider p 2 = p 3 = 0. Therefore, the coefficient p 1 is expressed as . (4.22) In addition, the amplitudes of the stable modes can be expressed as Substituting equation (4.23) into equation (4.15), the Stuart-Landau equation describing the amplitude dynamics of the fundamental mode can be obtained as where the coefficient σ is expressed as The equilibrium amplitude of the fundamental mode, e.g. alternate bars, can be obtained by setting the time derivative in equation (4.24) equal to zero. Therefore, the dimensionless amplitude of alternate bars can be expressed as In the above, the factor 2 is involved due to the complex conjugation.

Results and discussion
In the numerical computations, we consider the mass density of fluid ρ f = 1000 kg m −3 , mass density of particles ρ p = 2650 kg m −3 , porosity of particles ρ 0 = 0.3, transverse slope coefficient ϑ = 0.56 and acceleration due to gravity g = 9.81 m s −2 .
To present the results from the linear perspective, we first examine how the mega riverbedpatterns evolve with different transverse modes. To this end, we consider a transitional flow regime and set the characteristic value of the shear Reynolds number as R * = 50. In addition, the reference values of Shields number and the relative roughness number are taken as Θ = 0.2 and d r = 0.0001, respectively. This combination produces the flow Froude number Fr = 0.166. Figure 2 shows the aspect ratio β as a function of dimensionless wavenumber k for different transverse modes m (= 1, 2, 3 and 4). Note that each transverse mode corresponds to the specific rhythmic patterns in both longitudinal and transverse directions, as conceptually sketched in the insets of figure 2a-d. To be specific, m = 1 corresponds to one row of alternate bars, while m = 2 corresponds to two rows of alternate bars, i.e. central bars. In addition, m = 3 and 4 comprising three and four rows of alternate bars correspond to the channelized bars. In each panel of figure 2, the stability map for a given m is shown on the β(k) plane. The marginal stability curve (solid line) makes a distinction between the unstable zone (shaded portion) and the stable zone (white portion). It defines the locus of the vanishing growth rate on the β(k) plane, i.e. Re(ω m ) = 0. Therefore, the unstable and stable zones correspond to Re(ω m ) > 0 and Re(ω m ) < 0, respectively. From a physical viewpoint, in the unstable zone (zone confined to the marginal stability curve), the incipient riverbed-patterns with a given transverse mode grow with time. On the other hand, in the stable zone (zone outside the marginal stability curve), the incipient riverbed-patterns with a given transverse mode decay yielding no patterns. Figure 2 shows that the patterns are unable to grow below a threshold aspect ratio β c for which the slope of the β(k) curve disappears, i.e. dβ/dk = 0. The β c therefore sets the limit for a stable alluvial channel having no riverbed-patterns. This indicates that the patterns are likely to form if the aspect ratio is larger than its threshold value. It is evident that the threshold aspect ratio increases with an increase in transverse mode (figure 2). For a given β β c , the unstable zone captures longer wavenumbers as the transverse  mode increases. Note that in the subsequent sections, we consider only the fundamental mode (m = 1). Figure 2 does not offer an understanding of how the marginal stability curve for a given set (R * , m) is sensitive to the Shields number and the relative roughness number. To explore this aspect, figure 3 shows the marginal stability curves for different Shields numbers and relative roughness numbers. The shear Reynolds number and the transverse mode are considered to be R * = 50 and m = 1 (fundamental mode), respectively. In figure 3a, the β(k) curves are shown for a given relative roughness number (d r = 0.0001) and different Shields numbers (Θ = 0.2, 0.3 and 0.4 yielding different flow Froude numbers Fr = 0.166, 0.204 and 0.235, respectively). As the shear Reynolds number is considered to be constant (R * = 50), an increase in Shields number in figure 3a refers to a reduction in particle size, as evident from equation (3.4). The unstable zone for a given Shields number increases as the Shields number increases (figure 3a). This is credited to the presence of significant sediment suspension that causes a destabilizing effect. The threshold aspect ratio reduces with an increase in Shields number. However, the longer wavenumbers at a large aspect ratio are stabilized as the Shields number increases. On the other hand, in figure   (R * = 50 and Θ = 0.2) in figure 3b, the particle size becomes a constant. It follows that an increase in relative roughness number produces a reduction in the flow depth. This causes an increase in the friction factor. The unstable zone for a given relative roughness number increases as the relative roughness number increases (figure 3b), lowering the threshold aspect ratio. This causes the threshold aspect ratio to reduce. It is interesting to figure out how the growth rate of patterns evolves with the Shields number and the relative roughness number. To this end, we consider the fundamental mode (m = 1) and set the shear Reynolds number as R * = 50 (transitional flow regime). Figure 4a,b displays the dimensionless growth rate Re(ω 1 ) as a function of dimensionless wavenumber k for different Shields numbers and relative roughness numbers, respectively. In each panel, the growth rate is examined for three different cases, β < β c , β = β c and β > β c . In figure 4a, the Re(ω 1 )(k) curves are shown for a given relative roughness number (d r = 0.0001) and different Shields numbers (Θ = 0.2, 0.3 and 0.4 yielding different flow Froude numbers Fr = 0.166, 0.204 and 0.235, respectively). The growth rate curves corresponding to β < β c , β = β c and β > β c are shown by solid, dotted and dashed lines, respectively. As the Re(ω 1 ) = 0 sets the limit for the stable patterns, the shaded zone above Re(ω 1 ) = 0 represents the unstable zone. For β < β c (solid lines), the growth rate is expected to be negative. It is apparent that the growth rate for a given Shields number increases (decrease in negative magnitude) with a marginal increase in wavenumber depicting a peak at a shorter wavenumber, and thereafter it reduces (increase in negative magnitude) with a further increase in wavenumber (figure 4a). In general, the growth rate for a longer wavenumber dampens (increase in negative magnitude) as the Shields number increases. For β = β c (dotted lines), the growth rate for a given Shields number follows a similar trend with the wavenumber. Note that the growth rate for a given Shields number is negative except at the critical wavenumber, where it is zero. The growth rate for a longer wavenumber diminishes (increase in negative magnitude) as the Shields number increases. For β > β c (dashed lines), the growth rate for a given Shields number changes from negative to positive at a shorter wavenumber, attaining a positive peak with an increase in wavenumber and then it reduces with a further increase in wavenumber. The growth rate makes another changeover from positive to negative at a longer wavenumber and reduces (increase in negative magnitude) as the wavenumber increases ( figure 4a). Importantly, for a given Shields number and β > β c , the crossings (where the growth rate changes its sign) correspond to the two points on the marginal stability curve β(k). For a given intermediate wavenumber, the growth rate enhances with the Shields number. However, for a given shorter (or longer) wavenumber, it reduces (increase in negative magnitude) as the Shields number increases. In For β < β c (solid lines), the growth rate for a given relative roughness number increases (decrease in negative magnitude) with a small increase in wavenumber attaining a peak and thereafter reduces (increase in negative magnitude) with a further increase in wavenumber. The growth rate for a given wavenumber diminishes (increase in negative magnitude) with an increase in relative roughness. For β = β c (dotted lines), the growth rate for a given relative roughness number is negative everywhere except at the critical wavenumber, which produces a vanishing growth rate. The growth rate for a given wavenumber reduces (increase in negative magnitude) as the relative roughness number increases. Similar to figure 4a, for β > β c (dashed lines), the growth rate for a given relative roughness number changes its sign twice at shorter and longer wavenumbers (figure 4b). For a given intermediate wavenumber, the growth rate increases with an increase in relative roughness number. However, for a given shorter (or longer) wavenumber it decreases (increase in negative magnitude) as the relative roughness number increases. Figures 3 and 4 highlight the marginal stability curves and the growth rates in a transitional flow regime. However, riverbed-patterns are sensitive to flow regimes. Therefore, the evolution of patterns with shear Reynolds number is an important aspect. Figure 5a shows the marginal stability curves corresponding to the fundamental mode (m = 1) for different shear Reynolds numbers (R * = 50, 100 and 500). The Shields number and the relative roughness number are taken as Θ = 0.2 and d r = 0.0001, respectively. This combination produces different but nearly equal flow Froude numbers Fr = 0.166, 0.168 and 0.169. Note that R * = 100 and 500 correspond to a rough flow regime. As the Shields number is kept constant (Θ = 0.2), an increase in shear Reynolds number refers to an increase in particle size (see equation (3.4)). The unstable zone for a given shear Reynolds number reduces with an increase in shear Reynolds number due to an increase in particle size (figure 5a). Therefore, the threshold aspect ratio increases as the shear Reynolds number increases. The longer wavenumbers at a large aspect ratio are stabilized as the shear Reynolds number increases. It is worth emphasizing that for a given shear Reynolds number R * , the growth rate of patterns for a given aspect ratio β(> β c ) attains a maximum value at some characteristic wavenumber, called the maximum unstable wavenumber. For a given shear Reynolds number, the locus of the maximum growth rate is shown in figure 5a (see dotted lines). In fact, the loci of the maximum growth rate for different shear Reynolds numbers almost  coincide. Figure 5b shows the dimensionless growth rate Re(ω 1 ) as a function of wavenumber k corresponding to β < β c , β = β c and β > β c for different shear Reynolds numbers (R * = 50, 100 and 500). For β < β c (solid lines), the growth rate for a given shear Reynolds number increases (decrease in negative magnitude) with a small increase in wavenumber reaching a peak and thereafter decreases (increase in negative magnitude) with a further increase in wavenumber. The growth rate for a given longer wavenumber increases (decrease in negative magnitude) as the shear Reynolds number increases. For β = β c (dotted lines), the growth rate for a given longer wavenumber also increases (decrease in negative magnitude) as the shear Reynolds number increases. However, for β > β c (dashed lines), the growth rate for a given intermediate wavenumber lessens as the shear Reynolds number increases. Note that the growth rate for a given shear Reynolds number changes its sign twice at shorter and longer wavenumbers (figure 5b). Figure 6 offers comparison of the computed dimensionless wavelength of patterns with the available experimental data [10][11][12][13]17]. The comparison is performed for the fundamental mode of riverbed-patterns, e.g. alternate bars (see the conceptual sketch in the inset of figure 6) by making the pertinent variables on a par with the experimental data. A brief summary of the experimental conditions used for the comparison is given in table 1. Note that some of the experimental data in figure 6 primarily correspond to the bedload-dominated channels. Therefore, the theoretical results for them in figure 6 are obtained by neglecting the role of sediment suspension. In this context, it is worth mentioning that an accurate estimation of the transverse slope coefficient ϑ from the experimental data is a difficult proposition. Therefore, the ϑ is kept as a free parameter. It is found that for ϑ ≈ 0.2, the theoretical predictions are in satisfactory agreement with the experimental data, as most of the data are confined to the ±20% error band ( figure 6). In fact, the choice of ϑ = 0.2 is close to the previously proposed value (= 0.56) of the transverse slope coefficient [52].
To present the results from the weakly nonlinear perspective, we first explore the amplitude of patterns for a given set (β, k, R * , Θ, d r ). Figure 7 shows the contours of dimensionless amplitude on the β(k) plane. The flow regime is considered to be transitional (R * = 50). The Shields number and the relative roughness number are taken as Θ = 0.2 and d r = 0.0001, respectively. For a given aspect ratio β(> β c ), the amplitude varies with the wavenumber over a certain range. For a given  Kinoshita [10] Yoshino [11] Ikeda [12] Lanzoni [13] Boraey [17]   aspect ratio, the amplitude amplifies towards the shorter wavenumbers. On the other hand, for a given wavenumber, the amplitude varies insignificantly with the aspect ratio ( figure 7). The sensitivity of the amplitude of patterns to the key variables remains a key aspect. Therefore, it is further interesting to envision how the amplitude of patterns varies with the pertinent variables, such as Θ, d r and R * . Figure 8a-c shows the maximum unstable wavenumber k m and the dimensionless amplitude A as a function of the Shields number Θ, relative roughness number d r and shear Reynolds number R * , respectively. In figure 8a, the k m (Θ) and A(Θ) curves are shown for d r = 0.0001 and R * = 50. The maximum unstable wavenumber reduces with an increase in Shields number, while the amplitude increases. In figure 8b, the k m (d r ) and A(d r ) curves are shown for Θ = 0.2 and R * = 50 on a semi-logarithmic scale. The maximum unstable wavenumber increases marginally as the relative roughness number increases, while the amplitude remains nearly invariant with the relative roughness number. In addition, in figure 8c, the k m (R * ) and A(R * ) curves are shown for Θ = 0.2 and d r = 0.0001. The maximum unstable wavenumber appears to remain constant with an increase in shear Reynolds number, while the amplitude reduces attaining a constant value for large shear Reynolds numbers (R * > 300). Figure 9 shows comparison of the computed dimensionless amplitude of patterns with the experimental data corresponding to the bedload-dominated channels [10][11][12][13]17]. The comparison is performed for the fundamental mode of riverbed-patterns, e.g. alternate bars (see the  conceptual sketch in the inset of figure 9). The experimental conditions used for the comparison are furnished in table 1. The ϑ is set equal to 0.2, as was also considered in figure 6. The theoretical predictions match well with the experimental data, as most of the data excepting a few fall within the ±20% error band (figure 9).

Conclusion
We analyse the sensitivity of mega riverbed-patterns of finite amplitude without an external forcing from both linear and weakly nonlinear perspectives. The depth-averaged flow model and the particle transport model are coupled by means of the Exner equation, which governs the dynamics of the physical system that describes a steady-uniform flow in a straight channel with an erodible channel-bed. The stability analyses include the standard linearization technique and the centre-manifold-projection technique from linear and weakly nonlinear perspectives, respectively. The stability of mega riverbed-patterns is found to be largely sensitive to the channel Yoshino [11] Ikeda [12] Lanzoni [13] Boraey [17]  aspect ratio, wavenumber of patterns, shear Reynolds number, Shields number and relative roughness number. The marginal stability, growth rate and amplitude of patterns are explored.
The unstable zone for a given Shields number increases with an increase in Shields number due to the destabilizing effect caused by the sediment suspension. In addition, the unstable zone for a given relative roughness number increases, as the relative roughness number increases. However, the unstable zone for a given shear Reynolds number reduces as the shear Reynolds number increases. For a given shear Reynolds number, Shields number and relative roughness number, the threshold aspect ratio increases with an increase in transverse mode.
Within the unstable zone, the growth rate for a given intermediate wavenumber increases with an increase in Shields number. Further, the growth rate for a given intermediate wavenumber increases as the relative roughness number increases. However, the growth rate for a given intermediate wavenumber reduces as the shear Reynolds number increases. The loci of the maximum growth rate for different shear Reynolds numbers approximately coincide.
Within the unstable zone, the amplitude of patterns intensifies towards shorter wavenumbers. The amplitude increases with an increase in Shields number, while it remains nearly constant with an increase in relative roughness number. However, the amplitude reduces as the shear Reynolds number increases.
From the linear perspective, this study reveals that among the chosen range of characteristic parameters, the marginal stability curve is more sensitive to the Shields number and the shear Reynolds number than to the relative roughness number. However, the growth rate is largely sensitive to the relative roughness number. From the weakly nonlinear perspective, the finite amplitude of patterns is most sensitive to the Shields number.
It is worth emphasizing that the limitation of the linear stability analysis is its inefficacy to predict the finite amplitude of patterns. This can be resolved by performing a weakly nonlinear stability analysis. However, the weakly nonlinear stability analysis is not capable of addressing the complete nonlinear interactions of the physical system. This may require a sophisticated numerical treatment.
Essentially, this study provides an insight into the stability of mega riverbed-patterns from both linear and weakly nonlinear perspectives. However, the theoretical results in the presence