Modelling Intense Combined Load Transport in Open Channel

Abstract

Granular flow is modelled under the following conditions: Steady-state uniform turbulent open-channel solid–liquid flow carrying combined load at high solids concentration above a plane mobile bed. In the combined load, a portion of transported particles is transported as collisional bed load and the rest as suspended load supported by carrier turbulence. In our modelling approach, we consider one-dimensional flow and take into account a layered structure of the flow with the intense combined load. Principles of kinetic theory of granular flow are employed together with the mixing-length theory of flow turbulence in order to predict distributions of solids concentration and velocity in sediment-water flow of the given flow depth and longitudinal slope in an open channel. Components of the model are tested and calibrated by results of our laboratory experiments with lightweight sediment in a recirculating tilting flume.

1. Introduction

In an open channel with a mobile bed composed of loose sediment particles, flow of water interacts with the top of the bed and the interaction causes transport of sediment provided that a bed shear criterion for transport is satisfied [1]. In steep streams and high flow rates, bed shearing can become so intense that the induced flow rate of transported sediment becomes a considerable portion of the total flow rate of mixture of water and sediment in the channel. Depending on flow conditions and properties of sediment particles, the particles are transported either as bed load, or as suspended load or as combined load [2]. In transport of bed load (contact load), sediment particles are supported predominantly by mutual contacts, while particles transported as suspended load interact with flowing carrier instead of other particles and carrier turbulence keeps them suspended in the flow. In transport of combined load, both support mechanisms are effective, a part of transported particles is supported by contacts (often dominantly collisions) and the rest is suspended by turbulence. Due to variability of local conditions along the flow depth, the combined load transport can occur even if particles of the same size and density are transported in the flow.

Flow with intense transport of sediment tends to develop a layered structure across the flow depth as observed in laboratory experiments [3,4]. In such a flow of high-concentrated mixture, sediment particles are non-uniformly distributed across the flow depth above a typically plane top of the bed. Dominating mechanisms of particle support are different in the different layers and affect their thickness. For engineering purposes, a modelling approach to transport of sediment has been primarily empirical and produced transport correlations relating the sediment flow rate with the flow depth and the channel longitudinal slope without considering the internal structure of the flow. However, information about the structure and its exploration in modelling of the sediment transport is important to make model results more complete and physically sound.

Detailed experimental studies which include information on the internal structure of mixture flow (basically information on distributions of sediment concentration and velocity across the flow depth) with intense sediment transport over an erodible plane bed are scarce. More information is available for pressurized flows than for flows in open channels. In a pressurized pipe, the high bed shear condition at the top of stationary deposit is easy to reach due to the high hydraulic gradient i (frictional head loss over a unit length of a pipe), $i={\mathsf{\tau }}_0/\left({\mathsf{\rho }}_{f}g{R}_{ho}\right)$, where τ₀ is the bed shear stress, ρ_f is the fluid density, g is the gravitational acceleration and R_h0 is the hydraulic radius of the bed associated area of a flow cross section [5]. An ability of flow to transport a solid particle is evaluated using the bed Shields parameter θ₀, which the ratio of the shear force exerted by fluid on the particle and the submerged weight of the particle, ${\mathsf{\theta }}_0={\mathsf{\tau }}_0/\left({\mathsf{\rho }}_s-{\mathsf{\rho }}_f\right)gd$, where ρ_s is the solids density and d is the particle size. In pressurized flows, transport of solids is typically very intense, with the maximum volumetric flow rate of solids exceeding 25 per cent of the total flow rate of mixture [6]. For laboratory pipes, various measuring techniques to sense distributions in a pipe cross section have been available, for instance, a radiometric profiler for a concentration distribution [5,6,7]. Results of pipe experiments were used primarily to evaluate the friction condition at the top of the erodible bed at high bed shear and resulted in relationships between the bed friction coefficient and the bed roughness modified to include the interaction of transported sediment with the bed surface [5,8,9]. Furthermore, the pipe experiments, which typically cover a broad range of bed shear conditions, were exploited to define different modes of sediment transport. The widely used approximate criterion, based on a pressurized-pipe experiment, requires that the suspension ratio, i.e., the ratio of the bed shear velocity, $u_{*0}=\sqrt{{\mathsf{\tau }}_0/{\mathsf{\rho }}_f}$, and the terminal settling velocity of a particle, exceeds 1 (or 1.25) to enable transport of particles in the turbulent suspension mode [8]. For lower values of the suspension ratio, the non-suspension mode (basically contact-load transport) dominates.

Flume experiments on the internal structure of flow with intense sediment transport are currently more common as appropriate measuring techniques have become available recently for open channel flows of relatively high concentrated solid–liquid mixtures. Camera-based measuring techniques rely on tracking of particles in a near-wall domain of the flow to extract their local velocities [10,11] and concentrations [12] from a sequence of flow images recorded through a glass wall of a flume. So far, experiments exploiting these techniques targeted primarily purely collisional (bed-load) granular flows [3,13,14,15] covering a relatively broad range of bed shear conditions.

An alternative measuring principle, based on the acoustic Doppler effect, is employed by a measuring instrument ACVP to measure simultaneously profiles of velocity and concentration of particles across the depth of sediment-laden flow [4]. The instrument’s high resolution and frequency enables to quantify local turbulent characteristics and an analysis of measured local parameters allows for an identification of layers in flow with intense sediment transport [4]. So far, few test runs have been reported for very similar bed shear conditions.

Modelling efforts taking the internal structure of intense sediment transport into account have intensified recently. Works focus on CFD numerical two-phase models, e.g., [16,17], and alternatively on less computationally challenging algebraic models considering one-dimensional steady uniform flow. The latter models exploit principles of kinetic theory of granular flows and are therefore very suitable for a description of collision dominated transport [3,13,15,18,19]. The kinetic theory provides constitutive formulae which can be exploited to mutually relate distributions of local stresses, concentrations and velocities of colliding particles. In the models, the constitutive relations are combined with momentum balance equations for solids and mixture to provide relations among global quantities describing flow and transport (flow rates of mixture and its phases, the flow depth and the longitudinal slope). Current kinetic-theory based contact-load models assume linear distributions of velocity and concentration across the collisional transport layer and the uniform distribution of velocity in the water layer above the transport layer [3,15,19]. If the additional dense sliding layer is considered below the collisional transport layer, then it exhibits uniform distributions of concentration and velocity [18,19]. Basically, the existing models differ in chosen forms of constitutive relations and in assumptions made for conditions at flow interfaces.

