Two-Dimensional Two-Phase Depth-Integrated Model for Transients over Mobile Bed

: Fast geomorphic transients may involve complex scenarios of sediment transport, occurring near the bottom as bed load (i.e., saltating, sliding, and rolling) or as suspended load in the upper portion of the flow. The two sediment transport modalities may even coexist or alternate each other during the same event, especially when the shear stress varies considerably. Modeling these processes is therefore a challenging task, for which the usual representation of the flow as a mixture may result in being unsatisfactory. In the present paper, a new two-phase depth-averaged model is presented that accounts for variable sediment concentration in both bed and suspended loads. Distinct phase velocities are considered for bed load, whereas the slip velocity between the two phases is neglected in the suspended load. It is shown that the resulting two-phase model is hyperbolic, and the analytical expression of the eigenvalues is provided. The entrainment/ deposition of sediment between the bottom and the bed load layer is based on a modified van Rijn transport parameter, whereas for the suspended sediment a first-order exchange law is considered. A numerical finite-volume method is used for the simulation of three dam break experiments found in the literature, which are effectively reproduced in terms of both free surface elevation and bottom deformation, confirming the key role played by the solid concentration variability even for two-phase models. DOI: 10.1061/(ASCE)HY.1943-7900 .0001024. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://


Introduction
Morphological evolution in river, estuarine, and tidal environments involves the interplay of fluid flow, sediment transport, and loose bed deformation. During extreme events, such as flash floods, avalanche-induced flood waves, debris flows, or dam collapses, the above processes may evolve with comparable time scales. The resulting morphological evolution may lead to dramatic consequences in terms of damages and loss of human lives (Brooks and Lawrence 1999). Analysis and prediction of these fast morphological transients are therefore mandatory for hazard assessment (Sturm 2013). The present paper aims to contribute to this field by presenting a two-phase depth-integrated model suitable for fast unsteady flows, involving sediment transport and bed deformation.
During unsteady morphological processes, the sediment entrained from the bed is transported through bed load and suspended load. The former occurs under moderate bottom shear stress, whereas the latter pertains to higher bottom shear stress.
Bed load motion is strongly affected by particle-bottom and particle-particle collisions and by the drag received by the fluid. The suspended load is mainly characterized by the convection by the carrying fluid, often with negligible slip velocity and particle contact. In the presence of a strong spatial and/or temporal variability of the bed shear stress, the two transport modalities may coexist or alternate each other.
Experimental modeling of fast geomorphic transients encounters strong difficulties. In fact, high-resolution measurements in both time and space of flow field, sediment transport, and bottom deformation are tremendously expensive and beyond the capabilities of most laboratories. With the growing availability of computational resources, the mathematical modeling of these processes is becoming an increasingly interesting alternative for practitioners and researchers.
The present study follows a deterministic approach, describing the main features of the sediment transport in terms of timeaveraged flow properties only. This approach has the great advantage that the sediment dynamics may be analyzed without a detailed knowledge of the whole process at the price of losing some information concerning the turbulence dynamics. Although this approach is the most used in engineering applications, different analyses have been alternatively developed accounting for the turbulence effect on the sediment transport. For instance, starting from experimental evidence and following a stochastic approach, Papanicolaou et al. (2002a) developed a theoretical model for the inception of sediment motion, accounting for near-bed turbulent structures and bed microtopography. Wu and Chou (2003), incorporating the probabilistic features of the turbulent fluctuations and of the bed-grain geometry, investigated the probability of rolling and lifting for the sediment entrainment. Cheng (2006) showed that the mobility probability of a bed particle may be either enhanced or weakened by an increase of the shear stress fluctuation. In the case of low sediment entrainment, the mobility probability is increased by the turbulence, whereas it is reduced by the shear stress fluctuation if the average bed shear stress becomes relatively high. Wong et al. (2007) designed a detailed experiment to predict the probability density function for the particle virtual velocity and the thickness of the active layer, showing that the statistics of tracer displacements can be related to the macroscopic aspects of the bed load. Furbish et al. (2012) provided a probabilistic definition of the bed load sediment flux. Their formulation is shown to be consistent with experimental measurements and simulations of particle motion. Additionally, either numerical solution of the Reynoldsaveraged Navier-Stokes (e.g., Duran et al. 2012;Marsooli and Wu 2015) equations or of the direct and large eddy simulations (e.g., Keylock et al. 2005;Soldati and Marchioli 2012) of the turbulent flows coupled with sediment particle motion provided useful insights into the role of the coherent structures on erosion/ deposition dynamics.
In the following only depth-integrated models are considered and discussed. These models do not explicitly account for the dynamics of the very near-bed zone, i.e., the roughness layer. In such a layer, since the flow around sediment particles is strongly three-dimensional and influenced by wakes shed by grains, the velocity profile can significantly deviate from the logarithmic one (Byrd and Furbish 2000;Wohl and Thompson 2000). Considering that the mixing from wakes shed by particles induces a change in the eddy viscosity (Lopez and Garcia 1996;Nikora and Goring 2000;Defina and Bixio 2005). Lamb et al. (2008) assumed a mixing length proportional to the roughness height and derived a parabolic velocity profile, instead of a logarithmic one, in the layer.
Depth-integrated models may be further distinguished between coupled and decoupled ones. In the coupled models it is assumed that the sediment transport and the bottom evolution develop synchronously (Cao and Carling 2002). On the other hand, decoupled models assume a time-scale hierarchy by which hydrodynamics is usually considered to be faster than the sediment transport and the bottom evolution.
Common examples of decoupled models are those built up by supplementing a proper fixed-bed hydrodynamic model with a sediment continuity equation (the so-called Exner equation). In the simplest formulation (Graf 1998), the solid discharge is further assumed to instantaneously adapt to the transport capacity, which is estimated by means of empirical relationships proposed for uniform flow conditions (Graf 1998;Wang and Wu 2005). In many real situations this hierarchy is not acceptable, and the application of these models is questionable. Limitations of the decoupled approach have been discussed in the literature Garegnani et al. 2011) along with the drawbacks of models based on immediate adaptation of the solid discharge to the transport capacity (Simpson and Castelltort 2006;Di Cristo et al. 2006;Singh et al. 2004;Xia et al. 2010).
Among the existing coupled (i.e., nonequilibrium) morphological models, a further distinction arises from the representation of the fluid-sediment motion. They may be classified either as mixture or two-phase models, which is the type used herein. To highlight the features of two-phase models, it is useful to first discuss the more popular mixture models. For relatively low solid concentrations, the rheological behavior of the mixture may be represented through clear-water friction law (Wu 2007;Wu and Wang 2007;Sabbagh-Yazdi and Jamshidi 2013). As far as hyperconcentrated mud flows are considered, non-Newtonian constitutive relations able to describe the shear-thinning behavior of the flow are used in the model based on full (Ancey 2012) or simplified (Di Cristo et al. 2014b, c, 2015 wave dynamics. The description of a stratified flow with clear water above the mixture leads to the two-layer models, with equal (Fraccarollo and Capart 2002) or distinct (Capart and Young 2002;Li et al. 2013) velocities in the layers. However, within the transport layer no distinction is made between the motion regime of sediments and water. The interaction between mixture and clear-water layers is expressed through an interface shear stress based on the analogy with the multilayer shallow water models. Furthermore, most of these models (Capart and Young 2002;Savary and Zech 2007;Swartenbroekx et al. 2013) assume constant sediment concentration in the transport layer. These models are effective in the analysis of fast morphological transients (Spinewine et al. 2007;Chen and Peng 2006), but the assumption of constant concentration under highly unsteady conditions has been recently questioned. Li et al. (2013) suggested that sediment concentration has to be considered as one of the unknowns of the numerical model, proposing an enhanced two-layer formulation through the application of the fundamental mass conservation law for sediment. Their numerical tests support the conclusion that bed load concentration variability has to be taken into account, if a detailed description of the sediment routing is sought. The mixture models lack any explicit representation of the features of different transport regimes, i.e., bed load and suspended load, which are comprehensively lumped with the behavior of the mixture layer. Furthermore, in these models a hyperbolicity loss may occur in both subcritical and supercritical flow regimes (Savary and Zech 2007;Greco et al. 2008b;Savary and Zech 2008).
Two-phase modeling is an effective alternative for analyzing the morphohydrodynamics of rivers, debris flows, and snow avalanches (Armanini 2013). Usually, these models are deduced by averaging the conservation principles of mass and momentum for the liquid-solid mixture considered as an equivalent continuous fluid characterized by unique physical characteristics and a unique velocity value, obtaining a phase-averaged system of equations with an unknown variable concentration (e.g., Dewals et al. 2011;Canelas et al. 2013). The system of partial differential equations is hyperbolic, and it may be solved through standard finite volume schemes (Garegnani et al. 2011;Rosatti and Begnudelli 2013). Alternatively, Greco et al. (2012a) proposed a two-phase model that separately considers the liquid and solid phases, accounting for the difference between their velocities and preserving the hyperbolic nature of the system (Evangelista et al. 2013). However, in Greco et al. (2012a) the hypothesis of a constant bed load concentration has been assumed and the suspended load has not been considered. Recent research suggests that these two assumptions should be reconsidered. Indeed, the results by Li et al. (2013), even if referred to mixture models, suggest that the hypothesis of constant bed load concentration may represent a strong limitation. On the other hand, Zhang et al. (2013) recommend that the simulation of both bed load and suspended load may be required to analyze transients with a wide range of shear stress.
In the present paper a two-phase depth-integrated model is proposed, which is an extension of the preliminary version presented at the River Flow International Conference (Di Cristo et al. 2014a). The model accounts for both the bed and suspended load. As far as the former is concerned, both the liquid-solid velocities difference and the concentration variability are considered. The suspended load is still described assuming the concentration variability, but neglecting the slip velocity between the two phases. The entrain-ment\deposition of sediments between the bottom and the bed load is evaluated by a formula based on the modified van Rijn mobility parameter, whereas a diffusive vertical flux is assumed to drive the sediments toward the upper region of flow, where the suspended sediment transport occurs. The model is numerically integrated using a finite volume method, and its performance is tested against literature experimental test cases, reporting also the comparison with other existing models.
The paper is structured in the following way: the proposed model is presented in the next section. In the first subsection the governing equations are given, whereas the closures, the model mathematical characterization, i.e., its hyperbolic nature and a concise presentation of the numerical model are reported in the last two subsections. Then, the results of the model in reproducing experimental data are presented, along with a comparison with other literature models. Finally, the conclusions are drawn.

Governing Equations
In the proposed two-phase model the following hypotheses are assumed: • The liquid (ρ l ) and solid (ρ s ) densities are constant; • The sediment is uniformly graded (with diameter d) and noncohesive; • There is no inflow/outflow from sidewalls and free-surface; and • Standing bed is saturated with a porosity p.
In the depth-integrated framework, the following shallow-water assumptions are also considered: the vertical components of both acceleration and velocity are neglected; the hydrostatic pressure distribution along the vertical axis is assumed. While these conditions are not strictly verified in the near field of fast geomorphic transients (e.g., during the first instants and in the tip region of a dam break), shallow-water depth-integrated models are widely applied for simulating such events (e.g., Soares-Frazão et al. 2012;Li et al. 2013). In addition, it is supposed that the volume concentration, C s;b , along the vertical axis of the bed load region is constant and that the suspended sediment passively follows the motion of the fluid phase (Greco et al. 2012b).
The bed load dynamics is described considering separately the liquid and solid phase, with distinct velocities and accounting for the momentum exchange between them, instead of assuming an equivalent homogeneous fluid with a unique velocity value, i.e., as a water-sediment mixture. Similarly to most of the geophysical flow models (e.g., Pitman and Le 2005;Pudasaini et al. 2005;Pelanti et al. 2008), the lift and virtual (added) mass forces are neglected. As far as the latter force is concerned, Pudasaini (2012) has shown that its introduction in a two-phase model produces a strong coupling (in both time and space) between the streamwise and cross-stream velocity components in the differential terms. However, the inclusion of this force allows only a slight improvement of the model performance in predicting fast processes. On the other side, it has been shown that this additional term, modifying the differential structure of the model, may cause a loss of hyperbolicity and therefore the mathematical well posedness of the system equations is not guaranteed.
The governing equations, reported in the following, derive from the mass and momentum conservation for the liquid phase [Eqs.
(2) and (5)]. Eq. (3) represents the mass conservation of sediment moving as suspended load. Since it is assumed that the sediment velocity is equal to the liquid one in the region where suspended transport occurs, there is no drag between the two phases and therefore the momentum conservation equation for the suspended sediment is not needed. Finally, Eq. (6) is the equation for predicting bed deformation. The complete set of equations reads in which t = time; g = gravity acceleration; r ¼ ðρ s − ρ l Þ=ρ l ; and h ¼ z w − z B , where z w and z B are the free surface and bottom elevation, respectively. In Eqs.
(1)-(5), δ l denotes the liquid-phase volume for unit bottom surface, δ s;b (resp. δ s;s ) is the solid-phase volume transported as bed (resp. suspended) load for unit bottom surface so that h ¼ δ l þ δ s;b þ δ s;s . U l (resp. U s ) is the phase-averaged water (resp. solid) velocity vector, e B is the bottom erosion/ deposition rate, and e s;b−s is the sediment mass exchange between bed and suspended load. The second-order tensor U l U l (resp. U s U s ) represents the diadic product of the phaseaveraged water (resp. solid) velocity with itself. Finally, denoting with D the stress due to drag exchanged between the two phases, the source terms of momentum equations S l and S s;b are in which τ B;l and τ B;s are the bottom shear stresses on the liquid and the solid phases, respectively. The drag force of the water on the solid particles, D, is evaluated as where C D = bulk drag coefficient. The shear stress acting on the solid phase τ B;s is expressed as in which μ d is the dynamic friction coefficient. Eq. (10) accounts for both frictional, expressed through Mohr-Coulomb law, and interparticle collisional (Bagnold 1956) stresses. Following Seminara et al. (2002), the shear stress on the liquid phase is evaluated by the following relation: where s B = bottom slope. The first term is evaluated by means of the Chezy uniform flow formula, C Ch being the dimensionless Chezy coefficient. The bottom entrainment/deposition is expressed through the following formula proposed by Pontillo et al. (2010): in which w s denotes the sediment settling velocity and C s;b is the bed load concentration. The dimensionless mobility parameter T accounts for the excess of the mobilizing stresses onto the bottom surface with respect to the resisting ones (van Rijn 1984). A large number of experiments have shown that the settling velocity reduces as the particle concentration increases. The following semiempirical formula (Richardson and Zaki 1954) is therefore considered to evaluate the sediment settling velocity: in which w t = terminal settling velocity of a single particle in an indefinite fluid. According to Baldock et al. (2004), the exponent n is about 2.5 for particles with diameter of 1 mm, whereas it increases up to 5 for smaller sediments.
The mobility parameter T is herein defined as where τ c = threshold shear stress for sediment motion and jτ B j ¼ μ s rgδ s;b is the Mohr-Coulomb stress at the bottom, with μ s the static friction coefficient. Under clear-water conditions, Eq. (12) states that the erosion rate scales with the 3=2 power of the van Rijn transport parameter, which is consistent with van Rijn findings (van Rijn 1984). The solid exchange between the bed and suspended load is modeled through a first-order kinetic law (Wu et al. 2000) in which C s;s represents the depth-averaged suspended sediment concentration, C Ã s;s is the corresponding capacity value. The exchange is modulated by β and ω coefficients: the former relates the depth-averaged values to the local ones; the latter expresses the adaptation of suspended load and it is usually assumed as the sediment settling velocity (i.e., ω ¼ w s ), as it is also done herein. The expression proposed by Armanini and Di Silvio (1988) is used to evaluate β.
The capacity value for suspended sediment concentration is estimated through the following formula proposed by Wu et al. (2000) and Wu (2007): where θ 0 ¼ τ 0 =ðρ l gdrÞ is the Shields parameter computed through the modulus of the shear stress τ 0 at the bed without considering the transport layer, and θ c is the corresponding threshold value for the sediment transport initiation.

Model Closures
The α and C D coefficients may be estimated from existing empirical formulas (e.g., Maude and Whitmore 1958), which however introduce other parameters. As an alternative, in the present paper both coefficients are evaluated based on the analysis of uniform flow conditions. To this aim, the model is first applied to a uniform flow characterized by a bottom slope s B . In such a condition, the two-phase conservation Eqs.
(1)-(6) reduce to the following set of relations: Similarly to Parker et al. (2003), the following scaling law for the bed load volume for unit bottom area is assumed: with k 1 a dimensionless coefficient. Although Eq. (21) was deduced only for the low Shields parameter, i.e., θ 0 ≤ 0.1 (Fernandez-Luque and Van Beek 1976), recent experiments (Lajeunesse et al. 2010) have confirmed its validity up to θ 0 ≈ 0.2. In the present analysis, Eq. (21) is therefore applied even for a higher Shields number.
The peculiarities of the solid particles motion in the bed load, through saltation, rolling, and sliding have been thoroughly investigated through experimental studies, which have suggested that sediment velocity is different from that of the carrying fluid. Several formulas have been proposed for its evaluation, witnessing the importance of its correct computation for bed load modeling. In particular, Meland and Norrman (1966) deduced an empirical expression of the sediment average transport velocity in terms of shear velocity, roughness size, and particle diameter based on a series of experiments with glass beads rolling on a bed of homogenously sized particles. The dimensional nature of this formula limits its validity to the range of the investigated experimental conditions. Fernandez-Luque and van Beek (1976), starting from experiments carried out with a loose bed, proposed the following expression of the particles average transport velocity U p : in which u Ã = shear velocity; u Ãc = corresponding value in the Shields critical condition; and c a = dimensionless constant approximately equal to 11.5. A theoretical consideration about the dynamics of the bed load sediment transport led Bridge and Dominic (1984) to deduce the following expression for the bed grain velocity: Sekine and Kikkawa (1992), presenting a deterministicprobabilistic model to investigate the nature of the bed load motion, proposed the following expression for the bed load layer averaged mean velocity of saltation: The effectiveness of the dimensionless parameters of Sekine and Kikkawa (1992) for describing the motion of sediment particles over transitionally rough beds has been successively confirmed by Papanicolaou et al. (2002b) and Ramesh et al. (2011). Seminara et al. (2002), in deriving an entrainment-based model of sediment transport that neither satisfies nor suffers from the drawbacks of the Bagnold constraint, proposed a slight modification of the Fernandez Luque and van Beek (1976) formula, which reads with the dimensionless coefficient c 0 a ranging between 8 and 9. Recently Julien and Bounvilay (2013), based on a dimensional and regression analysis carried out considering bed load particles on smooth and rough rigid plane surfaces, proposed a simple singleparameter relation, which expresses the bed load particle velocity in terms of the shear velocity and of the logarithm of the Shields parameter of the boundary roughness.
In what follows, following Seminara et al. (2002), the solidphase average velocity in the bed load layer is assumed to be with k 2 an experimental dimensionless coefficient. By postulating the validity of Eqs. (21) and (26), the following expression of the bed load solid discharge is deduced: Eq. (27) has the same structure of the well-known Meyer-Peter and Müller (1948) formula, which is exactly reproduced provided that the k 1 k 2 product is set equal to the Meyer-Peter and Müller coefficient (K MPM ). K MPM ranges from about 4, as indicated in the reanalysis of original Meyer-Peter and Müller's dataset described in Wong and Parker (2006), to 12, used in the numerical simulations reported in El Kadi Abderrezzak and Paquier (2011). The original and most used value of 8 (Meyer-Peter and Müller 1948) is adopted in what follows. Assuming the classical value K MPM ¼ 8, the two empirical parameters k 1 and k 2 are fixed by considering the bounds deriving by the consistency of the model, as shown in the following.
The water velocity may be computed through Chezy's law Finally, it is postulated that the shear stress acting on the liquid phase may be represented as follows: with c 1 a nonnegative parameter smaller than unity, i.e., 0 ≤ c 1 ≤ 1.
In fact, the case c 1 ¼ 0 corresponds to the Bagnold's hypothesis, i.e., the shear between fluid and bottom reduces to the critical value (Bagnold 1956). On the other hand, the condition c 1 ¼ 1 implies that the shear stress acting on the liquid phase equals the corresponding value in absence of sediment transport, i.e., no momentum is transferred to the solid phase. However-as it will be shown later-a more restrictive upper bound may be specified for it. While clear indications may be found in the literature for estimating the C Ch and K MPM coefficients in their well-defined variability ranges, the dimensionless nonnegative coefficient c 1 represents a free model parameter. In the "Results" section, classical literature values are assumed for C Ch and K MPM , while the c 1 coefficient is allowed to vary in order to investigate its influence on the model predictions.
Substituting the relations (21), (26), (28), and (29) into Eq. (17), the following expression of the drag coefficient may be easily obtained: The substitution of Eqs. (21) and (26) into the momentum equation of the solid phase in the bed load layer, Eq. (18), gives the following expression for α: Expressions (30) and (31), strictly valid only in uniform flow, are used herein also in nonuniform conditions considering the local and instantaneous values of s B and τ 0 for a fixed value of c 1 . As far as the c 1 value is concerned, Eq. (31) shows that the positivity of the α coefficient imposes the following upper bound: The considered closures suggest a way to select the value for the k 1 coefficient, which has been experimentally found to vary between 0.66 (Seminara et al. 2002) and 2.51 (Lajeunesse et al. 2010). Indeed, rewriting the transport stage parameter T as and the concentration C s;b as with K s as the ratio of the bed load layer thickness to sediment diameter, the bottom entrainment/deposition condition (19) leads to the following expression for K s : Moreover, accounting for Eqs. (21) and (35) may be equivalently rewritten in terms of the bed load volume for unit bottom area as follows: Eqs. (35) or (36) indicates that the positiveness of K s implies the following condition on k 1 : Furthermore, for sufficiently large values of the shear stress, i.e., ðθ 0 − θ c Þ ≫ θ c , as those corresponding to sheet-flow regime, Eq. (35) can be approximated as and therefore the bed load concentration asymptotically approaches the value Since the asymptotic concentration Eq. (39) cannot exceed the sediment concentration in the erodible bottom, an additional condition for the k 1 value has to be respected In what follows, the value of k 1 is evaluated as the average between the lower Eq. (37) and upper Eq. (40) bounds It is easy to verify that for common values of the porosity (p) and of the static friction coefficient (μ s ), Eq. (41) provides values for the k 1 coefficient within the range of empirical values mentioned above. Furthermore, assuming the validity of the Meyer-Peter and Müller formula, the k 2 coefficient is determined as In Fig. 1, the consistency of the above set of closures is verified by comparing the prediction of the dimensionless saltation height provided by Eq. (35), with available experimental (Lee and Hsu 1994;Nino et al. 1994;Nino and Garcia 1998;Lee et al. 2000) and numerical (Wiberg and Smith 1985) results. Since unfortunately the considered references do not specify the values of porosity and of the static friction coefficient, Eqs. (35) and (41) have been applied considering two reasonable pairs of (μ s , p), namely, (0.5, 0.6) and (1.0, 0.4). On the other hand, accordingly with the values provided for the dimensionless threshold shear stress in the reference data, θ c has been assumed equal to 0.03 [ Fig. 1(a)] in the comparison with data of Lee and Hsu (1994) and Wiberg and Smith (1985), and equal to 0.06 in the comparison with data of Lee et al. (2000), Nino and Garcia (1998) and Nino et al. (1994), [ Fig. 1(b)]. Fig. 1 shows that Eqs. (35) and (41) provide relatively accurate predictions of the bed load layer thickness up to values of the Shields parameter order of unity. The fairly good agreement justifies the use of the relation (21) for the sediment volume for unit bottom area in combination with the entrainment formulation proposed by Pontillo et al. (2010) up to θ 0 ≈ 1.

Model Properties and Numerical Method
In order to show the hyperbolic character of the presented flow model, system Eqs. (1)-(6) is rewritten in quasi-linear form. Accounting for Eqs. (34) and (36) and without considering the source terms, it reads in which, denoting with U and V the x and y components of velocity vector for both phases, the unknowns' vector W is ð44Þ and the C, A, B, matrices may be easily deduced from Eqs. (1)-(6), through standard algebra. Following Courant and Hilbert (1961), the mathematical character of system (43) is investigated by looking for the eigenvalues of the matrix with n x and n y the director cosines of an arbitrary direction in the (x; y) plane of the unitary vector n. The eigenvalues read in which the derivative of the dimensionless bed load layer thickness with respect to δ s;b has the following expression: Accounting for (47) eigenvalues λ 7,8 may be equivalently rewritten as follows: From Eqs. (46) and (48) it follows that, independently of the n unitary vector, the matrix M possesses only real eigenvalues. Therefore, the present two-phase model is always hyperbolic, and the characteristics theory allows to define the correct number of conditions on each boundary of the computational domain.
The model represented by Eqs.
(1)-(6) may be equivalently rewritten in a compact form as follows: and Vector N represents the nonconservative terms in the partial differential system, arising from the bed slope source term.
The system (49) can be solved with any of the numerical schemes commonly used for hyperbolic partial differential equations (PDEs). The finite volume solver used in (Leopardi et al. 2002;Greco et al. 2012a) has been adapted to solve the PDEs of the two-phase model, along with an appropriate treatment of the bed slope source term N (Valiani and Begnudelli 2006;Greco et al. 2008a). To this aim, with reference to a structured rectangular mesh Eq. (49) is written in the following semidiscrete conservative form In Eq. (52), the overbar denotes the averaging over the computational cell of area A 0 ; l k is the length of the kth side of the cell, n k is the normal vector and H k is the average value of the flux on the same side, defined as being F 0 and G 0 the vectors of the numerical fluxes, modified as follows to include the slope terms: z is the bed elevation at the side of the cell opposite the one on which flux has to be evaluated; the terms in the square bracket are considered null if negative (Greco et al. 2008a).

Time integration of Eq. (52) is performed with a predictorcorrector (McCormack) schemē
The numerical fluxes at the interfaces are computed by a threepoint parabolic interpolation of the conserved variables values. In the predictor stage, two cells on a side of the interface and one on the opposite side are considered, vice versa in the corrector stage. The numerical stability of the proposed method is guaranteed provided that the Courant-Friedrichs-Lewy condition is satisfied for the largest eigenvalue [Eq. (46)].

Test Cases and Results
In the next two subsections the proposed model is tested against two laboratory experiments: a one-dimensional dam break, over a dry erodible bed (Capart and Young 1998), and a two-dimensional dam break, over both dry and wet bed (Soares-Frazão et al. 2012). Finally, in the last section of this paragraph, the present model is compared to four existing non-equilibrium models.

One-Dimensional Dam Break
The first test case is the fast geomorphic transient experimentally investigated by Capart and Young (1998). The experiments were carried out at National Taiwan University, and they consist of small-scale laboratory dam break of initial water depth h 0 ¼ 10 cm over an erodible bed in a prismatic rectangular channel. Notably, a very light sediment was used (density ρ s ¼ 1,048 kg m −3 ) with d ¼ 6.1 mm. Scouring propagates both upstream and downstream of the dam, where intense erosion occurs. Apart from the near-field evolution soon after the dam removal, the flood wave exhibits a rather regular shape characterized by a steep sediment-laden bore, at the front of the wave, and an enduring weak hydraulic jump at the center of the wave.
As indicated by the experimenters, the bottom porosity p has been fixed equal to 0.6, whereas the sediment free-fall velocity w t in Eq. (13) is assumed equal to 0.067 m=s. The settling velocity w s is computed through Eq. (13) at each point and time accordingly to the actual concentration value and with the n value fixed equal to 2.5. The values of the static and dynamic friction coefficients are μ s ¼ 0.52 and μ d ¼ 0.32, respectively. The dimensionless Chezy coefficient has been evaluated by Griffiths' (1981) formula for a value of the h=d ratio of about 12. The threshold Shields number was fixed at the classical value of θ c ¼ 0.047 and the Meyer-Peter and Müller coefficient (K MPM ) has been assumed equal to 8. The k 1 and k 2 coefficients have been evaluated through Eqs. (41) and (42), respectively, and their values are k 1 ¼ 1.05 and k 2 ¼ 7.62. Finally, the upper bound value of the free parameter c 1 , deduced by Eq. (32) is 0.44.
Simulations have been carried out with a grid size Δx ¼ 0.010 m and Δt ¼ 1=4; 096 s. The computational domain was sufficiently long to exclude any influence of the boundary conditions. Three different values of the c 1 parameter, namely c 1 ¼ 0, c 1 ¼ 0.2 and c 1 ¼ 0.4, have been considered. In Fig. 2 two snapshots of the experimental results from Fraccarollo and Capart (2002), corresponding to t ¼ 0.4 s and t ¼ 0.5 s after dam removal, are compared with the computed results. The numerical results show a very limited sensitivity to the c 1 value and moreover they indicate that the model predictions closely agree with the main features of the process, i.e., the celerity of the downstream tail, the free surface profile upstream and downstream the dam, and the scour of the bottom. The shape of the scour strongly resembles the experimental one, with a steep adverse slope just downstream the original dam location (x ¼ 0), followed by a nearly horizontal scoured bed. A general slight underestimation of the maximum scour occurring just upstream the bore is however observed at t ¼ 0.4 s. The observed weak hydraulic jump is also qualitatively reproduced in the simulations, with a bore appearing more upstream than in the experiments and with a sharper front.
As far as the sediment transport reproduction is concerned, Fig. 3(a) depicts in the space-time plane the suspended sediment discharge values q s;s ¼ δ s;s U l divided by the total solid discharge q s;tot ¼ δ s;s U l þ δ s;b U s , while Fig. 3(b) reports the space-time evolution of the ratio K s d=h. Even if in a large portion of the plane the suspended transport represents a small percentage, about 2%, of the total solid discharge, the map shows that there are some areas in which it increases up to 20%. The suspended solid discharge represents an appreciable contribution to the solid discharge only in a limited portion of the (x; t) plane, while it is absent in most of the region downstream to the original dam (i.e., x > 0), although in this region the Rouse number is less than 1 (results not shown herein). Such a result may be explained accounting for that, downstream the original dam position, the bed load thickness saturates the full flow depth [ Fig. 3(b)] and therefore the solid discharge is entirely conveyed as bed load.

Two-Dimensional Dam Break
An example of a two-dimensional fast geomorphic transient involving a wide range of the Shields parameter values is provided by the experiments carried out within the NSF-Pire project (Soares-Frazão et al. 2012).
The tests concern dam break waves expanding over a flat mobile bed, in a 3.6 m wide, 36 m long flume, whose geometry is reported in Fig. 4. The breached dam is represented by two impervious blocks and a 1.0 m wide gate located between the blocks. The sudden rise of the gate induces a flood wave expanding along both longitudinal and transversal directions. An initial 85-mm thick layer of coarse sand was put down upon the fixed bed, from 1 m upstream to 9 m downstream the gate. Sediments were constituted of a uniformly graded sand with d ¼ 1.6110 −3 m with relative density r ¼ 1.63, with a bottom porosity p ¼ 0.42. The sediment free-fall velocity w t is 0.18 m=s. Also in this test case the settling velocity w s has been computed through Eq. (13) with n ¼ 2.5 and considering the actual concentration value. The following values of friction coefficients have been assumed μ s ¼ 0.73 and μ d ¼ 0.63. The value of the k 1 coefficient through Eq. (41) is k 1 ¼ 1.09. The threshold Shields parameter and the Meyer-Peter and Müller coefficients have been fixed equal to θ c ¼ 0.047 and K MPM ¼ 8, as in the previous test, so that k 2 ¼ 7.34. The dimensionless Chezy coefficient has been similarly evaluated using the Griffiths' formula. Here the ratio h=d is about 200. The upper bound of the c 1 parameter is 0.29.
Two configurations were experimentally investigated: (1) an initial water level of 47 cm in the upstream reservoir and no water downstream (dry-bed test); (2) an initial water level of 51 cm in the upstream reservoir and a water level of 15 cm downstream (wet-bed test). The time evolution of the water level was measured at eight gauges by means of ultrasonic probes (Fig. 4), whose location is indicated in Tables 1 and 2 for dry and wet bed test, respectively. The final topography was measured by a bottom profiler with 5 cm resolution along y. Further details about the experimental procedure may be found in the paper by Soares-Frazão et al. (2012).
Both the dry-and wet-bed experiments have been simulated by means of a non-uniform mesh of about 35,000 cells, with variable size in x and y directions. The smallest cells, used to discretize the erodible floodplain, have size Δx ¼ Δy ¼ 2.5 · 10 −2 m. The adopted time step was Δt ¼ 1=2, 048 s. Freefall has been considered at the outlet section of the flume, whereas impervious boundaries have been considered for the flume sidewalls.
With reference to Test Case 1 (dry-bed), Fig. 5 compares measured and computed time series of free-surface elevation at the gauge points, obtained with three different values of the c 1 parameter, namely, 0, 0.1 and 0.2. Measures from symmetrical gauge points are grouped on the same plot.
An estimate of the experiment reproducibility has been provided by Soares-Frazão et al. (2012) resulting in mean observed standard deviation between σ mean ¼ 0.006 ÷ 0.016 m with maximum values being between σ max ¼ 0.018 ÷ 0.032 m, depending on the considered gauge. It is noticed that in all the gauges the arrival time of the surge caused by the dam failure is well captured, along with the general trend of the free-surface elevation decay after the surge transition.
The experimental and simulated final bottom topographies for three values of the y coordinate (y ¼ 0.2 m, y ¼ 0.7 m and y ¼ 1.45 m) are compared in Fig. 6, still considering the same three different c 1 values of Fig. 5. A slight but systematic underprediction of the deposition is observed in the simulated profile. This performance appears satisfactory if the scattering between the results of different repeated experimental runs is accounted for. Indeed, Soares-Frazão et al. (2012) estimated mean and maximum standard deviation of σ mean ¼ 0.008 m and σ max ¼ 0.029 m, respectively, with the latter value referring to the most intensely scoured zone. Moreover, the results depicted in both Figs. 5 and 6 confirm the limited influence of the c 1 parameter on the results quality. Fig. 7 reports the vector plot of both water and sediment velocities at different times (t ¼ 2 s; t ¼ 5 s, t ¼ 20 s), showing the differences between the velocity fields of the two phases. In particular, the different alignment of the velocities vectors of the two phases is  evident for t ¼ 5 s, after that the flood wave impacted the sidewall and it was reflected toward the channel axis. The fluid flow is more responsive than the sediment to the impact of the wave. As far as the far-field t ¼ 20 s snapshot is considered, the sediment transport has ceased in the recirculation zone past the rigid blocks. Moreover, the symmetry of the velocity vectors respect to the longitudinal axis confirms the ability of the adopted numerical scheme to predict symmetric results. With reference to the same instants, the wide range of the Shields parameter of this flow is witnessed in Fig. 8. Finally, Fig. 9 represents the instantaneous values of C s;b for the same times of Fig. 7. At all times, a steep transversal gradient of the concentration is observed in the narrow channel between the   With reference to Test Case 2 (wet bed), Figs. 10 and 11 report the time series of the free-surface elevation at the different gauge points and of the final topography for the three longitudinal sections y ¼ 0.2, 0.7 and 1.45 m, respectively. The sensitivity respect to the c 1 parameter is also represented. The results show that the present model is able to reproduce satisfactorily even in this test the wave propagation process (Fig. 10), independently of the c 1 value. Moreover, the computed bed profile (Fig. 11) is characterized by bedforms in the scour hole with a comparable length than in the experiments, whereas the remaining of the profile is less wavy compared than the experimental one.
The vector plot of both water and sediment velocities at different instants (t ¼ 2 s, t ¼ 5 s, t ¼ 20 s) are represented in Fig. 12. As far as the direction of the liquid and solid velocity is concerned, the presence of the water downstream the dam tends to dampen the differences. On the other hand, the initial quiescent water downstream the dam obstacles the momentum diffusion, which leads to a significantly different shear stress distribution with respect to the dry-bed test case. Indeed, whereas the range of the shear stress values encountered by the flow is comparable with that of the previous test case, the spatial distribution is characterized by a more pronounced shear stress concentration in the region downstream the corner, as shown in Fig. 13.
Along with the different shear stress distribution, the wet-bed test case differs significantly from the dry bed one also for the bed load concentration distribution. To enlighten such an aspect, the C s;b distribution is represented in Fig. 14 with reference to the same instants considered for the previous case. At the first snapshot (t ¼ 2 s), in fact, spatial gradients are more pronounced than in the The nonuniform distribution evolves in time toward a more homogeneous one. In the near field, however, the capability of the present model to account for variable concentration seems fundamental for the bed load sediment routing.

Comparison with Models in the Literature
In this section, the results of present model are compared with those obtained with four different models discussed in the literature review.
The comparison concerns the main underlying assumptions of the different models, the evaluation of their specific parameters, the computational complexity (herein intended as the number of equations to be solved), along with the agreement with the experimental tests considered in the previous sections.
As detailed in the "Model Closures" section, the present model essentially contains three dimensionless parameters, i.e. C Ch ; K MPM , and c 1 . The parameters C Ch and K MPM may be evaluated based on extensive literature indications, whereas for c 1 lower and upper bounds can be estimated. As far as the computational complexity is concerned, the one-dimensional (resp. two-dimensional) form of the proposed model needs the solution of five (resp. seven) differential equations expressing conservation principles of mass and momentum. Additionally, the bed evolution equation [Eq. (6)] has to be solved, which is however computationally less expensive than the other ones.
As far as the one-dimensional test-case is concerned, the singlephase model of Wu and Wang (2007) and the two-phase model of Greco et al. (2012a) have been considered for comparison. The one-dimensional model by Wu and Wang (2007) is a single-phase mixture model, which considers both the suspended and bed load and accounts for variable bed load concentration. It is slightly less computationally expensive than the presented model, since it requires the solution of four differential equations, plus the bed evolution one. The inertia of the bed load sediment is considered through an empirical spatial lag between the actual bed load solid transport rate and the capacity value. As a consequence, in addition to the Manning coefficient, two empirical parameters defining the nonequilibrium adaptation length of total load sediment transport have to be defined. Moreover, a correction factor for the transport stage number in the van Rijn (1984) formula (k t ) is introduced. It has been shown by the authors that while the results' sensitivity to the adaptation length value was limited, the correction factor k t significantly affected the predicted erosion magnitude. The two-phase model of Greco et al. (2012a) is constituted by four conservation laws plus the bed deformation equation. The suspended sediment motion is not accounted for and the sediment concentration in the bed load is assumed to be constant. The concrete model application needs the estimation of the Chezy coefficient and of the bed load concentration. The latter has been assumed to be equal to the bed concentration (Greco et al. 2012a). Fig. 15 compares the results for the one-dimensional test of the proposed model and of the two considered literature models. Fig. 15 indicates an evident improvement of the present model with respect to that by Greco et al. (2012a). In particular, the latter model fails to reproduce the observed weak hydraulic jump, with a gradual variation of the free surface and a very different position of the downstream waterfront. A significant underestimation of the bed scour is also noted. Present results support the consideration formulated by Li et al. (2013) that the assumption of a constant bed load concentration may fail during highly unsteady flows. Conversely, the present model performs similarly to the mixture model by Wu and Wang (2007), both in terms of   (Fig. 15). Although the mixture model may appear more attractive for the lower computational complexity, the agreement in the bed erosion significantly depends on the calibrated value of the correction factor k t .
For the two-dimensional test cases, the comparison involves the single-phase model of Canelas et al. (2013) and the two-layer model of Swartenbroekx et al. (2013). The mixture two-dimensional model of Canelas et al. (2013) exhibits a much smaller computational complexity than the present one, being constituted by four conservation type laws plus the bed evolution one. Similar to the Wu and Wang (2007) model, a spatial lag between the actual bed load discharge and the equilibrium value is introduced to mimic the effects of the bed load inertia in the layer. The spatial lag is computed through an ad hoc formula that includes three additional calibration parameters fixed through a heuristic adjustment process. The computational complexity of the two-layer model of Swartenbroekx et al. (2013) is slightly smaller than that of the present model. Indeed, it is composed of six conservation equations plus the bed evolution one. Similarly to the two-phase model of Greco et al. (2012a, b), it does not account for the suspended load and the sediment concentration in the bed load is assumed constant. The sediment inertia in the bed load layer is fully described through the balance equation for the mixture momentum in the transport layer. The shear stresses between the layers are expressed through two constant friction factors, which have been determined through calibration against experimental results. Fig. 16 (resp. Fig. 18) compares the results of the present model for the two dimensional Test Case 1 (resp. Case 2) in terms of free-surface elevation with those of Canelas et al. (2013) and Swartenbroekx et al. (2013). Fig. 17 (resp. Fig. 19) is the counterpart of Fig. 16 (resp. Fig. 18) in terms of final topography. Both free-surface elevation history (Figs. 16 and 18) and final bottom  (Figs. 17 and 19) are reproduced with an accuracy comparable to that of the model by Swartenbroekx et al. (2013) and with a slight improvement with respect to the mixture model of Canelas et al. (2013), despite the proper calibration of the three additional parameters. However, all models exhibit a slight but systematic under-prediction of the experimentally observed deposition.

Conclusions
A two-phase depth-averaged model able to deal with both bed load and suspended sediment transport has been proposed. The mathematical model, based on mass and momentum conservation equations for liquid and sediment phases, accounts for variable concentration both in the bed load and in the suspended load region. The entrainment/deposition of sediments from the bed toward the bed load layer is evaluated by a formula based on a modified van Rijn mobility parameter, whereas for the exchange between bed and suspended load a first-order exchange law is considered. The adopted set of closure relations is shown to comply, under uniform conditions of flow, with several empirical scaling laws for sediment transport and to allow for relatively accurate evaluation of the bed load layer thickness up to values of the Shields parameter order of unity. Two of the three dimensionless parameters of the model, the Chezy and the Meyer-Peter and Müller formula coefficients, may be evaluated based on extensive literature indications. The third one, c 1 , is allowed to vary in a range limited by theoretically deduced lower and upper bounds.
It has been proved that the proposed model is hyperbolic and the analytical expression of the eigenvalues has been provided. A numerical method based on a finite-volume approach has been used for the simulation of three experiments concerning three different dam breaks, showing a good agreement between simulated and experimental results. The results show that accounting for the variability concentration in the two-phase formulation leads to a neat improvement of the model performance. Finally, for all test, it has been demonstrated that the value of the free parameter c 1 has only a marginal influence on the results' quality. A further confirmation of this conclusion could be obtained through future application of the model to a wider class of morphodynamic transients.