Stability, resonance and role of turbulent stresses in 1D alluvial flows

Linear stability analysis is used to investigate the behavior of small perturbations of a uniform flow in a straight channel with an erodible bed composed by a unisize sediment. A shallow-water flow model is employed and bedload sediment transport is assumed. The mathematical structure of the linear problem, in terms of the eigenvalues and their associated eigenvectors is explored in detail and information is gathered on the wavespeed and growth rate of the perturbations as a function of their wavelength and of the relevant flow and sediment parameters. Several aspects of the solution are discussed, with particular focus on the behaviour in the transcritical region where the Froude number approaches unity. An approximate solution for the roots of the eigenrelationship is presented, which is not uniformly valid in the transcritical region, leading to the appearance of an unphysical instability. A regular perturbation expansion is then introduced that allows for the elimination of this singularity.


Introduction
Linear stability analysis has long been recognized as a powerful tool to investigate pattern formation in alluvial flows. Indeed, as Kennedy first pointed out [1], bed forms can be treated as sand waves that propagate (upstream or downstream) and grow (or decay) depending on a set of flow and sediment parameters. The sediment transport model used for the analysis is usually based on a mass conservation statement for the sediment phase, the Exner equation, coupled to a suitable empirical relationship for the sediment load that creates the necessary feedback between the bed wave, the secondary flow generated by the undulation of the bed, the resulting bed shear stress distribution, the sediment load and, eventually, the erosion-deposition processes that can either amplify or damp the bed wave according to the classic mechanism of instability. As for the flow, models of different complexity have been used, starting from the simplest 1 3 1D depth-averaged de Saint Venant system up to 3D direct numerical simulations, mostly depending upon the required detail for the prediction of the flow field and, in particular, of the bed shear stress [2].
Several bed patterns have been successfully studied in term of stability analysis, starting from bars [3,4], to dunes and antidunes [5][6][7][8], to ripples [9,10], to transverse or oblique bed forms like sand ridges [11] or diagonal dunes [12][13][14], just to mention the most common bed forms that are encountered in river environments. Although the topic of morphodynamic instability may appear as fairly settled, these analyses are periodically revisited and new details are continuously added to the overall picture [15][16][17].
In the following, the linear stability of a one-dimensional shallow-water flow over an erodible single-grain bed will be formulated and the nature of the resulting eigenvalues and eigenvectors will be discussed in detail. By this choice, we focus on possibly the simplest model for the flow and the sediment, the de Saint Venant-Exner (SVE) model. As it will become clear later on, none of the longitudinal bed forms previously mentioned are found unstable in the framework of the 1D-SVE model. The only instability about which information is gathered is that of roll-waves, which does not even fully qualify as a morphodynamic instability, since the presence of an erodible bed is not a necessary condition for their formation. Roll-waves are, in fact, easily observed in supercritical flows over a flat unerodible bed, being essentially the result of a free-surface instability [18,19]. On the other hand, at the cheap prize of a relatively simple flow and sediment transport model, the analysis will provide insightful information on the interactions between flow and bed waves and, consequently, on the dynamics of the flow-bed system. In particular, we aim at investigating: (i) the roll-wave hydrodynamic and morphodynamic instabilities; (ii) the effect of including a diffusive depth-averaged normal stress term in the SV model; (iii) the hydrodynamic resonance observed on a sinusoidal fixed bed in the transcritical ( Fr ∼ 1 ) region; (iv) the effect of coupling/decoupling flow and sediment dynamics.

