A finite element method for partially erodible bed evolution coupled with multiphase flows

This paper introduces a continuous finite element model for two phase flows coupled with evolving bedforms. The method extends the positive definite non-oscillatory finite element algorithm (NFEM) capabilities to predict sediment transport, multiphase flow, and evolution of the resulting interfaces limiting fluids, erodible sediment layers and non-erodible strata. Free surface hydrodynamics is simulated by the integration of the governing equations for incompressible multiphase flows, including the advection and reinitialization of a phase function to update the interface location. Sign preservation of the algorithm is essential both for fluid interface tracking procedure and for the landform tracking. In the first case, the property avoids unphysical overshoots along the free surface. In the second case, a positive definite thickness of the erodible layer of sediment is mandatory to account for interaction between evolving cohesionless sediment layers and rigid beds. To compute sediment layer thickness, the method incorporates a conservation law to balance evolution of the bed position with sediment mass flux, including a general template for saturated flux calculation. The framework permits the simulation from bed load dominant transport to saltation dominant transport. Fluid/terrain interface is explicitly captured by a surface tracking methodology whereby the mesh is adapted to the bedform. Numerical experiments explore several saturated sediment flow formulae for interactive evolution of dune/atmosphere system, as well as stringent dam-break type problems on channels with dry beds and partially and totally erodible beds.


Introduction
The numerical simulation of cohesionless sediment transport past evolutionary landforms is frequently needed for planetary boundary layer flow research and applications, as well as for free surface flow studies in rivers, estuaries and coastal regions.These solutions of one phase or two phase flows are partially confined by an evolving fluid/terrain interface, thus the complete method demands some specifics about interface tracking.Two families of numerical procedures are widely recognized to model fluid interfaces, those methods based on surface tracking, and those methods based on surface capturing (comprehensive reviews can be found, among profuse literature, in Refs.[1,2]).An efficient handling of the fluids/sediment interface must also include a straightforward treatment of the erodible/non-erodible interface.This new interface acts when rigid strata appear due to erosion or when the problem involves the presence of structures and obstacles.For instance, alluvial areas stemmed from deposition can compose evolving sediment layers, while concrete structures, bed rocks, armored beds or coastal defenses can compose non-erodible layers.Furthermore, interaction between rigid beds and margins and sedimentary regions often is significant in the design of structures in riverbanks such as bridges and viaducts.
The main goal of this work is the simulation of sediment transport and partially erodible bed evolution coupled with free surface hydrodynamics.The fluid flow model incorporates a fluid interface capturing method, establishing the interface implicitly by the advection of a phase function φ.A surface tracking method follows the interface between fluids and sediment bed.The latter course of action includes a mesh adaptive algorithm such that the grid evolves and fits with the bottom boundary of the flow domain (see a simple illustration in Fig. 1).Present approach only takes action on those elements containing the interface and, unlike similar procedures developed for two-phase flows (see e.g.[3,4]), neither special boundary conditions nor enriched finite element spaces are required because solid phase is deactivated.
A major difficulty tracking the erodible/non-erodible front evolution is the choice of an algorithm able to circumvent spurious residuals arising during propagation, with minimal mass and momentum imbalances.Positivity preservation property is essential for an algorithm to avoid nonphysical computed values of the layer thickness.Extensive numerical experiments prove that this property is also beneficial to decrease the oscillatory behavior of the scheme once boundary motion takes place [5].We present a continuous numerical model based on the Finite Element Method (FEM) to address current multiphase flow and partially erodible bed evolution.Equations of motion for two fluids are solved by the Non-oscillatory Finite Element Method denominated NFEM.For multiphase flow problems, the method follows some ideas of the conservative level set method [6] by reconstructing interface sharpness in a reinitialization step [7].The NFEM preserves sign by limiting element flux contributions derived from a high order algorithm and an independent low order solution.High order algorithm is the Characteristic Based Scheme (CBS), incorporating second order accuracy extension for transient advective fields [8].Low order algorithm is a first order upwind FEM, following the idea of selecting a nearly minimum added diffusion to maintain positivity.For features of the background of NFEM, see Refs.[8,9].Two main benefits of the method are of fundamental relevance for this research.First, concerning fluids interface capturing, to preserve positivity of phase function after advection and reinitialization steps, while preventing spurious momentum transfer between phases.Second, concerning erodible/non-erodible interfaces, to maintain physical bounds for the thickness of erodible stratum by a conservative procedure.Moreover, the embedded adaptive grid algorithm is optimized for a minimum number of loops each timestep to update bottom position.
To compute the evolution of the solid/fluid interface a conservation law is integrated with the fluid flow solution.The conservation law balances temporal evolution of the bed position with spatial variation of the sediment mass flux.Sediment flux is mainly ascribed to saltation of grains, hence a proper definition of a saltation flux law is needed.For many years a template to saltation flux computation has been a formula where mass flux is proportional to the cube of the friction velocity.At first proposed by Bagnold [10], the formula assumes a saturated state determined by the equilibrated total momentum transfer from the flow to the grains.Several alternate models of saltation were presented over the years.In Section 2.1 we introduce a general format for saturated flux computation, and some well-known historical and more recent laws are exemplified and referenced in this framework.An adjusted Bagnold's law to account for the critical dependence of the saltation transport on a threshold friction velocity [11], as well as variable flow direction, has given satisfactory results for aeolian transport in several relevant morphodynamics problems (see e.g.[12,13]).In recent years, advances in experimental, theoretical and numerical results show important insights to determine the dependence between saltation sediment mass flux and wall shear stress.In particular, new functional relationships can detail the characteristics of saturated transport in the limiting cases of bed load (usually in water flows) and saltation dominant process (usually in air flows) [14,15].The present model is applied to sediment transport by water and air flows, thus we scrutinize different sediment flux laws incorporated in the numerical model.Further features on the selected options are outlined in Section 2.1.
Section 2 reports the continuous governing equations for the fluids and the sediment (2.1), and the formulation of the numerical method (2.2).Section 2.3 treats relevant characteristics of the method regarding flux limiting, interface tracking, and mesh adaptation procedure.Design of numerical experiments and the corresponding results are discussed in Section 3. Remarks in Section 4 concludes the paper.

