## Abstract

We derive a two-layer depth-averaged model of sediment transport and morphological evolution for application to bedload-dominated problems. The near-bed transport region is represented by the lower (bedload) layer which has an arbitrarily constant, vanishing thickness (of approx. 10 times the sediment particle diameter), and whose average sediment concentration is free to vary. Sediment is allowed to enter the upper layer, and hence the total load may also be simulated, provided that concentrations of suspended sediment remain low. The model conforms with established theories of bedload, and is validated satisfactorily against empirical expressions for sediment transport rates and the morphodynamic experiment of a migrating mining pit by Lee *et al.* (1993 *J. Hydraul. Eng.* **119**, 64–80 (doi:10.1061/(ASCE)0733-9429(1993)119:1(64))). Investigation into the effect of a local bed gradient on bedload leads to derivation of an analytical, physically meaningful expression for morphological diffusion induced by a non-zero local bed slope. Incorporation of the proposed morphological diffusion into a conventional morphodynamic model (defined as a coupling between the shallow water equations, Exner equation and an empirical formula for bedload) improves model predictions when applied to the evolution of a mining pit, without the need either to resort to special numerical treatment of the equations or to use additional tuning parameters.

## 1. Introduction

The majority of coastal and river morphodynamic models employed in practice consist of coupled systems of depth-averaged flow mass and momentum equations, a bed-update equation, and a sediment transport formula. We refer to this type of model as a conventional morphodynamic (CM) model. A feature common to almost all CM models is the use of empirical formulae for estimation of sediment transport rates, necessary for closure of the morphodynamic model. However, the vast number of such formulae available in the literature makes selection of the most appropriate expression difficult, leading to considerable uncertainty. For example, sediment transport in open channel flows is typically divided into bedload and suspended load; but, although different mechanisms govern these two modes of transport, a reliable method to distinguish one from the other has yet to be provided [1]. Therefore, ambiguity in identification of the mode of transport present (i.e. bedload versus suspended load) is reflected in the selection of the closure formula, thus adding to uncertainty. Further uncertainty in CM models also arises from our lack of understanding of the fundamental mechanics behind sediment transport. In particular, the effect of bed slope on sediment transport should be included in reliable morphodynamic models of any kind [2,3]. However, the bed slope effect is often neglected, or accounted for through additional tuning parameters, adding empiricism to CM models and thus uncertainty when no data are available for calibration.

Alternative morphodynamic models include two-phase and two-layer models, which may be scientifically more insightful than CM models, but at the cost of increased mathematical complexity and computational demand. Thus, such alternatives tend to be more appealing to a scientific audience than to the practitioner community. Among two-layer models, it is worth noting the early work by Fraccarollo & Capart [4], who introduced the idea of an erosion rate estimated using simple concepts from open channel hydraulics and soil mechanics, replacing part of the empiricism inherent to sediment transport formulae with physical mechanisms that drive bed erosion. The model comprised clear water flowing on top of a constant-density sediment–water mixture, which in turn had the same average density as the non-mobile bed underneath; both fluid layers moved at the same speed. Later, Spinewine [5] extended the work by Fraccarollo & Capart [4] to account for the different velocities and concentrations in the two fluid layers, while assuming that the average density of the transport layer was constant. The latter restriction was removed by Li *et al.* [6], who considered a variable-density lower layer; the variability in density was in turn estimated via an empirical formula for sediment transport. All the aforementioned models simulate clear water over a transport layer, and track the evolution of a distinct physical interface dividing both layers. (For a comprehensive review of depth-integrated two-layer models, see [7].) Two-phase models are also strong candidates for simulating morphological evolution and sediment transport, yielding interesting insights [2,8,9]. However, we restrict our present study to depth-averaged two-layer approaches.

In an attempt to reconcile the scientific and practitioner communities interested in morphological modelling, we propose a simplified two-layer morphodynamic model, which may also be used to enhance CM approaches. The model idealizes shallow water–sediment flow as being divided into two layers with temporally and spatially varying densities in the plane parallel to the mean bed. To deal with the inherent ambiguity in the distinction between different modes of transport, the thickness of the lower layer is fixed at a small value, distinguishing the present model from previous two-layer approaches. Then, in the spirit of Fraccarollo & Capart [4], simple constitutive equations are used to estimate the driving mechanism of bed erosion and other closure terms, such that selection of a particular sediment transport formula is not required. The model is primarily designed for bedload-dominated sedimentary processes, although it is demonstrated that suspended load can also be simulated provided the suspended sediment concentration is low. The model is compared against an experimental study of bedload-driven migration of a mining pit, and then employed to derive an analytical, phenomenologically meaningful expression (free from tuning parameters) for the bed slope effect on bedload. This expression can be included in CM models, enhancing their accuracy without increased empiricism, as verified after comparison against the migrating pit experiment.

The body of the paper is organized in three parts. The first part (§2) deals with the description of the proposed model and its underlying assumptions. Analysis, as well as experimental and mathematical validation of the model are included in the second part (§3). The third part (§4) deals with the effect of bed slope on bedload, presenting the derivation and application of a slope-related term that can be employed within CM models; comparison against morphodynamic experimental data is also presented. Section 5 summarizes the key findings.

## 2. The model

### 2.1. Description

