A simple non-equilibrium bedload transport equation for the formation of dune in a shallow-water flow over an erodible bed

In this work, we consider the long-standing problem of capturing dune formation in an erodible-bed channel at subcritical speed by using a reduced order model of depth-averaged equations. The pioneering study by Reynolds [1] showed that the standard Saint-Venant-Exner equations are unconditionally stable at subcritical Froude number. Hence, the use of depthaveraged flow equations, which are commonly used by the hydraulic community, prevents the formation of bedforms as dunes. Recently, Cañada-Pereira & Bohorquez [2] have proposed a simple sediment transport formulation able to capture the formation of dune when coupled with the Saint-Venant equations. We replace the standard Exner equation with a non-equilibrium sediment transport equation that includes the following necessary ingredients: first, a phase shift in the particle entrainment rate; second, a particle diffusivity and an eddy viscosity. Subsequently, we solve the linear stability problem of an erodiblebed channel and show that the neutral curve properly captures the bed instability both in subcritical regime (i.e. dune) and supercritical flow (i.e. antidune and roll wave). Finally, we corroborate the capabilities of the model by means of non-linear numerical simulations which reproduce the growth of dune and antidune in agreement with experiments.


Introduction
The prediction of the occurrence of bedforms and its morphology has been a case of study of great interest due to its influence in major engineering problems such as floods and the effects on the landscape in general.The main theoretical frameworks used to study this problem can be classified, according to the flow equations, into rotational formulations using Navier-Stokes equations, and its depth-integrated form known as the shallow water equations (also referred to as Saint-Venant equations).Morphodynamic models add supplementary equations that describe sediment transport processes [3].
The Navier-Stokes approach has proven to be a good option to predict with accuracy the evolution of these phenomena, whose standardised numerical scheme and closure equations provide reliable results checked with experimental data.Even though a rotational formulation seems to be a great solution, one of its major drawbacks lies in the huge computational cost that confines its use rather in supercomputers or for a simplified and low-resolution computational domain.Within this strategy, many studies have been carried out, such as Niemann et al. [4] and Van Duin et al. [5].Alternative morphodynamic models for dune formation based on the Saint-Venant equations are scarce and require an ad-hoc bed slope correction term that triggers the instability of the bed [6,7].Hence, existing morphodynamic models predict the formation of bedforms with precision at the expense of relatively high computational cost or approximately using some ad-hoc formulation of physical processes to trigger the instabilities.
An alternative approach was first introduced by Kennedy [8] and continued by other authors such as Charru et al. [9] who proposed the incorporation of a phase shift to the transport rate and the shear stress, respectively.This study is motivated by the need to set a simple and standard model that accurately and efficiently describes the instability of an erodible bed at the inception of sediment motion.Hence, we extend the one-dimensional morphodynamical model of shallow-water flow over an inclined erodible bed by Bohorquez & Ancey [10,11], which couples the Saint-Venant equation of hydraulic engineering to a bed-load model advocated by Ancey & Heyman [12], to predict the dune/antidune instability of the erodible bed as well as the roll wave instability of the free surface.The key point is the inclusion of a phase lag between the sediment entrainment rate and the depth-averaged flow velocity, as originally proposed by Cañada-Pereira & Bohorquez [2].