A kinetic-theory based transport model is applicable to flow with combined load transport if modifications are introduced to account for the ability of carrying liquid to suspend transported sediment particles. A principle of the turbulent support of particles assumes that the submerged weight of a burden of sediment at a certain position within a transport layer is either entirely or partially balanced by a lift force produced by a diffusive action of turbulent eddies on the particles in the layer. The lifting effect is described by a turbulent diffusion equation based on a balance of granular fluxes due to particle settling and particle dispersion as in the theory by Rouse and Schmidt (e.g., [20,21]). An evaluation of the turbulent diffusive action requires a local characterization of carrier turbulence and is done by using the Prandtl mixing-length theory in the Rouse–Schmidt concept.

So far, kinetic-theory based combined-load transport models similar to those for contact-load transport have been proposed only for pressurized flows [22,23]. They consider the transport layer divided into two parts—the upper suspended-load layer, in which particles interact exclusively with the eddies, and the lower turbulent collisional layer, in which a portion of particles is supported by turbulent eddies and the rest by collisions. To account for the lift effect, the turbulent-diffusion term is included in the momentum balance equation for the local granular normal stress at the bottom of the collisional layer. The model by Berzi and Fraccarollo [23] applies the suspension-ratio criterion from [8] to identify the top of the suspended-load layer (called fully turbulent layer by Berzi and Fraccarollo) at which the fluid shear stress has the minimum value required to suspend transported particles. The fluid shear stress exceeds this threshold stress throughout the fully turbulent layer and reaches a value equal to the flow induced total shear stress (a value of Shields parameter) at the bottom of the layer. The same total shear stress is assumed to act at all elevations between this interface and the top of the bed. In the model, the mixing length l is assumed to increase linearly with the distance from the top of stationary deposit y (l = κy, where κ is von Karman constant). The distribution of velocity is linear across the dense layer and the adjacent turbulent collisional layer and it is uniform at elevations above the top of the turbulent collisional layer. The distribution of concentration is linear in the dense layer and in the turbulent collisional layer and it is non-linear, following the Rouse–Schmidt solution, above the top of the turbulent collisional layer.

Our previous work was concerned with intense transport of contact load dominated by intergranular collisions. We built our experimental database of results for contact-load transport of lightweight sediments in a laboratory tilting flume [24,25,26] and used it to test constitutive relations of the classical kinetic theory and the extended kinetic theory at conditions typical for intense contact-load transport [25]. Selected constitutive relations based on the classical kinetic theory were employed in a collisional transport model for intense bed load in open-channel flow [27]. For given slope and depth of the flow, the model predicts the total flow rate and the flow rate of sediment while considering the internal structure of the flow described by linear vertical distributions of velocity and concentration across the collisional layer. In [28], we tested which criteria are suitable to assess whether collisional contact is a dominating particle-support mechanism at various elevations within a transport layer of our observed flows. We exploited experimental data at flow conditions where the transition between predominantly contact load and predominantly suspended load were most likely to occur. For selected tests, experimentally observed distributions of velocity and concentration enabled to evaluate relevant turbulence-associated quantities, as the mixing length, at various vertical positions in the transport layer. A distribution of the mixing length was extracted from the experimental data using two alternative approaches. One approach isolated the local mixing length from the Rouse–Schmidt equation and the other from a Karman-constant containing formula suitable for a transport-layer condition [4]. The latter equation includes the turbulence damping effect quantified using Richardson number. It modifies the local mixing length in flow exhibiting steep gradients of local velocity and concentration, which is a condition typical for a transport layer. The mixing length from the turbulent-diffusion equation is interpreted as the local mixing length required to maintain particles suspended at the local position. It is assumed that this length must be bigger than the available mixing length at the position. This mixing length is related to the local flow conditions and determined by the equation with the Richardson number. Furthermore, the local mixing length should be bigger than the particle size to ensure that turbulent eddies are able to suspend the particles. It appeared that the assumptions were satisfied in the upper part of a transport layer if the bed Shields parameter exceeded say 1.5 (corresponding with the suspension ratio of about 1.2) and so showed that the threshold at which a relevant part of the transport layer exhibited turbulent suspension of particles corresponded well with the suspension-ratio criterion from [8]. This analysis provided information which could be further exploited in a development of a combined-load transport model for open-channel flow.

The main aim of the work reported in this paper is to develop a transport model for open-channel flow with intense combined load which includes predictions of profiles of velocity and concentration across the flow depth, and predictions of other relevant quantities at all interfaces of the layered structure of the sediment laden flow. The model should be simple enough to be easily applicable for engineering purposes. Moreover, the model should contribute to a further refinement of combined-load modelling options offered by the previous models designed for pressurized flows, e.g., by modifications of some simplifying assumptions regarding distributions of relevant quantities (shear stress, mixing length, velocity).

The additional major objective of this work is to calibrate model components and to validate the model by new original experimental results produced in our laboratory flume. So far, no suitable experimental data covering a required sufficiently broad range of the combined-load transport conditions and including the distributions of concentration and velocity have been made available.

2. Materials and Methods

A modelling approach is employed which enables to predict characteristics of fully developed steady uniform turbulent open-channel flow with intense transport of combined load. The approach is based on appropriate theories and considers a layered structure of the sediment-laden flow. The resulting transport model considers conditions at identified layer interfaces and conditions within individual layers to simulate mutual relations among the longitudinal slope, flow depth and distributions of velocity and concentration of sediment.

The transport model is calibrated and validated by own experiments which produced sufficiently detailed information about tested flow conditions including measured profiles of sediment concentration and velocity.

2.1. Combined-Load Model

2.1.1. Modelled Conditions