A Cartesian frame of reference (*x*,*z*) is adopted, where *x* and *z* are the stream-wise and vertical coordinates, respectively (the transverse coordinate is not considered in this paper); and time is denoted by *t*. The water–sediment mixture is divided into two layers: the lower concerned with bedload transport, and the upper representing sediment in suspension. Although the model is primarily designed to deal with bedload transport, relaxation of the assumption—common to most two-layer models—of an upper water layer allows suspended load to be considered at dilute concentrations. The water–sediment mixture is assumed to be an incompressible continuum, with each layer experiencing zero vertical acceleration so that the flow is hydrostatic and parallel to the mean bed. As with most models for hydrostatic flows, we assume a small bed slope and refer to the flow parallel to the bed as being ‘horizontal’—non-negligible bed slopes and their influence on morphological evolution are treated in detail in §4. A single sediment size and uniform bed porosity are considered. In each layer, the horizontal velocity, *u*(*x*,*t*), and sediment concentration, *c*(*x*,*t*), are assumed uniform with depth (see figure 1), but varying in the stream-wise direction and time. The bedload layer is assumed to have constant, arbitrary, vanishing thickness and variable density. This permits simulation of the (often ambiguously defined) bedload layer as a near-bed transport zone, whose sediment concentration can vary from zero for no sediment transport, to a maximum or saturation value. While typical two-layer models (e.g. [4,5,10]) aim to track the evolution of an interface dividing two layers with different but homogeneous densities (usually those of water and water–sediment mixture), in the present model such an interface is instead thought of as an imaginary line whose sole role is to delimit the near-bed region; in other words, it sets a tracking volume near the bed. The idea behind this assumption is to translate into the mathematical model the inherent ambiguity in the distinction between different modes of transport by fixing an arbitrary, near-bed layer concerned with bedload, and to compensate for this arguably restricting condition by relaxing the assumption of an upper layer composed of clear water (i.e. allow sediment to enter the upper layer). This key feature distinguishes the present model from other two-layer models, which is why we refer to it as a quasi-two-layer (Q2L) model.

The model may simulate three different modes of transport:

— Mode 0: no sediment transport. For flow conditions below the threshold of sediment motion, both layers consist of pure water (figure 1

*b*).— Mode 1: bedload only. Flow conditions are such that only the lower layer carries sediment and the upper layer consists of pure water (figure 1

*c*). The model is primarily concerned with this type of transport.— Mode 2: total load. We relax the assumption of an upper layer consisting of clear water, such that at higher flow conditions, the bedload layer has reached a saturation point and sediment entrains into the upper layer, where it is treated as suspended load (figure 1

*d*). Only low concentrations of suspended sediment are considered.

It is assumed that the fluid is always in motion, encompassing hydrodynamic conditions leading to the modes described above. Cases not modelled include: still fluid; sediment being present in the upper layer when the lower layer has not reached saturation; sudden entrainment caused by intense turbulence or a lateral source of sediment; and fine cohesive sediments (which may enter suspension mode without necessarily going through the bedload stage). Sheet-flow transport, where a distinct interface occurs between the lower transport layer and the upper pure water layer, may be simulated as Mode 1. Flows carrying highly concentrated suspended loads are outside the scope of the present model, noting that different approaches (e.g. [2,11,12]) may be required in order to simulate such flows more accurately. Bedload-dominated problems are the main objective of the present model. However, the option of modelling low-concentration suspended load has been incorporated as a ‘safety valve’ that allows the user to study a problem without *a priori* certainty that only bedload transport will take place. To model entrainment and deposition, the bed erosion rate is estimated from the conservation of horizontal momentum at the bed interface, following Fraccarollo & Capart [4].

### 2.2. Governing equations

#### 2.2.1. Derivation

The governing equations are derived by (i) applying conservation laws to overall water–sediment mass in the upper layer, sediment mass in the upper layer, sediment mass in the lower layer, overall water–sediment horizontal momentum in the upper layer, overall water–sediment horizontal momentum in the lower layer and sediment mass in the bed; (ii) considering mass and momentum exchanges between all three layers (bed included); and (iii) assuming thickness of the lower layer to be constant in space and time. Contribution of (a small) bed slope to momentum fluxes is considered. The following set of equations is obtained:
*a*
*b*
*c*
*d*
*e*
*f*where subscripts ‘0’ and ‘1’ refer to the lower and upper layers, *L*_{0} and *L*_{1}, respectively; *ρ*, *h* and *u* are the layer density, depth and horizontal velocity, respectively; *g* is gravitational acceleration; *z*_{b} is the bed level with respect to a fixed horizontal datum (subscript b refers to the bed layer, *L*_{b}, assumed to be static; i.e. *u*_{b}=0); *i* and *j* are net mass and momentum exchanges between layers (taken as positive in the upward direction) through interfaces denoted by superscripts (*i*) and (*b*) (figure 2); *e*^{(b)} denotes net water–sediment volumetric exchange between *L*_{b} and *L*_{0} (constant value of the bed bulk density, *ρ*_{b}, has been assumed). Note that the assumption of uniform velocity profiles implies the Boussinesq profile coefficient (commonly found in depth-averaged models) is equal to unity in the momentum balance equations, and so is not included here. A similar remark applies to the concentration profiles; i.e. a uniform, fully mixed profile, in conjunction with uniform velocity, yields a profile factor equal to unity in the mass conservation equations. The key assumption of ∂*h*_{0}/∂*t*=∂*h*_{0}/∂*x*=0 is implicit in (2.1). Observe that (2.1b) and (2.1c) imply that sediment particles are transported by the fluid at the same stream-wise speed as the whole water–sediment mixture. In other words, *u*_{sk}=*u*_{k} is assumed, where *u*_{sk} represents the stream-wise velocity of sediment particles in layer *L*_{k} (*k*=0,1). This is a sensible assumption for the upper layer, given that, for suspended load, sediment is expected to be transported at about the same speed of the flow [13]. As for bedload, in appendix A, a Lagrangian model for particle saltation is used to test the sensibility of this hypothesis, revealing that sediment transport rates predicted by the model remain virtually unaffected. Therefore, equivalence is pragmatically assumed between the sediment particle stream-wise velocity, *u*_{s0}, and the velocity of the corresponding water–sediment mixture, *u*_{0}. (Such equivalence is also backed by experimental evidence, see [14].) In this first stage of model development, we try to retain simplicity whenever possible, acknowledging that additional modifications may be needed eventually to render the model more applicable to real, complex scenarios (see §§3.6 and 5).

