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

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 time-averaged 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 Reynolds-averaged 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 (Cao et al. 2002; 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 entrainment\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.

The $\alpha$ 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., ${\theta }_0\le 0.1$ (Fernandez-Luque and Van Beek 1976), recent experiments (Lajeunesse et al. 2010) have confirmed its validity up to ${\theta }_0\approx 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:with $c_b=\sqrt{\tan {\mu }_d}{w}_s/u_{*c}$.

Moreover, Sekine and Kikkawa (1992), presenting a deterministic-probabilistic 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 readswith the dimensionless coefficient $c_a^{\prime }$ 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 single-parameter 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 solid-phase average velocity in the bed load layer is assumed to bewith $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_{\text{MPM}}$). $K_{\text{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_{\text{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\le c_1\le 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_{\text{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_{\text{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 $\alpha$:

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 ${\tau }_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 $\alpha$ 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$ asand the concentration $C_{s,b}$ aswith $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., $({\theta }_0-{\theta }_c)\gg {\theta }_c$, as those corresponding to sheet-flow regime, Eq. (35) can be approximated asand 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 (${\mu }_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 (${\mu }_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, ${\theta }_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 ${\theta }_0\approx 1$.

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 readsin which, denoting with $U$ and $V$ the $x$ and $y$ components of velocity vector for both phases, the unknowns’ vector $\mathbf{W}$ isand the $\mathbf{C}$, $\mathbf{A}$, $\mathbf{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 matrixwith $n_x$ and $n_y$ the director cosines of an arbitrary direction in the ($x,y$) plane of the unitary vector $\mathit{n}$. The eigenvalues readin which the derivative of the dimensionless bed load layer thickness with respect to ${\delta }_{s,b}$ has the following expression:

Accounting for (47) eigenvalues ${\lambda }_{\mathrm{7,8}}$ may be equivalently rewritten as follows:

From Eqs. (46) and (48) it follows that, independently of the $\mathbf{n}$ unitary vector, the matrix $\mathbf{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:in whichand

Vector $\mathbf{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 $\mathbf{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 $k$th side of the cell, ${\mathbf{n}}_k$ is the normal vector and ${\mathbf{H}}_k$ is the average value of the flux on the same side, defined asbeing ${\mathbf{F}}^{\prime }$ and ${\mathbf{G}}^{\prime }$ the vectors of the numerical fluxes, modified as follows to include the slope terms:$\tilde{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 predictor-corrector (McCormack) scheme

The numerical fluxes at the interfaces are computed by a three-point 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)].

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.

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.