Linear stability analysis of plane beds under ﬂows with suspended load

. Plane beds develop under ﬂows in ﬂuvial and marine environments; they are recorded as parallel lamination in sandstone beds, such as those found in turbidites. However, whereas turbidites typically exhibit parallel lamination, they rarely feature dune-scale cross lamination. Although the reason for the scarcity of dune-scale cross-lamination in turbidites is still debated, the formation of dunes may be dampened by suspended load. Here, we perform, for the ﬁrst time, linear stability analysis to show that ﬂows with suspended load facilitate the formation of plane beds. For a ﬁne-grained bed, suspended 5 load can promote the formation of plane beds and dampen the formation of dunes. These results of theoretical analysis were veriﬁed with observational data of plane beds under open-channel ﬂows. Our theoretical analysis found that suspended load promotes the formation of plane beds, which suggests that the development of dunes under turbidity currents is suppressed by the presence of suspended load.


Introduction
The interactions between fluids and erodible surfaces generate small-scale topographic features called bedforms both on terrestrial surfaces (e.g., riverbeds, deserts, and deep-sea floors) and on extra-terrestrial surfaces (Bourke et al., 2010;Gao et al., 2015;Hage et al., 2018;Cisneros et al., 2020).Such bedforms are preserved in sedimentary rocks as sedimentary structures such as cross-and parallel lamination (Harms, 1979).The types of sedimentary structures observed vary among different types of rocks.Turbidites typically exhibit parallel lamination (Bouma, 1962), whereas they rarely feature dune-scale crosslamination (Talling et al., 2012).However, the opposite is true for fluvial deposits; i.e., dune-scale cross laminae are often observed in riverine sandstone (Miall, 2010).
Although the reason for the paucity of dune-scale cross-lamination in turbidites is still debated (Lowe, 1988;Arnott, 2012;Schindler et al., 2015;Tilston et al., 2015), it could be attributed to the presence of suspended load.For example, in the case of open-channel flows, nearly flat bed waves and low-angle dunes have been observed in suspension-dominated rivers (Smith and McLean, 1977;Kostaschuk and Villard, 1996;Bradley et al., 2013;Ma et al., 2017).Additionally, flume experiments have suggested that dune height decreases with increasing suspended load flux (Bridge and Best, 1988;Naqshband et al., 2017).
Therefore, the influence of suspended load on the suppression of dune development and the formation of plane beds is worth investigating.
The relationships between sediment transport modes and the formation of plane beds have received little attention in theoretical works that performed linear stability analyses.The reason could be because previous studies have succeeded in predicting the wavelength of dunes and antidunes without considering suspended load (Colombini, 2004;Di Cristo et al., 2006;Colombini and Stocchino, 2008;Vesipa et al., 2012;Bohorquez et al., 2019).However, this assumption is not appropriate for analyzing open-channel flows where the suspended load is not negligible, such as flows in rivers with a fine sediment bed (de Almeida et al., 2016;Sambrook Smith et al., 2016).Moreover, although some research has considered both bed-and suspended load (Engelund, 1970;Nakasato and Izumi, 2008;Bose and Dey, 2009), the hydraulic conditions of these analyses were limited, and the results were tested using only observational data of dunes and antidunes.
Therefore, in order to investigate the effect of sediment transport mode on the formation of plane beds, we performed a linear stability analysis of bedforms under open-channel flows carrying suspended load.The model introduced in Nakasato and Izumi (2008) was extended in this study to evaluate plane bed formation under various conditions of sediment diameter and flow depth.To evaluate the suspended load effect, linear stability analyses were performed on flows both with and without suspended load.Further, we tested our stability diagrams against observational data of plane beds.Our theoretical analysis reveals for the first time that suspended load promotes the formation of plane beds, which has implications for interpreting sedimentary structures in turbidites.

Methods
Linear stability analysis of fluvial bedforms can provide the wavelengths of perturbations (i.e., bed waves) that grow over time (Colombini, 2004;Bohorquez et al., 2019).We employ the two-dimensional Reynolds-averaged Navier-Stokes equations as the governing equations for flows and the quasi-steady assumption to neglect the unsteady terms in the flow equations.The eddy viscosity is evaluated using a mixing-length approach.In this study, bed-load discharge is estimated using the Meyer-Peter and Müller formula modified as described in Wong and Parker (2006).The entrainment rate of suspended load is estimated using the relationship proposed in de Leeuw et al. (2020).See the following section for details.To test the results of linear stability analyses against the observational data of plane beds, we plotted stability diagrams in the parametric space of hydraulic parameters.