#### 2.2.2. Erosion rate and shear stresses

Following Fraccarollo & Capart [4], the bed erosion rate is estimated from conservation of longitudinal momentum at the bed interface. Across interface (*b*), momentum flux has to be single-valued, and so *j*^{(b)} ought to be computed from variables at either side of the interface, yielding *u*_{b}=0), where *e*^{(b)}=*i*^{(b)}/*ρ*_{b}, can then be estimated as

As stated previously, the present idealization assumes the fluid is always in motion; however, should *u*_{0}=0 occur at some point in the domain at a certain time, the condition *e*^{(b)}=0 for *u*_{0}=0 is imposed in order to avoid a mathematical error being introduced.

Herein, *u*_{0}, and the bedload layer average density, *ρ*_{0}; namely
*C*^{(b)} is a friction coefficient, one of the main calibration parameters within the present model. By making *ρ*_{0} (as opposed to the water density, *ρ*_{w}), the idea is to incorporate, even if crudely, the influence of the sediment transport contribution to the total bed shear stress.

Following Spinewine [5] and Zech *et al.* [10], the bed interface is treated as a failure plane, such that the shear stress *ρ*_{w} is the density of water, *τ*_{c} is the critical yield stress, obtained from the Shields curve, and *φ* is the soil friction angle, taken as equal to the angle of repose. Note that both *τ*_{c} and *φ* depend on sediment characteristics, such as the type of sediment, its density and its particle diameter. The term |*u*_{0}|/*u*_{0} is included to ensure that

Flux of horizontal momentum at (*i*) is also required to be single-valued regardless of whether variables from the top or bottom side of the interface are invoked. However, such a flux does not evolve freely, but instead depends on the bed erosion rate, *e*^{(b)}, thus requiring only one of the two shear stresses at interface (*i*) (i.e. *i*), *C*^{(i)} is a second calibration coefficient, the first being *C*^{(b)} in equation (2.3).

Closure relationships for shear stresses herein adopted prioritize simplicity and ought to be considered as exploratory measures, which should be revised in the future by comparison against high-quality experimental data.

### 2.3. Inter-layer mass and momentum fluxes

Starting from Mode 0, when bed material is initially eroded, *L*_{b} and *L*_{0} exchange *water–sediment* mass. Conservation of volume within *L*_{0} (i.e. *h*_{0}= constant) requires a compensating flux of *water* between *L*_{0} and *L*_{1}. This is true until *L*_{0} gets saturated with sediment (i.e. *c*_{0} has attained its maximum permitted value, *c*_{0 }), in which case the *water–sediment* mass mixture eroded/deposited from/onto *L*_{b} is compensated by an equal amount of *water–sediment* mass exchanged between *L*_{0} and *L*_{1}. Therefore, the net mass fluxes through interfaces (*b*) and (*i*) are expressed as

The above expressions underpin the logical requirement that no sediment may be deposited onto the bed when operating as Mode 0 (where *L*_{0}).

The corresponding *sediment* mass fluxes are the sediment components of the total water–sediment mass exchanges and are expressed as
*ρ*_{s} is the density of sediment.

Exchange of horizontal momentum between layers takes place when mass crosses the interface from one layer to an adjacent layer. This occurs at both interfaces (*i*) and (*b*) when *i*^{(b)}≠0 (and hence *i*^{(i)}≠0). As mentioned before, momentum fluxes at such interfaces are to be single-valued. However, at the bed interface, we compute the horizontal momentum flux from variables at the upper part of interface (*b*); namely

This is because (2.10) implies the physically meaningful condition that *L*_{0} faces solely a resistive bed friction proportional to the square of its velocity when the model is operating as Mode 0. Moreover, analysis of equation (2.10) enables correct interpretation of the momentum flux, *j*^{(b)}. This is important because, as remarked upon by Iverson & Ouyang [7], such a flux has been occasionally reported in the literature to yield the seemingly illogical conclusion that *L*_{0} gains momentum due to mass crossing the interface (*b*) from a state originally at rest (and hence with no initial momentum). Under no sediment transport conditions, the total resistance encountered by the flow is *L*_{0} is reduced by a factor of −|*i*^{(b)}*u*_{0}|, and thus any apparent gain of momentum of *L*_{0} is in actuality a reduction of the diffusive momentum or basal friction [7].

As with the bed, the net flux of momentum at (*i*) can be evaluated from variables at either side of the interface. However, unlike (*b*), for (*i*) there is no preferential candidate based on physical significance, and so to ensure consistency with the previous arbitrary parametrization of *L*_{0} are invoked, yielding

Equations (2.3)–(2.11) close the set of governing equations given by (2.1).

## 3. Model validation

We devote this section to analysis of the model and comparison against established theory and experiments on sediment transport rates and morphodynamics. Special emphasis is given to bedload-dominated problems because this is the main aim of the present model, as previously remarked upon. Further details of the mathematical treatments, proofs and derivations that follow can be found in [16].

### 3.1. Analytical solution for bedload

We consider steady uniform flow over an erodible bed with bedload transport exclusively present (*c*_{0}≤*c*_{0 mx} and *c*_{1}=0). This case permits derivation of an analytical solution to the Q2L model, which can then be used to compare the present model against bedload theory, including validation against empirical formulae. The volumetric bedload transport rate, *q*_{b}, is evaluated as *q*_{b}=*h*_{0}*c*_{0}*u*_{0}. Note that the *sediment* bedload rate should strictly be computed as *q*_{b}=*h*_{0}*c*_{0}*u*_{s0}. However, this would require an additional equation relating *u*_{s0} to the model output *u*_{0}. Appendix A proves that the assumption *u*_{0}≈*u*_{s0} appears sensible from quantitative and pragmatic perspectives and is thus adopted herein.

For steady uniform flow that is initially above the threshold of sediment motion, equilibrium conditions for sediment transport are expected to be reached eventually. Given that *h*_{0} is a constant within the present model, equilibrium-state values have to be found solely for *c*_{0} and *u*_{0} in order to compute *q*_{b}. Once bed erosion has initiated, steady sediment transport conditions can only occur once *e*^{(b)} decreases to zero. This happens when both *c*_{0} and *u*_{0} have reached certain values, *c*_{0 eq} and *u*_{0eq}, respectively, such that the force exerted by the water–sediment flow on the bed surface equals its resistance to erosion. In other words, *e*^{(b)}=0 to occur. Hence, *ρ*_{0 eq}, and thus *c*_{0 eq}, can be found. Then, an expression for *u*_{0 eq} simply follows from (2.3), allowing calculation of the bedload transport rate, equal to *h*_{0} *c*_{0 eq} *u*_{0 eq}, from

Note that bedload is independent of the second calibration parameter, *C*^{(i)}.

### 3.2. Mathematical agreement with bedload formulae

Inspection of (3.1) reveals that the bedload transport rate predicted by the present model follows the general form
*q*_{b}=*F*(*τ*−*τ*_{c})*τ*^{1/2}, where *F* is an expression often taken as a constant obtained from a best-fit curve to laboratory data, and *τ* (*q*_{b}=*h*_{0}*c*_{0}*u*_{s0}, and relating *u*_{s0} to *u*_{0} by means of a Lagrangian saltating particle model, *q*_{b} instead follows the form

Also note that for sheet-flow conditions (a special regime of bedload), where

### 3.3. Validation (empirical bedload formulae)

Figure 3 compares model predictions (equation (3.1)) against corresponding values from popular empirical bedload formulae in [22] (MP & M), [19] (Y), [21] (A & M), [23] (W), [18] (N) and [24] (FL & vB). Two particle diameters, *D*, are considered: 0.5 and 2.0 mm. These particle sizes are chosen because bedload is likely to be the main mode of transport for *D*>0.3 mm [25]. Parameter values are *s*≡*ρ*_{s}/*ρ*_{w}=2.65, *φ*=32.1^{°} and *h*_{0}=10*D* (more on *h*_{0} in §3.4). Three values of *C*^{(b)} are investigated; namely, 0.01, 0.03 and 0.06.

The present model predictions fall within the band of estimates delimited by the empirical formulae considered; this band is representative of the well-known uncertainty in the quantification of bedload. Here, the model predictions are truncated where *c*_{0}=*c*_{0 mx} is reached. Experiments by Spinewine [5] show that values of sediment concentration within the bedload layer tend to be confined to the range [0.21,0.25], and so a slightly larger value of *c*_{0 mx}=0.3 is selected in order to extend the model prediction curves, noting that such a selection purely acts as the Mode 1 limit; in other words, the larger the value of *c*_{0 mx}, the longer the model can operate as bedload only. Comparison between figure 3*a* and figure 3*b* demonstrates that the overall behaviour of the Q2L model in relation to established empirical formulae is independent of particle size, within the range of parameters considered. The recommended value of *C*^{(b)} depends on the reference formula, but the overall behaviour of the model matches empirical predictions, confirming the mathematical agreement discussed above.

### 3.4. Bedload layer thickness, *h*_{0}

Except for certain specific cases, such as sheet flow, a reliable, general method does not exist by which to determine the thickness of the bedload layer. Here, we prescribe an arbitrary, yet realistic, thickness of the bedload layer (i.e. *h*_{0}) that represents the *vicinity* of the bed where bedload occurs. Figure 4 investigates the sensitivity of the bedload predicted by the model to values of *h*_{0} in the range of [2*D*,20*D*]. This range is selected noting estimates of bedload layer thickness in [26,27]. The particle diameter is 1.0 mm. The influence of the arbitrary *h*_{0} on the predicted *q*_{b}, for a given *h*_{0} is the variation in the range of model validity when operating as Mode 1. The curves are plotted up to the point where *c*_{0}=*c*_{0 mx}=0.3, beyond which sediment is considered (within the framework of the present model) to be transported as suspended load (Mode 2). A larger value of *h*_{0} allows the model to operate as Mode 1 over a wider range of *c*_{0 mx}, such as

From figure 4, it can be observed that *h*_{0}≈10*D* yields a relatively wide range of validity for Mode 1, up to *et al.* [2] arrive at the conclusion that the thickness of the bedload layer is ∼10*D*. Therefore, the value *h*_{0}=10*D* is selected as the default unless otherwise stated.

### 3.5. No transport and total load (Modes 0 and 2)

Although the Q2L model is primarily designed for bedload-dominated problems, analysis of the other two possible modes of transport (i.e. no transport and total load) can lead to further useful insights. For example, the simple case of steady uniform flow below the threshold of motion can be invoked to prove that, in general, *C*^{(b)} is not equal to (in fact, is larger than) typical bed friction coefficients derived from Chézy or Manning formulae, which are employed in standard one-layer hydrodynamic models. Moreover, under these conditions, values for the ratio *C*^{(b)}/*C*^{(i)} can be proposed using either an assumed or a measured vertical velocity profile. For instance, if the flow velocity profile follows a power law as suggested by Soulsby [13], an expression for *C*^{(b)}/*C*^{(i)} as a function of *h*_{1} and *h*_{0} can be obtained, which may be useful (see below).

Of more interest is the analysis of Mode 2 (total load), which arises from relaxing the assumption of an upper layer consisting always of pure water, as a form of compensation for the arguably limiting constraint of a fixed bedload layer thickness. The user does not need to be concerned about violating the bedload-only condition during the simulation, but sediment may enter the upper layer at times and locations where flow is sufficiently fast, so long as near-bed transport dominates. In this way, Mode 2 operates as a ‘safety valve’ that adds flexibility to the model. To demonstrate the potential of the model to deal with total load, an analytical solution for Mode 2 has been derived (details not presented here, for brevity), and its results compared against empirical expressions for total transport proposed by Engelund & Hansen [28] and van Rijn [29]. Let *h*_{T}≡*h*_{0}+*h*_{1}=10 m, *D*=0.2 mm (uniform sediment), and *h*_{0}=10*D*. Unlike the bedload-only case, the analytical solution for total load requires both tuning parameters *C*^{(b)} and *C*^{(i)}. The additional degree of freedom for calibration (with respect to Mode 1) can be removed by proposing a value for the ratio *C*^{(b)}/*C*^{(i)} based on hydrodynamic considerations, as mentioned at the beginning of this section. For values of *h*_{1} and *h*_{0} considered, *C*^{(b)}/*C*^{(i)}≈5 can be prescribed. Figure 5 compares predictions by the selected empirical formulae and the present model for *C*^{(b)}= 0.05, 0.056, 0.06 and constant *C*^{(b)}/*C*^{(i)}=5. The Q2L model fits the formula by Engelund & Hansen [28] better than that of van Rijn [29]. For *C*^{(b)}=0.056, the agreement achieved between the model prediction and that by Engelund & Hansen [28] is outstanding over the range of parameters studied. The model predictions shown in figure 5 are for a bedload layer saturation value of *c*_{0 }=0.25, following Spinewine [5]. Further studies (not included here) demonstrate that the model shows little sensitivity to the selection of *c*_{0 }, over the realistic range *c*_{1} for curves shown in figure 5 are ∼0.001, hence *c*_{0}≫*c*_{1}, thus verifying the condition that bedload is the predominant mode of transport.

It is not intended that figure 5 be used to promote use of the model for problems dominated by suspended transport. Instead, figure 5 merely indicates the potential of the model to cope with scenarios where complex transport patterns occur, where the distinction between bedload and suspended load is unclear. Such potential is worth further exploration, but it is outside the scope of the present paper, which is on bedload-governed cases.

### 3.6. Validation (morphodynamic experiment)

To test the model, we select for comparison an experiment carried out by Lee *et al.* [30]. This experiment studies the migration of a mining pit due to bedload driven by a steady current. For reasons that become evident in §4, we focus on the second stage of the experiment, referred to by Lee *et al.* [30] as the ‘diffusion period’. Figure 6 illustrates this comparison. For reference, predictions by Chen *et al.* [31] are also included. In [31], a CM model is employed (i.e. hydrodynamic and morphological models coupled via a sediment transport empirical formula) with a bed-update equation modified by inclusion of an adaptation length (the distance it takes bedload to adjust from a non-equilibrium state to an equilibrium one). Figure 6 includes predictions in [31] for two values of this adaptation length; namely, 1 cm (labelled ‘Chen *et al.* 1’) and 2 cm (‘Chen *et al.* 2’—their reported best fit). Two predictions by the present model are also shown. In the first case (‘Q2L model 1’), values of the calibration parameters *C*^{(b)}=6.55×10^{−3} and *C*^{(i)}=4.50×10^{−2} are used, which yield the correct migration speed of the pit. However, the value reported by Lee *et al.* [30] of the upstream transport stage (ratio of bed shear stress to critical shear stress) is not replicated. In fact, when we aim to reproduce such a value (≈1.77), the predicted migration of the pit is significantly faster than the one reported. This has motivated the introduction of a further calibration parameter, similar to that of [31], in the bed-evolution equation (2.1f), leading to: ∂*z*_{b}/∂*t*=−*η*_{e}*e*^{(b)}. (See discussion in §4.3 regarding the abnormally low value of *F*_{*} required to reproduce the reported experimental settings in the context of a CM model.) Here, *η*_{e} deals with the irregularity in shape and size of the bed material present in the experiment, not accounted for in the model derivation (that assumes perfectly uniform sediment), which impact the bed’s packing fraction and thus the vertical distribution of its bulk density, *ρ*_{b}. (See appendix B for further details.) No detailed information on the bed composition is given in [31] that could aid the development of a sophisticated representation of *η*_{e}, and hence a constant value is assumed as a first step. A value of *η*_{e}=0.15 leads to the curve labelled ‘Q2L model 2’ in figure 6, where not only the correct pit migration speed is achieved, but also the reported value of the transport stage is replicated.

The present model without any modification (Q2L model 1) is able to predict the correct migration speed and final depth of the pit. However, from a qualitative perspective, significant discrepancies occur, with the overall bed level higher than observed, precisely because the low transport stage upstream ensures that no erosion takes place downstream from the pit (flow is from left to right in the plot). The second prediction by the model (Q2L model 2) gives better results, especially from a qualitative viewpoint. Note that this curve properly predicts the observed inflection point in the bed profile upstream of the pit, and the steepness of the upstream face of the pit. The latter condition is not predicted in [31]. Instead, the ‘Chen *et al.* 1’ results develop an unrealistic peak upstream of the pit, whereas the ‘Chen *et al.* 2’ results introduce diffusion that avoids development of the unrealistic peak, but at the cost of over-diffusing the pit profile.

The present experimental comparison indicates the potential of the model to deal with real-world morphodynamic problems, achieving predictions of quality at least comparable to previous studies. It, however, also highlights the need for future investigation into the idealized assumptions underpinning the derivation of the model. In the next section, the model is used to derive an expression for morphological diffusivity that can be incorporated into CM models for application to cases similar to the one studied here.

## 4. Bed slope influence on bedload: morphological diffusion

Under scenarios involving steep local bed slopes, such as mountain streams and certain beaches, the gradient of the bed elevation may play an important role in sediment transport processes and morphological evolution [1,32] as confirmed through laboratory experiments (e.g. [33–35]). Additionally, inclusion of the effect of bed slope on bedload may prevent the evolution of unrealistic morphological oscillations without the need to resort to specialized numerical techniques, while enhancing the physical significance of the model—bed gradients help diffuse out spurious bed features [36], thus the term ‘morphological diffusivity’. Methods used to account for the influence of bed slope include: addition of a slope-related diffusivity term to a sediment transport formula, which typically translates into an additional calibration parameter [36]; semi-empirical models based on Bagnold’s ideas (e.g. [17,37,38]); formulae explicitly derived for sloping beds, which often imply a significant degree of empiricism or complexity (e.g. [33,39,40]); and modification of the threshold of motion for sloping beds by inclusion of the weight of the particle at rest. In this section, we propose an analytically derived expression for morphological diffusivity that can be incorporated into a CM model for use in bedload-dominated problems.

### 4.1. Quantification of the bed slope influence on bedload

Here, we assess the bed slope influence on bedload using the ratio of bedload transport on a sloping bed to the bedload that would occur in a horizontal channel for the same bed shear stress. This is expressed as
*a*_{0} is a given value of *τ*, *β* is the bed slope angle, *τ*_{cβ} represents the threshold of particle motion for a bed of arbitrary slope and *τ*_{ch} is the threshold for a horizontal bed. Both quantities are related by incorporating the effect of gravity on particles at rest on a given slope, such that

Hence, using the analytical solution derived for bedload under steady uniform flow, (3.1), the bed slope influence, *Π*_{β}, can be rewritten (replacing *τ*) as
*Π*_{2}≡(*τ*−*τ*_{cβ})/(*τ*−*τ*_{ch}) and *Π*_{1}≡*ρ*_{0h}/*ρ*_{0β}. (Subscripts *β* and h denote sloping and horizontal beds, respectively.) Invoking the principle of equilibrium sediment transport conditions, such that *τ*=*τ*^{(b)}_{b}⇒*e*^{(b)}=0, as described in §3, expressions for *ρ*_{0 h} and *ρ*_{0β} can be obtained, yielding

The above equation depends on other variables besides the bed slope angle, *β*; namely: *τ*_{ch}, *τ*, *φ* and *h*_{0}. Hence, a sensitivity analysis is undertaken using the following data (typical of a bedload-dominated scenario): *D*=0.5 and 2.0 mm; *s*=2.63; *h*_{0}=5*D* and 15*D*; *τ*/*τ*_{ch} (not to be confused with *τ*/*τ*_{cβ})=2, 5 and 20; and *φ*=31 and 37^{°}. The results are plotted in figure 7, where it can be seen that the bed slope influence seems to be governed by the bed shear stress and the angle of repose, but insensitive to the selection of *D* and *h*_{0}/*D* within the ranges of values considered. The bed slope influence also vanishes (*et al.* [34]. This is confirmed mathematically by taking the limit of equation (4.4) as

Further study on the sensitivity of *Π*_{β} to *D* and *h*_{0} (not included here for brevity, but see [16]) demonstrates that, for all combinations of values considered,

The above equation is solely dependent on *β*, *φ* (through equation (4.2)) and the bed shear stress, *τ*. It is also worth remarking that both the exact and approximate expressions for quantifying the bed slope influence on bedload are independent of the calibration parameter *C*^{(b)}.

### 4.2. The morphological diffusion

A common way of accounting for the influence of bed slope is to modify the bedload formulae originally derived for horizontal channels by adding a term that promotes(inhibits) sediment transport in down(up)-sloping beds [36,37,41]. Such a term is proportional to the bed slope, and can be added as follows:
*ε*_{β} is a proportionality parameter related to morphological diffusion; *ε*_{β} is often taken as an additional tuning parameter in morphodynamic models [36,41]. Assuming, for convenience, positive unidirectional flow (i.e.

The ratio *q*_{bβ}/*q*_{bh} (consistent with our definition of bed slope influence, *Π*_{β}; see equation (4.1)) varies linearly with *S*_{b}; the parameter *ε*_{β} is the slope of the line (with *y*-intercept equal to 1). From figure 7, it can be observed that the nonlinear expression derived for *Π*_{β}≡*q*_{bβ}/*q*_{bh} exhibits quasi-linear behaviour consistent with (4.7) over a relatively wide range of *β*. Thus, an expression for *ε*_{β} can be proposed based on the bed slope influence predicted by the present model.

We extract analytically a value of *ε*_{β} by obtaining the slope of the line tangent to the curve *Π*_{β} versus *S*_{b} (as given by the approximation (4.5)) at the origin

By rewriting (4.5) as
*τ*_{cβ}/*τ*_{ch} depends on

Invoking (4.8), the morphological diffusivity is thus given by

The proposed diffusivity parameter depends on sediment characteristics (through *τ*_{ch} and *φ*) and the bed shear stress, vanishing for large values of the latter (i.e. *ε*_{β}→0 for *Π*_{β} (equation (4.4)) against linear fits to *Π*_{β} with line slopes given by (4.9). Agreement of the linear fit with the original expression for *Π*_{β} is very good for all negative bed slopes, mild adverse slopes and large bed shear stress. Discrepancies increase for low *τ* and steep adverse slopes. The bed slope influence is well described by the linear approximation for *Π*_{β} given by (4.4); however, additional tests carried out by the authors indicated that the expression yielded is significantly more complicated and depends on more variables than (4.9), with negligible quantitative improvement.

### 4.3. Incorporation of morphological diffusion into a conventional morphodynamic model

For this test, we define a CM model as the coupling between the one-dimensional nonlinear shallow water equations and the Exner equation via an empirical formula for bedload transport. We invoke a bedload formula of the form of Meyer-Peter & Müller’s [22]; namely: *Φ*=*F*_{*}(*θ*−*θ*_{c})^{3/2}; where *Φ* is the non-dimensional bedload transport; *F*_{*} is a non-dimensional constant obtained from model calibration; and *θ* and *θ*_{c} are non-dimensional forms of bed shear stress and critical shear stress, respectively. Morphological diffusivity is incorporated into the CM model by use of (4.6). For comparison, the model by Bailard & Inman [37] is also investigated. This model is of interest because Bailard & Inman [37] proposed a formula for bedload on a sloping beach of the form of (4.6), from which an analytical expression for *ε*_{β} may be deduced, provided that

We have derived the above expression by inserting the definition of bed slope influence, (4.1), in the formula by Bailard & Inman [37], and so *ε*_{B&I} is directly comparable to *ε*_{β}. Note that, although they have different roots, both *ε*_{B&I} and *ε*_{β} predict an inverse proportionality with *ε*_{β}, *ε*_{B&I} depends solely on the angle of repose and hence is independent of bed shear stress (i.e. it does not vanish at large *τ*). A similar expression to *ε*_{B&I} has been used by other authors in order to include the effect of the bed slope (e.g. [42]).

We return to the mining pit investigated previously. The objectives are to test the effect that inclusion of morphological diffusion has on a CM model and the influence of bed shear stress on morphological diffusion (i.e. compare *ε*_{β} versus *ε*_{B&I}). Figure 9 illustrates this comparison. It is worth mentioning that in order to replicate faithfully the experimental set-up reported in [30], coefficient *F*_{*} in the bedload equation takes a value of 2.3. This value is much lower than commonly used values proposed in [22] (i.e. *F*_{*}=8) and [24] (i.e. *F*_{*}=5.7). This implies some uncertainty in the reported set-up and appears to justify the use of a modified bed-update equation in the Q2L model (i.e. addition of *η*_{e}), as discussed in §3.6.

The CM model without morphological diffusion predicts the evolution of a very steep front and an unrealistic peak similar to one of the predictions by Chen *et al.* [31]. Inclusion of morphological diffusion prevents development of the aforementioned peak, and both diffusivity parameters *ε*_{β} and *ε*_{B&I} yield similar results, with the main difference between them being the slope of the predicted propagating front—which is steeper in the latter case. Without more detailed experimental data it is difficult to reach a definite verdict regarding the proposed *ε*_{β}, but two points can be made with certainty: (i) use of *ε*_{β} yields realistic results whose agreement with experiments is at least of comparable quality to previous studies and (ii) *ε*_{β} and *ε*_{B&I} lead to similar results, although the former has the advantage of translating correctly hydrodynamic effects into the influence of bed slope on bedload (i.e. *ε*_{β}→0 for

The same numerical scheme and set-up have been used for all predictions presented here. The CM model is discretized using second-order central differences in space, with fourth-order Runge–Kutta integration in time. All computations were undertaken on a uniform grid with Δ*x*=0.01 m and a time step Δ*t*=0.01 *s*. Discharge is imposed at both upstream and downstream boundaries; water depth is prescribed solely at the upstream boundary, and a transmissive condition is invoked at the downstream end; for the bed level, transmissive conditions are used in both boundaries. A ramp function is employed to bring the model from its initial condition of rest to the desired discharge, and when hydrodynamic steady state is reached, the sediment transport module is activated allowing the bed to evolve. However, it is important to highlight that even though numerical techniques have been widely used (e.g. [36]) to prevent evolution of spurious high-frequency oscillations such as the peak predicted (in certain cases) upstream of the pit, the morphological diffusion discussed here is underpinned by phenomenological observations (e.g. [3]). In other words, morphological diffusion is not merely a remedy for numerical instability but is also necessary for a realistic prediction of the morphodynamic model. Consider, for example, the case of the migrating hump studied in [43] within the framework of the present CM model. Under certain assumptions, the problem has an analytical solution by means of the method of characteristics, which, despite avoiding numerical oscillations, predicts a vertical wall (

## 5. Conclusion

A simplified two-layer model has been introduced for the prediction of sediment transport rates and morphological evolution for bedload-dominated scenarios. The model differs from previous two-layer models primarily through its treatment of the lower layer, which is modelled here as a thin layer of arbitrary but realistically small thickness. Thus, the inherent ambiguity in the definition of a near-bed transport layer is incorporated within the mathematical framework. Model results are found to be weakly dependent on the arbitrary selection of the lower layer thickness, within the range of values studied; a value of *h*_{0}≈10*D*, in agreement with previous works on bedload, is recommended. Although the model is devised specifically for bedload-dominated problems (such that *c*_{1}→0), relaxation of the assumption of an upper layer consisting of pure water allows the total load to be modelled to satisfactory accuracy (see figure 5) provided the concentration of suspended sediment is low. The Q2L model, successfully validated against empirical expressions for bedload, has then been compared against the experiment of a migrating mining pit in [30], with satisfactory agreement that rivals predictions in [31]. An analytical expression has been derived for morphological diffusion (equation (4.9)) which permits easy modification of bedload empirical formulae, originally derived for nearly horizontal flumes, in order to render them applicable to steep stream-wise slopes, provided the bed slope angle, *C*^{(b)}, *C*^{(i)}, *τ*^{(b)}_{b}, *τ*_{cβ}, *u*_{s0}=*u*_{0} or *ρ*_{b}=const.

## Data accessibility

This article does not include new experimental data.

## Authors' contributions

S.M. and A.G.L.B. planned the research. S.M. developed and programmed the model, carried out the simulations and interpreted the results. Both the authors wrote the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

S.M. was supported by the Mexican National Council for Science and Technology (Conacyt) through scholarship no. 310043 and the Mexican Ministry of Education (SEP).

## Acknowledgements

We thank the anonymous referees, whose insightful comments have led to a much improved paper. S.M. thanks the University of Edinburgh and Stanford University, where he was previously based and where much of this work was developed, and Joanna A. Zielińska for her help in reviewing some of the derivations.

## Appendix A

To study the relationship between *u*_{s0} and *u*_{0}, an expression of the form *D*_{*}≡*D*[(*s*−1)*g*/*ν*^{2}]^{1/3} is the non-dimensional particle diameter, with *s*≡*ρ*_{s}/*ρ*_{w} and *ν* being the kinematic viscosity of water; *u*_{*} is the bed friction velocity and *c*_{0}, and using the substitution *q*_{b}=*h*_{0}*c*_{0}*u*_{s0} as

Note that the above equation follows the form of bedload expressions proposed in [20,21]. It should be noted that, unlike (3.1), the above equation does not depend on *C*^{(b)} nor *h*_{0}; instead, calibration values of *C*^{(b)}≈0.01, both expressions yield very similar results. However, this is only the case for the values of calibration coefficients *u*_{0}=*u*_{s0}, appear valid from a practical viewpoint.

## Appendix B

Consider an ‘active layer’ of erodible bed material, of thickness, *z*′_{b}, and depth-variable density, *ρ*′_{b}(*z*); the layer is located such that its upper face corresponds with the bed interface at *z*_{b} (i.e. the layer is an upper sublayer of *L*_{b}). The layer’s average density, *z*_{κ}, between *z*=*z*_{b}−*z*′_{b} and *z*=*z*_{b}; i.e. *z*_{(κ)}≡*z*_{b}−*κz*′_{b}. Defining the mass exchange through the bed surface interface as *z*_{b}/∂*t*=∂*z*′_{b}/∂*t*, it follows that
*η*_{e} relates to the vertical profile of the bed average density, *z*′_{b} of an active erodible layer, which is expected to be different from zero in real experiments, such as the one considered here.

### Remark

If *ρ*′_{b}(*z*) varies linearly over *z*′_{b}, then

## Appendix C. Notation list

c_{k}(x,t) | sediment concentration of layer L_{k} (k=0,1) |

c_{b} | bed sediment concentration |

c_{0 mx} | maximum allowed value for c_{0} (saturation concentration for L_{0}) |

(C^{(b)},C^{(i)}) | Q2L model calibration parameters |

D | sediment particle diameter |

e^{(b)} | net water–sediment volumetric exchange between L_{b} and L_{0} (the erosion rate) |

g | gravitational acceleration |

h_{k} | thickness (depth) of layer L_{k} (k=0,1) |

i^{n} | net water–sediment mass exchange between layers through interface n (n=(b),(i)) |

i^{n}_{s} | net sediment mass exchange between layers through interface n (n=(b),(i)) |

j^{n} | net horizontal momentum exchange between layers through interface n (n=(b),(i)) |

L_{k} (k=b,0,1) | bed, lower and upper layer, respectively |

q_{b} | volumetric bedload transport rate per unit width |

q_{T}=total (bedload+suspended) | sediment transport rate per unit width |

whole-depth-averaged stream-wise velocity | |

u_{k}(x,t) | parallel-to-bed stream-wise velocity of layer L_{k} (k=b,0,1) |

u_{sk} | stream-wise sediment particle velocity within layer L_{k} (k=0,1) |

bed slope | |

t | time |

(x,z) | Cartesian frame of reference with stream-wise and vertical coordinates, respectively |

z_{b} | bed level with respect to a datum |

β | bed slope angle |

ε_{β} | morphological diffusivity parameter |

ε_{B&I} | morphological diffusivity derived from Bailard & Inman [37] |

η_{e} | additional Q2L model tuning parameter related to non-uniform bed material and packing fraction |

Π_{β} | bed slope influence as defined in (4.1) |

ρ_{k} | m density of layer L_{k} (k=b,0,1) |

ρ_{s} | density of sediment |

ρ_{w} | density of water |

τ_{c} | critical bed shear stress |

τ | bed shear stress |

shear stress exerted by the fluid on the bed surface | |

τ^{(b)}_{b} | bed resistance to erosion |

shear stress at the bottom and top of interface (i), respectively | |

φ | angle of repose |

subscripts β and h | denote sloping and horizontal bed, respectively. |

- Received November 27, 2017.
- Accepted January 19, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.