The modelled solid–liquid flow is a one-dimensional inclined gravity-driven open-channel flow with intense transport of sediment over an erodible bed with a plane surface. Above the bed, the flow is composed of layers in which different dominating mechanisms of particle support act and it results in different distributions of velocity and concentration across the different layers. Based on these assumptions, different theories and equations are applied in the individual layers and at their interfaces.

The assumed layers and interfaces are (Figure 1):

• The stationary bed (granular deposit with the plane surface expressed as the 0-boundary);
• the permanent-contact (dense sliding) layer (its top is the d-boundary);
• the combined-load (collision and turbulent suspension) layer (its top is the c-boundary);
• the suspended-load (turbulent suspension) layer (its top is the water surface, the w-boundary).

2.1.2. Model Assumptions

Specific assumptions are proposed for the individual layers and their boundaries to emphasize dominating conditions and mechanisms and at the same time to simplify the conditions so that an analytical solution of the model is feasible. Some assumptions originate from experimental experience with the modelled type of solid–liquid flow.

In the permanent-contact layer (the d-layer), the local volumetric concentration of solids c is higher than the maximum concentration to which constitutive relations of the classical kinetic theory apply. The submerged weight of transported particles is assumed to be supported solely by interparticle contact which is virtually permanent. The granular body is sheared and hence there is a velocity gradient throughout the layer which is maximum at the top of the layer and tends to zero at the bottom. The solids distribution is not uniform across the layer either with the local concentration at the bottom of the layer (the bed concentration c₀, here and below an index indicates the position at which the quantity is evaluated) being larger than the concentration at the top of the layer c_d. This lower concentration at the top is responsible for a predominantly collisional contact among particles at this interface. Hence, the support of particles is assumed to be purely collisional at the d-boundary and this boundary is the lowest elevation at which the collisional contact applies.

Above this interface, i.e., in the regions occupied by the combined-load layer and the suspended-load layer, turbulent lift contributes to balance the submerged weight of the particles in the direction perpendicular to the interface. In the combined-load layer (the c-layer), both particle support mechanisms (collision and turbulent lift) act at each elevation with the increasing contribution of turbulent suspension at higher elevations in the layer. At the top of the combined layer the collisional support vanishes. Our experiments (discussed below) suggest that distributions of both concentration and velocity are virtually linear across the layer as is the case with purely collisional transport layers investigated previously [24,25,26].

In the suspended-load layer (the s-layer) adjacent to the top of the combined-load layer, the distributions are no longer linear, and their shape is governed by turbulent diffusion which is no longer affected by colliding particles [29] as is the case in the combined-load layer. All particles are suspended by the lifting action of turbulent fluctuating velocities of carrying fluid.

Throughout all layers, the local velocities u and their gradients γ are assumed to be equal for liquid and sediment. Previous investigations justify this assumption through experimental evidence [26] and theoretical analysis [30]. The analysis showed that even a very small difference in local velocities of solids and fluid produces a sufficiently large drag force on particles in a transport layer of inertial flow.

2.1.3. Model Principles and Applied Theoretical Concepts

The modelled flow conditions require to simulate a flow behavior of solid–liquid mixture in which interactions between the two phases (solid–liquid interaction) and interactions among transported solid particles (particle–particle interaction) are equally important. The particle–particle interactions are either sporadic although intense (collisions) or permanent (sliding contact). Particle-fluid interactions include buoyancy, drag and turbulent lift. The last is active only if appropriate conditions are satisfied.

The particle–particle interactions are described using principles of the kinetic theory of granular flows. It exploits constitutive relations for local particle stresses (normal and shear) and for a balance of particle fluctuation energy resulting from interparticle collisions in a sheared granular body. The classical kinetic theory considers that particles interact exclusively by mutual binary collisions which covers collision-driven transport of particles at low to moderate local concentrations and where particles exchange momentum through uncorrelated collisions, e.g., [3,13,15]. The extended kinetic theory covers correlated motion of particles in the environment of high local concentration of particles and modifies the constitutive relations for the condition called the dense limit, e.g., [18,19]. The threshold concentration at which relations of the classical theory cease holding because particle velocity fluctuations due to collisions begin to be correlated is called the freeze concentration (c_f).

The mixing-length concept based on the Prandtl theory is exploited to solve the distributions of velocity and fluid shear stress in the suspended-load layer. In the model, the effect of the gradients of concentration and velocity on the local mixing length using Richardson number is applied as in [4]. Furthermore, the distribution of the mixing length is used to describe interactions of suspended particles with turbulent eddies for purposes of a determination of distribution of sediment concentration in the suspended-load layer using the turbulent diffusion concept with the Rouse–Schmidt formula. In high concentrated flows like flows with intense transport of sediment, the local terminal settling velocity of a particle should include a hindered settling term by Richardson and Zaki [22,23] and it is included in the model.

Furthermore, the mixing-length concept is exploited to express fluid turbulence suppression due to the presence of colliding particles [29] inside the combined-load layer. In the flow with colliding particles, the local mixing length depends on the mean distance between particles at the location and hence decreases with the increasing local concentration of solids towards the bottom of the layer. It results in diminishing of the local effective viscosity relating the turbulent fluid stress with the fluid shear rate. Furthermore, the local fluid effective viscosity has an additional granular-like component related to the exchange of momentum between the fluid shells surrounding the fluctuating particles and it can be expressed as added-mass effect [29] using the added mass term introduced to constitutive relations by [13]. In the model calculations, the suppression of turbulence diminishes the fluid shear stress at the d-boundary and affects the thickness of the d-layer.

2.1.4. Model Features

The model exploits momentum balances to relate stresses at different interfaces. Furthermore, the model exploits constitutive relations by Garzó and Dufty [31] simplified to conditions at dense limit with the freeze concentration at the d-boundary [32]. The constitutive relations take care of local stresses and other relevant parameters at the d-boundary and assist to determine the elevation of the d-boundary, y_d (y is the distance from the top of the bed in the direction perpendicular to the bed). Above the d-boundary, the applied assumptions (including experiment-based ones) enable to abandon the constitutive relations and to simulate the flow exclusively using a combination of momentum balance and particle-fluid interactions.

