Morphodynamic model of the lower Yellow River : flux or entrainment form for sediment mass conservation ?

Sediment mass conservation is a key factor that constrains river morphodynamic processes. In most models of river morphodynamics, sediment mass conservation is described by the Exner equation, which may take various forms depending on the problem in question. One of the most widely used forms of the Exner equation is the flux-based formulation, in which the conservation of bed material is related to the stream-wise gradient of the sediment transport rate. An alternative form of the Exner equation, however, is the entrainmentbased formulation, in which the conservation of bed material is related to the difference between the entrainment rate of bed sediment into suspension and the deposition rate of suspended sediment onto the bed. Here we represent the flux form in terms of the local capacity sediment transport rate and the entrainment form in terms of the local capacity entrainment rate. In the flux form, sediment transport is a function of local hydraulic conditions. However, the entrainment form does not require this constraint: only the rate of entrainment into suspension is in local equilibrium with hydraulic conditions, and the sediment transport rate itself may lag in space and time behind the changing flow conditions. In modeling the fine-grained lower Yellow River, it is usual to treat sediment conservation in terms of an entrainment (nonequilibrium) form rather than a flux (equilibrium) form, in consideration of the condition that fine-grained sediment may be entrained at one place but deposited only at some distant location downstream. However, the differences in prediction between the two formulations have not been comprehensively studied to date. Here we study this problem by comparing the results predicted by both the flux form and the entrainment form of the Exner equation under conditions simplified from the lower Yellow River (i.e., a significant reduction of sediment supply after the closure of the Xiaolangdi Dam). We use a one-dimensional morphodynamic model and sediment transport equations specifically adapted for the lower Yellow River. We find that in a treatment of a 200 km reach using a single characteristic bed sediment size, there is little difference between the two forms since the corresponding adaptation length is relatively small. However, a consideration of sediment mixtures shows that the two forms give very different patterns of grain sorting: clear kinematic waves occur in the flux form but are diffused out in the entrainment form. Both numerical simulation and mathematical analysis show that the morphodynamic processes predicted by the entrainment form are sensitive to sediment fall velocity. We suggest that the entrainment form of the Exner equation might be required when the sorting process of fine-grained sediment is studied, especially when considering relatively short timescales. Published by Copernicus Publications on behalf of the European Geosciences Union. 990 C. An et al.: Morphodynamic model of the lower Yellow River


Introduction
Models of river morphodynamics often consist of three elements: (1) a treatment of flow hydraulics; (2) a formulation relating sediment transport to flow hydraulics; and (3) a description of sediment conservation.In the case of unidirectional river flow, the Exner equation of sediment conservation has usually been described in terms of a flux-based form in which temporal bed elevation change is related to the stream-wise gradient of the sediment transport rate.That is, bed elevation change is related to ∂q s /∂x, where q s is the total volumetric sediment transport rate per unit width and x is the stream-wise coordinate (Exner, 1920;Parker et al., 2004).This formulation is also referred to as the equilibrium formulation, since it considers sediment transport to be at local equilibrium; i.e., q s equals its sediment transport capacity q se , as defined by the sediment transport rate associated with local hydraulic conditions (e.g., bed shear stress, flow velocity, stream power, etc.) regardless of the variation in flow conditions.Under this assumption, sediment transport relations developed under equilibrium flow conditions (e.g., Meyer-Peter and Müller, 1948;Engelund and Hansen, 1967;Brownlie, 1981) can be incorporated directly in such a formulation to calculate q s , which is related to one or more flow parameters such as bed shear stress.
An alternative formulation, however, is available in terms of an entrainment-based form of the Exner equation, in which bed elevation variation is related to the difference between the entrainment rate of bed sediment into the flow and the deposition rate of sediment on the bed (Parker, 2004).The basic idea of the entrainment formulation can be traced back to Einstein's (1937) pioneering work on bed load transport and has been developed since then by numerous researchers so as to treat either bed load or suspended load (Tsujimoto, 1978;Armanini and Di Silvio, 1988;Parker et al., 2000;Wu and Wang, 2008;Guan et al., 2015).Such a formulation differs from the flux formulation in that the flux formulation is based on the local capacity sediment transport rate, whereas the entrainment formulation is based on the local capacity entrainment rate into suspension.In the entrainment form, the difference between the local entrainment rate from the bed and the local deposition rate onto the bed determines the rate of bed aggradation-degradation and concomitantly the rate of loss-gain of sediment in motion in the water column.Therefore, the sediment transport rate is no longer assumed to be in an equilibrium transport state, but may exhibit lags in space and time after changing flow conditions.The entrainment formulation is also referred to as the nonequilibrium formulation (Armanini and Di Silvio, 1988;Wu and Wang, 2008;Zhang et al., 2013).
To describe the lag effects between sediment transport and flow conditions, the concept of an adaptation length-time is widely applied.This length-time characterizes the distance-time for sediment transport to reach its equilibrium state (i.e., transport capacity).Using the concept of the adaptation length, the entrainment form of the Exner equation can be recast into a first-order "reaction" equation, in which the deformation term is related to the difference between the actual and equilibrium sediment transport rates, as mediated by an adaptation length (which can also be recast as an adaptation time) (Bell and Sutherland, 1983;Armanini and Di Silvio, 1988;Wu and Wang, 2008;Minh Duc and Rodi, 2008;El kadi Abderrezzak and Paquier, 2009).The adaptation length is thus an important parameter for bed evolution under nonequilibrium sediment transport conditions, and various estimates have been proposed.For suspended load, the adaptation length is typically calculated as a function of flow depth, flow velocity, and sediment fall velocity (Armanini and Di Silvio, 1988;Wu et al., 2004;Wu and Wang, 2008;Dorrell and Hogg, 2012;Zhang et al., 2013).The adaptation length of bed load, on the other hand, has been related to a wide range of parameters, including the sediment grain size (Armanini and Di Silvio, 1988), the saltation step length (Phillips and Sutherland, 1989), the dimensions of particle diffusivity (Bohorquez and Ancey, 2016), the length of dunes (Wu et al., 2004), and the magnitude of a scour hole formed downstream of an inerodible reach (Bell and Sutherland, 1983).For simplicity, the adaptation length can also be specified as a calibration parameter in river morphodynamic models (El kadi Abderrezzak and Paquier, 2009;Zhang and Duan, 2011).Nonetheless, no comprehensive definition of adaptation length exists.
In this paper we apply the two forms of the Exner equation mentioned above to the lower Yellow River (LYR) in China.The LYR describes the river section between Tiexie and the river mouth and has a total length of about 800 km. Figure 1a shows a sketch of the LYR along with six major gauging stations and the Xiaolangdi Dam, which is 26 km upstream of Tiexie.The LYR has an exceptionally high sediment concentration (Ma et al., 2017), historically exporting more than 1 Gt of sediment per year with only 49 billion tons of water, leading to a sediment concentration an order of magnitude higher than most other large lowland rivers worldwide (Milliman and Meade, 1983;Ma et al., 2017;Naito et al., 2018).However, the LYR has seen a substantial reduction in its sediment load in recent decades, especially since the operation of the Xiaolangdi Dam beginning in 1999 (Fig. 1b), because most of its sediment load is derived from the Loess Plateau, which is upstream of the reservoir (Wang et al., 2016;Naito et al., 2018).Finally, the bed surface material of the LYR is very fine, as low as 15 µm.This is much finer than the conventional cutoff of wash load (62.5 µm) employed for sediment transport in most sand-bed rivers (National Research Council, 2007;Ma et al., 2017).
When modeling the high-concentration and fine-grained LYR, it is common to treat sediment conservation in terms of an entrainment-based rather than a flux-based formulation.This is because many Chinese researchers view the entrainment formulation as more physically based, as it is capable of describing the behavior of fine-grained sediment, which when entrained at one place may be deposited at some distant location downstream (Zhang et al., 2001;Ni et al., 2004;Cao et al., 2006;He et al., 2012;Guo et al., 2008).However, the entrainment formulation is more computationally expensive and more complex to implement.Because the differences in prediction between the two formulations do not appear to have been studied in a systematic way, here we pose our central questions.Under what conditions is it valid to use the entrainment form of the Exner equation, and under what conditions may the flux form be used?Or more specifically, which form of the Exner equation is most suitable for the LYR?
Here we study this problem by comparing the results of flux-based and entrainment-based morphodynamics under conditions typical of the LYR.The organization of this paper is as follows.The numerical model is described in Sect. 2. In Sect.3, the model is implemented to predict the morphodynamics of the LYR with a sudden reduction of sediment supply, which serves to mimic the effect of the Xiaolangdi Dam.We find that the two forms of the Exner equation give similar predictions in the case of uniform sediment, but show different sorting patterns in the case of sediment mixtures.In Sect.4, we conduct a mathematical analysis to explain the results in Sect.3; more specifically, we quantify the effects of varied sediment fall velocity in the simulations.Finally, we summarize our conclusions in Sect. 5.