Governing Parameters
The instability of a system is illustrated as a contour diagram of the perturbation growth rate ω i (Fig. 1).Generally, theoretical studies of bedforms based on linear stability analyses describe the transition of bedform phases in the parametric space of wavenumber k and Froude number Fr, which are given by: where λ denotes the perturbation wavelength, Ũ0 is the depth-averaged flow velocity of the uniform flow, g is the gravitational acceleration (= 9.81 m 2 /s), and h0 is the flow depth of the uniform flow.Hereafter, we denote dimensional variables using a tilde (˜).
Stability diagrams described on the k-Fr plane have been commonly used to predict the development of dunes and antidunes (Kennedy, 1963).A few studies have used other combinations of dimensionless numbers such as the friction coefficient C versus Fr (Colombini and Stocchino, 2008) and the relative roughness D/ h0 on the k-Fr plane (Bohorquez et al., 2019).
Although the classic k-Fr diagrams are widely accepted, we cannot use this approach to evaluate whether plane bed formation can be predicted reliably because plane beds have extremely small wavenumber or have an infinite wavelength (i.e., they are flat).Therefore, we illustrate stability diagrams as contour maps of the maximum growth rate ω i,max of the instability on the Re p -Fr plane with fixed D and D-Fr plane with fixed Re p to investigate the impact of suspended load on the formation of plane beds.Here, D denotes the dimensionless particle diameter D/ h0 .
The instability of a system is illustrated as a contour diagram of the perturbation growth rate ω i (Fig. 1).We can rewrite Eq.
(A30) as: Thus, we can obtain the growth rate ω i as a function of k for a given combination of Fr, D, h0 .In this study, we assume that the system is stable if ω i is not positive for all k within the domain [k min , k max ] for a given Fr, D, h0 combination.In contrast, the system is assumed to be unstable if ω i is positive for some k (Fig. 1).We describe stability diagrams as contour maps of the maximum growth rate of the instability in the parametric space of (Re p , Fr) (Fig. 2) and (D, Fr) (Fig. 3).
The particle Reynolds number Re p ranges from 5.62 to 15.9 ( D = 0.125-0.25 mm).For the D-Fr diagram, we employed Re p = 5.62 and 15.9 as the fixed value of particle diameter.The dimensionless particle diameter D ranges from 5.0 × 10 −2 to 5.0 × 10 −5 in D-Fr diagram (Fig. 3).

Linear Stability Analysis
Here we present the formulation of the problem and the method used to solve the differential equations.

Flow equations
The governing equations for flows are the two-dimensional Reynolds-averaged Navier-Stokes equations.On erodible beds, the flow adjustments occur immediately relative to the bed adjustments (Fourrière et al., 2010).Therefore, we employ the quasi-steady assumption to neglect the unsteady terms in the flow equations (Colombini, 2004;Yokokawa et al., 2016).
Under the quasi-steady assumption, the dimensionless forms of the Reynolds-averaged Navier-Stokes equations and continuity equation for incompressible flow are described as: where u and w are the flow velocities in xand zdirection, respectively; p denotes the pressure; S is the bed slope; and T ij (i, j = x, z) is the Reynolds stress tensor.
We employ a Boussinesq-type assumption to close the flow equations: Then, the eddy viscosity ν T is evaluated using a mixing-length approach: where l is the mixing length, κ is the Kármán coefficient (= 0.4), h is the flow depth, Z denotes the bed height, and R is the height of the reference level at which the flow velocity is assumed to vanish in a logarithmic profile (Fig. A1).
In the above equations, the system is nondimensionalized as follows: (u, w) = (ũ, w)/ũ f0 (12) where D is the non-dimensional diameter of a bed particle, ũf0 denotes the shear velocity in the basic flat-bed state, and ρ is the water density (= 1000 kg/m 3 ).The shear velocity in the basic flat-bed state ũf0 is obtained as: As the flow is continuous, the system can be rewritten using the stream function ψ defined as: Then, Eqs. ( 4) and ( 5) are rearranged to: Eliminating p from Eqs. ( 18) and ( 19), we obtain:

Advection-diffusion equations for suspended sediment
We also assume a quasi-steady state for the advection-diffusion equation for suspended sediment, which is formulated as: Here, F x and F z are the normalized fluxes of suspended sediment in xand zdirections, respectively, given by: where c denotes the concentration of suspended sediment and w s is the settling velocity of sediment.We assume that the diffusion coefficient of suspended sediment is equal to the eddy viscosity ν T .Based on Eqs. ( 22) and ( 23), Eq. ( 21) is reformulated as: The settling velocity of sediment w s is calculated using a relationship given in Ferguson and Church (2004): where the constants C 1 and C 2 are set to the values for natural sand: C 1 = 18 and C 2 = 1.0.
The particle Reynolds number Re p is defined as: where R s is the submerged specific density and ν is the kinematic viscosity of the fluid (= 1.0 × 10 −6 m 2 /s).The submerged specific density R s is defined as: where ρs denotes the density of the bed particles (= 2650 kg/m 3 ).

Transformation of variables
We employ the following transformation of variables to apply the boundary condition at the bed and flow surfaces: The derivatives with respect to x and z are described as follows: where ∂ x denotes the partial derivative with respect to x.Using the above transformation of variables approach, the height of the water surface and the reference level correspond to η = 1 and η = 0, respectively.

Boundary condition
The boundary conditions include a vanishing flow component normal to the water surface, and vanishing stresses normal and tangential to the water surface as follows: where u = (u, w) is the velocity vector, e denotes the unit vector, and T is the stress tensor.The subscripts ns and ts denote directions normal and tangential to the water surface, respectively.
At the bed, the boundary conditions include the vanishing flow components normal and tangential to the bed.
where the subscripts nb and tb denote directions normal and tangential to the bed, respectively.The vectors e ns , e ts , e nb , and e tb , and the tensor T are defined as: The boundary conditions for the suspended sediment flux at the flow surface and bed are as follows: where F = (F x , F z ) is the flux vector of suspended sediment and Ẽs is the entrainment rate of the sediment calculated as Ẽs = ws E s .In this study, the dimensionless coefficient E s is estimated using the relationship proposed in de Leeuw et al. (2020): where C 3 was set to 5.73 × 10 −3 and coefficients e 1 , e 2 , and e 3 were set to 1.31, 1.59, and −0.86, respectively.

Basic state
The basic flow state for linear stability analysis is a uniform flow over a flat bed.Under this condition, the hydraulic parameters u, w, h, Z, R, and c are described as: where the subscript 0 denotes a parameter in the basic state.The governing equations of flows can be simplified as: with the boundary conditions: With Eqs. ( 46)-( 50), we can obtain the following logarithmic law for the flow velocity: Then, the friction coefficient C z is obtained by the direct integration of Eq. ( 51) from η = 0 to η = 1: Now, we consider the logarithmic law of the open-channel flows as: with z 0 = D/12 (Colombini, 2004).It should be noted here that the bed roughness can be modified by the sediment transport (Dietrich and Whiting, 1989).Additionally, we set the origin of z-axis at a distance of D/6 below the top of the bed particles (Fig. A2).By setting the top of the bed particles as z = D/6, the reference level R 0 is positioned below the top of bed particles.
Therefore, the domain in which the mixing-length approach cannot be applied is restricted near the bed.
Under the above uniform flow condition over a flat bed, Eq. ( 24) can be rewritten as: with the following boundary conditions: Here, c b is the near-bed concentration of suspended sediment.Under the basic state, the entrainment and deposition rates of the suspended sediment are balanced.Thus, c b is described as: where C 3 was set to 5.73 × 10 −3 and coefficients e 1 , e 2 , and e 3 were set to 1.31, 1.59, and −0.86, respectively.
By integrating Eq. ( 54), we obtain the suspended sediment distribution in the basic state as follows:

Temporal development of bed configurations
The development of the bed configuration can be described by the Exner equation considering the suspended load as follows: where λ p denotes the sediment porosity, B denotes the height of the bed-load layer, t is time, and qB denotes the bed-load discharge per unit width.In the case without suspension, the development of the bed configuration associated with suspended load is ignored by setting the coefficient α s in Eq. ( 60) to 0. In the case of the stability analysis with suspension, the coefficient α s take a value of 0 or 1 depending on the sediment transport regime (Eq.( 71)).
Equation ( 60) is nondimensionalized as: In this study, dimensionless bed-load discharge per unit width is estimated using the Meyer-Peter and Müller formula modified as described in Wong and Parker (2006); this equation is given as: where C 4 and e 4 were set to 3.97 and 1.5, respectively.Here, θ b is the Shields stress at the top of bed-load layer and θ c is the critical Shields stress for particle motion.These variables can be expressed as follows: where θ 0 is the Shields stress of the base flow, τ b denotes the shear stress at the top of the bed-load layer, θ ch denotes the critical Shields stress under the flat-bed conditions, and µ is a constant set to 0.1 (Fredsøe, 1974).The shear stress τ b is described as: where η b is the dimensionless thickness of the bed-load layer and is obtained as: where B 0 and R 0 denote the height of the top of the bed-load layer and the reference level in the basic state, respectively.
According to Colombini (2004), the thickness of the bed-load layer h b is estimated as follows: where l b denotes the relative saltation height, τ r is the shear stress at the reference level, and τ c is the critical shear stress.
In this study, the sediment transport regimes are classified using the threshold conditions of sediment motion in Brownlie (1981) as follows: The coefficients α b and α s in Eq. ( 60) were set to 0 when θ 0 < θ ch and set to 1 when θ ch ≤ θ 0 .

Linear Analysis
We impose an infinitesimal perturbation on the basic state.Then, with the use of boundary conditions, we can solve the differential equations to get the growth rate of the perturbation.Please see the appendix for details of linear analysis.

Compilation of published data
The stability diagrams were assessed using an observational dataset pertaining to open-channel flows compiled from the literature, as summarized in Tables A1-A5.We compiled from the literature a total of 56 sets of data for Fig. 2 and 59 sets of data for Fig. 3.The flow depth, the flow velocity, and particle diameter ranges from 0.0209 to 1.11 m, 0.349 to 1.66 m/s, and 0.138 to 0.32 mm, respectively.
We used the data of plane beds in which the sediment transport mode could be identified, i.e., plane bed with suspension.We identified whether sediment particles were transported as suspended load or not based on the suspended sediment concentration.
For comparison with the theoretical analysis results, we used the data of dunes and antidunes with wavenumbers with the range where T represents the water temperature in degrees Celsius.A value of 20 • C was assumed for data when T was not reported.

Re p -Fr diagram
The contour maps of ω i,max on Re p -Fr plane show that the stable region, which denotes that hydraulic conditions where the plane bed appear larger in the diagram with suspension than in that without suspension (Fig. 2).In the case of the stability analysis without suspension, a stable region do not appear when D/H = 10 −4 and the growth rate decreases with increasing Re p (Fig. 2a).For the phase diagram with D = 10 −3 and D = 10 −2 , a stable region appears at 0.8 < Fr < 1.2 (with D = 10 −3 ) and 0.7 < Fr < 0.9 (with D = 10 −2 ) (Fig. 2c, e), and the growth rate increases with increasing Re p .
The phase diagrams for the case of the stability analysis with suspension show that a stable region appear at 0.4 < Fr < 0.9-1.0(Figure 2b, d and f).Also, the growth rate at Fr > 1 of the diagram with suspension is higher than that without suspension (Figure 2a-d).In the case of the shallow flow (D = 10 −2 ), the value of growth rate do not much differ between the diagrams with and without suspension (Figure 2e and f).
Comparing the results of theoretical analysis and the observational data, all the plane bed data are within unstable region in the case without suspension (Fig. 2a and c).The analysis with suspension shows that all the plane bed data plot in the stable region when D = 10 −3 (Fig. 2d), whereas two data points out of 10 points plot in the stable region when D = 10 −4 (Fig. 2b).
Table 1 also shows that the error rate, which denotes the number of plane bed data plotted in the unstable region, is smaller in the case with suspension than that in the case without suspension.
As expected, most dune and antidune data plot in the unstable region, whereas several data points of dunes and antidunes plot in the stable region in both cases with and without suspension (Fig. 2; Table 1).