Formulation of the 1D flow and sediment transport models
The starting point of our analysis is the 1D dimensionless form of the governing equations of hydrodynamics in a straight channel, written in the sloping Cartesian coordinate system sketched in Fig. 1. The triplet composed by the fluid density and by the uniform area velocity U * and depth D * is used for nondimensionalization (a star is used to denote dimensional variables). Moreover, a different scaling is chosen for the longitudinal coordinate [20], which is appropriate for the context of the shallow-water flow model adopted and reflects on the choice of the time scaling: where unstarred variables are dimensionless and a Chézy conductance coefficient is introduced here expressed through Keulegan relationship [21] as a function of the uniform flow depth and of the constant roughness K * , which in turn is set equal to 2.5 times the sediment grain size d * [22]. Typical values of C for alluvial river flows fall in the range 10-20.
The flow domain is bounded by the curve z = B(x, t) , which represents the actual bed, and by the free where D is the local flow depth. These two quantities, together with the local value of the depth-averaged velocity U are chosen as unknown functions. It is customary to assume the level B to coincide with the so-called 'reference level', the distance above the bed at which the logarithmic velocity profile vanishes and the bed shear stress is evaluated. The hydrodynamics is then governed by the following equations, which express mass conservation and the momentum balance [23]: where T t and T n are the bed shear stress and the depth-averaged Reynolds normal stress, respectively. Moreover, the uniform-flow law yields: where Fr is the Froude number of the undisturbed uniform flow.
Note that the normal stress term on the right hand side of (4) is usually disregarded, leading to the well-known de Saint Venant equation. However, as it will be shown in the following, this diffusive term provides a crucial stabilizing effect in the large wavenumber range and will be therefore retained in the analysis.
Several closures are available in the literature [23] for the turbulent stresses, which can be formally derived by assuming suitable profiles for the velocity and the eddy viscosity along the vertical. In particular, we adopt the formulation of [24], whereby both the eddy viscosity and the longitudinal velocity are assumed to depend on U and D, yielding: where the Von Kármán constant is taken as 0.4 and n is a binary switch used to turn the normal stress term on (1) and off (0) in the analysis.
The system (4-3) is complemented by the Exner equation imposing mass conservation of sediments, which takes the dimensionless form: where Φ is the sediment load function, q * s is the volumetric sediment discharge per unit width, p s is sediment porosity, s is the relative density of the sediment and is a small O(10 −3 ) parameter, which represents the ratio between the flux of sediments and the flow discharge, often interpreted as the ratio of the characteristic timescales of the flow and of the sediment transport.
Lastly, assuming dominant bedload, the sediment load function is related to the local value of the Shields parameter through the classical formula: where the values of C and m have been chosen as 0.0495 and 3.97, respectively [25].

Linearization
In the present section, the stability of a uniform flow in an infinitely wide channel with active sediment transport is explored, looking for the temporal evolution of a perturbation of this base state. The perturbed system can be thought as a set of traveling waves of given wavelength, the amplitude of which can either grow or decay in time, meaning that the base state is either unstable or stable. Linearization, which involves the assumption of small (strictly infinitesimal) amplitude of the perturbations, leads at the leading order to the formulation of an eigenvalue problem, the solution of which yields the complex wave speed (the eigenvalue) as a function of the wavenumber of the perturbation and of all the relevant flow and sediment parameters. No information on the amplitude of the perturbations is gathered at a linear level, but the eigenvectors corresponding to each eigenvalue can be used to assess, for example, the amplitude of the free-surface wave with respect to the bed wave.
To this end, we expand variables as: where is an arbitrarily small dummy parameter, i is the imaginary unit and c.c. stands for the complex conjugate of the preceding quantity. Moreover, the complex wavespeed w is related to the celerity and the growth rate of the perturbations by: and the wavenumber compares the dimensional wavelength of the perturbation with the longitudinal scale C 2 D * , which has been chosen for nondimensionalization. Introducing the wavenumber k, scaled with the uniform flow depth D * , we derive: Accordingly with the shallow-water framework considered herein, a value of ∼ O(1) represents waves that are much longer than the average flow depth, C being typically of O(10) for alluvial flows. Indeed, this provides a new physical interpretation of the long ( ≪ 1 ) and short ( ≫ 1 ) wave ranges, which will be widely investigated in the following. In particular, taking the limit → ∞ is not in contrast with k ∼ O(1) , which corresponds to the shortest wavelengths that can be resolved in a shallow-water framework, provided that this limit is thought of as the inviscid limit (C → ∞) of the solution, whereas friction dominates as → 0. Substituting the expansion (11) into the governing equations and collecting terms with the same power of , we obtain a series of algebraic problems that can be solved in cascade. In particular, at O( ) the following system is obtained: The three matrices in (15) represent the variable fluxes, the source terms and the diffusive terms appearing in the governing equations, respectively.
Moreover, we have: where Hence, the parameter tends to zero as 0 → C and is of order O(C −3 ) for active sediment transport ( 0 > C ). All the results presented in the following have been obtained with a value of C equal to 20, typical of alluvial flows, which corresponds to a grain-depth ratio d of about 0.0015 and to a value of that varies in the range 0.0002-0.003 depending on the Froude number. Equation (14) reveals itself as a classic eigenvalue problem, whereby the system admits a non-trivial solution only for those specific values of w, the eigenvalues, for which: where is the identity matrix and det(⋅) stands for the determinant of the array. The eigenvalues are the roots of the cubic eigenrelation: Although the roots of (19) can be obtained in closed form, their explicit expressions are not particularly instructive. Hence, before discussing the solution of the complete morphodynamic problem, we have chosen to introduce a series of simplified cases, each of interest on its own account, which will progressively disclose both the physics and the mathematics hidden beneath the solution of the linear eigenvalue problem (14), gradually unraveling the complex flow-bed interactions that ultimately control the growth and propagation of small perturbations in 1D alluvial flows.