Model formulation
In this paper, we present a one-dimensional morphodynamic model for the lower Yellow River.The fully unsteady Saint-Venant equations are implemented for the hydraulic calculation.Both the flux form and the entrainment form of the Exner equation are implemented in the model for sediment mass conservation.For each form of the Exner equation, we consider both the cases of uniform sediment (bed material characterized by a single grain size) and sediment mixtures.Since the sediment is very fine in the LYR, the component of the load that is bed load is likely negligible (e.g., Ma et al., 2017), so we consider only the transport of suspended load.Considering the fact that most accepted sediment transport relations (e.g., the Engelund and Hansen, 1967, relation) underpredict the sediment transport rate of the LYR by an order of magnitude or more (Ma et al., 2017), in our model we implement two recently developed generalized versions of the Engelund-Hansen relation, which are based on data from the LYR.These are the version of Ma et al. (2017) for uniform sediment and the version of Naito et al. (2018) for sediment mixtures.In cases considering sediment mixtures, we also implement the method of Viparelli et al. (2010) to store and access bed stratigraphy as the bed aggrades and degrades.
Since the aim of this paper is to compare the two formulations of the Exner equation in the context of the LYR rather than reproduce site-specific morphodynamic processes of the LYR, some additional simplifications are introduced to the model to facilitate comparison.The channel is simplified to be a constant-width rectangular channel, and bank (sidewall) effects and floodplain interactions are not considered.The channel bed is assumed to be an infinitely deep supplier of erodible sediment with no exposed bedrock, which is justifiable because the LYR is fully alluvial and has been aggrading for thousands of years, as copiously documented in Chinese history.Finally, water and sediment (of each grain size range) are fed into the upstream boundary at a specified rate, and at the downstream end of the channel we specify a fixed bed elevation along with a normal flow depth.These restrictions could be easily relaxed so as to incorporate site-specific complexities of the Yellow River.Because of the severe aggradation of the LYR developed before the Xiaolangdi Dam operation, the LYR is famous for its hanging bed (i.e., bed elevated well above the floodplain) and no major tributaries need be considered in the simulation.

Flow hydraulics
Flow hydraulics in a rectangular channel are described by the following 1-D Saint-Venant equations, which consider fluid mass and momentum conservation, where t is time, h is water depth, q w is flow discharge per unit width, g is gravitational acceleration, S is bed slope, u is depth-averaged flow velocity, C f is dimensionless bed resistance coefficient, and C z is the dimensionless Chézy resistance coefficient.In our model, the fully unsteady 1-D Saint-Venant equations are solved using a Godunov-type scheme with the HLL (Harten-Lax-van Leer) approximate Riemann solver (Harten et al., 1983;Toro, 2001), which can effectively capture discontinuities in unsteady and nonuniform open channel flows.In this paper, the full flood hydrograph of the LYR is replaced by a flood intermittency factor I f (Paola et al., 1992;Parker, 2004).According to this definition, the river is assumed to be at low flow and not transporting significant amounts of sediment for time fraction 1−I f and is in flood at constant discharge and active morphodynamically for time fraction I f .In the long term, the relation between the flood timescale t f and the actual timescale t is t f = I f t.With the consideration that a river is in flood only for a fraction of time, here we introduce I f into the time derivative of all governing equations so that the flood timescale t f is implemented in the simulation.This notwithstanding, the results we exhibit later in this paper are all cast in terms of actual timescale t.Full hydrographs are considered in the Supplement.

Flux form of the Exner equation
When dealing with uniform sediment, the flux form of the Exner equation can be written as where λ p is the porosity of the bed deposit, and z b is bed elevation.Sediment transport is regarded to be in a quasiequilibrium state so that the sediment transport rate per unit width q s equals the equilibrium (capacity) sediment transport rate per unit width q se .When considering sediment mixtures, an active layer formulation (Hirano, 1971;Parker, 2004) is incorporated in the flux-based Exner equation so that the evolution of both bed elevation and surface grain size distribution can be considered.In this formulation, the riverbed is divided into a wellmixed upper active layer and a lower substrate with vertical stratigraphic variations.The upper active layer therefore represents the volume of sediment that interacts directly with suspended load transport and also exchanges with the substrate as the bed aggrades and degrades.Discretizing the grain size distribution into n ranges, the mass conservation relation for each grain size range can be written as where q si is volumetric sediment transport rate per unit width of the ith grain size range (taken to be equal to its equilibrium value q sei in the flux formulation), F i is the volumetric fraction of surface material in the ith grain size range, f I i is the volumetric fraction of material in the ith grain size range exchanged across the surface-substrate interface as the bed aggrades or degrades, and L a is the thickness of the active layer.
For bedform-dominated sand-bed rivers, L a is often related to the height of dunes (Blom, 2008) so that the vertical sorting processes due to bedform migration can be considered.
In this paper, a constant value of L a is implemented in the simulation.Summing Eq. ( 5) over all grain size ranges, one can find that the governing equation for bed elevation in the case of sediment mixtures is the same as Eq. ( 4) upon replacing q s with q sT = q si , where q sT denotes the total sediment transport rate per unit width summed over all size ranges.Reducing Eq. ( 5) with Eq. (4), we get Therefore, the flux formulation Eqs. ( 4) and ( 6) are implemented as governing equations for sediment mixtures, with Eq. ( 4) describing the evolution of bed elevation and Eq. ( 6) describing the evolution of surface grain size distribution.The exchange fractions f I i between the active layer and the substrate are calculated using the following closure relation.
That is, the substrate is transferred into the active layer during degradation, and a mixture of suspended load and active layer material is transferred into substrate during aggradation.In Eq. ( 7), f i | z b −L a is the volumetric fraction of substrate material just beneath the interface, p si = q si /q sT is the fraction of bed material load in the ith grain size range, and α is a specified parameter between 0 and 1.The formulation is adapted from Hoey and Ferguson (1994) and Toro-Escobar et al. (1996), who originally used it for bed load.In this paper, a value of 0.5 is specified for α.
The method of Viparelli et al. (2010) is applied in our model to store substrate stratigraphy and provide information for f i | z b −L a (i.e., the topmost sublayer in Viparelli et al., 2010).The reader can refer to the original reference of Viparelli et al. (2010) for more details or refer to An et al. (2017) for a concise description of how to implement this method in a morphodynamic model.When solving the flux form of the Exner equation, a first-order upwind scheme is implemented to discretize the spatial derivatives, and a first-order explicit scheme is implemented to discretize the temporal derivatives.