The model requires input information on properties of liquid (density ρ_f and dynamic viscosity μ_f) and sediment (density ρ_s, equivalent diameter d and terminal settling velocity w) and on flow (the depth H and the longitudinal slope ω). It employs constants related to constitutive relations (bed concentration c₀, concentration at d-boundary c_d, effective coefficient of restitution of particles in dry conditions ε, particulate friction coefficient at the bed β₀, constant for solids stress ratio at the d-boundary C_bd) and to turbulent-diffusivity relations (von Karman constant κ, constant for initial mixing-length position C_ini, particle diffusivity constant C_η). The model predicts elevations of the d-boundary, c-boundary and distributions of concentration and velocity as major outputs.

• Inputs: Solid/liquid properties: d, ρ_s, w, ρ_f, μ_f
• Flow: H, ω
• Constants: c₀, c_d, ε, β₀, κ, C_bd, C_ini, C_η
• Major outputs: Positions of interfaces between layers: y_d, y_c
• Distributions of velocity and solids concentration in layers: u(y), c(y)

2.1.5. Model Equations and Computational Procedure

An analytical solution of the combined-load transport model is described below in a step-by-step procedure using equations in their final form based on applied conditions. The general assumptions have already been discussed; more specific assumptions are added below where appropriate. Some described equations are the same as in previously published collisional-transport models and they are derived elsewhere. References are added to appropriate sources.

In the equations below, all quantities are dimensionless. The conversion to their dimensional counterparts using g, d, ρ_s, ρ_f and the specific gravity S = ρ_s/ρ_f is summarized in Appendix A. The notation of dimensionless quantities and corresponding dimensional quantities is the same for sake of simplicity.

The procedure starts with conditions at the d-boundary as in [19,27].

The permanent-contact (dense sliding) layer (the d-layer and its boundaries: The d-boundary, the 0-boundary).

As the first step, a value of the coefficient of wet restitution e_d (subject to iteration in the procedure) is introduced. The coefficient e differs from ε by the additional damping effect of lubrication forces on collisions expressed through the fluid dynamic viscosity, μ_f, and the granular temperature, T (a measure of local particle velocity fluctuations due to intergranular collisions) [18]. In the procedure, the introduction of e_d enables to express T_d at the d-boundary as in [19]

$$T_d={\left(\frac{62.1}{S\mathsf{\cdot }R}\mathsf{\cdot }\frac{1+\mathsf{\epsilon }}{\mathsf{\epsilon }-e_d}\right)}^2={\left(62.1\mathsf{\cdot }{\mathsf{\mu }}_f\mathsf{\cdot }\frac{1+\mathsf{\epsilon }}{\mathsf{\epsilon }-e_d}\right)}^2,$$

(1)

where $R=\frac{{\mathsf{\rho }}_f\mathsf{\cdot }d}{{\mathsf{\mu }}_f}\mathsf{\cdot }\sqrt{g\mathsf{\cdot }\frac{S-1}{S}\mathsf{\cdot }d}$ is particle Reynolds number.

It is assumed that turbulence does not suspend any particles at the d-boundary. All particles are supported by collisions and particle velocity fluctuations begin to be correlated at this boundary. Hence, constitutive relations of the extended kinetic theory apply and the one for the solids normal stress (s-index for solids) reads

$${\mathsf{\sigma }}_{sd}=2\mathsf{\cdot }{c}_d\mathsf{\cdot }{G}_d\mathsf{\cdot }\left(1+e_d\right)\mathsf{\cdot }{T}_d,$$

(2)

in which the radial distribution function at contact for dense limit [32]

$$G_d=\frac{2-c_f}{2\mathsf{\cdot }{\left(1-c_f\right)}^3}\mathsf{\cdot }\frac{c_{cr}-c_f}{c_{cr}-c_d}\mathsf{\cdot }{c}_d,$$

(3)

where c_cr is the critical concentration and c_f is the freeze concentration (we assume c_d = c_f).

The solid stress ratio β_d can be isolated from the constitutive relation expressing the balance of fluctuation energy of colliding particles for the assumption of negligible diffusion flux in the energy balance as in [19]

$${\mathsf{\beta }}_d=\sqrt{\frac{24\mathsf{\cdot }{J}_d}{5\mathsf{\cdot }\mathsf{\pi }\mathsf{\cdot }{L}_d}\mathsf{\cdot }\frac{\left(1-e_d\right)}{\left(1+e_d\right)}},$$

(4)

in which the concentration-related function by Garzó and Dufty [31] for dense limit [19]

$$J_d=\frac{1+e_d}{2}+\frac{\mathsf{\pi }}{4}\mathsf{\cdot }\frac{\left(3\mathsf{\cdot }{e}_d-1\right)\mathsf{\cdot }{\left(1+e_d\right)}^2}{24-\left(1-e_d\right)\mathsf{\cdot }\left(11-e_d\right)},$$

(5)

and the correlation length L_d = 1. For calculating the solids shear stress at the d-boundary using the stress ratio β_d, it appeared necessary to expand the relation with the empirical constant C_βd which in a simplistic way compensated for effects not included in the constitutive relation. This compensation was required to match a predicted velocity gradient with an experimentally determined velocity gradient across the combined-load layer which was considered constant across the entire layer. The effects not covered by the used constitutive relation may include the effect of the diffusion flux (neglected in the balance above) and/or the effect of inhomogeneity (non-uniform distribution of sediment and only partial collisional suspension above the interface) on the solids shear stress. Thus

$${\mathsf{\tau }}_{sd}=C_{\mathsf{\beta }d}\mathsf{\cdot }{\mathsf{\beta }}_d\mathsf{\cdot }{\mathsf{\sigma }}_{sd},$$

(6)

The local shear-induced solids shear stress, τ_s, is also related to local c and T in a constitutive relation for the shear stress which relates τ_s to the local solids shear rate γ_s (the gradient of longitudinal velocity u) [18]. This information enables to express the shear rate at the d-boundary

$${\mathsf{\gamma }}_{sd}=5\mathsf{\cdot }\sqrt{\mathsf{\pi }}\mathsf{\cdot }\frac{1+e_d}{4\mathsf{\cdot }{J}_d}\mathsf{\cdot }{C}_{\mathsf{\beta }d}\mathsf{\cdot }{\mathsf{\beta }}_d\mathsf{\cdot }\sqrt{T_d},$$