Governing equations
The continuous model for fluid flow solution is the set of Navier-Stokes equations for a Newtonian incompressible fluid with density ρ f and dynamic viscosity µ f , in Here, u(x, t) is the velocity field (x = (x l ), l = 1, D), D is the number of space dimensions, p is the pressure, g is the acceleration of gravity, τ is the viscous stress tensor and [t 0 , T ] is the time interval.The suitable boundary and initial conditions are discussed in the corresponding sections for the numerical solution and numerical experiments.
The interface between the two fluids is established by the value φ = 1/2 of a phase function φ ∈ [0, 1], such that it has zero value in the domain occupied by the fluid with ρ f = ρ 1 , and value of one in the domain occupied by the fluid with ρ f = ρ 2 .The reference densities ρ 1 and ρ 2 are taken as air and water densities, respectively.Phase function is computed by means of the advective transport equation ∂φ ∂t with boundary conditions and initial condition φ(x, t 0 ) = φ 0 (x) in Ω .
The domain Ω in R D is bounded by Γ = Γ − + Γ + , and φ and q φ are known (the latter a vector) functions (from now on overline designates known values).The inflow boundary is denoted by Γ − while Γ − q includes slip condition if suitable, Γ + = {x ∈ Γ : (u • n b ) > 0} is the outflow boundary and n b is the outward unit normal to the boundary.The evolution of the effective sediment layer is governed by the conservation law is the thickness of the erodible sediment layer.The bulk density of the sediment ρ s = ρ m (1 − λ), with ρ m and λ as the density of the grain material and the porosity (volume fraction of voids), respectively.The (vertically integrated over the sediment layer) sediment mass flux is denoted as q s , and . Flux ρ s K∇ H h represents avalanche flux.The resulting term also operates as a natural slope limiter and is computed as anisotropic diffusion, where diffusion coefficient , s C = tan α is the critical slope and α is the angle of friction; hence, K depends critically on the bottom slope.The diffusivity β is specified in terms of temporal and spatial resolution of the numerical model [12].
The sediment flux q s is essentially ascribed to saltation process.Realistic saltation flux models depend significantly on a detailed calculation of the wall shear stress field.In this work several models are used to compute bedload sediment fluxes.Among them, the simplest model is a potential formula of the type q s ∝ u ∥u∥ 4 customized to a particular experiment [16].This formula avoids a refined stress computation, but its usage is restricted to comparisons between numerics and laboratory experiments with precise specifications.
Potential formulae are frequently written in terms of the friction velocity.These formulae currently assume a saturation condition (see some relevant references in Table 1).To cover either established or novel potential formulae, we propose the following general form, where C ϵ is an empirical coefficient which in some cases also depends on ∥u * ∥, u * = u * u ∥u∥ −1 , the friction velocity u * = √ ρ −1 f τ w , while τ w denotes the wall shear stress.Normally the exponents ϵ 1 , ϵ 2 , and ϵ 3 are integers.The threshold value u t of the friction velocity, which accounts for the angle γ between local flow current direction and slope (see derivation in [12]), is given as where θ is the local slope angle, and u t0 is the threshold friction velocity for a flat bed.The format of Eq. ( 5) takes into account the critical dependence of the saltation transport on a threshold friction velocity as well as a variable flow direction.Table 1 illustrates some relevant historical and recent saturated flow transport laws adapted to the general framework introduced in Eq. (5).In case of the laws proposed in Refs.[17,18], the coefficient C ϵ depends on the shear velocity and on the threshold velocity, while α, z 1 are parameters of the model, and z m is a mean saltation height; further, d is the grain size and g = ∥g∥.
The second law implemented is the general adjusted Bagnold's law [11], indicated in Table 1 in its fourth row, where C is an empirical coefficient [11].A third option is explored.This is a quadratic law according to fifth row of Table 1 and is recommended in Ref. [14], where authors find that saturated flux scales asymptotically as u 2 * for saltation, while flux scales as u 3 * under water flows (sixth row of Table 1).Ref.