Entrainment form of the Exner equation
The entrainment-based Exner equation for uniform sediment is In Eq. ( 8), v s is the fall velocity of sediment particles; E is the dimensionless entrainment rate of sediment normalized by sediment fall velocity; C is the depth-flux-averaged volume sediment concentration; and r o = c b /C is the recovery coefficient of suspended load, which denotes the ratio between the near-bed sediment concentration c b and the fluxaveraged sediment concentration C. By definition, r 0 is related to the concentration profile of suspended load and is expected to be no less than unity in cases appropriate for a depth-averaged shallow-water treatment of flow and morphodynamics.Therefore, the first term on the right-hand side of Eq. ( 8), i.e., v s • E, denotes the sediment entrainment rate per unit area; the second term on the right-hand side of Eq. ( 8), i.e., v s • r 0 • C, denotes the sediment deposition rate per unit area.
For the sediment fall velocity v s , we compare two widely used relations: the relation of Dietrich (1982) and the relation of Ferguson and Church (2004).Results show that these two relations give almost the same fall velocity for the bed material load of the LYR, whose grain sizes typically fall in the range of 15 to 500 µm.Therefore, only the relation of Dietrich (1982) is implemented in our simulations in this paper.Readers can refer to Sect.S1 of the Supplement for more details.
In the entrainment formulation the sediment transport rate q s is not necessarily in its equilibrium state, but the dimensionless entrainment rate E is taken to be at capacity.The sediment transport rate q s is calculated according to the following continuity relation.
For the dimensionless entrainment rate E, we assume that sediment transport reaches its equilibrium state (q s = q se ) when the sediment deposition rate and the sediment entrainment rate balance each other (r 0 C = E).Therefore, E can be back-calculated from q se as For the depth-flux-averaged sediment concentration C, another equation is implemented describing the conservation of suspended sediment in the water column.
The entrainment-form Exner equation for sediment mixtures also uses the active layer formulation described in Sect.2.2.Mass conservation of each grain size range can be written as where the subscript i denotes the ith size range of sediment grain size.Summing Eq. ( 12) over all grain size ranges, we get the governing equation for bed elevation.
Reducing Eq. ( 12) with Eq. ( 14) we get the governing equation for surface fraction F i . 1 C. An et al.: Morphodynamic model of the lower Yellow River The governing equation for the sediment concentration of each grain size C i can be written as and the sediment transport rate per unit width for the ith size range q si obeys the following continuity relation: In the entrainment formulation, the closure relation for f I i is the same as that used in the flux formulation (i.e., Eq. 7), and the substrate stratigraphy is also stored and accessed using the method of Viparelli et al. (2010).When discretizing the entrainment form of the Exner equation, a first-order upwind scheme is implemented for the spatial derivatives, and a first-order explicit scheme is implemented for the temporal derivatives.

Uniform sediment
To close the Exner equations described in Sect.2.2 and 2.3, equations for equilibrium sediment transport rate q se (q sei ) are still needed.For the simulations using uniform sediment, we implement the generalized Engelund-Hansen relation proposed by Ma et al. (2017).This equation is based on data from the LYR and can be written in the following dimensionless form: where q * s is the dimensionless sediment transport rate per unit width (i.e., the Einstein number), and τ * is dimensionless shear stress (i.e., the Shields number).They are defined as where D is the characteristic grain size of the bed sediment (here approximated as uniform); τ b is bed shear stress; and R is the submerged specific gravity of sediment defined as (ρ s −ρ)/ρ, in which ρ s is the density of sediment, and ρ is the density of water.The sediment submerged specific gravity R is specified as 1.65 in this paper, which is an appropriate estimate for natural rivers and corresponds to quartz.
In the relation of Ma et al. (2017), the dimensionless coefficient α s = 0.9 and the dimensionless exponent n s = 1.68.These values are quite different from the original relation of Engelund and Hansen (1967), in which α s = 0.05 and n s = 2.5.Ma et al. (2017) demonstrated that such differences imply that the riverbed of the LYR is dominated by lowamplitude bedform features (dunes) approaching the upperregime plane bed.According to this finding, form drag is then neglected in our modeling, and all of the bed shear stress is used for sediment transport.

Sediment mixtures
We implement the relation of Naito et al. (2018) to calculate the equilibrium sediment transport rate of size mixtures.Using field data from the LYR, Naito et al. (2018) extended the Engelund and Hansen (1967) relation to a surface-based grain-size-specific form, in which the suspended load transport rate of the ith size range is tied to the availability of this size range on the bed surface: where N * i is the dimensionless sediment transport rate in the ith size range, and u * is shear velocity calculated from the bed shear stress τ b .
The transport relation itself takes the form in which D i is the characteristic grain size for sediment in the ith size range, D sg is the geometric mean grain size in the active layer, and τ * g is the dimensionless bed shear stress associated with D sg .The parameters τ * g , coefficient A i , and exponent B i are calculated as follows.
If A i and B i are specified as constant values in Eq. ( 24), then the sediment transport rate for each size range depends only on the flow shear stress and the characteristic grain size of this size range, without being affected by other size ranges.But according to Eqs. ( 26) and ( 27), the coarser the sediment the smaller the values of A i and B i will be, thus leading to reduced mobility for coarse sediment (and increased mobility for fine sediment) due to the presence of grains of other sizes.Thus the relations in Eqs. ( 26) and ( 27) serve as a hiding function that allows for grain sorting.We note that a form of the Engelund-Hansen equation for mixtures was introduced by Van der Scheer et al. (2002) and implemented by Blom et al. (2016).Blom et al. (2017) further extended this relation to a more general framework capable of including hiding effects.These forms, however, have not been calibrated to the LYR data and are thus not suitable for the LYR.