Formulation of the problem
Figure 1 shows a simplified representation of the problem we address and the main variables that affect the morphology of bedforms and water waves.The depth-averaged velocity is denoted by v [m s -1 ] and the flow depth is given by h(x, t) = y s − y b with y b (x, t) [m] and y s (x, t) [m] being the bed and free surface elevations, respectively.As usual, t [s] is time, x [m] is the streamwise coordinate, g [m s -2 ] represents the gravitational acceleration, and the bed slope is defined as tan θ = −∂ x y b .The governing equations of the water depth, speed and the height of the bed can be written in the dimensional form as [11,12]: The proposed equations correspond to the standard mass (1) and momentum (2) conservation equations by Saint-Venant, and a modified Exner equation (3), where the term of the divergence of sediment flux has been substituted by a non-equilibrium entrainment-deposition rate expression (−∂ x q s = κ γ − λ).The head loss induced by hydraulic resistance, τ b , is evaluated as τ b /ρ = f /8 v |v| where f is the dimensionless Darcy-Weisbach friction factor.The extra term, ∂ x (ν h ∂ x v), represents a depth-averaged Reynolds stress where ν is the eddy viscosity [m 2 s -1 ].In (3), ζ b is the dimensionless bed porosity, κ [s -1 ] and λ [m s -1 ] represent the deposition and entrainment rates, respectively.
The Saint-Venant-Exner equations ( 1)-( 3) are supplemented by an advection-diffusion equation for the particle activity, γ(x, t) [m] (solid volume of particles in motion per unit streambed area), that reads [10,13] where D u [m 2 s -1 ] denotes the particle diffusivity, and u s = β v [m s -1 ] is the mean velocity of moving particles, with the constant coefficient in the range of β ≤ 1.
The closure equation for the Darcy-Weisbach friction factor, f , is evaluated from Colebrook equation considering fully developed turbulent flow and rough regime, in which d = d/h is the grain roughness relative to flow depth.An approximate estimation of eddy viscosity is given by the mixing length equation ν = ν t h ( f v 2 /8) 1/2 with the dimensionless parameter ν t ∼ O (1).The precise dependence of particle diffusivity D u on the flow variables remains unknown but can be related to the eddy viscosity using the turbulent Schmidt number, Sc = ν/D u , which has been set up to Sc = 0.5.For the sake of simplicity, we set ν t = 1 and β = 1.Furthermore, κ is assumed constant, and λ is evaluated as established in Bohorquez & Ancey [10], Lajeunesse et al. [14] and also Charru et al. [9].

Linear stability theory
With a suitable scaling with respect to a length scale h 0 , a velocity scale v 0 , and a characteristic particle activity γ 0 , corresponding to a uniform flow, the dimensionless governing equations have been obtained by making the following change of variables: Note that the subscript 0 denotes the value of the base flow from now on.Furthermore, to simplify the notation, we introduce the equivalent bed slope of the base flow as S 0 = f 0 v 2 0 /(8 g h 0 ) and recall that γ 0 = λ 0 /κ.These base flow conditions are not limited necessarily to a sloped terrain but can be easily used in other phenomena as dam-break and glacial outburst floods under quasi-uniform flow and quasi-steady regime.
The linear perturbation equations are obtained by the standard procedure of asymptotic expansion.Hence, we decompose the perturbed variable as the sum of the background flow plus a small amplitude disturbance, i.e. (z, φ, η, u) = (− x S 0 , 1, 1, 1) + (z , φ , η , u ).Substituting the decomposition into (1)-( 4), expanding in Taylor series the friction factor as a function of the dimensionless grain roughness d0 (see (5)), and the dimensionless saturated (or equilibrium) particle activity φ eq = λ/λ 0 about the base flow, retaining only the terms of the order O( ) and removing the prime decoration, we end up with the linear perturbation equations: The derivative term resulting from the Taylor series expansion of f was neglected due to d0 1.The dimensionless parameters in ( 7)- (10), denoted below by a circumflex, are related with the dimensional quantities by means of the following scales: According to Section 2, the dimensionless deposition rate reads κ = c d (s − 1) 1/2 /(Fr d1/2 0 ), with c d = 0.1 (coefficient inferred from the particle flight time).The relative density mismatch for silt in water is typically s = 2.65.Also, the proposed entrainment rate function is given by κe = 14 κ d0 (Sh − Sh cr ).Here, Sh denotes the Shields number which represents the dimensionless shear stress: The threshold value of motion for sand and gravel in hydraulically rough turbulent flow is typically Sh cr = 0.02 (from Julien [15]).A suitable empirical relation for particle activity in uniform flow arises from the relations used for the entrainment and deposition rates, being Sh 0 the value of the Shields number corresponding to the base flow.Therefore, the derivative term resulting from Taylor expansion of φ eq in ( 9) and ( 10) yields: Finally, the solution to ( 7)-( 10) can be written as (z, φ, η, u) , where the eigenvector is denoted by T ≡ (ζ, Φ, Γ, U) T , the dimensionless wavenumber by k, and the dimensionless complex frequency by ω = ωr + i ωi (Schmidt & Henningson [16]).
A key part of the present model is the phase shift or lag in the particle entrainment rate, φ eq , with respect to the flow variables.This physical phenomenon was first introduced by Kennedy [8,17] and can be described as a lag between the flow speed used to compute the entrainment rate, φ eq (u( x − φ/ k, t )), and the actual water speed of that point, u( x, t ).In accordance with Kennedy's [8] criterion, we set −π/2 ≤ φ ≤ 0 or 3π/2 ≤ φ ≤ 2π due to the periodic property of the wave.This lag, therefore, affects the term of the Taylor Red diamonds correspond to dunes (D) collected by Cheng [18] and Bradley & Venditti [19], green circles correspond to antidunes (AD) collected by Recking et al. [20] and Carling & Shvidchenko [21], and blue squares correspond to roll-wave experiments (RW) by Brock [22].Also, the black diamond and asterisk respectively highlight the dune and antidune experiments by Naqshband et al. [23] and Simons [20] simulated in figure 3. (b) Contour lines of the constant growth rate in the {Sh 0 , k}plane for d0 = 3 × 10 −3 and experiments with similar characteristic grain size.Thin contour lines: ωi = 0.001, 0.005, 0.01, 0.02, 0.04, 0.08 for antidunes and one fourth of these values for dunes.
Evaluating the derivatives in time and space of the differential equations ( 7)- (10) for the wave term, exp[i ( k x − ω t )], we obtain the eigenproblem A • T = 0: where the matrix A can be reduced to second order by means of the relation established between Γ = f (U) and ζ = f (Φ, U) described in previous works, see Bohorquez & Ancey [10].Note that the linear solutions to the problem share the term exp{i [ k x − ( ωr + i ωi ) t ]} = exp[i( k x − ωr t )] exp( ωi t ) from which the amplification factor of the wave can be easily obtained, i.e. exp( ωi t ).The base flow is unstable (stable) when the growth rate is positive, i.e. ωi > 0 (negative ωi < 0).