Hydrodynamic stability: roll-wave formation
We start form the simple hydrodynamic problem, whereby the bed is assumed to be flat and unerodible. The latter conditions can be easily imposed on the eigensystem (18), by setting B 1 to zero (flat bed) and dropping the Exner equation (unerodible bed). By doing so, the system reduces to The hydrodynamic eigenvalues are the roots of the quadratic eigenrelation: where Let us first consider the case N = 0 , whereby the effect of the normal stress is neglected and the flow equations reduce to the classic de Saint Venant formulation. The wavespeed and the growth rate of the perturbations, are plotted in Fig. 2 as a function of the wavenumber, for several values of the Froude number. Note that the two eigenvalues are sometimes labeled 'fast' (h+) and 'slow' (h−) , because their wavespeed is either larger or smaller than unity (i.e. of the depth-averaged, uniform flow velocity used for nondimensionalization).
Both the growth rate and the wavespeed clearly tend to well defined limits for ≪ 1 and ≫ 1 . In order to formally derive these limits, we extract the real and imaginary parts of the characteristic polynomial (25), obtaining the following nonlinear algebraic system: In the long-wave range, the limit of (28) for → 0 yields: Hence, the longest waves propagate downstream, the slow one (h−) being always stable, the fast one (h+) being marginal (i.e. of vanishing growth rate). Note that the function f in (29) is positive for Fr > 2 , so that for small but finite values of the wavenumber, the growth rate of the fast eigenvalue is positive, leading to instability.
In the short-wave range, taking the limit of (28) for → ∞ provides: Hence, the shortest waves travel with the characteristic speed * h± = U * ± √ gD * . Moreover, for Fr > 2 the fast eigenvalue is unstable. Indeed, as shown in Fig. 3b, this is an example of a short-wave instability, whereby the maximum growth rate is found in the limit ≫ 1 . Note that the solution provided by (30) is regular [26], since the growth rate is bounded, though there is no wavelength selection mechanism because, as soon as the Froude number exceeds the value of 2, the shortest modes are all equally unstable. We recognize this instability as associated with roll-wave formation, a phenomenon which has been studied at length in the framework of the de Saint Venant model by means of linear [18,19,27], weakly nonlinear [28,29] and nonlinear [30] stability analyses.
Secondly, we switch on the depth-averaged normal stress ( N = 1 ) and extract the real and imaginary parts of the characteristic polynomial (25) Taking the limit of (31) in the long-wave range yields: A comparison with (29) shows that the normal stress term provides, at most, an O( 2 ) contribution, which does not alter the behaviour of the longest waves.
The solution in the short-wave range, being not uniformly valid in the limit, requires some care, but we eventually obtain: Hence, independently of the Froude number, the growth rate of both eigenvalues becomes negative (i.e. the uniform flow is stable). Moreover, the celerity of both waves tends to unity, showing how diffusion dominates for large values of . Note also that in the inviscid limit ( C → ∞ ), N tends to infinity as well as the group N 2 ∕C 4 , even if k ∼ O(1) . This confirms that the short-wave limit remains significant in the shallow-water framework considered. Considering the growth rate of the unstable fast (h+) eigenvalue, we eventually obtain: Therefore, for Fr > 2 (f > 0) the growth rate is positive for small but finite wavenumbers, whereas is negative for the shortest waves, necessarily reaching a maximum for intermediate wavenumbers, as shown in Fig. 3, where the growth rate of the hydrodynamic eigenvalues is presented, with (a) and without (b) the inclusion of the depth-averaged normal stress, for a value of C equal to 20. A comparison between the two pictures makes evident the role of the diffusive term. The wavenumber of maximum amplification, is found for intermediate wavenumbers and not in the short-wave limit anymore. Indeed, this 'most unstable' perturbation is the one that grows faster and, thus, it is expected to prevail on the others, at least in a linear framework. It is important to point out that the inclusion of a diffusive term into the 1D momentum equation has long been recognized [27,30,31] to provide a cut-off in the short-wave range. However, this term was mainly introduced in order to regularize the discontinuous shocks that appear in the solution when the classic de Saint Venant equation is used and the parameter N was considered as a dummy parameter controlling its intensity. In the present approach, the formal derivation of N borrowed from [24] and, in particular, its direct proportionality with the non dimensional Chézy coefficient C, provides a new insight on the role of this additional 'viscous' term in the determination of the wavelength of maximum amplification of roll-wave unstable perturbations.