Numerical modeling of the LYR using the two forms of the Exner equation
In this section, we conduct numerical simulations using both the flux form and the entrainment form of the Exner equation, with the aim of studying under what circumstances the two forms give different predictions.Numerical simulations are conducted in the setting of the LYR.We specify a 200 km long channel reach for our simulations, along with a constant channel width of 300 m and an initial longitudinal slope of 0.0001.Bed porosity λ p is specified as 0.4.Based on field measurements of the LYR available to us, we implemented a dimensionless Chézy resistance coefficient C z of 30, which corresponds to a dimensionless bed resistance coefficient C f of 0.0011.For the entrainment form of the Exner equation, we specify the ratio of near-bed sediment concentration to flux-averaged sediment concentration r 0 (r 0i ) = 1.Such a value of r 0 (r 0i ) corresponds to a vertically uniform profile of sediment concentration and will thus give a maximum difference between the prediction of the entrainment form and the prediction of the flux form.More discussion about the effects of r 0 is presented in Sect.4.3.
A constant flow discharge of 2000 m 3 s −1 (corresponding to a flow discharge per unit width q w of 6.67 m 2 s −1 ) is introduced at the inlet of the channel with the flood intermittency factor I f estimated as 0.14 (Naito et al., 2018).The downstream end is specified far from the river mouth to neglect the effects of backwater.Therefore, the bed elevation is held constant and the water depth is specified as the normal flow depth at the downstream end of the calculational domain.The above flow discharge per unit width q w combined with the bed slope S as well as the bed resistance coefficient C f leads to a normal flow depth of 3.69 m.In our simulation, we use the height of bedforms in the LYR to determine the thickness of the active layer (Blom, 2008).According to the field survey of Ma et al. (2017), the characteristic height of bedforms in the LYR is about 20 % of the normal flow depth, which can fall in the range suggested by the data analysis of Bradley and Venditti (2017).This eventually leads to an estimate of active layer thickness of L a = 0.738 m.The sublayer in the substrate to store the vertical stratigraphy is specified with a thickness of 0.5 m.
Two cases are considered here.In the first case, the sediment grain size distribution of the LYR is simplified to a uniform grain size of 65 µm.This is based on the measured grain size distribution of bed material at the Lijin gauging station, which has a median grain size of D 50 = 66.6 µm, a geometric mean grain size of D g = 65.5 µm, and a geo-metric standard deviation σ g = 2.0, as shown in Fig. 1c.In the second case, we consider the effects of sediment mixtures.The grain size distribution of the initial bed is based on the bed material at the Lijin gauging station, as shown in Fig. 1c, but we renormalize the measured grain size distribution with a cutoff for wash load at 15 µm as suggested by Ma et al. (2017).The renormalized grain size distribution for the initial bed as implemented in the case of sediment mixtures is shown in Fig. 2, with a total number of grain size fractions of 5.In both cases, simulations start with an equilibrium state in which sediment supply rate, sediment transport rate, and equilibrium sediment transport rate are the same so that the initial state of the channel is in equilibrium.Then we cut the sediment supply rate (of each size range) to only 10 % of the equilibrium sediment transport rate and keep this sediment supply rate.This is to mimic the reduction of sediment load in the LYR in recent years, as shown in Fig. 1b.The grain size distribution of sediment supply in the case of sediment mixtures is shown in Fig. 2.
The 200 km channel reach is discretized into 401 cells, with cell size x of 500 m.In the case of uniform sediment, we specify a time step for morphologic calculation t m = 10 −4 years and a time step for hydraulic calculation t h = 10 −6 years.In the case of sediment mixtures, we specify a time step for morphologic calculation t m = 10 −5 years and a time step for hydraulic calculation t h = 10 −6 years.Computational conditions are briefly summarized in Table 1.The computational conditions we implement are much simpler than the rather complicated conditions of the actual LYR.But it should be noted that the aim of this paper is not to reproduce specific aspects of the morphodynamic processes of the LYR, but to compare the flux form and entrainment form of the Exner equation in the context of conditions typical of the LYR.

Case of uniform sediment
In this case, we implement a uniform grain size of 65 µm for both the bed material and sediment supply.Such a grain size is nearly equal to the observed median grain size (or geometric mean grain size) of bed material at the Lijin gauging station.The relation of Ma et al. (2017) is implemented to calculate the transport rate of bed material suspended load.This relation provides an equilibrium sediment transport rate per unit width q se of 0.0136 m 2 s −1 under the given flow discharge, bed slope, and sediment grain size.With a flood intermittency factor I f of 0.14, this further gives a mean annual bed material load of 47.8 Mt a −1 .Adding in wash load according to the estimate of Naito et al. (2018), the total mean annual load is 86.9 Mt a −1 , a value that is of the same order of magnitude as averages over the period 2000-2016 (89-126 Mt a −1 depending on the site), i.e., since the operation of the Xiaolangdi Dam beginning in 1999 (Fig. 1b).The sediment supply rate q sf we specify at the upstream end of the channel is only 10 % of the equilibrium sediment transport www.earth-surf-dynam.net/6/989/2018/Earth Surf.Dynam., 6, 989-1010, 2018  rate (i.e., the sediment supply rate is cut by 90 % from the equilibrium state) such that q sf = 0.00136 m 2 s −1 .
Figure 3 shows the modeling results using the flux form of the Exner equation.As we can see in the figure, the bed degrades and the sediment load decreases in response to the cutoff of sediment supply.Such adjustments start from the upstream end of the channel and gradually migrate downstream.Figure 4 shows the modeling results using the entrainment form of the Exner equation.A comparison between Figs. 4 and 3 shows that the entrainment form and the flux form give very similar predictions in this case.The entrainment form provides a somewhat slower degradation (at the upstream end the flux form predicts a 3 m degradation, whereas the entrainment form predicts a 2.3 m degradation) and a more diffusive sediment load reduction.Such more diffusive predictions of sediment load variation can be ascribed to the condition of nonequilibrium transport that is embedded in the entrainment form.This issue will be studied analytically in Sect. 4.Here we present the results for only 0.2 years after the cutoff of sediment supply, since the differences between the predictions of the two forms tend to be the most evident shortly after the disruption but gradually diminish as the river approaches the new equilibrium (El kadi Abderrezzak and Paquier, 2009).Modeling results over a longer timescale will be discussed in Sect.4.3.
To further quantify the differences between the predictions of the two forms, we propose the following normalized parameter: where y denotes an arbitrary variable calculated by the morphodynamic model, and subscripts F and E denote results using the flux form and the entrainment form, respectively.Therefore, δ (y) denotes the difference between the prediction of the two forms y F and y E normalized by the prediction of the flux form y F .Table 2 gives a summary of the maximum values of δ along the channel at different times in the case of uniform sediment.The values of δ for both z b and q s are presented.As we can see from the table, the maximum value of δ(z b ) along the calculational domain stays within 4 % in the first 0.2 years after the cutoff of sediment supply.This indicates that the flux form and the entrainment form can indeed give almost the same prediction in terms of bed elevation in this case.But in the case of the sediment load per unit width q s , the maximum value of δ(q s ) can be as high as 20 %, indicating that even though the two forms give qualitatively similar patterns of evolution in terms of sediment load as shown in Figs. 3 and 4, a quantitative difference is clearly evident due to the more diffusive nature of the predictions of the entrainment form.The value of δ(q s ) is largest at the beginning of the simulation and is then gradually reduced with time.It should be noted that the values of δ(z b ) depend on the choice of elevation datum.In this paper bed elevation at the downstream end is fixed as 0 m, which serves as the elevation datum.In the simulation of this paper, the maximum value of δ(z b ) almost always occurs at the upstream end where bed elevation does not deviate far from the initial value of 20 m.
The above results show that the flux form and the entrainment form can provide similar predictions of the LYR when the bed sediment grain size distribution is simplified to a uniform value of 65 µm.To understand under what conditions the two forms will lead to more different results, we conduct an idealized run using the entrainment form in which the sediment fall velocity v s is arbitrarily multiplied by a factor of 0.05.That is to say, we keep the sediment grain size at 65 µm in the computation of the Shields number, but let the sediment fall velocity in Eqs. ( 8) and ( 10) equal only 1/20 of the value calculated by the relation of Dietrich (1982) from this grain size.With a much smaller, and indeed intentionally unrealistic, sediment fall velocity the entrainment form predicts very different results as shown in Fig. 5.The adjustment of the sediment load becomes even more diffusive in space: it takes almost the entire 200 km reach for the sediment load to adjust from the upstream disruption to the equilibrium transport rate.Meanwhile, there is barely any bed degradation at the upstream end after 0.2 years, in correspondence with the fact that the spatial gradient of q s becomes quite small.In Table 2 we also exhibit the δ values for this idealized run.It is no surprise that both δ(z b ) and δ(q s ) are high, as the entrainment form and flux form predict very different patterns with such an arbitrarily reduced sediment fall velocity.
In Sect.S2, we also conduct numerical simulations with hydrographs.Results indicate that our conclusions based on constant flow discharge also hold when hydrographs are considered: the flux form and the entrainment form (with the sediment fall velocity not adjusted) of the Exner equation give very similar predictions using a characteristic grain size of 65 µm.