(7)

The local fluid shear stress τ_f (f-index for fluid) is related to the local fluid shear rate through the local fluid effective viscosity. Note that the shear rates are assumed to be equal for liquid and solids at any elevation in the flow (γ = γ_f = γ_s). Following [29], the local fluid effective viscosity at the d-boundary is split into two components

$${\mathsf{\mu }}_{fd}={\mathsf{\mu }}_{fd,turb}+{\mathsf{\mu }}_{fd,gran},$$

(8)

The granular-like component is proportional to the particle viscosity,${\mathsf{\tau }}_{sd}/{\mathsf{\gamma }}_d$, and the added-mass term a_d,

$${\mathsf{\mu }}_{fd,gran}=\frac{a_d}{S}\mathsf{\cdot }\frac{{\mathsf{\tau }}_{sd}}{{\mathsf{\gamma }}_d},$$

(9)

where from [13]

$$a_d=\frac{1+2\mathsf{\cdot }{c}_d}{2\mathsf{\cdot }\left(1-c_d\right)},$$

(10)

It is assumed that the turbulent component of the local effective viscosity, associated with the local mixing length, is negligible at the d-boundary, ${\mathsf{\mu }}_{fd,turb}=0$, and therefore

$${\mathsf{\tau }}_{fd}={\mathsf{\mu }}_{fd,gran}\mathsf{\cdot }{\mathsf{\gamma }}_d$$

(11)

A relationship determining the elevation of the top of the dense sliding layer (i.e., also the thickness of the dense sliding layer) is obtained from the stress balance at the 0-boundary expressed by a combination of momentum balances and constitutive relations and it reads

$$y_d=\frac{{\mathsf{\beta }}_0\mathsf{\cdot }{\mathsf{\sigma }}_{sd}-{\mathsf{\tau }}_{sd}-{\mathsf{\tau }}_{fd}}{\frac{2+\left(c_0+c_d\right)\mathsf{\cdot }\left(S-1\right)}{2\mathsf{\cdot }\left(S-1\right)}\mathsf{\cdot }\sin \mathsf{\omega }-{\mathsf{\beta }}_0\mathsf{\cdot }\frac{c_0+c_d}{2}\mathsf{\cdot }\cos \mathsf{\omega }}$$

(12)

This thickness is required together with the information about the distribution of the shear rate across the d-layer (assumingly linear) and the shear rate at the bottom of the layer (γ₀ = 0) to determine the velocity at the top of the layer

$$u_d=\frac{{\mathsf{\gamma }}_d+{\mathsf{\gamma }}_0}{2}\mathsf{\cdot }{y}_d=\frac{{\mathsf{\gamma }}_d}{2}\mathsf{\cdot }{y}_d,$$

(13)

At the 0-boundary, velocity is neglected (u₀ = 0) and so is the fluid shear stress (τ_f0 = 0). The solids shear stress is related to the solids normal stress through the yielding value given by the coefficient β₀. The normal stress is obtained from the momentum balance

$${\mathsf{\sigma }}_{s0}={\mathsf{\sigma }}_{sd}+\frac{c_0+c_d}{2}\mathsf{\cdot }{y}_d\mathsf{\cdot }\cos \mathsf{\omega },$$

(14)

and thus

$${\mathsf{\tau }}_{s0}={\mathsf{\beta }}_0\mathsf{\cdot }{\mathsf{\sigma }}_{s0},$$

(15)

The resulting bed Shields parameter

$${\mathsf{\theta }}_0={\mathsf{\tau }}_{s0}+{\mathsf{\tau }}_{f0}={\mathsf{\tau }}_{s0},$$

(16)

is compared with a value of the bed Shields parameter obtained from the flow quantities which are inputs to the model (the flow depth H and the slope ω)

$${\mathsf{\theta }}_{0,flow}=\frac{H\mathsf{\cdot }\sin \mathsf{\omega }}{S-1},$$

(17)

in which the effect of mixture density on the Shields parameter is neglected. This is justified by the low value of S of the lightweight sediment which produces low values of the specific gravity of mixture. In the procedure, the condition ${\mathsf{\theta }}_{0,flow}={\mathsf{\theta }}_0$ selects a value of e_d for which values of other quantities at the d-boundary are determined.

Finally, the profiles of velocity and concentration are calculated for the conditions of a linear distribution of the shear rate and of a uniform distribution of the concentration gradient throughout the d-layer

$$u=\frac{{\mathsf{\gamma }}_{sd}}{2\mathsf{\cdot }{y}_d}\mathsf{\cdot }{y}^2,$$

(18)

$$c=c_0+\frac{c_d-c_0}{y_d}\mathsf{\cdot }y,$$

(19)

The combined-load (collisional and suspended) layer (the c-layer, the c-boundary).

In the c-layer, a certain part of transported particles exchanges momentum through uncorrelated collisions. At the top of the c-layer, the concentration of colliding particles and the particle stresses vanish. The gradients are assumed constant across the c-layer. At the top of the c-layer, the velocity gradient ${\mathsf{\gamma }}_c={\mathsf{\gamma }}_d$, and the concentration gradient

$${\mathsf{\delta }}_c=\frac{c_c-c_d}{y_c-y_d},$$

(20)

This determines a value of Richardson number at the top of the c-layer

$$R{i}_c=-\frac{S\mathsf{\cdot }\cos \mathsf{\omega }\mathsf{\cdot }{\mathsf{\delta }}_c}{{\mathsf{\gamma }}_c^2},$$

(21)

and it is further used to identify the mixing length at the c-boundary [4]

$$l_c=\sqrt{1-c_c}\mathsf{\cdot }\left(1-C_l\mathsf{\cdot }R{i}_c\right)\mathsf{\cdot }\mathsf{\kappa }\mathsf{\cdot }\left(y_c-y_{ini}\right),$$

(22)

in which the initial position y_ini is estimated as a multiple of the thickness of the d-layer, $y_{ini}=C_{ini}\mathsf{\cdot }{y}_d$, and a value of C_ini is calibrated by experiments. The liquid shear stress at the c-boundary