Results
There are only two independent parameters that characterise the stability of the base flow when solving for the eigenvalue problem ( 15  select Sh 0 instead of Fr using (12) with u = 1.The neutral curves have been obtained by varying Sh 0 and d0 , and solving the dispersion relation (both real and imaginary part) for ω and k.The dispersion relation was obtained by setting the determinant of A to zero, i.e.D( k, ω) ≡ |A| ≡ 0. It links the complex frequency, ω, to the real wavenumber, k, and provides the bedform growth rate ( ωi ) and, consequently, the neutral curves ( ωi = 0) which determine the occurrence of roll wave, dune or antidune.
The solutions to the neutral curves lie in the parameter space {Fr, k, d0 } or {Sh 0 , k, d0 } but, for the sake of simplicity, they have been projected onto the planes {Sh 0 , d0 } and {Sh 0 , k} in accordance with the bedforms diagrams by Carling & Shvidchenko [21].On the one hand, figure 2(a) shows the values of Sh 0 for the range of wavenumbers k where the bulk of the experiments occur (i.e. 3 × 10 −1 ≤ k ≤ 5 for dunes and antidunes, and 4 × 10 −2 ≤ k ≤ 7 × 10 −2 for roll waves).On the other hand, figure 2(b) draws the neutral curve and the growth rate in the {Sh, k}-plane at the specific value of d0 = 3 × 10 −3 , where the three types of instabilities co-exist.The distinction between both bedforms has been made regarding the sign of the wave speed (ĉ = ωr / k), which is positive for downstream migrating dunes and negative for upstream migrating antidunes.The neutral curve for roll waves shares the positive wave speed but can be distinguished from the dune instability of the erodible bed due to the supercritical regime (Fr > 1) of the base flow.
In addition, the maximum and minimum values of Sh 0 in the neutral curve of the dune and antidune is represented in figure 2(a) with solid lines as function of d0 .Both curves nearly overlap for d0 > 10 −3 .Also, the minimum value of Sh 0 required for the development of roll wave is shown in solid line.Note that the curves Sh 0 ( d0 ) in figure 2(a) define four regions.The gray region corresponds to the absence of sediment motion (i.e.Sh 0 < Sh cr ).Just above the threshold of sediment motion, dune or antidune occur depending on the value of Sh 0 , as depicted by the areas in red and green, respectively.Finally, the critical value for the free surface instability defines the blue region which lies above the antidune.In the three cases, a good agreement is observed with existing experiments (symbols).
At this point, it is worth mentioning that the wavenumber and Shields number observed in experiments with dunes and antidunes correspond to the most dangerous mode predicted by the linear stability theory.As a matter of fact, figure 2(b) shows that experimental points overlap with the regions of highest growth rate.It is fair to affirm that the stability theory correctly captures the wavelength selection mechanism observed in experiments.
In order to check the reliability of the proposed results, we benchmarked the linear stability results against non-linear simulations in the two following cases: first, we replicate the conditions in one of the antidune experiments carried out by Simons [20]; second, we mimic the experiment on dune by Naqshband et al. [23].We followed the same procedure as described at length in Bohorquez & Ancey [10].We refer the reader to [10] for further details on the numerical scheme, boundary conditions and simulation set up.In both simulations we evaluated the dimensional growth rate, ωi h 0 /v 0 [s -1 ], by best fit of the bed elevation amplitude as a function of time, see figures 3(a) and 3(c).The numerical growth rate only differs from the linear theory by 11.35% for the dune experiment and 2.32% for the antidune.Therefore, the introduction of a phase lag in the model allows us to capture the bed instabilities in both subcritical and supercritical regimes with an accurate prediction of the linear growth rate.The agreement is higher for lower values of φ because the first-order term of the Taylor expansions of φ ϕ eq increases with φ.To conclude, figure 3(b) represents the evolution of the bed perturbation elevation in time for the dune, verifying the downstream migration direction as commonly observed in dunes.A different behaviour is obtained for antidunes, see figure 3(d), where the migration of the bedform occurs opposite to flow.These results have been corroborated by the linear stability theory, which shows a positive and negative value of ĉ, respectively.

Conclusions
In this paper, we have shown that the morphodynamic model by Bohorquez & Ancey [10,11] can be readily generalised to predict the formation of dunes at subcritical flow speed.To achieve this goal, one only needs to include a phase lag between the sediment entrainment rate and the depth-averaged flow velocity, as proposed by Cañada-Pereira & Bohorquez [2].We have formalised this great result by means of a simple linear stability theory.At the same time, the use of non-linear numerical simulations of the original model equations allowed us to verify the accuracy of the linear stability theory in realistic experimental conditions at the laboratory scale.Overall, for the first time, a morphodynamic model based on the shallowwater equations can predict the dune and antidune formation observed in one-dimensional flow experiments and natural phenomena.
The nondimensional groups controlling the erodible bed instability are the Shields (or Froude) numbers, Sh (Fr), the grain roughness relative to the flow depth, d0 , and dimensionless wavenumber, k.At subcritical speed, the particle diffusivity and the phase lag between erosion rate and flow velocity trigger destabilisation of the flat bed and the growth of dune while, for supercritical flow, the particle diffusivity provokes a wavelength selection of antidune in-phase with water wave.These arguments contrast with the common belief that the only model able to predict the dune-antidune instability is based on rotational or Navier-Stokes flow equation.
We conclude that erosion-deposition morphodynamic models based on the advectiondiffusion equation (Balmforth & Vakil [24]; Ancey & Heyman [12]; Bohorquez & Ancey [10]) outperform the traditional Exner equation to study pattern formation.Taking into account the analogy of the morphodynamic model ( 1)-( 4) and well-known non-equilibrium sediment transport equations, see Bohorquez & Ancey [11], we point out that the inclusion of particle diffusivity and phase lag effects can be used to readily extend the capabilities of existing models [3].Furthermore, depth-averaged equations offer the new possibility of simulating bedforms in large-scale geophysical flows as real rivers, which is particularly useful for practitioners of hydraulic engineering.

E3SFigure 1 .
Figure 1.Sketch of an open channel flow over an erodible bed.

Figure 2 .
Figure 2. (a) Unstable region in the plane {Sh 0 , d0 } for values in the range of k observed in experiments.Red diamonds correspond to dunes (D) collected by Cheng[18] and Bradley & Venditti[19], green circles correspond to antidunes (AD) collected by Recking et al.[20] and Carling & Shvidchenko[21], and blue squares correspond to roll-wave experiments (RW) by Brock[22].Also, the black diamond and asterisk respectively highlight the dune and antidune experiments by Naqshband et al.[23] and Simons[20] simulated in figure3.(b) Contour lines of the constant growth rate in the {Sh 0 , k}plane for d0 = 3 × 10 −3 and experiments with similar characteristic grain size.Thin contour lines: ωi = 0.001, 0.005, 0.01, 0.02, 0.04, 0.08 for antidunes and one fourth of these values for dunes.