Case of sediment mixtures
In this section we consider the morphodynamics of sediment mixtures rather than the case of a uniform bed grain size implemented in Sect.3.1.The grain size distribution of the initial bed is based on field data at the Lijin gauging station and is shown in Fig. 2. Using the sediment transport relation of Naito et al. (2018) for mixtures, such a grain size distribution combined with the given bed slope and flow discharge leads to a total equilibrium sediment transport rate per unit width q seT of 0.0272 m 2 s −1 .With a flood intermittency factor I f of 0.14, this further gives a mean annual bed material load of 95.5 Mt a −1 .Adding in wash load according to the estimate of Naito et al. (2018), the total mean annual load is 173.7 Mt a −1 , a value that is of the same order of magnitude as averages over the period 2000-2016 (89-126 Mt a −1 depending on the site), i.e., since the operation of the Xiaolangdi Dam beginning in 1999 (Fig. 1b).The sediment supply rate of each grain size range is set at 10 % of its equilibrium sediment transport rate.This results in a total sediment supply rate of q sf = 0.00272 m 2 s −1 and a grain size distribution of the sediment supply (shown in Fig. 2) that is identical to the grain size distribution of the equilibrium sediment load before the cutoff.That is, the grain size distribution of sediment supply does not change; only the total sediment supply is reduced by 90 %.Again we exhibit simulation results for only 0.2 years here, a value that is enough to show the differences between the two forms, flux and entrainment, as applied to mixtures.Modeling results over a longer timescale are presented in Sect.4.3.
Figure 6 shows the simulation results using the flux form of the Exner equation.As a result of the reduced sediment supply at the inlet, bed degradation occurs first at the upstream end and then gradually migrates downstream.The total sediment transport rate per unit width q sT is also reduced as a response to the cutoff of sediment supply.More specifically, the evolution of q sT shows marked evidence of advection, with at least two kinematic waves being observed within 0.2 years.Actually, as illustrated by Stecca et al. (2014Stecca et al. ( , 2016)), each grain size fraction should induce a migrating wave.As shown in Fig. 6b, the fastest kinematic wave migrates beyond the 200 km reach within 0.06 years, and the second fastest kinematic wave migrates for a distance of about 60 km in 0.2 years.Figure 6c and d show the results for the surface geometric mean grain size D sg and geometric mean grain size of suspended load D lg , respectively.As can be seen therein, both the bed surface and the suspended load coarsen as a result of the cutoff of sediment supply.This represents armoring mediated by the hiding functions of Eqs. ( 26) and ( 27).Such coarsening is not evident near the upstream end, possibly due to the inverse slope visible in Fig. 6a.Similarly to the variation of q sT , the patterns of time variation of both D sg and D lg also exhibit very clear kinematic waves, with migration rates about the same as those of q sT .
Figure 7 shows the simulation results obtained using the entrainment form of the Exner equation.In general, the patterns of variation predicted by the entrainment form have similar trends and magnitudes to those predicted by the flux form: the bed degrades near the upstream end, the suspended load transport rate is reduced in time, and both the bed surface and the suspended load coarsen as a result of the cutoff of sediment supply.But the results based on the two forms exhibit very evident differences when multiple grain sizes are included.That is, the results predicted by the entrainment form are sufficiently diffusive so that the variations of q sT , D sg , and D lg (Fig. 7b, c, and d) do not show the advective character seen in Fig. 6. Figure 7c, however, shows the same armoring as in the case of calculations with the flux form.No clear kinematic waves can be observed in Fig. 7. Table 3 gives a summary of the values of δ in the case of sediment mixtures.The prediction of bed elevation is not affected much when multiple grain sizes are considered, with δ(z b ) being no more than 3.5 % within 0.2 years.The δ values of q sT , D sg , and D lg are, however, relatively large since the two forms predict quite different patterns of variations, as shown in Figs. 6 and 7.
The results shown in Fig. 8 have also been calculated using the entrainment form of the Exner equation, but here the sediment fall velocities v si used in Eqs. ( 14)-( 16) are arbitrarily multiplied by a factor of 20.That is, we still apply the grain size distribution in Fig. 2, but the sediment fall velocities implemented in the simulation are 20 times the corresponding fall velocities calculated by the relation of Dietrich (1982).In the case of uniform sediment in Sect.3.1, we arbitrarily reduce the sediment fall velocity to force a difference between the predictions from the entrainment form and those from the flux form.Here we arbitrarily increase the sediment fall velocity with the aim of determining under what conditions the sorting patterns predicted by the two forms converge.As we can see in Fig. 8, with such a larger and intentionally unrealistic sediment fall velocity, the general trend of variations predicted by the entrainment form does not change, but the results show a notably less diffusive pattern.The variations of q sT , D sg , and D lg show more advection compared with Fig. 7, and at least two kinematic waves appear within 0.2 years.It should be noted that even though these kinematic waves appear after we arbitrarily increase the sediment fall velocity, they are more diffusive than those obtained from the flux formulation and also migrate with a slower celerity compared with those predicted by the flux form, especially for the fastest kinematic wave in the modeling results.
Table 3 summarizes the δ values for this run.The values of δ(z b ) become smaller with arbitrarily increased sediment fall velocities except for t = 0.06 years.A relatively large value of δ(z b ) at t = 0.06 years occurs near the downstream end of the channel, where the entrainment form predicts some slight degradation.Also, δ(q sT ) is quite large at t = 0.01 years and 0.03 years, even though the results for the case of increased fall velocities become qualitatively more similar to the prediction of the flux form.This is because the flux form and the entrainment form with arbitrarily increased sediment fall velocities predict different celerities for the fastest kinematic wave.The error δ(q sT ) becomes smaller from t = 0.06 years as the fastest kinematic wave migrates beyond the channel reach.The error δ(D lg ) behaves similarly to δ(q sT ), with δ(D lg ) being quite large at t = 0.01 years and 0.03 years near the fastest kinematic wave, but gradually becoming smaller as time passes.The error δ(D sg ) stays low within the whole 0.2-year period, possibly because the fastest kinematic wave of D sg has a small magnitude, as shown in Fig. 8c.
In Sect.S3, we present additional numerical cases that are similar to the cases in this section, except that hydrographs are implemented instead of constant discharge.Results indicate that our conclusions based on constant flow discharge also hold when hydrographs are considered.The flux form and the entrainment form (with the sediment fall velocity not adjusted) of the Exner equation predict quite different patterns of grain sorting, with the flux form exhibiting a more advective character than the entrainment form.