$${\mathsf{\tau }}_{fc}={\mathsf{\mu }}_{fc}\mathsf{\cdot }{\mathsf{\gamma }}_c,$$

(23)

where the effective viscosity equals to the local turbulent viscosity (the granular viscosity does not apply due to the absence of colliding particles) [29],

$${\mathsf{\mu }}_{fc,turb}=\frac{1-c_c}{S}\mathsf{\cdot }{l}_c^2\mathsf{\cdot }{\mathsf{\gamma }}_c,$$

(24)

Moreover, the fluid shear stress can be expressed by the momentum balance at the c-boundary

$${\mathsf{\tau }}_{fc}=\frac{H-y_c}{S-1}\mathsf{\cdot }\sin \mathsf{\omega },$$

(25)

and a combination of the two equations for the local fluid shear stress determines the position of the top of the combined-load layer, y_c.

The local concentration at this boundary, c_c, is obtained from the turbulent diffusion balance (Schmidt–Rouse equation). Its use is justified by the assumption that all particles are suspended by turbulence at this elevation,

$$c_c=-C_{\mathsf{\eta }}\mathsf{\cdot }\frac{\sqrt{S\mathsf{\cdot }{\mathsf{\tau }}_{fc}}\mathsf{\cdot }{l}_c}{w\mathsf{\cdot }{\left(1-c_c\right)}^n}\mathsf{\cdot }{\mathsf{\delta }}_c,$$

(26)

The equation is solved iteratively for c_c. The velocity at the top of the c-layer is determined simply from the assumption of the linear distribution of velocity across the c-layer

$$u_c=u_d+{\mathsf{\gamma }}_c\mathsf{\cdot }\left(y_c-y_d\right),$$

(27)

The profiles of velocity and concentration throughout the c-layer,

$$u=u_d+{\mathsf{\gamma }}_{sd}\mathsf{\cdot }\left(y-y_d\right),$$

(28)

$$c=c_d+{\mathsf{\delta }}_c\mathsf{\cdot }\left(y-y_d\right)$$

(29)

The suspended-load (turbulent suspension) layer (the s-layer, the w-boundary).

Near the top of the s-layer, particles are often absent, and the solids concentration tends to zero. At the top of the layer, i.e., at the water surface, the liquid shear stress vanishes. The profile of velocity is obtained from the logarithmic equation using the boundary condition at the c-boundary,

$$u=y_c\mathsf{\cdot }{\mathsf{\gamma }}_c\mathsf{\cdot }\ln \left(\frac{\exp \left(\frac{u_c}{y_c\mathsf{\cdot }{\mathsf{\gamma }}_c}\right)}{y_c}\mathsf{\cdot }y\right),$$

(30)

The concentration profile is obtained from the set of equations below:

$${\mathsf{\tau }}_f=\frac{H-y}{H-y_c}\mathsf{\cdot }{\mathsf{\tau }}_{fc},$$

(31)

in which the effect of mixture density is neglected, $\mathsf{\gamma }=\frac{y_c}{y}\mathsf{\cdot }{\mathsf{\gamma }}_c$, $Ri=-\frac{S\mathsf{\cdot }\cos \mathsf{\omega }\mathsf{\cdot }\mathsf{\delta }}{{\mathsf{\gamma }}^2}$, $\mathsf{\delta }=-c\mathsf{\cdot }\frac{w\mathsf{\cdot }{\left(1-c\right)}^n}{C_{\mathsf{\eta }}\mathsf{\cdot }\sqrt{S\mathsf{\cdot }{\mathsf{\tau }}_f}\mathsf{\cdot }l}$ and $l=\sqrt{1-c}\mathsf{\cdot }\left(1-C_l\mathsf{\cdot }Ri\right)\mathsf{\cdot }\mathsf{\kappa }\mathsf{\cdot }\left(y-y_{ini}\right)$.

If these equations are combined, then they produce one implicit relationship

$$-c\mathsf{\cdot }w\mathsf{\cdot }{\left(1-c\right)}^{n-0.5}=\mathsf{\delta }\mathsf{\cdot }{C}_{\mathsf{\eta }}\mathsf{\cdot }\sqrt{S\mathsf{\cdot }{\mathsf{\tau }}_f}\mathsf{\cdot }\left(1+C_l\mathsf{\cdot }\frac{S\mathsf{\cdot }\cos \mathsf{\omega }\mathsf{\cdot }\mathsf{\delta }}{{\mathsf{\gamma }}^2}\right)\mathsf{\cdot }\mathsf{\kappa }\mathsf{\cdot }\left(y-y_{ini}\right)$$

(32)

Equation (22) is solved iteratively to get c at any elevation y inside the s-layer for estimated gradient $\mathsf{\delta }=\frac{c-c_{-1}}{y-y_{-1}}$ (c₋₁, y₋₁ being values at the previous elevation).

This completes the calculation procedure and the resulting layered structure of flow can be plotted including the positions of the boundaries of the individual layers and the distributions of velocity and solids concentrations across all layers. Those are further compared with experimentally determined profiles of velocity and concentration obtained from laboratory experiments described below.

2.2. Experimental Work

2.2.1. Experimental Set-Up

Experiments for calibration and validation of the combined-load model were carried out in a recirculating titling flume (Figure 2) at Water Engineering Laboratory of Czech Technical University in Prague. The set-up was used in several recently reported sediment-transport experiments and it was described in detail elsewhere [24,26]. The flume is 6-m long, 0.2-m wide and allows for various longitudinal slopes including steep ones. In a typical sediment-transport test, an appropriate combination of the longitudinal slope and the mixture flow rate is installed to set required flow conditions for an experimental observation. The mixture flow rate is controlled by the speed of a centrifugal pump (No. 3 in Figure 2) equipped with a variable frequency drive.

2.2.2. Measuring Techniques

The measuring techniques include a magnetic flow meter (No. 5 in Figure 2) for mixture flow rate, the set of two differential pressure transmitters for the mixture density (DPTs are installed in the vertical pipes of the circuit) and a set of ultrasonic water level sensors along the length of the flume to check the uniformity of mixture flow in the flume and to determine the flow depth depending on the thickness of the sediment bed at the bottom of the flume. In a measuring cross section perpendicular to the flow, distributions are measured of flow velocity and sediment concentration between the top of the sediment bed and the water surface.