Numerical solution
The numerical solution of Eqs. ( 1)-( 4) is performed by the non-oscillatory finite element method designated as NFEM.The method corrects a conservative and sign-preserving predictor algorithm (typically with large diffusion error) with anti-diffusive contributions.Corrections are calculated by limiting the difference between the contributions of a high order (HO) method, and those of a predictor method (here called low order solution (LO)).
To report the flux correction procedure in a compact manner, we denote as B = { p h , u h i , φ h , h h } the finite element solution of Eqs. ( 1)-( 4).The high order corrected solution is designated as Bi at node i, of the finite element mesh Ω = ⋃ Ω j , ( j = 1, E), at time (n + 1)∆t, and is built as High order corrected solution results from updating the low order solution b at time (n + 1)∆t by the sum of corrected anti-diffusive element contributions denoted as Ã.In Eq. ( 8), Ã j is the corrected anti-diffusive contribution of element j, while A j is the difference between the contribution of the high order scheme A H O and the contribution of the low order scheme A L O corresponding to element j.Sum extends over e, the number of elements j surrounding node i; summation convention is not used.To compute the anti-diffusive corrections it is necessary to define the elementwise correcting functions c j 's.The correcting functions depend on nodal high order solution, on nodal low order solution, and on the element contribution to the node of the m variables of the problem.Range of c j 's is: 0 ≤ c jm ≤ 1, c jm ∈ c j , (m = 1, D + 1) for Eqs. ( 1) and (2), and m = 1 for Eqs.(3) and (4).For c jm = 1 the corrected solution equals the high order solution while for c jm = 0 solution equals the low order solution [9].
To impose the positivity of the corrected solution, we introduce at node i the bounds B min i , B max i (here written for a scalar equation to simplify notation).Thus if To specify B min i and B max i , bounds given by Ref. [20] (see equations (17') and (18') therein) are appropriate where all variables in Eqs. ( 9) and ( 10) are positive definite.
To introduce specifics about the computation of the correcting functions, we rewrite in compact form the corrected solution (8) at node i (see Eq. ( 15) in [8]) for a scalar quantity as Here c j is recast as A small number ς is integrated to avoid the vanishing of the denominator.A straightforward derivation of these correcting functions is reported in Ref. [8], Appendix B. Direct extension of scalar correction (m = 1) for the complete set of equations of motion (m = 1, D + 1) could result in inappropriate constraints on some of the variables (see Ref. [21]).Besides, in the computation of two fluids flow, proper synchronization of anti-diffusion is relevant to shrink oscillations in the vicinity of the contact discontinuity, in particular for small values of density ratio ρ 1 /ρ 2 .The action only on velocity components [21] is an effective approach to coordinate corrections.Hence correcting functions for a node i are defined as c i = min(c im ), where m (m = 1, D) indicates the velocity component.Whereas correcting functions are established at nodes, Eq. ( 8) requires the calculation of the functions by element j; so the corresponding c j is c j = min(c i ), ∀ nodes i ∈ j.
To complete the formulation of the method, we examine the basis of the low order and high order algorithms.The proposed low order approach is an upwind finite element monotone scheme independent of the high order procedure.For an efficient computation a median dual cell-based procedure is recommended, gathering contributions to nodes in an elementwise manner, and then preserving the finite element structure defined in the high order algorithm [8].Results in the computation of the LO solution are essentially the same by median dual finite volume method or Galerkin linear finite elements with lumped mass approximation.Hence the low order predictor method is written via the Galerkin formulation for the sake of simplicity.Now, consider the advection transport equation of a scalar B with a source term Q, and the same boundary conditions as Eq. ( 3).After time discretization, the LO scheme written as an equivalent differential equation is where k is an artificial diffusivity tensor such that k is a scalar depending on an element characteristic length in the streamline direction, and superscript n + 1/2 is defined as ( We assume a standard first order split approach to compute values of b n+1/2 , u n+1/2 , and Q n+1/2 ; thus the final low order computation of the unknown can be written as where The advection-source transport solution after time discretization for the high order solution is determined by the Characteristic Galerkin method [8].The method for the advection of a scalar B, gives Operator (a • ∇) (2) for a vector field a is Prior to the discussion on significant details about the extension of the limiting procedure for two fluids flow, phase function, reinitialization of the phase function, and sediment mass conservation (Section 2.3), it is necessary to establish the low order and the high order finite element discretization.

Finite element method
The two solutions of the equations of motion follow the same steps as a standard fractional step algorithm.After time discretization, velocity solution is where velocity at time level (n + 1)∆t is computed as the sum of a predictor velocity u * = u n + ∆u * , and a pressure correction velocity increment ∆u * * .However, LO algorithm differs fundamentally from the HO algorithm in the predictor velocity calculation.Henceforth, we introduce two steps to report the formulation for both methods.First, starting from the transport solution, we state the finite element discretization for velocity predictor step and complete the solution with the pressure solution and correction step.Second, based again on transport solution, we formulate the finite element method for phase function advection transport and for the sediment conservation law.
From now on, we introduce the notation LO( ) comprising upwind fluxes computation such that for a scalar field ψ and a vector field a.
Low and high order finite element discretizations of Eqs. ( 1)-( 2) are equipped with the finite element spaces V h i , U h i , W h , and P h , (i = 1, D) defined as where Γ q is the portion of boundary with prescribed velocity, denoted as q, and Γ p specifies the portion of the boundary with prescribed pressure denoted as p.
Low order method is defined as: where v h i L corresponds to the lumped mass matrix, L O subindex signifies low order solution, and (v, w) Ω = ∫ Ω v w dΩ .Intermediate velocity increment is calculated as where Pressure solution is nearly identical to its high order formulation (introduced in Eq. ( 23)) and is not reproduced for brevity.We refer the reader to the full description of the pressure step of the high order method for details.Furthermore, pressure step is nearly always avoided in LO calculation.Among other practical aspects of the complete method, Section 2.3 reports particulars about this cost effective circumvention.Likewise, velocity correction is not presented for the low order approach.
The discrete low order solution for phase function transport (3) is formulated by recalling final form ( 14) and upwind operator LO.We introduce the following finite element spaces for this solution and for high order solution: where r h L corresponds to lumped mass matrix.Phase updating can be computed efficiently during velocity predictor step, although limiting process for interface tracking is performed once high order solution and low order solution are attained.
To formulate the discrete low order form for the bedform's evolution (4), low order and high order finite element discretizations are stated in terms of the spaces S h and H h , defined as where Γ h is specified after Eq. ( 4), and corresponds with the solid/fluid interface surface.Solution of bedform evolution is formulated as: where s h L corresponds to lumped mass matrix, the advective velocity u s ≡ q s hρ s , and Standard first order split has been assumed to compute diffusion term coming from avalanche transport.Computation of advective velocity u s needs previous knowledge of the stress field.This step is performed before low and high order computation of the thickness h n+1 and its details are reported after the definition of the high order algorithm.High order solution for the set of Eqs. ( 1)-( 2) is attained by extending the continuous characteristic based split FEM (CBS) to incorporate preservation of second order accuracy for transient advective field condition (see Appendix A in Ref. [8]).
Predictor velocity calculation is determined by the Characteristic Galerkin method [8].First step solution is formulated with Eq. ( 15) as basis: where for vector fields v, w, c, and Ω I is the domain without elements with sides belonging to the boundary.Supplemental integration by parts of viscous term is omitted for brevity.Second step involves pressure computation from the continuity equation.Solution is formulated as: Find ( p h ) n+1 ∈ P h for all t ∈ [t 0 , T ], such that where Final pressure field is p n+1 = p n + ∆ p.Last step updates velocity field through Eq. ( 16), where ∆u * * is given by ) The solution of the discrete form for the advective transport of phase function (3) is formulated by taking into account Eq. ( 15); then the method is: Find and Equations ( 22) and ( 25) demand a priori computation of velocities at n + 1/2.For that purpose we propose a one level second order predictor-corrector procedure.To outline the procedure, Eq. ( 22) is condensed by introducing the operator G such that A predictor velocity û is then computed as where values at (n + 1/2)∆t were assumed as those corresponding at n∆t.Then, values of u n+1/2 in Eq. ( 26) are calculated as u n+1/2 = 1 2 (u n + û).Computation of u n+1/2 is performed before the phase function advection step, hence intermediate velocity field is ready to be used by Eq. ( 26) and by Eq. (25).
To attain a high order solution for the erodible interface tracking, we extend time discretization given by Eq. ( 15) for the advective-diffusion transport law (4).Solution of bedform evolution is formulated as: Find h h ∈ H h for all t ∈ [t 0 , T ], such that introducing ∂Γ + h as the portion of the bedform boundary with a prescribed sediment mass flux q s .Time integration of the anisotropic diffusion term coming from avalanche transport follows the same procedure as source Q in Eq. (15).Second order terms coming from this time integration of avalanche transport term have been omitted in Eq. ( 27) for conciseness.
Solution of the sediment/fluid interface tracking requires recovering of stress field from the updated velocity field to establish velocity u s .Friction velocity is recovered as where u ∆ = u(z ∆ ) is the velocity value taken at height z ∆ over the bedform, κ = 0.41 is the Von Karman constant, and z 0 is the equivalent roughness length.The height z ∆ is estimated in terms of mesh resolution and its value in Eq. ( 28) is assumed at time n∆t.Previous to solution (27), computation of q n+1 s and u )/2 is achieved by elementwise usage of law ( 5) and an a posteriori nodal recovery.
To finish the finite element formulation, we inspect its global conservation properties.Compact matrix forms of the high order and low order finite element discretizations are written now as respectively, where C, R H are the left hand side matrix and the right hand side for the HO model, and M ′ L , R L are a conservative diagonal matrix and the right hand side for the LO model.By combining Eqs.(29), taking into account the definition of the method specified in Eq. ( 8), global solution for the total number of elements E can be put as The finite element method produces a conservative matrix (M ′ L − C) (thus verifies the conservation condition 1 ), while all non-conservative terms of R H − R L coming from Eqs. ( 15) and ( 14) cancel each other once calculation of correction (30) is performed.Hence a conservative corrected solution is reached (see particulars in Ref. [7]).However, high-order global non-conservative residuals can still be created in the calculation of R H − R L due to practical reasons.This consequence comes from the computation of low order contribution by a median dual cell-based procedure and its subsequent balance with typical elementwise manner of finite element structure of the high order algorithm.A simple intermediate step removes completely these high order residuals.Its rationale is as follows.
Calculation of total positive and negative element contributions, M + and M − , respectively, can be straightforwardly done as where N is the total number of nodes.Then the corrected solution (11) is modified as [( min This small-scale modification removes high order mass residuals from the correction procedure, while bounds in algorithm (31) preserve nodal bounds B max i and B min i .The redistribution of anti-diffusion arising from the global correction functions is inconsequential in terms of errors in standard norms, as well as in terms of interface resolution errors (see experiments in Section 3 of Ref. [7]).The resulting mass errors are of the order of round-off values.

On flux limiting and interfaces tracking
Solution of two fluids flow equations by the correcting method requires, strictly speaking, the implicit solution of Eq. ( 23) and the implicit solution of the pressure equation of the low order procedure to build in the upcoming correction.To avoid this computational charge, a reduced procedure is recommended by performing correction only for the predictor velocity field.Given that predictor velocity field is arbitrary, ensuing steps of the split algorithm ensure conservation.Usage of the reduced method has not noticeable differences in compare with the complete algorithm, also retaining the capability to preserve interface resolution [22].Otherwise, the average extra computational cost of the abridged NFEM respect to the procedure without corrections is about 1% for incompressible flows (serial).
The limiting procedure discussed in preceding sections is not able to take action to control artificial velocity jumps across the fluids interface because customary bounds do not identify the presence of the interface.An assessment to distinguish interface location is integrated into the limiting procedure and it is fully reported in Subsection 2.5 of Ref. [22].Besides, a reinitialization step complements the limiting step to maintain resolution of the contact discontinuity [23], and to add a reduced amount of diffusion to avoid wiggles on the interface.Standard reinitialization of phase function combines artificial compression and diffusion in a non-linear advection diffusion equation (see e.g.[6]).In this work a non-linear anisotropic diffusion equation merges both effects into an anisotropic diffusivity along the lines of that proposed by the authors in Ref. [22].The resulting iterative reconstructing scheme provides a very fast convergence while positivity of solution is preserved by integrating an ad-hoc correction of diffusive element contributions [7].
Lastly, sign preservation supplied by the NFEM is essential in the computation of partially erodible bed evolution, where rigid strata can appear due to complete removing of the sediment layer or the presence of non-erodible obstacles.Once limiting process is performed after finite element computation (27), the method ensures conservation and positivity of the thickness of the sediment layer.Nonetheless, position of the interface fluid-sediment must be tracked and updated to achieve a modified grid for flow calculation.Next subsection describes this issue. 1The conservation condition for a symmetrical matrix J is defined as

Fluid-bedform interface. Mesh adaptation algorithm
We report concisely a practical grid adaptation method to treat the fluid-bedform interface tracking.Grid adaptation requires the following original mesh data: (a) coordinates of all nodes and their connectivity; (b) all edges and the pair of nodes defining each one; and (c) edges defining each element of the mesh.Data can be efficiently stored using linked lists.Grid adaptation is conducted after computation of h n+1 and bedform position is updated.The method is arranged in an efficient manner by the following three steps: 1. Loop over nodes: h i is calculated at coordinates (x i , y i ) of every node i (with vertical coordinate z i of the original mesh).Nodes are labelled as • Under the interface: z i + ς < h i , Node flag is 0 (inactive).
• Over the interface: z i − ς > h i .Node flag is 1 (active).
• Coincident with the interface: The number ς is of the order of round-off value.
2. Loop over edges: pair of nodes delimiting every edge are labelled as • Both nodes have flag 0: edge is under the bedform; edge flag is 0 (inactive).
• Both nodes have flag 1: edge is over the bedform; edge flag is 1 (active).
• Both nodes have flag 2: edge is coincident with the interface; edge label is 2 (coincident).
• One node has flag 2 and the other node has flag 0 or 1; edge flag is 0 or 1 depending on the flag of the second node.• One node has flag 0 and the other node has flag 1; edge is cut by the bedform and its flag is 3 (intersected).In this situation, coordinates of the cut-off node are calculated and stored.Nevertheless, if the intersection is close to one of the edge's nodes, we could get too small elements in comparison with those of the original mesh.To avoid this inconvenience, we define a tolerance ξ such that, if 2), we change the flag of Node 1 to label 2. Practical values for the tolerance are ξ ∈ [0.05, 0.10], giving meshes that fit the bedform without too small elements.If those special cases are found, the step 2 (loop over edges) must be run again to avoid indeterminacy of other edges containing modified nodes.

Loop over elements:
To simplify, we consider the case of two-dimensional meshes of triangles.
• The three edges have flag 0: element is not considered for the adapted mesh.
• The three edges have flag 1: element is active and its geometry is identical to that of the original mesh.
• At least one edge has flag 2: the element is active only if any node of the element has label 1.
• One edge has flag 3: the three nodes of the element must have label 0, 1 and 2, respectively (see Fig. 3).
Hence, a new triangle is formed with the node with flag 2, the node with flag 1 and the cut-off node of the edge with label 3 (blue triangle in Fig. 3).• Two edges have label 3: we can distinguish two different cases.The first one occurs when just one node is label 1; in this case a new triangle is formed with both cut-off nodes of flag 3 edges and the node of flag 1 (see Fig. 4).The second case takes place when two nodes have label 1.Then we firstly calculate the diagonals' length of the quadrilateral formed by the two nodes of label 1 and both cut-off nodes of label 3 edges (see Fig. 5).Then, two new triangles are formed by splitting the quadrilateral by its shortest diagonal, as can be seen in Fig. 5.
Finally, once activated elements are defined, the subsequent steps are the assembling of the new mesh and the prescription of the boundary condition on those edges in contact with the bedform.Pressure values at new nodes on edges with label 3 are taken from previous nodes belonging to the same edge.Errors introduced in this process are negligible as the fluid-bedform interface movement is very slow, hence new nodes are very close to the former ones.
When the time scale of flow variability is much shorter than that of the bedform evolution, grid adaptation algorithm can be done every certain number of timesteps.When "severe flow scenario" is adopted, defined either by similar time scales of fluid flow and saltation flux or by the saltation flux amplified few orders of magnitude [24], mesh adaptation in the experiments should be carried out every time step.

Channel with partially erodible bed
We reproduce the experiment reported by Struiksma [16].Here the model must consider four interfaces: air/water, water/cohesionless sediment, cohesionless sediment/non-erodible bed and water/non-erodible bed.The lab test was performed in a straight flume with an effective length of 11.5 m, an uniform width of 0.2 m, and a vertical side rigid wall of 0.5 m.Layers of uniform sand (d 90 = 0.6 mm) and non-erodible bed were prepared according to two layouts (sketches superimposed in Fig. 6).In the first arrangement the flume bottom is completely filled by the sediment, while in the second arrangement the same layer is partially replaced by a medium gravel (8-16 mm) obstacle (blue line in Fig. 6).In both cases the initial bottom level of sand has a trench of depth 0.04 m and length of 2 m (red line in Fig. 6).
The experiment imposes a steady water discharge (from left to right in Fig. 6) of 9.2 l/s with an initial water depth such that free surface is placed at z = 0.362 m.The structured two dimensional grid of triangles of dimension 17.25×0.237m has an average element size2 δ = 0.005 m.No-slip condition is imposed at the bottom; left boundary has a prescribed influx condition.Right boundary reproduces a water discharge over a weir to keep a subcritical flow in the flume, and is implemented by imposing no-slip condition up to a height of 0.323 m and outlet conditions to the remaining part of the boundary.For air ρ 1 = 1.205 kg/m 3 , µ 1 = 1.808 • 10 −5 kg/(ms), and for water ρ 2 = 998 kg/m 3 , µ 2 = 1.0 • 10 −3 kg/(ms).The ad-hoc first potential formula discussed in Section 2.1 applies to compute bedload sediment flux, q s /ρ s = C ϵ u ∥u∥ 4 , where C ϵ = 0.00049 s 4 /m 3 and u is the depth-averaged water velocity.Flux formula is scaled by a factor of 60 to transform 4 min of simulated time into 4 h of real time.Numerical tests reveal very small deviations from the strict severe flow scenario results.Simulations last 255 s and timestep ∆t = 0.001 s.During the first 15 s sediment transport is deactivated to reach an initial steady-state hydrodynamics; subsequently bedload transport law is active.
Figs. 7 and 8 summarize the series of experimental and numerical results by representing bottom level at (a) t = 1 h, (b) t = 2 h, (c) t = 3 h, and (d) t = 4 h for first and second configuration, respectively.Results are shown superimposed on the initial bed.Numerics for the first layout match satisfactorily with laboratory results (Figs. 7(a)-7(d)), particularly for t = 1 h and t = 2 h with L 2 -norm errors e(L 2 ) 1h = 6.05 • 10 −3 m and e(L 2 ) 2h = 6.86 • 10 −3 m, respectively.The sharp advancing front of the trench is accurately captured at t = 1, 2 and 3 h, whereas front at t = 4 h is not evinced by the lab experiment (Fig. 7(d)), increasing the error to e(L 2 ) 4h = 1.22 • 10 −2 m.Persistence of the sharp front is also detected in numerical results of Ref. [5].The two small hollows at t = 4 s, appearing in lab results placed at x ≈ 3 m and x ≈ 5 m, are not shown by the present model or by the model in Ref. [5].Interestingly, results for trench evolution (Fig. 7(d)) are close with those achieved in Ref. [5] by a finite element Fig. 7. Struiksma's experiment: numerical results with layout 1. Black dashed line is the initial bed profile, black solid line indicates numerical results and red dots show laboratory outputs [16].(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)model of the depth integrated equations.Here, front position is moderately better tracked mainly due to the detailed calculation of the vertical vortex dynamics in the trench.Also, position of the bedform downstream the trench is closer to laboratory outputs for t = 4 h, although somewhat overpredicted.This slight overprediction is also present in the second layout (Figs.8(c) and 8(d)), and can be ascribed to deviations in the sediment outflow through the right boundary respect to the laboratory set-up.A good agreement with measurements is also obtained for the second configuration (Fig. 8), even when the non-erodible bed becomes uncovered (Figs.8(c) and 8(d)).Errors range from e(L 2 ) 1h = 6.94 • 10 −3 m to e(L 2 ) 4h = 1.40 • 10 −2 m, where the increase comes from the aforementioned overprediction of the bedform.However, sharp front position and the silting process of the trench are properly reproduced for the second arrangement.Despite the fact that the sediment flux law has not a threshold condition for sediment onset of motion, erosion excess is apparent only in a reduced region close to the upper boundary.Similar effect can be recognized in Refs.[5,16,25].

Dam break problem for erodible bed
We examine a dam break flow in a straight channel with horizontal sedimentary bed according to the specifications of the laboratory test conducted by Spinewine and Zech [26].The problem includes three interfaces: air/water, water/erodible sediment and air/erodible sediment.Initial layout is depicted in Fig. 9. Thickness of erodible layer has the value of 0.1 m, and water height in the reservoir h u = 0.35 m.Bed material is composed by sand with d 50 = 1.82 mm, angle of friction α = 30 o , density ρ m = 2683 kg/m 3 , and porosity λ = 0.5.Simulations are performed in a two dimensional triangular mesh with average element length of δ = 0.01 m.Slip condition is imposed in all boundaries, including the bed-fluid interface to avoid the artificial intrusion of air bubbles into the water phase.Fluids parameters are those used in the previous experiment.
The vertical gate located at position x = 0 is completely removed at t = 0; final time of the tests is t = 1.25 s, and time increment ∆t = 1.0 • 10 −5 s.Two series of simulations are performed with the same set-up, but using two saturated flow transport laws.The first law is the adjusted Bagnold law, generalized for variable flow direction and critical dependence on a threshold friction velocity.This law corresponds to that specified in the fourth row of Table 1 and by Eq. (7).The second law is a cubic law proposed in Ref. [14], and is specified in the sixth row of Table 1.In both cases, aeolian sediment flux is neglected.Roughness length z 0 = d 50 /10 = 1.82 • 10 −4 m and the distance to compute the friction velocity (Eq.( 28)) is assumed as z ∆ = 10z 0 ; flat bed threshold friction velocity t0 is calculated from Ref. [14] as: u t0 = ( 0.12 ( = 0.06 m/s.Figs. 10 and 11 show numerical results by using the first and second law, respectively, for t = 0.50, 0.75, 1.00 and 1.25 s.In these figures, green and purple lines represent the free surface position obtained in laboratory and in the simulation, respectively, and the red line indicates numerical results of the bed position.The blue line indicates the interface position measured in the lab between clean water and the layer of water mixed with densely packed grains, and orange line indicates the interface position measured in the lab between the layer of water-moving grains mixture and motionless sediment grains. Erosion at weir location is somewhat overpredicted when law (7) is applied (see Figs. 10(a)-10(d)).This fact also affects the free surface position in the numerical simulation, that differs moderately from laboratory outputs.However, for x/ h u > 2, bed position is fairly well captured by the numerical model, where its position can be found between the two aforementioned intermediate interfaces.Mean L 2 -norm errors of the free surface and bedform positions are e(L 2 ) fs = 6.08 • 10 −2 and e(L 2 ) b = 2.80 • 10 −2 , respectively.More accurate results are achieved for the sediment-water location by using the cubic law (see Figs. 11(a)-11(d)).In this case, erosion at weir position is not as pronounced as with the previous law, and the interface is mainly placed in an intermediate position between the bed/mixed water and the mixed/clean water interfaces measured in laboratory.Moreover, the free surface is better captured, except in the front part.These enhanced results are evinced by measured errors, e(L 2 ) fs = 4.51 • 10 −2 and e(L 2 ) b = 1.81 • 10 −2 .Slight discrepancies observed could be caused by the artificial compression effect in the reconstruction of the interface, reducing excessively diffusion originated in the transport step of the algorithm.Even so, velocity of the dambreak wave is well estimated, as can be recognized in Figs.11(a [26] and numerical outputs by Eq. ( 7).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 11.DamBreak test, comparison between laboratory data [26] and numerical results by law of Table 1, 6th row.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 12. Layout of the experiment and dune's shape.

Evolution of an isolated dune
This experiment scrutinizes the shape evolution and velocity migration of an isolated dune.The stable shape of a dune, also called equilibrium shape, recurrently has a short life and suffers a progressive elongation leading ultimately to a blow-out [24].If the meta-stable stage has sufficient span, an approximately constant migration velocity can be computed (see e.g.[27]).The series of simulations studies a medium to severe ambient wind flow over five evolving isolated dunes, starting from sandpiles with heights h 0 = 7.5, 6.0, 3.0, 1.5, and 1.0 m.The bottom of the computational domain (Fig. 12) is defined by a flat layer (either erodible with constant thickness h b , or not removable by the flow) and by the initial sandpile, which has a symmetry axis with its origin at the middle x = x 0 Sandpile is composed by well sorted quartz sand with sediment particles of mean diameter d 50 = 0.4 mm.The shape of the initial sandpile h(r ) is defined as where λ s = h 0 /0.075, and r = |x − x 0 |.In the simulations we assume h b = 0. On the left boundary an ambient wind is prescribed by a constant velocity profile u e = (u e , 0), while on the right boundary zero pressure condition is imposed.On the upper boundary and lower (fluid/sediment interface) boundary we define slip and no-slip conditions, respectively.To reach an approximately fully developed flow to compute a realistic saturated sediment flux q s (Eq.( 5)), a non-erodible bottom is placed from the left boundary with a length of 0.67h 0 (see light-grey shaded part in Fig. 12).The numerical 2D domain has linear triangular finite elements with average element length δ ≈ h 0 /38.Resolution of the final mesh comes from a convergence analysis to give accurate results with an affordable computational cost.The remaining simulation parameters are u e = 11 m/s, ρ f = 1.205 kg/m 3 , µ f = 1.808•10 −5 kg/(ms), α = 32 • , ρ m = 2650 kg/m 3 , and λ = 0.5; therefore ρ s = 1325 kg/m 3 .The parameters related to the flow-sediment interaction are: z 0 = 0.001 m, assumed in the range of 1/30 to 1/10 of the grain nominal diameter d 50 [28], and u t = 0.22 m/s.In this test we examine two laws for the computation of sediment fluxes.Both are shown in fourth (Lettau law) and fifth (Duran law for air) rows of Table 1.To fast the computation of sediment flux, it is scaled by a time multiplier 1440 (number of minutes per day), assuming the time scale of flow variability much shorter than that for the dune evolution.
Figs. [13][14][15] show a comparison between both laws for the shape evolution of isolated dunes with heights h 0 = 7.5, 3.0 and 1.5 m.Fig. 16 displays the history h/ h 0 (t/T 0 ) for various h 0 and for both laws, where characteristic time T 0 = λ s /(2.0u e ).Fig. 17 shows the migration velocity U and length of the dune L in function of the initial height h 0 , also for both laws.
For initial heights h 0 equal or greater than 6.0 m the dune remains stable for a long period of time (Fig. 16).In these cases, results computed by Lettau law show a slight amplification followed by a reduction of the same order, and the resulting dune ends up having approximately the same initial height.However, results computed by the Duran law have a more sensible amplification that remains for the complete meta-stable stage.For h 0 smaller than 3.0 m the dune disappears at around half of the simulation time while this process occurs earlier for the   Lettau law.In particular, for h 0 < 1.5 m the dune becomes very unstable and vanishes suddenly.Therefore, shape stabilization is reached for sandpiles of initial heights bigger than the range between 3.0 and 6.0 m, depending on the law for the computation of the sediment fluxes.Application of the Duran law exhibits solutions closer to those obtained previously by 3D simulations in Ref. [24], while height of unstable dunes are overpredicted by both laws.Vortices with axes parallel to the main flow direction, main source of the formation of the barchan dune horns, can add stability to the complete landform.Distribution of these vortices cannot be captured by the two-dimensional simulation, and this is one relevant reason to the overprediction of unstable shapes.[10] where c is a constant, although this constant value is considerably different for each law.Value of 35 m 2 /day is the result of the application of the Duran law.Instead, Lettau law produces values very close to those reported in Ref. [24] (see Figure 5); celerity of the largest dune, U = 1.92 m/d for h 0 = 7.5 m, is around half of that reported in [29].While data in Ref. [27] shows that small and large dunes may not be homothetic, the self-similar initial conditions assumed in this work provide a reasonable explanation for the disparity between results and observations (see a detailed discussion in Ref. [24]).The relation L(h 0 ) is linear for both laws (red solid lines in Fig. 17).These results are consistent with data in [30] (see Figs. 3 and 4) and with Ref. [24].Fig. 18 shows the shape evolution of the initial sandpile of h 0 = 7.5 m each 100 h, along the time interval [0, 500] h.Once the sandpile begins to migrate, it losses its vertical symmetry and starts to have similar shape characteristics as a typical barchan dune (Fig. 18(b)).Fig. 18(b) also shows small transitory depositions of sediment at the bottom of the lee side of the dune (see also Fig. 15).Sediment accumulation is a result of the combination of the sudden reduction of flow velocity and initiation of the downstream vortex distribution.These ripples are strongly attenuated once landform shape reaches a meta-stable state; this state can be appreciated at t ≈ 200 h (Fig. 18(c)).

Conclusions
The model proposed combines an advanced multiphase flow solution with differential bedform evolution laws and with a grid adaptation algorithm.The coupling permits simulations of one phase and free surface flows over erodible beds as well as over beds constituted by a combination of cohesionless and non-removable materials.The formation of multiple fronts and layers of fluid and/or sediment during the computation could result in spurious residuals and false mass and momentum exchanges.These unphysical residuals were drastically shrank by the application of a conservative positivity preserving scheme (NFEM) for fluids interface tracking, for fluid/sediment interface, or for erodible/non-erodible front evolution.Although usage of NFEM in recent years to simulate flows over partially erodible beds has been remarkably successful, previous implementation was focused on the depth integrated fluid flow equations coupled with depth integrated sediment transport conservation laws.Here, integration of the complete equations of motion permits a detailed computation of vertical dynamics.Hence formation of local vortexes and more realistic distribution of stresses facilitate a more elaborated representation of the sediment and fluid layers dynamics.
Numerical experiments explore extreme conditions of interfaces and fronts evolutions.In the case of a channel with partially erodible bed, the complex dynamics produced by the existence of a deep trench and of a medium

Fig. 2 .
Fig. 2. Edge (black line) cut by the bedform (red line).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .Fig. 4 .
Fig. 3. Case with one edge of label 3. Red line indicates the bedform position.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .Fig. 6 .
Fig. 5. Case with two edges with label 3 and two nodes with label 1. Red line indicates the bedform position.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .Fig. 9 .
Fig. 8. Struiksma's experiment: numerical results with layout 2. Black dashed line and blue line are the initial erodible and non-erodible bed profiles, respectively, black solid line indicates numerical results and red dots show laboratory outputs [16].(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) )-11(d) by comparing experimental and numerical results of the surge front position.

Fig. 10 .
Fig.10.DamBreak test, comparison between laboratory data[26] and numerical outputs by Eq. (7).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .
Fig. 16.Dune height in function of dimensionless time t/T 0 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) . The computed values of U with both laws fits the classic relation U = cH −1

Fig. 17 .
Fig. 17.Isolated dune length L and velocity U as function of its initial height h 0 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .
Fig. 18.Dune: 7.5 height (Lettau law).Velocity field.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Some saturated sediment flow formulae.