Adjustment of sediment load and the adaptation length
In Sect.3.1, our simulation shows that in the case of uniform sediment, the flux form and the entrainment form of the Exner equation give very similar predictions for a given sediment size of 65 µm.However, if we arbitrarily reduce the sediment fall velocity by a multiplicative factor of 0.05, the prediction given by the entrainment form will become much more diffusive in terms of both z b and q s .The diffusive nature of the entrainment form and the important role played by the sediment fall velocity can be explained in terms of the governing equation.
In the entrainment form, the equation governing suspended sediment concentration is i.e., the same as Eq. ( 11).The sediment transport rate per unit width q s = huC = q w C, and the dimensionless entrainment rate E = r 0 q se /q w .In order to simplify the mathematical analysis, here we consider only the adjustment of sediment concentration in space and neglect the temporal derivative in Eq. ( 29) so that we get www.earth-surf-dynam.net/6/989/2018/where L ad can be identified as the adaptation length for suspended sediment to reach equilibrium.This definition of adaptation length is similar to those in Wu and Wang (2008) and Ganti et al. (2014).
If we consider the spatial adjustment of sediment load shortly after the cutoff of sediment supply, we can further neglect the nonuniformity of the capacity (equilibrium) transport rate q se along the channel, and Eq. ( 30) can be solved with a given upstream boundary condition.That is, with the boundary condition Eq. ( 30) can be solved to yield Here q sf is the sediment supply rate per unit width at the upstream end.According to Eq. ( 33), q s adjusts exponentially in space from q sf to q se , which also coincides with our simu-Earth Surf.Dynam., 6, 989-1010, 2018 www.earth-surf-dynam.net/6/989/2018/lation results in Sect.3.1, as shown in Figs.3-6.The adaptation length L ad is the key parameter that controls the distance for q s to approach the equilibrium sediment transport rate q se .More specifically, q s attains 1 − 1/e (i.e., 63.2 %) of its adjustment from q sf to q se over a distance L ad .Therefore, the larger the adaptation length, the slower q s adjusts in space so that the more evident lag effects and diffusivity are exhibited in the entrainment form.In the flux form, however, the sediment load responds simultaneously with the flow conditions so that L ad = 0 and q s = q se along the entire channel reach.
For the case of uniform sediment in Sect.3.1, q w = 6.67 m 2 s −1 and r o is specified as unity.Therefore, the value of L ad is determined only by the sediment fall velocity v s .Figure 9 shows the value of the adaptation length L ad for various sediment grain sizes, with the sediment fall velocity v s calculated by the relation of Dietrich (1982).From the figure we can see that L ad decreases sharply with the increase in grain size, indicating that the lag effects between sediment transport and flow conditions are evident for very fine sediment but gradually disappear when sediment is sufficiently coarse.For the sediment grain size of 65 µm implemented in Sect.3.1, the corresponding L ad = 1.88 km, which is much smaller than the 200 km reach of the computational domain.In this case and in general, the predictions of the flux form and the entrainment form show little difference when L ad /L 1, where L is domain length.However, if we arbitrarily multiply the sediment fall velocity by a factor of 0.05, then L ad becomes 37.60 km.With such a large adaptation length, it is no surprise that the entrainment form gives very different predictions from the flux form.
The evolution of bed elevation z b can also be affected by the value of L ad .For example, in the case of uniform sediment in Sect.3.1, the flux form corresponds to an adaption length of zero.As a result, the flux form yields a spatial derivative of q s near the upstream end that is relatively large, thus leading to fast degradation from the upstream end.In the case of the entrainment form, however, the spatial derivative of q s is small with a large L ad , thus leading to a slower and more diffusive bed degradation.This is especially evident when we arbitrarily reduce the sediment fall velocity by a factor of 0.05 while keeping grain size invariant.
The above analysis also holds for sediment mixtures, except that each grain size range will have its own adaptation length.Here we neglect the temporal derivative in Eq. ( 29) and analyze only the spatial adjustment of sediment load.If we neglect the spatial derivative in Eq. ( 29) and conduct a similar analysis for sediment concentration, we would find that the temporal adjustment of sediment concentration is www.earth-surf-dynam.net/6/989/2018/Earth Surf.Dynam., 6, 989-1010, 2018 also described by an exponential function of time, in analogy to Eq. (33).

Patterns of grain sorting: advection vs. diffusion
In Sect.3.2 we find that the flux form and entrainment form of the Exner equation provide very different patterns of grain sorting for sediment mixtures: kinematic sorting waves are evident in the flux form but are diffused out in the entrainment form.The diffusivity of grain sorting becomes smaller and the kinematic waves appear, however, if we arbitrarily increase the sediment fall velocity by a factor of 20.In this section, we explain this behavior by analyzing the governing equations.
First we rewrite the sediment transport relation of Naito et al. (2018) in the following form.q sei = F i q ri (34) Substituting Eq. (34) into Eq.( 6), which is the governing equation for surface fraction F i in the flux form, we get Equation ( 36) can be written in the form of a kinematic wave equation with source terms as follows.
Earth Surf.Dynam., 6, 989-1010, 2018 www.earth-surf-dynam.net/6/989/2018/Here c Fi is the ith celerity of a kinematic wave and SF i denotes source terms.Since the surface geometric mean grain size D sg , the total sediment load per unit width q sT (which equals the equilibrium sediment transport rate q seT ), and the geometric mean grain size of sediment load D lg are all closely related to the surface grain size fractions F i , the evolution of these three parameters shows marked advective behavior when simulated by the flux form of the Exner equation.However, the evolution of bed elevation z b is related to ∂q sT /∂x, which is dominated by diffusion if q sT is predominantly slope dependent (as is the case here).The advectiondiffusion character of the flux form of the Exner equation for sediment mixtures has been documented thoroughly in a series of papers (e.g., Stecca et al., 2014Stecca et al., , 2016;;An et al., 2017).
The reader can reference these papers for more details.Now we turn to the entrainment form of the Exner equation.Combined with the sediment transport rate per unit width q si = huC i = q w C i and the dimensionless entrainment rate E i = r 0i q sei /q w , Eqs. ( 16) and ( 15) can be written as where Eq. ( 40) denotes the conservation of suspended sediment and Eq. ( 41) denotes the conservation of bed material.
If we rewrite Eq. ( 40) in the form then q si can be solved iteratively.With an initial guess of q si = q sei and neglecting the temporal derivatives, we obtain the second-order solution of q si as Details on the iteration scheme are given in Sect.S4.Substituting Eqs. ( 43) and (34) into Eq.( 41), we find that Expanding out the last two terms in Eq. ( 44) using the chain rule, after some work the relation for the conservation of bed material can be expressed as follows.
Here c Ei is the celerity of a kinematic wave, ν i is the diffusivity coefficient, and SE i denotes source terms.From Eq. ( 45) we can see that the governing equation for F i in the entrainment form is an advection-diffusion equation rather than the kinematic wave equation of the flux form.The surface geometric mean grain size D sg is governed by Eq. ( 45), which describes the variation of the surface fractions F i from which it is computed.The equilibrium sediment transport rate q sei is governed by Eq. ( 45) because we implement a surface-based sediment transport relation as shown in Eq. ( 34).According to Eq. ( 43), the total sediment load per unit width q sT and the geometric mean grain size of sediment load D lg must also be closely related to the surface grain size fractions F i .Therefore, the diffusion terms in Eq. ( 45) can lead to dissipation of the kinematic waves in Fig. 7b, c, and d.
From Eq. ( 47), we can also see that the diffusivity coefficient ν i is related to the sediment fall velocity v si : the larger www.earth-surf-dynam.net/6/989/2018/Earth Surf.Dynam., 6, 989-1010, 2018 the sediment fall velocity, the smaller the diffusivity coefficient.Thus when we increase the sediment fall velocity arbitrarily by a factor of 20 in Sect.3.2, the kinematic waves become more evident as a result of the reduction of diffusivity.
Moreover, if we compare the celerity of kinematic waves in both the flux form and the entrainment form, we have where L adi is the adaptation length for the ith size range as defined by Eq. ( 31).More specifically, the value of r ci depends on ∂q ri /∂x.For our numerical simulation in Sect.3.2, ∂q ri /∂x > 0 as a result of bed degradation progressing from the upstream end, thus leading to a positive value of r ci and an entrainment celerity c Ei that is smaller than the corresponding flux celerity c Fi .This is consistent with our numerical results: the kinematic waves in Fig. 8 predicted by the entrainment form are somewhat smaller than the kinematic waves in Fig. 6 predicted by the flux form.