In the experiments reported here, the distribution of sediment velocities above the bed was measured by the stereoscopic method using a set of two high-speed cameras (Figure 3). The method and the equipment are described in [26]. In some tests, the longitudinal velocity of sediment particles was also measured at various elevations in the flow cross section by the horizontal-oriented UVP sensor (the technique and method described in [24,33]). Longitudinal velocities of fluid were measured along the flow depth by a Prandtl tube [24]. The distribution of volumetric concentration of transported particles across the flow depth was measured by the laser stripe technique using a camera and a stripe-producing laser (Figure 3) [12,26,34].

2.2.3. Tested Solids

Flow of mixture of water and lightweight sediment above a plane bed composed of the same sediment was tested in the flume. The sediment was a fraction of plastic and virtually mono-size particles (fraction code SUN25) of an approximately prismatic shape (red particles in Figure 1). The experimentally determined properties of the tested solids are summarized in

Table 1. Properties of tested solids: Particle equivalent diameter d, specific gravity S, terminal settling velocity w, particle Reynolds number R.

	SUN25
d [mm]	2.8
S [-]	1.28
w [mm/s]	76.5
w [-]	0.987
R [-]	217

.

2.2.4. Experimental Flow Conditions

Tests were carried out for flow conditions which ensured a plane surface of the erodible bed, intense transport of sediment and combined load regime. Basically, the first two conditions were controlled by a value the bed Shields parameter, while the last condition required to stay within a certain range of values of the suspension ratio as suggested in [8]. These requirements, together with operability limits of the measuring set-up, constrained ranges of flow quantities attainable in the tests. Nevertheless, a sufficiently broad range of conditions was reached and included the conditions very near the limits of the combined-load regime. In the entire data set, the longitudinal slope varied between 0.25 deg and 2.04 deg and the mixture flow rate between 5.1 L/s and 12.0 L/s. The resulting flow depth varied from 42 mm to 83 mm and the bed Shields parameter from 0.45 to 2.0.

2.2.5. Experimental Data Set

Two experimental campaigns were carried out. The first campaign provided experimental data at different longitudinal slopes and flow depths to cover flow conditions in a sufficiently broad range of the bed Shields parameter to calibrate some of the model constants. The tests included measurements of velocity distributions using three independent techniques (camera, UVP, Prandtl tube) and measurements of concentration profiles using the laser-stripe method. Later, a validation test campaign was executed with the same sediment at slightly different flow conditions, primarily at bigger flow depths and mixture flow rates and milder longitudinal slopes than in the calibration campaign. The conditions captured as broad as possible range of values of bed Shields parameter including those near the limits of the upper plane bed regime. The validation tests did not include UVP measurements due to time constraints. However, the calibration tests confirmed a previously observed good agreement among results of all three methods in the particle rich parts of the flow (the combined-load layer and the lower part of the suspended-load layer) and thus the absence of UVP data in the validation data set did not pose a problem for model validation. Selected results of the validation tests are discussed below.

3. Results

3.1. Experimental Results

In the following, results of selected validation test runs are summarized and discussed. The selection covers a wide range of bed shear conditions from those at the maximum Shields parameter attainable in the experiments to those at so low value of the Shields parameter that the presence of suspended load in the transported sediment is negligible (the suspension ratio of the bed shear velocity $u_{*0}=\sqrt{{\mathsf{\theta }}_0\mathsf{\cdot }S}$ and the terminal settling velocity w is significantly lower than 1).

Table 2. Combined-load experiment in tilting flume: Experimental conditions for Tests 1 to 5 (longitudinal slope ω, flow depth H, flow Reynolds number Re, bed Shields parameter θ_0,flow, suspension ratio u_*0,flow/w; all quantities are dimensionless).

	Test 1	Test 2	Test 3	Test 4	Test 5
ω [-]	0.0257	0.0142	0.0154	0.0098	0.0070
H [-]	20.4	25.5	20.0	25.1	22.7
Re [-]	4.0 × 104	6.1 × 104	4.1 × 104	5.1 × 104	4.0 × 104
θ0,flow [-]	1.88	1.30	1.10	0.88	0.57
u*0,flow/w [-]	1.57	1.31	1.20	1.08	0.86

summarizes experimental conditions for five selected test runs from the validation data set. In Table 2, θ_0,flow is obtained from the measured H and ω using Equation (17) and u_*0,flow is based on θ_0,flow.

3.1.1. Concentration Profiles

The measured shapes of concentration profiles exhibit systematic local changes in trends which indicate an existence of the layered structure of the tested flow and suggest positions of the layer interfaces (the right-hand-side panels in Figure 4). In the lowest part of a profile, there is a characteristic sharp change in the concentration gradient which becomes considerably steeper at elevations where the local concentration c drops below say 0.45. While the resolution of the measurement is not sufficient to unambiguously detect the elevation and the associated value of the local concentration at which the change actually occurs, it can be hypothesized that the change in the concentration gradient is associated with the change between the predominantly permanent contact regime and predominantly collisional regime of particle interaction in the flow. In the model, this change is associated with the interface between the combined-load layer and the permanent-contact layer. The observed region of the approximately linear profile above the assumed top of the permanent-contact layer justifies the model assumption of a linear profile in the combined-load layer.

Another systematic and detectable feature of the measured concentration profiles is that the profiles can be no longer considered linear at elevations with c lower than say 0.25. This is different from concentration profiles in purely collisional flows, which can be approximated by a line from the bottom to the top of the transport layer [3,26]. Therefore, the non-linear shape of the profile can be attributed to the dominating turbulent suspension of particles and the lower boundary of the non-linear profile can be considered the interface between the combined-load layer and the suspended-load layer. The region occupied by the non-linear profile is the largest in the flow with the highest value of u_*0,flow/w (Test 1: u_*0,flow/w = 1.57, see Table 2) where it spans a distance between the assumed top of a thin combined-load layer and some position very near the water surface, i.e., the interval 5 < y < 20 approximately (Figure 4). If the u_*0,flow/w value decreases (Tests 2 to 5), then the region of full suspension narrows, its top departs from the water surface and leaves an increasingly thick particle-free zone below the water surface. At the lowest u_*0,flow/w (equal to 0.86 in Test 5), the full-suspension region is marginal and can be neglected. Note also that the permanent-contact layer tends to become thinner if the applied bed shear stress (the bed Shields parameter) decreases. In general, this layer is considerably thinner than the other layers at all bed shear conditions covered by Tests 1–5.