Forced solution: the resonance
The formation of roll waves examined in Sect. 4 is a typical 'free' stability problem, whereby, for any given set of the flow parameters, disturbances of different wavelength are found to propagate either upstream or downstream while growing or decaying in amplitude. Let us now consider the 'forced' problem of flow over a rigid sinusoidal bottom of small amplitude (see [4,32] for an introduction to the topic of free-forced interactions and resonance in morphodynamic problems). Beside giving rise to a variety of interesting stability problems related to multiple steady states in supercritical flows [33], we want to investigate here the resonance between bed and free-surface waves that takes place in the transcritical Fr ∼ 1 region, often invoked in the past as the mechanism responsible for bed form instability [34,35].
The forced problem is easily solved by extracting the first two rows from (15), which can be rewritten as the non-homogeneous algebraic system: that is easily solved yielding A comparison of (35) with (23) immediately reveals a potential resonant behaviour of the system. Indeed, if one of the eigenvalues of the free problem vanishes, the system (35) becomes indeterminate: the external forcing (the rigid sinusoidal bed) is exciting one natural frequency of the system (the eigenvalue). Hence, we expect to observe resonance if perturbations of the flow do exist that are both neutral (vanishing growth rate) and stationary (vanishing wave speed), or, correspondingly, if Moreover, the solution of the forced problem (36) provides the flow response to a given sinusoidal bed perturbation of amplitude B 1 . Recalling that det(A h ) is a complex quantity, the relative amplitude of the free surface oscillation H 1 is eventually obtained: Neglecting the role of the normal stress ( N = 0 ), the limit of (38) for → ∞ is found to be: and is shown as a thick black line in Fig. 4a where the relative amplitude of the free surface oscillation is presented for several values of . A resonance is detected in the transcritical region ( Fr ≃ 1 ) in the short-wave 'inviscid' range ( ≫ 1 ), which is progressively damped as is reduced. Moreover, the amplitude of the free surface oscillation: (i) tends to vanish (flat surface) in the limit Fr → 0 ; (ii) is smaller (larger) than the amplitude of the bed wave for Fr smaller (larger) than 1∕ √ 2 ; (iii) tends to unity (constant depth) for Fr → ∞. The maximum amplitude H 1M , which is obtained from (37)(38) for Fr = 1 , reads: and is shown in Fig. 4b as a function of with and without the effect of the normal stress T n . As expected, the 'viscous' term associated with the depth-averaged normal stress acts to damp the resonance, but it must be pointed out that, for a value of C equal to 20, its effect starts to become noticeable for values of of O(10 2 ) , corresponding to k of O(1) , or, in other words, to wavelengths scaling with the flow depth, which represent the shortest waves that can be resolved in the framework of a depth-averaged shallow-water flow model.

Morphodynamic stability: the quasi-steady case
A third instructive simplification of the complete morphodynamic problem is related to the so-called quasi-steady approximation. Most of the existing linear stability theories on bed form formation are based on this procedure, which dates back to [1]. The rationale behind the approach is that sediment transport evolves on a time scale much longer than that of the flow, so that, as bed perturbations grow, the flow can be thought as instantaneously adapting to bed changes. The quasi-steady approximation can be formally enforced by neglecting time derivatives in the flow equations, thus imposing steady flow over a slowly evolving sinusoidal bed. Setting N = 0 , the following degenerate eigenvalue problem is obtained: which provides a single complex eigenvalue w b , where the suffix b stands for 'bed' since w b is representative of the bed dynamics alone. It is interesting to note that the procedure for the determination of this eigenvalue can be formally split in two parts: (i) using the first two lines of (41) U 1 is obtained as a function of B 1 , as in the forced solution (36); (ii) the latter is inserted into the third row of (41), yielding: In the short-wave inviscid limit we eventually obtain Both the growth rate and the wavespeed are proportional to the small parameter , thus implying that the timescale on which bed waves evolve and propagate is much longer than that of the flow perturbations, as expected. The growth rate is consistently negative and so the uniform base flow is stable, independently of the Froude number. Conversely, the celerity changes sign at Fr = 1 , so that the bed wave propagates downstream (upstream) in the subcritical (supercritical) regime.
Finally, resonance clearly affects both the solutions since the intensity of the eigenvalue becomes unbounded for Fr = 1 . This behaviour is due to the fact that, under the quasisteady approximation, the bed acts as an external forcing for the flow. Indeed, this does not come unexpectedly, since the resonant forced solution (36) has been used in the determination of the bed eigenvalue. Note, however, that resonance does not produce any amplification of the bed wave. On the contrary, short-wave perturbations are strongly damped in the transcritical region.