Modeling implications and limitations
In Sect.3, two numerical cases are presented to compare the flux form and the entrainment form of the Exner equation, but only within 0.2 years after the cutoff of sediment supply.Here we run both numerical cases for a longer time (5 years).
Table 4 shows the results of the case of uniform sediment (as described in Sect.3.1) within 5 years, and Table 5 shows the results of the case of sediment mixtures (as described in Sect.3.2) within 5 years.For both cases, the δ values corresponding to relative deviation between the flux and entrainment forms become quite small after 1 year, thus validating our assumption that the predictions of the two forms tend to be most evident shortly after disruption, but gradually diminish over a longer timescale.Moreover, if the water and sediment supply are kept constant for a sufficiently long time, the flux form and entrainment form of the Exner equation predict exactly the same equilibrium in terms of both the channel slope and the bed surface texture.Under such conditions, the sediment transport rate (of each size range) equals the equilibrium sediment transport rate (of each size range) and also equals the sediment supply rate (of each size range).
Based on the numerical modeling and mathematical analysis in this paper, we suggest that the entrainment form of the Exner equation be used when studying the river morphodynamics of fine-grained sediment (or, more specifically, sediment with small fall velocity).This is because the adaptation length L a and the diffusivity coefficient ν i are large for fine sediment, but the flux form of the Exner equation does not account for lag effects or the diffusivity of individual size fractions, thus leading to unrealistic simulation results.Such unrealistic simulation results can include an It should be noted that in the morphodynamic models of this paper, we implement the mass and momentum conservation equations for clear water (i.e., Eqs. 1 and 2) to calculate flow hydraulics instead of the mass and momentum equations for water-sediment mixture as suggested by Cao et al. (2004Cao et al. ( , 2006)).More specifically, Cui et al. (2005) have pointed out that when sediment concentration in the water is sufficiently small, bed elevation can be taken to be unchanging over characteristic hydraulic timescales, and the effects of flow-bed exchange on flow hydraulics can be neglected.For the two simulation cases in this paper, the volume of sediment concentration C drops from about 2 × 10 −3 to about 2 × 10 −4 in the case of uniform sediment and from about 4 × 10 −3 to about 4 × 10 −4 in the case of sediment mixtures due to the cutoff of sediment supply at the upstream end.These dilute concentrations validate our implementation of mass and momentum conservation equations for clear water.Our assumption is not necessarily correct for the entire Yellow River.Upstream of our study reach, especially upstream of Sanmenxia Dam, the flow is often hyperconcentrated (Xu, 1999).
Considering the fact that in our numerical simulations a constant inflow discharge (along with a flood intermittency factor) is implemented, and also considering that the morphodynamic timescale is much larger than the hydraulic timescale in our case, the quasi-steady approximation or even the normal flow approximation can be introduced to further save computational efforts (Parker, 2004).But one thing that should be noted is that in our simulation results in Sect.3, the bed exhibits an inverse slope near the upstream end.The normal flow assumption becomes invalid under such circumstances, thus requiring a full unsteady shallow water model.
By definition, the recovery coefficient r o is the ratio of the near-bed to the flux-depth-averaged concentration of suspended load and is thus related to the concentration profile.In our simulation r 0 is specified as unity.That is, density stratification effects of suspended sediment are neglected, and the vertical profile of sediment concentration is regarded as uniform.However, in natural rivers, the value of r 0 can vary significantly under different circumstances (Cao et al., 2004;Duan and Nanda, 2006;Zhang and Duan, 2011;Zhang et al., 2013).In general, the value of r 0 is no less than unity and can be as large as 12 (Zhang and Duan, 2011).Therefore, according to our mathematical analysis in Sect.4.1 and 4.2, r 0 = 1 corresponds to a maximum adaptation length L ad , a maximum diffusivity coefficient ν i , and a minimum ratio of celerities c Ei /c Fi , thus leading to the largest difference between the flux form and the entrainment form.When sediment concentration is sufficiently high, hindered settling effects reduce the sediment fall velocity.Considering the fact that the sediment concentrations considered in our simulation are fairly small, hindered settling effects are not likely significant.More study on stratification and hindered settling effects would be useful in the case of the LYR.
In this paper, a one-dimensional morphodynamic model with several simplifications is implemented to compare the flux-based Exner equation and the entrainment-based Exner equation in the context of the LYR.However, a site-specific model of the morphodynamics of the LYR without these simplifications would be much more complex.For example, in our 1-D simulation we observe bed degradation after the closure of the Xiaolangdi Dam, but we cannot resolve its structure in the lateral direction.In natural rivers, bed degradation is generally not uniform across the channel width, but may be concentrated in the thalweg.Moreover, the spatial variation of channel width and initial slope, which are not considered in this paper, are also important when considering applied problems.The abovementioned issues, even though not the aim of this paper, merit future research (e.g., He et al., 2012).Chavarrias et al. (2018) have reported that morphodynamic models considering mixed grain sizes may be subject to instabilities that result from complex eigenvalues of the system of equations.No such instabilities were encountered in the present work.