3.1.2. Velocity Profiles

The velocity profiles measured in Tests 1–5 all exhibit a region of an approximately linear profile which corresponds with the region delimited by the boundaries of the combined-load layer and which exhibits also a linear concentration profile (the left-hand-side panels of Figure 4). In this region, there is a very good agreement between local velocities measured by two different measuring techniques: The Prandtl tube (sensing local fluid velocity) and the stereoscopic technique (sensing local solids velocity). It indicates that a difference in local velocities of fluid and particles can be neglected. Additional measurements of the sediment velocity by UVP used in the calibration tests confirmed a good accuracy of the stereoscopic method and this was consistent with our previous observations for different lightweight sediment in the flume [26]. However, the agreement between results from the Prandtl tube and from the stereoscopic method is considerably weaker in the particle-lean region above the top of the linear profile (i.e., in the suspended-load layer) where a measured velocity profile is no longer linear. The local velocities by the stereoscopic method are systematically lower than those by the Prandtl tube and the difference seems to increase with the decreasing local concentration of sediment at positions closer to the water surface. It can be attributed to lower accuracy of the stereoscopic method in a flow domain with a low population of particles. The observed deviation is again consistent with our previous observations [26].

3.2. Model Predictions

Distributions of concentration and velocity are predicted by the model for the conditions of Tests 1–5 (see Table 2) following the procedure described in Section 2.1.5. The predictions include the elevations of interfaces in the sediment-water flow and they are presented in Figure 4 together with the experimentally determined distributions. Parts of the predicted profiles of velocity and concentrations are plotted in different colors to distinguish among predicted layers across the flow depth.

In the model calculations, the effective coefficient of restitution in dry conditions ε = 0.60 and the particulate friction coefficient β₀ = 0.65 are employed for the SUN25 sediment.

3.2.1. Concentration Profiles

Based on the model assumptions, the predicted distributions of solids are linear across the permanent-contact layer and the combined-load layer. The slopes of the profiles result from the predicted elevations of the layer interfaces (y_d, y_c) above the bed (y₀ = 0), from the chosen values of the local concentrations at the interfaces (c₀, c_d) and from the calculated c_c at y_c. The values of c₀ = c_cr (c_cr = 0.62) and c_d = 0.45 were chosen for all test runs to be consistent with the measured values of c at the positions of the abrupt change in a concentration gradient as discussed above. At the d-boundary, c_d was considered the freeze concentration (c_d = c_f) to satisfy the kinetic-theory based condition for dense limit (L_d = 1). The c_f is sensitive to the shape of particles and this justifies its chosen value which is slightly lower for our tested particles than for particles of simple shapes (e.g., 0.49 for cylindrical particles). In the suspended-load layer, a predicted shape of the concentration profile is sensitive to the particle diffusivity constant, C_η = 2 was used for all tested flow conditions.

3.2.2. Velocity Profiles

The model assumptions for velocity gradients define shapes of velocity profiles in the permanent-contact layer (the shear rate γ linearly distributed) and in the combined-load layer (the shear rate constant). The shear rate is assumed equal to zero at the bed (γ₀ = 0) and it is calculated by a constitutive relation at the d-boundary (Equation (7)). At the c-boundary, γ_c = γ_d. The predicted elevations y_d and y_c determine u_d and u_c. In the suspended-load layer, a predicted shape of the velocity profile is logarithmic and affected by the local mixing length (κ = 0.41) and by the initial position determined empirically using the constant C_ini, which is the same for all tested flow conditions (C_ini = 1.65).

4. Discussion

Comparisons of the predicted and measured distributions of velocity and concentration show very reasonable agreement and indicate that the model is capable of capturing the layered character of the flow with intense transport of combined load (Figure 4). The distributions inside the individual layers are captured well too, although more detailed modelling of processes at the local level may provide more precise information about the distributions, particularly in the suspended load layer. If the suspended-load layer is thick (Tests 1 to 3 in Figure 4), it tends to exhibit a complex shape of the concentration profile composed of a convex part and a concave part. This change in the profile shape in the central region of the layer cannot be reproduced by the current model.

The transport model has been calibrated and validated by results of our own experiments which produced sufficiently detailed information about the tested flows including information on the local level from the channel bed to the water surface. Those experimental results are quite unique and enabled to test model abilities in a broad range of conditions including those close to the limits of combined-load transport.

A prediction of distributions of velocity and concentration provides information that enables to determine flow rates of mixture and its phases for the modelled flow. It makes the model a useful tool from an engineering perspective as it is capable of simulating mutual relations among the longitudinal slope, the flow depth, the flow rate of sediment and the total flow rate of sediment-water mixture. In this sense, the model is a more physically sound alternative to currently widely used empirical transport models.

Note also that the proposed transport model does not require the bed friction law which the above-mentioned empirical transport models need. In such models, the friction formula for the top of the erodible bed interacting with particles transported above the bed is a significant source of uncertainty affecting a determination of the relation between the flow rate, the flow depth and the water surface elevation. The uncertainty is associated primarily with difficulties in the interpretation of the roughness of the erodible bed at the condition of intense transport of sediment above the bed. This problem is avoided in the proposed model and hence it solves the bed friction problem in a more elegant way.

However, future work is required to further refine the model. Modelling of the quite specific conditions at the bottom of the combined-load layer, which seem to differ from conditions at the bottom of a purely collisional layer and which are therefore more challenging to describe using the kinetic-theory based constitutive relations, must be addressed so that the reduction constant C_bd can be removed. Furthermore, modelling of interactions of suspended particles with turbulent flow in the sediment-rich region just above the top of the combined-load layer requires a refinement to reproduce more closely the characteristic shape of a concentration profile observed experimentally in this region.