D-Fr diagram
The contour maps of the maximum growth rate on D-Fr plane also show that the stable region is larger in the diagram with suspension than in that without suspension (Fig. 3).For fine sediment, the upper limit of the stable region is smaller in the diagram with suspension than in that without suspension (Fig. 3a, b), whereas that does not much differ in the medium-sand case (Fig. 3c, d).Comparing with the observational data, most plane bed data plot in the stable region in the case of the stability analysis with suspension (Fig. 3).Also, most dune and antidune data plot in the unstable region in both cases with and without suspension (Fig. 3).The error rate for plane bed data decreased from 1 to 0.6 (Re p = 5.62) and 0.45 to 0.18 (Re p = 15.9) by adding the term of suspended load (Table 2).For dunes and antidunes, the error rate do not differ between the cases with and without suspension, except for the antidune in the case of fine sediment where r e decreases from 0.6 to 0.2 (Table 2).

Discussion
The role of suspended load in the formation of plane beds and suppressing dune-scale instabilities is quantitatively illustrated as the broadening of the stable regions (Figs. 2 and 3).The stability diagrams show a good agreement with the observational data of plane beds under flows with suspension.The transition from dunes to plane beds has been explained by the spatial lag δ between the bed topography and the local sediment transport rate (Naqshband et al., 2014;van Duin et al., 2017).If the bed topography and sediment transport rate are entirely in-phase (δ = 0), dunes migrate downstream without growth or decay.
The dune height increases and decreases when the maximum sediment transport rate occurs upstream (δ < 0) and downstream (δ > 0) of the dune crest, respectively.Kennedy (1963) introduced the spatial lag in his flow model to account for the bedform growth and decay, and subsequent research has investigated the effect of spatial lag on the bedform development (McLean, 1990;van Duin et al., 2017).Recently, Naqshband et al. (2017) quantitatively observed the positive spatial lag under suspended load dominated flows in their flume experiments.Our analyses confirm that suspended load dampens the development of bed waves, thereby facilitating the formation of plane beds, and thus cannot be neglected in theoretical analyses for realistic predictions of bedforms.
We found that dunes are deformed under flows with suspended load, although further work is needed to investigate the amplitudes of dunes under such conditions.Field surveys have indicated the existence of low-angle dunes in suspendedload dominated rivers (Smith and McLean, 1977;Kostaschuk and Villard, 1996;Hendershot et al., 2016); moreover, flume experiments have indicated that dune height decreases with increasing suspended load flux (Naqshband et al., 2017;Bradley and Venditti, 2019).Theoretical analyses in Fredsøe (1981) have also predicted a decrease of dune steepness under unsteady flows with suspension where the flow discharges were being increased.In future works, nonlinear analyses should be done to obtain the amplitudes of dunes under flows with suspended load.
Ultimately, our linear analyses provide a possible explanation for the absence of dunes in turbidites: suspended load suppresses dune formation and facilitates plane-bed formation.Previous research has suggested that the formation of dunes is suppressed due to the insufficient time for dune development (Walker, 1965), the hysteresis effect under waning flow conditions (Endo and Masuda, 1997), the turbulence suppression by high suspended-sediment concentrations (Lowe, 1988), the lack of a sharp near-bed density gradient (Arnott, 2012), and the effect of clay-sized sediment on bed rheology (Schindler et al., 2015).Although these interpretations could explain the absence of dune-scale cross-lamination in turbidites, we show that dune formation is suppressed without considering the above conditions.Although the above conditions may contribute to the deformation of dunes, instead, we propose that the development of dune-scale bed waves under turbidity currents is restricted by the presence of suspended load.The model can be improved by the inclusion of such effect in the future studies.