Conclusions
In this paper, we compare two formulations for sediment mass conservation in the context of the lower Yellow River, i.e., the flux form of the Exner equation and the entrainment form of the Exner equation.We represent the flux form in terms of the local capacity sediment transport rate and the entrainment form in terms of the local capacity entrainment rate.In the flux form of the Exner equation, the conservation of bed material is related to the stream-wise gradient of sediment transport rate, which is in turn computed based on the quasi-equilibrium assumption according to which the local sediment transport rate equals the capacity rate.In the entrainment form of the Exner equation, on the other hand, the conservation of bed material is related to the difference between the entrainment rate of sediment from the bed into the flow and the deposition rate of sediment from the flow onto the bed.A nonequilibrium sediment transport formulation is applied here so that the sediment transport rate can lag in space and time behind changing flow conditions.Despite the fact that the entrainment form is usually recommended for the morphodynamic modeling of the LYR due to its finegrained sediment, there has been little discussion of the differences in predictions between the two forms.
Here we implement a 1-D morphodynamic model for this problem.The fully unsteady Saint-Venant equations are implemented for the hydraulic calculation.Both the flux form and the entrainment form of the Exner equation are implemented for sediment conservation.For each formulation, we include the options of both uniform sediment and sediment mixtures.Two generalized versions of the Engelund-Hansen relation specifically designed for the LYR are implemented to calculate the quasi-equilibrium sediment transport rate (i.e., sediment transport capacity).They are the version of Ma et al. (2017) for uniform sediment and the version of Naito et al. (2018) for sediment mixtures.The method of Viparelli et al. (2010) is implemented to store and access bed stratigraphy as the bed aggrades and degrades.We apply the morphodynamic model to two cases with conditions typical of the LYR.
In the first case, a uniform bed material grain size of 65 µm is implemented.We study the effect of the cutoff of sediment supply, as occurred after the operation of the Xiaolangdi Dam beginning in 1999.We find that the flux form and the entrainment form give very similar predictions for this case.Through quantification of the difference between the two www.earth-surf-dynam.net/6/989/2018/Earth Surf.Dynam., 6, 989-1010, 2018 forms with a normalized measure of relative difference, we find that the difference in the prediction of bed elevation is quite small (< 4 %), but the difference in the prediction of sediment load can be relatively large (about 20 %) shortly after the cutoff of sediment supply.
The results for the case of uniform sediment can be explained by analyzing the governing equation of sediment load q s .In the flux form, the volume sediment transport rate per unit width q s equals the local equilibrium (capacity) value q se .In the entrainment form, however, we find that the difference between q s and q se decays exponentially in space.The adaptation length L ad = q w /(v s r 0 ) is the key parameter that controls the distance for q s to approach its equilibrium value q se .The larger the adaptation length, the more different the predictions of the two forms will be.For computational conditions in this case, the adaption length is relatively small (L ad = 1.88 km).
In the second case the bed material consists of mixtures ranging from 15 to 500 µm.We find that the flux form and the entrainment form give very different patterns of grain sorting.Evident kinematic waves occur at various timescales in the flux form, but no evident kinematic waves can be observed in the entrainment form.The different sorting patterns are reflected in the evolution of surface geometric mean grain size D sg , total sediment load q sT , and geometric mean grain size of sediment load D lg , but are not reflected in the evolution of bed elevation z b .
The different sorting patterns exhibited in the case of sediment mixtures can be explained by analyzing the governing equation for bed surface fractions F i , i.e., the grainsize-specific conservation of bed material.We find that in the flux form, the governing equation for F i can be written in the form of a kinematic wave equation.In the entrainment form, however, the governing equation for F i is an advection-diffusion equation.It is the diffusion term that leads to the dissipation of kinematic waves.Moreover, in the advection-diffusion equation arising from the entrainment form, the coefficient of diffusivity is inversely proportional to the sediment fall velocity.In addition, under the condition of bed degradation the wave celerity is smaller than that arising from the flux form.
Overall, our results indicate that the more complex entrainment form of the Exner equation might be required when the sorting processes of fine-grained sediment (or sediment with small fall velocity) is studied, especially at relatively short timescales.Under such circumstances, the flux form of the Exner equation might overestimate advection in sorting processes and the aggradation-degradation rate due to the fact that it cannot account for the relatively large adaptation length or diffusivity of fine particles.water depth I f flood intermittency factor L a thickness active layer L ad adaptation length of suspended load p si volumetric fraction of bed material load in the ith size range q ri normalized sediment transport rate per unit width for the ith size range, defined by Eq. ( 34) q s volumetric sediment transport rate per unit width q se equilibrium volumetric sediment transport rate (capacity) per unit width q sf sediment supply rate per unit width q w flow discharge per unit width R submerged specific gravity of sediment r 0 user-specified parameter denoting the ratio between the near-bed sediment concentration and the flux-averaged sediment concentration

Figure 1 .
Figure 1.(a) Sketch of the lower Yellow River showing six major gauging stations and the Xiaolangdi Dam.(b) Annual sediment load of the LYR measured at three gauging stations since 1950.(c) Grain size distributions of both bed surface material and suspended load measured at six gauging stations of the LYR.

Figure 2 .
Figure2.Grain size distributions of both the initial bed and the sediment supply in the case of sediment mixtures.For the initial bed, the surface and substrate grain size distributions are the same.The grain size distribution of the initial bed is renormalized based on the field data at the Lijin gauging station.The grain size distribution of the sediment supply equals the grain size distribution of bed material load at equilibrium.Grain sizes in the range of wash load have been removed from both distributions.

Figure 3 .
Figure 3.The 0.2-year results for the case of uniform sediment using the flux form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) sediment load per unit width q s of the LYR in response to the cutoff of sediment supply.The inset shows detailed results near the upstream end.

Figure 4 .
Figure 4.The 0.2-year results for the case of uniform sediment using the entrainment form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) sediment load per unit width q s of the LYR in response to the cutoff of sediment supply.The inset shows detailed results near the upstream end.

Figure 5 .
Figure 5.The 0.2-year results for the case of uniform sediment using the entrainment form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) sediment load per unit width q s of the LYR in response to the cutoff of sediment supply.Sediment fall velocity v s is arbitrarily multiplied by a factor of 0.05 while holding bed grain size constant in this run.The inset shows detailed results near the upstream end.

Figure 6 .
Figure 6.The 0.2-year results for the case of sediment mixtures using the flux form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) total sediment load q sT , (c) surface geometric mean grain size D sg , and (d) geometric mean grain size of sediment load of the LYR in response to the cutoff of sediment supply.The inset shows detailed results near the upstream end.

Figure 7 .
Figure 7.The 0.2-year results for the case of sediment mixtures using the entrainment form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) total sediment load q sT , (c) surface geometric mean grain size D sg , and (d) geometric mean grain size of sediment load of the LYR in response to the cutoff of sediment supply.The inset shows detailed results near the upstream end.

Figure 8 .
Figure 8.The 0.2-year results for the case of sediment mixtures using the entrainment form of the Exner equation: time variation of (a) bed elevation z b and water surface (WS), (b) total sediment load q sT , (c) surface geometric mean grain size D sg , and (d) geometric mean grain size of sediment load of the LYR in response to the cutoff of sediment supply.Sediment fall velocities v si are arbitrarily multiplied by a factor of 20 in this run while keeping the grain sizes invariant.The inset shows detailed results near the upstream end.

Figure 9 .
Figure 9. Relation between adaptation length L ad and grain size D. The values of flow discharge per unit width q w and recovery coefficient r 0 are the same as those in Sect.3.1.The relation of Dietrich (1982) is implemented for sediment fall velocity.
kinematic wave corresponding to F i in the entrainment form c Fi celerity of the kinematic wave corresponding to F i in the flux form D sediment grain size E dimensionless entrainment rate of sediment F i volumetric fraction of surface material in the ith size range f I i volumetric fraction of sediment in the ith size range exchanged across the surface-substrate interface g gravitational acceleration h

Table 1 .
Summary of computational conditions for numerical modeling of the LYR.

Table 4 .
Quantification of the difference between predictions of the flux form and the entrainment form in the case of uniform sediment.The maximum δ values in the calculational domain are presented for 5 years.

Table 5 .
Quantification of the difference between predictions of the flux form and the entrainment form in the case of sediment mixtures.The maximum δ values in the calculational domain are presented for 5 years.