Morphodynamic stability: the fully coupled case
It may be useful at this point to remark that the resonant solutions we have discussed so far appear only because the flow and the bed system are considered as two distinct systems, which interact only mildly. Indeed, two conditions are required for the onset of a full resonance at Froude equal to unity: (i) the bed must be fixed or, equivalently, the dynamics of the bed is ignored; (ii) the flow must be steady, or, equivalently, the dynamics of the flow is ignored. Neither of these conditions apply in the coupled case, hence no resonance is expected to take place. This notwithstanding, moved by the very same results we have obtained in the previous sections, several researchers [34][35][36] have in the past investigated resonance in great detail and, pushing this concept to the limit, invoked resonance as the only real mechanism driving ripple and dune instability: Instability of flow systems of rippled or dune-covered equilibrium beds is shown to occur at finite growth rates for a range of wavelengths via a resonance mechanism occurring for surface waves and bed waves traveling at the same celerity [35].
A simple reasoning suggests however that if the flow-bed system is considered as a whole and the flow and the bed dynamics are fully coupled by correctly modeling the flow-bed interactions, no external forcing is present and, thus, no resonance can possibly appear [37]. Hence, the instability associated with transcritical resonance must be unphysical and artificially introduced at some point in the analysis. Solving this common misconception, which periodically reemerges in the study of bed forms [38], will be the task we will confront with in the present section, where the fully coupled case is considered.
To this end, we will ignore the role of the normal stress, which have been shown to damp resonance, and focus on the behaviour of the eigenvalues in the region (Fr → 1, → ∞) , where resonance, if any, is expected.
As before, in order to formally derive these limits, we extract the real and imaginary parts of the characteristic polynomial (19) and take their limits for → ∞ , yielding: In Fig. 5 the roots of the cubic equation (44) are shown. The three morphodynamic eigenvalues have been labeled (m1, 2, 3) in order of decreasing celerity, from the largest positive to the largest negative. By comparison with the hydrodynamic eigenvalues, showed by the dashed black lines, it appears that the (m1) eigenvalue is practically indistinguishable from (h+) , whereas (h−) overlaps (m2) in the supercritical regime and (m3) in the subcritical regime.
Although it could be tempting to discuss the behaviour of the morphodynamic eigenvalues in terms of the flow (h±) waves plus a new, slower, 'bed' wave, this would imply a change of sign of two roots at Fr = 1 . In fact, as shown in Fig. 5b, in the transcritical region none of the celerities vanishes and consequently there is no change of sign as in the hydrodynamic case: for any given Froude number, including unity, (m1) and (m2) propagate downstream and (m3) upstream [39]. Moreover, moving away from criticality, (m2) and (m3) tend to the quasi-steady (b) eigenvalue in the subcritical and supercritical regime, respectively.
It is worth noting that the cubic equation (44), which represents the short-wave limit of (19), coincides with the characteristic polynomial of the linearized system matrix of the Saint-Venant-Exner (SVE) model. Moreover, (44) has three distinct real roots, as it can be formally proved by considering its discriminant and making use of the fact that ≪ 1 . We obtain: Hence, the condition (46) ensures that the system of PDEs is hyperbolic for any value of the parameters. For this reason, the roots of (44), as well as their behaviour in the transcritical region, have been thoroughly discussed in the past due to the role this equation plays in the definition of the characteristic curves of the SVE hyperbolic system [39][40][41].
Conversely, this is probably the first time that a link is established with the stability of small amplitude morphodynamic waves in the short-wave range. This is ultimately due to the double interpretation of the limit → ∞ , which can be considered as either the inviscid limit for finite k or the short-wave limit for finite C. In this regard, the present contribution complements the works of [42], where wave hierarchies in alluvial flows were investigated, and of [20], where the stability of long waves in erodible channels was studied. In Fig. 6, the growth rates of the morphodynamic eigenvalues are shown, using the same color code as in Fig. 5. As expected, the red (m1) eigenvalue behaves as (h+) and, in particular, it becomes unstable ( Ω m1 > 0 ) as the Froude number exceeds the value of 2, leading to the formation of morphodynamic roll waves [24,31,[42][43][44]. Note that this is the only unstable mode of the morphodynamic problem, the other two eigenvalues being stable ( Ω m2,m3 < 0 ) and following the behaviour of (h−) in the supercritical and subcritical regimes, respectively. In the transcritical region, as shown in the zoom of Fig. 6b, (m2) and (m3) tend to (b) moving away from criticality.
Considering the right eigenvectors associated with the three morphodynamic eigenvalues, we obtain, in the limit → ∞: where the m are the roots of (44) and r is a generic constant. The eigenvectors express the relative amplitudes of the perturbations [24,45] and can be used to estimate, for each eigenvalue, the amplitude of the free-surface wave with respect to the bed wave, similarly to what we did in Sect. 5. We obtain: The relative amplitudes are plotted in Fig. 7, together with the forced solution given by (39).
We first remark that the eigenvector associated with the unstable (m1) eigenvalue takes values of O(10 3 −10 4 ) , thus confirming the hydrodynamic character of the roll-wave instability, whereby the amplitude of the bed perturbation is only a small fraction of the amplitude of the free-surface wave. This is ultimately due to the fact that the small parameter appears in the denominator of (48). Conversely, whenever (m2) and (m3) behave as (b) , m 2 , 1 (i.e. in the subcritical and supercritical regime, respectively) the relative amplitudes are of O(1) and follow closely the forced solution given by (39). This behaviour is open to a nice physical interpretation: the bed wave is so slow that is felt by the flow as a stationary wave, which produces O(1) free-surface oscillations as in the forced, fixed bed, case. Accordingly, the relative amplitude of the free surface is smaller (larger) than unity in the subcritical (supercritical) regime and tends to vanish as Fr ≪ 1 and to unity as Fr ≫ 1.
The most relevant information emerging from the analysis of the eigenvectors is, however, that no resonant behaviour is present in the transcritical region, as is evident from the zoom of Fig. 7b. This is ultimately due to the fact that, however small it may be, none of the wave celerities vanishes for Fr → 1 and so a full resonance cannot take place. This substantiates the previous intuition that the free system, whereby the bed and the flow dynamics fully interact, cannot possibly resonate the way a forced system does.
Unfortunately, having eliminated the possibility of a resonance-driven mechanism of instability, we are left with a very small catch in the net: to this point, the only longitudinal bed form we have been able to describe, is the result of the roll-wave instability, which drives the growth of a fast-traveling bed wave of amplitude much smaller than that of the free-surface wave. Any hope of describing in the present context the formation of longitudinal bed forms like dunes, antidunes and ripples clashes with the inability of the shallowwater flow model to provide the correct lag of the bed shear stress with respect to the bed perturbation for wavelengths that are of the order of the flow depth. A more refined, non hydrostatic, rotational flow model is required in order to fill this gap [8,10].

Morphodynamic stability: the approximate outer solution
Although the roots/celerities of the cubic equation (44) can be determined exactly, an approximate solution provides a valuable insight into the nature of the eigenvalues and eigenvector. The last simplification we want to introduce is similar to the one adopted (a) (b) Fig. 7 The amplitude of the free surface associated with the three morphodynamic eigenvalues for a unitary bed amplitude (a); plot (b) is an enlargement of the transcritical ( Fr ≃ 1 ) region. The dashed black line represents the forced solution given by (39) under the quasi-steady approximation and is related to the concept of decoupling flow and bed dynamics. Instead of neglecting the time derivatives in the flow equations as we did in Sect. 6, decoupling is here introduced by expanding the celerity and growth rate associated with each eigenvalue as: where we take advantage of the typically small value attained by the parameter . This approach was first introduced by [43], who investigated the wave propagation and the stability of shallow-water flows in erodible bed channels. The scientific community interested on the properties of the hyperbolic SVE system used the same approach to find approximate solutions for the slope of the characteristic curves and, more generally, to study the linear and nonlinear dynamics of the system [40,42,46]. More recently, the failure of the above expansion in the transcritical region has been investigated and an alternative solution valid in that region has been suggested [20,39]. In the following, the procedure outlined in the latter works will be revisited, extending the analysis to the growth rates and the eigenvectors to obtain a complete picture of the behaviour of small morphodynamic waves and of their stability. We will also show that in the transcritical region a spurious instability appears, which is related to the failure of the expansion (49) and, ultimately, to the unphysical resonance previously discussed.
Substituting (49) into the cubic polynomial (44) and equating likewise powers of yields the following equations.
At the leading O( 0 ) order, we obtain: which provides the simple solution: where h± are the celerities of the fast and slow eigenvalues of the hydrodynamic problem discussed in Sect. 4. Hence, the celerities of two eigenvalues reduce to those of the hydrodynamic case and the third one vanishes. For this reason, it seems natural to label the latter as the 'bed eigenvalue', as we did in Sect. 6. Indeed, it is not surprising that the bed eigenvalue vanishes at leading order: the information propagated through sediment transport, i.e. sediment waves, migrate at a speed that is typically much smaller than the speed of hydrodynamic information, which was the heart of the matter for the quasi-steady approximation. Moreover we note that, due to (16), the vanishing of implies C → ∞ and so the inviscid limit, whereby the flow does not respond to the drag action of the bed.
Moving to the next O( 1 ) order, we obtain: so that the wave speeds associated with the three morphodynamic eigenvalues can be approximated at O( ) by which were first presented by [40], with a different sign in the equation corresponding to the fastest eigenvalue (m1) that was corrected by [42]. Note that in order to remain coherent with the labeling scheme adopted for the exact solution (i.e. (m2) downstream propagation, (m3) upstream propagation), the two eigenvalues must be switched at criticality. Due to the smallness of the parameter , upon which the solution is built, (53) well represents the celerities of the complete morphodynamic problem (44). Indeed, outside of the transcritical region, the approximate solutions are almost indistinguishable from the exact ones presented in Fig. 5a. On the contrary, as Fr → 1 the approximate solution (53) fails for (m2) and (m3), as shown in Fig. 8a, where the celerities of the approximate eigenvalues (dashed) are compared with their exact counterparts (solid). We can easily recognize that the expansions for (m2) and (m3) are not regular in a neighbourhood of Fr = 1 , this also suggesting a way this singularity can be removed [47].
Similarly, substituting the expansion (49) into (45) and using (53) yields The growth rates associated with the three eigenvalues provide a similar picture with respect to the celerities previously discussed. In particular, outside of the transcritical region the exact solution, shown in Fig. 6a, is well approximated by (54). On the contrary, as Fr → 1 the approximate solution (54) fails for (m2) and (m3), as shown in Fig. 8c, where the approximate growth rates (dashed) are compared with their exact counterparts (solid).
In the transcritical region, which we recognize as the region where the bed wave forces a near resonant response of the flow, the growth rates of (m2) and (m3) are either positive, thus implying instability, or identical to the one obtained under the quasi-steady approximation, both the solutions being affected by the resonance as well.
Following the same procedure, we also obtain, from (48): which provide the approximate solutions shown in Fig. 8e. In the outer region, the approximate amplitude of the free-surface oscillation for a unitary bed perturbation is practically undistinguishable from the exact one given in Fig. 7a, whereas the presence of a resonance in the inner transcritical region is confirmed. Indeed, as compared to the decoupled quasi-steady solution, this approximate solution provides a weak coupling, whereby the flow and the bed dynamics mildly interact. In this regard, it is not surprising that the almost stationary bed (celerity is of O( ) ) forces a resonance of the slow flow eigenvalue, which, in the transcritical region, is itself almost stationary, being h2 ≃ 0 . However, resonance appears here as the result of a singular perturbation problem which makes the expansion (49) fail in the transcritical region. The main difference with respect to the quasi-steady case is that resonance affects here two eigenvalues instead of just one, and, most important, that this ultimately leads to the appearance of a new, though artificial, mode of instability. The presence in the transcritical region of a "very weak" linear instability, rapidly damped by nonlinear effects, was reported by [20]. In this regard, it is also worth mentioning the work of [35], who adopting a potential flow model showed that the free surface and the bed wave can resonate at criticality, generating an instability with a mechanism quite similar to the one we have just described. Although time derivatives were retained in the flow equations and the system of equations was formally fully coupled, this is again a case of weak coupling due to the peculiar flow model adopted. The potential flow model is in fact unable to provide a lag between the sediment transport and the bed, so that there is no real feedback between the flow and the bed dynamics and a resonance is possible.
The terminology 'inner' and 'outer', borrowed from boundary layer theory, refers here to the transcritical and the non-transcritical region, respectively.
In particular, let us consider a neighbourhood of Fr = 1 of the kind: where n defines the size of the set, n is a generic positive exponent and F is the O(1) inner variable. Substituting the expansion (56) into (44) yields which represents the inner equation. Taking the limit for → 0 we easily recognize that the undisturbed problem has the simple solution: where the isolated root 0 = 2 clearly represents the non-singular (m1) eigenvalue and can safely be ignored.
Following the lead of [40], an higher-order approximation for the vanishing roots is sought by setting where the the exponents n and m, making use of the so-called dominant balance method [47], are found both equal to one half.
Hence, at the leading order we obtain: which yields the inner-region approximate wavespeeds of the (m2) and (m3) eigenvalues, expressed in terms of the outer variable Fr: The inner approximate solutions (61) are compared with the exact solutions in Fig. 8b, showing a good agreement in the transcritical region. In particular, the inner solutions are almost indistinguishable from the exact ones in the subcritical (supercritical) branch of m2 ( m3 ), where the outer solution tends to b . Substituting (61) into (45), the leading-order contribution to the inner solution for the growth rate is found as: Following the same procedure, the leading-order contribution for the inner solution of the free-surface relative amplitude is obtained from (48) and reads (56) The comparison between the exact and the inner solutions for the growth rate and the freesurface relative amplitude is shown in Fig. 8d, f, respectively. They confirm the proper fit of the inner solutions and, in particular, the complete disappearance of the spurious instability associated with the outer solution, since the uniform flows results stable and the free surface oscillation remains limited. The weak coupling provided by the approximate solution is therefore sufficient to produce a well-behaved solution in the whole range of Froude, but care must be taken in the inner transcritical region, where the outer solution fails producing a resonant behaviour. Indeed, we stress once more that no real resonance is present in the morphodynamic problem at hand. It is only when one interferes with the strong coupling between flow and sediment dynamics, thus breaking or weakening this link, that resonance resurfaces again. This is particularly relevant for numerical methods that aim at accelerating the dynamics of the bed with respect to that of the flow [49,50], which should be carefully tested for the appearance of resonant behaviours that can affect the solution, in particular in the transcritical region.

Conclusions
The behaviour of small perturbations of a uniform flow in a straight channel with an erodible bed composed by a unisize sediment has been investigated by means of a linear stability analysis. The mathematical structure of the linear problem, in terms of the eigenvalues and their associated eigenvectors has been explored in detail and information has been gathered on the wavespeed and the growth rate of the perturbations as a function of their wavelength and of the relevant flow and sediment parameters. Several aspects of the solution are discussed, with particular focus and on the presence of a resonant behaviour in the transcritical region and on the damping role of the depth-average normal stress. In particular, a link between the eigenvalue problem emerging from the linear stability analysis in the short-wave limit and the eigenvalue problem that characterizes the classic problem of wave hierarchies in alluvial river flows is established.
Results show that any interference with the complex mechanisms that control the coupling between the flow and the bed dynamics is prone to produce unrealistic solutions that manifest themselves through the appearance of a spurious resonance in the transcritical region. This happens both in the quasi-steady case, whereby time derivatives are dropped in the flow equations, and in the case of weak coupling, whereby the solution is expanded in power series of the small parameter , which compares the timescales of flow and sediment transport and thus controls the strength of the coupling. The approximate solution found in the latter case is shown to be not uniformly valid in the transcritical region. A regular (non resonant) perturbation expansion has been presented that allows for the elimination of this singularity. Indeed, the availability of analytical approximate solutions for the eigenvalues and the eigenvectors can be of great importance in the formulation of numerical models for the solution of the nonlinear SVE system based on local linearization of the problem.
Despite the crude models used to describe the flow and the bed dynamics, a theoretical framework has been provided in which some classic results about stability, resonance and coupling have been revisited and find a new perspective. Possible future developments of this work could imply an extension of the approach to the case of a bed composed by a mixture of sediments of different sizes.
Funding Open access funding provided by Università degli Studi di Genova within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.