Conclusions
We investigated the influence of suspended load on the formation of plane beds under open-channel flows.The stability diagrams show that the stable region for finer sediments is wider in the diagram with suspension than that without suspension.
Further, the published data of plane beds with suspension coincide well with the stability diagrams where the suspension was considered.Our theoretical analysis found that suspended load promotes the formation of plane beds and suppresses the formation of dunes on the fine-grained bed.These results suggest that dune-scale cross lamination is absent in turbidites because the development of dunes in turbidity currents is restricted by the presence of suspended load.Additional theoretical work can be improved in the future studies by the inclusion of possible mechanisms for the absence of dunes in turbidites. accuracy: The above functions are substituted into Eq.(A2); then, we evaluate the equation at the Gauss-Labatte points, which are defined as: By combining the governing equations, boundary conditions, and closure assumptions, we obtain the following system of linear algebraic equations: a = (a 0 , a 1 , . . ., a N , D 1 ) (A14) M = 0, 0, 0, PR , Ľh , . . ., Ľh where a check mark (ˇ) denotes a linear operator associated with variable transformation from η to ζ.We obtain the following solution from Eq. (A12): Additionally, Eqs.(A9) and (A16) give: Similarly, we solve the eigenvalue problems for the sediment transport equations.By substituting Eq. (A1) into Eq.( 24), we obtain the following equations at the order of O(A): Table 1.Error rates for the case of fixed D. The parameter nc denotes the number of correctly classified data points and re is the error rate.

Figure 2 ,
D of the data plotted in the diagram ranges from D/3.16 to 3.16D.The data of which particle Reynolds number range from Re p /1.26 to 1.26Re p were chosen to plot Figures 3. To calculate the particle Reynolds number, the kinematic viscosity ν was assumed as follows (van den Berg and van Gelder, 1993): ν = 1.14 − 0.031 (T − 15) + 0.00068 (T − 15) 2 10 −6 (72)

Figure 1 .
Figure 1.Contour map of perturbation growth rate ωi without suspension.The particle Reynolds number and dimensionless particle diameter were set to Rep = 15.9 and D = 2.51×10 −4 , respectively.The dotted line denotes the threshold of sediment motion.The dashed lines denote the critical Froude numbers Fr cd and Frca for instabilities.The region where the growth rate is positive is highlighted in grey.

Figure 2 .
Figure 2. Contour maps of the maximum growth rate ωi,max of perturbations with a fixed dimensionless particle diameter D. Symbols are observational data.a, D = 10 −4 without suspension.b, D = 10 −4 with suspension.c, D = 10 −3 without suspension.d, D = 10 −3 with suspension.e, D = 10 −2 without suspension.f, D = 10 −2 with suspension.a and b, The range of D of observational data is from 3.16 × 10 −5 to 3.16 × 10 −4 .c and d, The range of D of observational data is from 3.16 × 10 −4 to 3.16 × 10 −3 .e and f, The range of D of observational data is from 3.16 × 10 −3 to 3.16 × 10 −2 .

Figure 3 .
Figure 3. Contour maps of the maximum growth rate ωi,max of perturbations with a fixed particle Reynolds number Rep. Symbols are observational data.a, Rep = 5.62 without suspension.b, Rep = 5.62 with suspension.c, Rep = 15.9 without suspension.d, Rep = 15.9 with suspension.a and b, The range of Rep of observational data is from 4.46 to 7.0749.c and d, The range of Rep of observational data is from 12.6 to 20.

Figure A1 .
Figure A1.Conceptual diagram of the flow.The dimensionless parameters u and w are the flow velocities in xand zdirections, respectively, h is the flow depth, Z denotes the bed height, and R is the height of the reference level at which the flow velocity is assumed to vanish in a logarithmic law.

Figure A2 .
Figure A2.Conceptual diagram of the sediment bed.The origin of z-direction is denoted by O.The parameter D is the dimensionless diameter of a bed particle, B0 is the height of the top of the bed-load layer in the basic state, and R0 is the height of reference level in the basic state.
Summary of data used for the stability diagram with D = 10 −4 Summary of data used for the stability diagram with D = 10 −3 Summary of data used for the stability diagram with D = 10 −2

Table 2 .
Error rates for the case of fixed Rep.The parameter nc denotes the number of correctly classified data points and re is the error rate.

Table A1 .
Summary of data used for the stability diagram with Rep = 5.62.

Table A2 .
Summary of data used for the stability diagram with Rep = 15.9.