sedFlow – a tool for simulating fractional bedload transport and longitudinal profile evolution in mountain streams

Especially in mountainous environments, the prediction of sediment dynamics is important for managing natural hazards, assessing in-stream habitats and understanding geomorphic evolution. We present the new modelling tool sedFlow for simulating fractional bedload transport dynamics in mountain streams. sedFlow is a one-dimensional model that aims to realistically reproduce the total transport volumes and overall morphodynamic changes resulting from sediment transport events such as major floods. The model is intended for temporal scales from the individual event (several hours to few days) up to longer-term evolution of stream channels (several years). The envisaged spatial scale covers complete catchments at a spatial discretisation of several tens of metres to a few hundreds of metres. sedFlow can deal with the effects of streambeds that slope uphill in a downstream direction and uses recently proposed and tested approaches for quantifying macro-roughness effects in steep channels. sedFlow offers different options for bedload transport equations, flow-resistance relationships and other elements which can be selected to fit the current application in a particular catchment. Local grain-size distributions are dynamically adjusted according to the transport dynamics of each grain-size fraction. sedFlow features fast calculations and straightforward preand postprocessing of simulation data. The high simulation speed allows for simulations of several years, which can be used, e.g., to assess the long-term impact of river engineering works or climate change effects. In combination with the straightforward preand postprocessing, the fast calculations facilitate efficient workflows for the simulation of individual flood events, because the modeller gets the immediate results as direct feedback to the selected parameter inputs. The model is provided together with its complete source code free of charge under the terms of the GNU General Public License (GPL) (www.wsl.ch/sedFlow). Examples of the application of sedFlow are given in a companion article by Heimann et al. (2015).


Introduction
Environmental models typically seek to predict the future state of a system, based on information about its current state and the mechanisms that regulate its evolution through time.In the case of sediment transport by flowing water in open channels, the temporal evolution of these variables is determined by a complex interaction of multiple processes includ-ing hydraulic water routing, sediment entrainment, erosion and deposition.In recent years many numerical models have been developed for simulating sediment transport in rivers.However, most of these models are intended for, and only applicable in, lowland rivers with gentle slopes.In mountain streams the effects of macro-roughness and shear stress partitioning have to be considered.Otherwise, sediment transport rates may be overestimated by several orders of magnitude F. U. M. Heimann et al.: The bedload transport model sedFlow (Rickenmann and Recking, 2011;Nitsche et al., 2011Nitsche et al., , 2012)).Based on these observations and contrasting the dominantly gradient-based definition used by, e.g., Wohl (2000), we define mountain streams as streams which are located within a mountainous region and in which the effects of macroroughness and shear stress partitioning play an important role in the sediment transport system.
Few sediment transport models have been specifically designed for mountain streams.Cui et al. (2006) developed the two Dam Removal Express Assessment Models (DREAM-1&2), based on the previous models of Cui and Parker (2005) and Cui and Wilcox (2008), to focus specifically on dam removal scenarios.The DREAM models therefore feature the simulation of (a) bank erosion during the downcutting of reservoir deposits, (b) transcritical flow conditions, (c) combined bedload and suspended load transport, (d) the details of gravel abrasion and (e) staged dam removal and partial dredging as options in the dam removal scenarios.Due to their specific focus, the wider applicability of the DREAM models is limited.Both the model of García-Martinez et al. (2006) and the model MIKE21C (DHI, 1999, with its modifications by Li and Millar, 2007) focus on a two-dimensional representation of hydraulic and sediment transport processes.Therefore, these models require more extensive input data and longer calculation times compared to one-dimensional model representations.As another example, the model of Papanicolaou et al. (2004) is intended for studying sediment transport under transcritical flow conditions by solving the unsteady form of the Saint-Venant equations, which results in long calculation times.
Other sediment transport models have been described by Mouri et al. (2011), Lopez and Falcon (1999) and Hoey and Ferguson (1994), all of which feature the one-dimensional simulation of fractional bedload transport using a simplified representation of the hydraulic processes.The model of Mouri et al. (2011) can represent a combination of debris flow, bedload and suspension load processes.In contrast, SEDROUT (Hoey and Ferguson, 1994) is designed to study the spatial and temporal evolution of local grain-size distributions.Therefore, it determines the composition of the sediment surface layer by a numeric iteration within each time step.In its latest version, SEDROUT has been also extended to deal with islands and other features of river bifurcation (Verhaar et al., 2008).However, neither the source code, the executable model binary nor a detailed description of the model implementation is available for any of the three models mentioned in this paragraph.
The model Tom Sed (formerly known as SEdiment TRansport model in Alpine Catchments (SETRAC)) was developed to study the influence of different shapes of channel cross sections on bedload transport in steep streams (Chiari et al., 2010;Chiari and Rickenmann, 2011).Therefore, the user can define cross sections with laterally varying bed elevations.The shape of a particular cross section stays the same during the complete simulation.Most published model applications used bedload transport calculations for a single grain size.In such a situation, all grain sizes and their spatial distribution are constant for the complete simulation.In this set-up, Tom Sed is slightly faster than real time in a typical application.A fractional transport approach with dynamic grain-size distributions is implemented in Tom Sed as well.However, it is rarely used due to the long calculation times.
The TOPographic Kinematic wave APproximation and Integration (Topkapi) model was originally developed as a rainfall-runoff model providing fast hydrologic simulations (Ciarapica and Todini, 2002).Later a sediment transport module was added, and this model version is called Topkapi ETH (Konz et al., 2011).The code is intended for the study of reach-scale sediment transport in the context of largescale hydrologic processes.Due to this scope that integrates different processes and scales, the model features a spatial as well as temporal sub-gridding approach.The hydrologic processes are simulated on a coarse two-dimensional grid with time steps that are an integer multiple of the time steps for the hydraulic and sediment transport processes.The latter two processes are simulated in a one-dimensional channel at a finer spatial resolution.This channel receives water from the hydrologic two-dimensional grid, but the morphodynamic changes due to bedload transport have no influence on the topography used for the hydrologic calculations.The channel cross section is represented by a rectangle and bedload transport is based on a single grain-size approach, in which local grain-size distributions do not change over time.In typical applications, a flood event of several days can be simulated within a few minutes of calculation time.
Currently, no model is available that combines short calculation times with easy use and up-to-date sediment transport equations for mountain catchments.The new model sedFlow (Figs. 1, A1, A2, A3) presented in this contribution has been developed to provide an efficient tool for the simulation of bedload transport in mountain streams.By efficient tool we mean a model that combines straightforward pre-and postprocessing of simulation data with fast calculation speeds.sedFlow is intended for temporal scales from the individual event (several hours to a few days) to longer-term channel evolution (several years).The envisaged spatial scale covers complete catchments at a spatial discretisation of several tens to a few hundreds of metres.The following elements were important for the development of sedFlow: 1. a sediment transport model provided together with its complete source code, open source and free of charge; 2. implementation of recently proposed and tested approaches for calculating bedload transport in steep channels accounting for macro-roughness effects; 3. individual calculations for several grain diameter fractions (fractional transport) resulting in more robust simulations as compared to single-grain approaches; 4. consideration of ponding effects of adverse slopes (uphill slopes in the downstream direction), e.g.due to sudden sediment deposition by debris flow inputs; 5. fast calculations for modelling entire catchments, and for automated calculation of multiple scenarios exploring a range of parameter space; 6. an object-oriented code design that facilitates flexibility in model development; 7. flexibility in model application featuring different options (Figs.A1, A2), e.g. for the bedload transport equations, which can be selected to fit the current application, as well as straightforward pre-and postprocessing of simulation data.
The model sedFlow thus fills a gap in the range of existing sediment transport models for mountain streams (Table 1) and the goals outlined above have led to the implementation described in the following sections.This implementation represents the current state of the model, and may be easily extended and adjusted in the future.
Examples of sedFlow application are given in a companion article by Heimann et al. (2015).

Hydraulic calculation
Hydraulic equations describe the temporal evolution of the three-dimensional flow field of the water continuum.A formalised description of the dynamics has been provided by Navier and Stokes (given in the form for incompressible flow).
where ρ is fluid density, v is flow velocity, t is time, p is pressure and µ is dynamic viscosity.f summarises other influencing body forces.If large-scale backwater effects are not present, as is often the case in steep channels of mountain streams, the energy slope can be approximated by the bed slope (i.e.assumption of kinematic wave propagation) and the complete Navier-Stokes equation can be reduced to the following simplified, cross-sectional averaged conservation of volume (e.g.Chow et al., 1988).
where Q is discharge, x is distance in flow direction, A is wetted cross-sectional area and Q lat is lateral water influx.

Flow routing
Within sedFlow a channel network joined by confluences can be simulated.At the upstream ends of the main channel and each of the user-defined tributaries, a discharge time series is input and has to be routed through the channel system.For the following discussion of hydraulic routing schemes we will differentiate between three cases: first, in ponding -i.e.water collecting behind sediment obstructions -the friction slope S f is approximately zero (S f ≈ 0).Second, in situations with parallel slopes, the friction slope approximately equals the channel bed slope S b (S f ≈ S b ), which is commonly true for steep S b .Third, the situations of moderate backwater effects cover all cases between the extremes of ponding on the one hand and situations with parallel slopes on the other.Especially in models not focussing on the details of the hydraulic routing, the kinematic wave approach (assuming the situation of parallel slopes) can be implemented using a temporally explicit Eulerian forward approach (van de Wiel et al., 2007;Chiari et al., 2010).Such an approach can be used in sedFlow as well (see Sect.A1).However, Eulerian forward approaches must assume that all parameters within one time step can be sufficiently approximated by their values at the beginning of the time step.To ensure this approximation, Eulerian forward approaches require very small time steps, especially for fast processes.This can be problematic when a relatively fast process, such as the routing of water, is combined with a relatively slow process, such as bedload transport including bed level adjustments.The water routing requires small time steps and thus calculation times that may be orders of magnitude too long and slow from the perspective of bedload dynamics.Therefore, as an alternative, an implicit discharge routing is implemented in sedFlow.The implicit routing is unconditionally stable, and thus has no requirements concerning the length of time steps.In sedFlow the approach of Liu and Todini (2002) is used, which omits time-consuming iterations and finds the solution for the kinematic wave analytically via a Taylor series approximation.However, the approach depends on a power-law representation of discharge as a function of water volume in a reach.This means that it can only be applied to the specific crosssectional shapes of infinitely deep rectangular or V-shaped channels in combination with a power-law flow-resistance equation.
The kinematic wave assumption of parallel slopes is usually valid for steep channel gradients, which are typical of mountain catchments.Nevertheless, the kinematic wave assumption can be problematic, especially in mountain streams, when tributaries deliver large amounts of sediment to the main channel within a relatively short time, e.g. during debris flow events.This may result in adverse channel slopes (uphill slopes in the downstream direction) and backwaters in the main channel, violating the assumptions of a kinematic wave.If configured for kinematic wave routing, sedFlow will abort simulations whenever adverse channel slopes occur.
When it is necessary to deal with adverse channel slopes, one has to drop the kinematic wave approximation and use a backwater calculation instead.Unfortunately, the backwater calculation is numerically intensive.Therefore, within sedFlow, a pragmatic approach can be selected to deal with adverse channel slopes: discharge is assumed to be uniform and thus equal along the entire channel for a given time step only increasing at confluences.This assumption of uniform discharge is reasonable because the temporal scale for water routing is orders of magnitude smaller than the temporal scale for morphodynamic adjustments.In the case of positive slopes, flow depth and velocity are commonly calculated using the bed slope as proxy for the friction slope (thus assuming the bed slope and water surface are parallel).However, the flow resistance may be adjusted such that a maximum Froude number is not exceeded.In cases of adverse channel slopes, the formation of ponding is simulated.That is, flow depth and velocity are selected to ensure a minimum gradient of hydraulic head, which is positive but close to zero, and thus ensures numeric stability, and approximately corresponds to the hydraulic gradient of ponding water.For bedload transport calculations the gradient of the hydraulic head is used, which by definition can only have positive slopes.Thus, the energy slope for bedload transport estimation is not the result of a backwater calculation, but it is the gradient between individual hydraulic head values, which under non-ponding conditions have been calculated independently from each other using the local bed slope as a proxy for friction slope.This approach is based on the assumption that the simulated system only consists of the two extreme cases of ponding on the one hand and parallel slopes on the other.At a spatial discretisation of several tens of metres, the assumption of the two extreme cases is valid for many mountain streams and it allows efficient simulation of ponding, by omitting numerically intensive backwater calculations.However, it has to be noted that this approach will produce large errors when intermediate cases of moderate backwater effects are part of the simulated system.In such systems, the first approach, which uses bed slope both as friction slope for Qualitative comparison of uniform discharge (left) and kinematic wave (right) flow routing approaches.The top row shows a hypothetical discharge time series at the upstream boundary, which can be used in both approaches.The point in time for the following rows is indicated by the dashed vertical line.In the uniform discharge approach, the current discharge value of the time series defines the discharge for all reaches of the simulated system at the current point in time (second row left).In contrast, for the kinematic wave approach, the temporal variability of discharge is reflected in a spatial variability as well (second row right).Therefore, in the uniform discharge approach, the spatial variation of flow depth (third row), as well as water surface (blue curve), is mainly a function of roughness and slope, which is determined by the river bed (black curve).For the kinematic wave approach, flow depth may also vary due to the spatial variation of discharge (third row right).In cases of uphill channel slopes in the downstream direction, the uniform discharge approach will reproduce the effects of ponding (fourth row left), while the kinematic wave approach cannot deal with such situations (fourth row right).
the hydraulic calculations and as energy slope for the sediment transport calculations, will produce better estimates of the transported sediment volumes, but it cannot accommodate adverse channel gradients.Heimann et al. (2015) have demonstrated that, despite their simplicity, the implemented hydraulic concepts (Fig. 2) appear to be sufficient for a realistic integrated representation of bedload transport processes.In that study, for the different hydraulic routing schemes, results have been obtained, which  are close to the observed morphodynamic changes and very similar among each other.

Flow resistance
The interaction of flowing water with the river bed and banks determines the relation between the average downstream velocity and the wetted cross section.This interaction is summarised as flow resistance, which can be described by the following physically based relation: where f is the Darcy-Weisbach friction factor, v is crosssectional mean flow velocity, g is gravitational acceleration and r h is hydraulic radius.The flow routing method of Liu and Todini (2002) requires a power-law relation between discharge and water volume within a reach.Therefore, the following flow-resistance law can be used in sedFlow.
Here j 1 , k and l are empirical constants and D x is the xth percentile diameter of the local grain-size distribution.Unless otherwise stated, all grain sizes refer to the surface layer.By selecting l = 1 6 , this formula represents a classic grain-size dependent Gauckler-Manning-Strickler relation.For the other variables, the values j 1 = 6.5, k = 1 and x = 84 (Eq.4b) have been found to reproduce observational data well for deeper flows with r h D 84 larger than about 7-10 (Rickenmann and Recking, 2011).If another flow routing is used, one can select the variable power-equation flow-resistance approach provided by Ferguson (2007) with the parameter values proposed by Rickenmann and Recking (2011), which was recommended also for applications in steep channels, including shallow flows with small relative flow depths r h D 84 : , where j 1 = 6.5 and j 2 = 2.5. (5) In general, flow resistance describes the effects of drag forces exerted on the bed and its structures.A part of this drag, namely skin drag, is responsible for the transport of sediment grains (e.g.Morvan et al., 2008).The drag created by larger-scale surface geometries, such as bed forms and channel shape features (e.g.bends and irregular channel width), may be summarised as macro-roughness, which reduces the energy available for the transport of sediment.
Here S red is the reduced slope that accounts for macroroughness effects, S is channel or hydraulic energy slope, f 0 is base-level flow resistance according to Eq. ( 4b), f tot is total flow resistance according to Eq. ( 5) and e is an exponent ranging from 1 to 2, with a typical value of e = 1.5.Within sedFlow, one can select to use S red based on Eq. ( 6) to account for macro-roughness.6), by comparison with bedload transport observations in steep mountain streams (Nitsche et al., 2011).The equations of Wilcock and Crowe (2003) were derived from fractional bedload transport data.The equation of Recking (2010) was developed for the estimation of total bedload transport rates.

2.2
In the same way as the equations of Meyer-Peter and Müller (1948), Fernandez Luque and van Beek (1976) and Soulsby and Damgaard (2005), the equation of Rickenmann (2001) (Sect.A4; especially its simplified version in Eq.A20) is a good example of the following generic type of bedload estimation methods: where b = q b √ (s−1)gD 3 is dimensionless bedload flux, θ = τ (s−1)ρgD is dimensionless bed shear stress, θ c is the dimensionless bed shear stress threshold for the initiation of motion, q b is bedload flux per unit flow width, D is grain diameter, a, b and d are empirical constants, τ = (ρ • g • r h • S) is bed shear stress and s = ρ s ρ is the density ratio of solids ρ s and fluids ρ.To account for macro-roughness, S can be replaced by S red in the calculation of the bed shear stress τ .In the case of the Rickenmann (2001) equation, a is the product of an empirical constant and the Froude number F r.In equations like Eq. ( 7), bedload transport is mainly a power law of the dimensionless bed shear stress that exceeds some threshold for the initiation of bedload motion.This threshold is known as the Shields criterion (Shields, 1936) with typical values ranging from 0.03 to 0.05.In natural channels, geometric complexity and thus energy losses increase at steep bed slopes S b .Therefore, Lamb et al. (2008) suggested the following empirical relation to account for increasing θ c values with increasing S b : The main objective of Lamb et al. ( 2008) is a theoretical explanation for the observed increase of θ c with increasing bed slopes S b .However, for very small values of S b , the power law of Eq. ( 8) predicts values of θ c that are approaching zero.This contrasts with the results of other studies (e.g.Recking, 2009) that observed roughly constant values of θ c larger than zero for small slopes.sedFlow can be configured to use either a constant threshold or a slope-dependent threshold according to Eq. ( 8) combined with a minimum value θ cMin .
The estimated bedload flux can be corrected for gravel abrasion according to the classic equation of Sternberg (1875), in which q b abr is bedload flux per unit flow width corrected for abrasion, λ is an empirical abrasion coefficient and X is the travel distance of the grains.Here the material loss due to abrasion is regarded as suspension throughput load: If grain-size fractions are treated individually, the calculated bedload capacity b needs to be normalised with F i , the relative volumetric portion of bed surface material of a grainsize fraction i, compared to the total surface material with D > 2 mm (e.g.Parker, 1990).Here F i can be interpreted as the availability of a certain grain-size fraction in the bed; for an example, see Eq. (A23) in Sect.A6 compared to Eq. (A20) in Sect.A4.Further details also have to be accounted for, such as the varying exposure of different grain-size fractions, grain-size-dependent grain-grain interactions, etc.This is commonly done using some sort of hiding function.Hiding functions not only focus mainly on grain exposure but also integrate all kinds of grain-size-dependent effects which are not covered by the capacity estimation methods.Within sed-Flow a relatively simple power-law hiding function can be used (Parker, 2008): as well as the hiding function of Wilcock and Crowe (2003): , where m wc = 0.67 Here θ ci is the θ c for the ith grain-size fraction, D i is the mean grain diameter for ith grain-size fraction, m is an empirical hiding exponent, D m is the geometric mean diameter of the local grain-size distribution and m wc is the hiding exponent according to Wilcock and Crowe (2003).The empirical exponent m ranges from 0 to −1, where m = −1 corresponds to the so-called "equal mobility" case in which all grains start moving at the same bed shear stress τ , and m = 0 corresponds to no influence by hiding at all.For x = 50, the values for m, which have been derived from various field observations, typically vary within a range from −0.60 to −1.00 (Recking, 2009), and unfortunately there are only few data points for D i > D 50 (Bathurst, 2013;Bunte et al., 2013).For consistency, the following θ ci,r is used in bedload transport calculations.
Within sedFlow two alternatives are implemented for the calculation of the correction factor γ : and In Eq. ( 13a), θ ci,r varies with discharge, as it depends on S red , which in turn is a function of the hydraulic radius r h .In Eq. ( 13b), suggested by Nitsche et al. (2011), θ ci,r is independent of discharge.The value of S c is calculated using Eq. ( 6), with the value of r h replaced by the critical hydraulic radius r h, c : Good arguments can be found for both approaches (Eqs.13a and 13b).Due to the lack of suitable data, it is unclear which approach is more plausible.

Evolution of channel bed elevation and slope
The temporal evolution of the longitudinal profile is simulated in sedFlow based on a finite-difference version of the general Exner equation (e.g.Parker, 2008).
Here η pore is pore volume fraction, z is elevation of the channel bed and q b lat is lateral bedload influx per unit flow width.Eq. ( 15) allows for the calculation of the new channel slope z x after each time step.Heretofore, infinitely deep rectangles are used within sedFlow as the shape of the crosssectional profiles, with the complete width defined as active width (i.e.sediment transport takes place over the complete width).
All three elements (the cross-sectional channel geometry, its alteration due to morphodynamics and the determination of the active width) are implemented as abstract classes.Thus, the presented realisations just represent the current state of the code and any programmer can easily extend the code to deal with more complex cross-sectional geometries.However, the implicit flow routing by Liu and Todini (2002), with its advantages in terms of simulation efficiency, requires www.earth-surf-dynam.net/3/15/2015/Earth Surf.Dynam., 3, 15-34, 2015  16) for the simulations and observations displayed in Fig. 3: obs denotes the reference data derived from observations, orig denotes the simulation results using the detailed channel geometry, and rect(Q) and rect(w) denote simulation results using a rectangular substitute channel based on averaged representative discharges (Q) and widths (w), respectively.The numerical values in the table are the rRMSDs between the two data sets determined by the row and column.The labels for the rows and columns (bold type face) are given in the diagonal of each panel.For example, at the Chirel the simulations using a detailed channel geometry and using a substitute rectangle based on averaged widths differ from each other by a rRMSD value of 0.28.infinitely deep rectangular or V-shaped channels (together with a simple power-equation flow-resistance law such as Eq. 4a).
Additionally, Stephan (2012) has studied the impact of the rectangular-shaped approximation and found that, at least during major flood events, it is negligible compared to the other uncertainties.Stephan (2012) recalculated bedload transport for the August 2005 transport event in the catchments of the Chiene, Chirel and Schwarze Lütschine.For details on the catchments and event characteristics see Chiari and Rickenmann (2011).For the simulations, Stephan (2012) used the one-dimensional model Tom Sed (Chiari et al., 2010), which allows for the definition of cross sections with laterally varying bed elevations.The simulations were repeated with a detailed channel geometry as presented in Chiari and Rickenmann (2011) and with two different rectangular substitute channels.The widths of the rectangular substitute channels w were determined based on a discharge Q rep that is representative of the simulation period.The channel width was selected to produce the same wetted cross-sectional area and hydraulic radius for the representative discharge as was simulated for the detailed channel geometry.In one approach the threshold discharge for the initiation of bedload motion Q c , as well as the maximum discharge of the simulation period Q max , were averaged to find Q rep , for which the representative channel width was determined.In the second approach, one width w was determined for both Q c and Q max and then the two widths were averaged to find the representative channel width.The detailed channel geometry produced results that were broadly similar to the rectangular substitutes when compared with the field observations on bedload transport (Fig. 3).To quantify the deviations between the different simulations and between the simulations and the observations, Table 2 summarises the relative root mean square deviations (rRMSD) which have been calculated according to the following equation: Here χ and ψ are two arbitrary data sets of length n, with χ and ψ being their average values, and i is a running index.
As can be seen from Table 2, the simulations on average differ from each other by a rRMSD value of only 0.23, while the simulations differ from the observations on average by a rRMSD value of 0.70.However, it has to be noted that these results are for a flood event with high discharge values.For low discharges close to the initiation of bedload motion, a detailed representation of channel geometry may play a more important role than the one suggested by Fig. 3 and Table 2.
For further details see Stephan (2012).
Finally, the introduction of more complex cross-sectional shapes raises the question of how these shapes are influenced and altered through morphodynamics.As far as we know, no generally accepted concepts are available for this problem.

Grain-size distribution changes
In sedFlow the alluvial substrate of the river is represented by a stack of horizontal layers with homogeneous grain-size characteristics.The topmost layer of the bed interacts with the flow and is typically called the active surface layer.The grain-size distribution of the active surface layer is used for the determination of the flow resistance, hiding processes and bedload transport capacity (Fig. 1).All deposited material is added to this layer; all eroded material is taken from it.The thickness of this layer determines the inertia of its evolving grain-size distribution.Therefore, the active surface layer and especially its thickness play an important role in the numeric representation of bedload transport systems (e.g.Belleudy and Sogreah, 2000).When the alluvium thickness is smaller than the expected usual active layer thickness, sedFlow makes use of the shape properties of bedrock to determine flow resistance and hiding.The thickness of the active surface layer may be set constant or dynamic as a multiple of some grain-size percentile.Three different approaches are available within sedFlow for the interaction between the active surface layer and the underlying subsurface alluvium.
The first method (Fig. A4) has been adapted from the one described by van de Wiel et al. (2007).Lower and upper thresholds are defined for the thickness of the active surface Earth Surf.Dynam., 3, 15-34, 2015 www.earth-surf-dynam.net/3/15/2015/layer.Whenever these thresholds are exceeded, sediment increments are incorporated from, or released to, the subsurface alluvium underneath until the active surface layer thickness again takes a value within the given thresholds.The sediment increments are stored as bed strata underneath the active surface layer.In this way the simulated river bed is able to remember its history.In contrast to the procedure of van de Wiel et al. ( 2007) the thickness of the sediment increments can be defined independently of the thickness of the active surface layer.Trivially, the thickness of the increments defines the minimum distance between the thresholds for the active surface layer thickness. 1The smaller this distance between the thresholds is, the more intense the interaction is between the active surface layer and the underlying subsurface alluvium.
The second approach (Fig. A5), which has been applied in various models (e.g.Hunziker, 1995), can be described as an extreme case of the first one, in which the two thresholds collapse to a single target thickness for the active surface layer.In this case, any addition or removal of material to or from the active surface layer is instantaneously balanced against the underlying subsurface alluvium.In this case of maximum interaction between active surface layer and subsurface alluvium, the subsurface alluvium is represented by one homogenised volume without any internal structure.The target thickness is usually determined at the start of a simulation 1 Here h s is the active surface layer thickness, within a time step, after erosion or deposition but before the interaction between the active surface layer and the subsurface alluvium; accordingly, h s,post is the active surface layer thickness after the layer interaction, h s is the thickness of sediment increments, |y| is the number of sediment increments that are incorporated from or released to the subsurface alluvium, thresh hs,low and thresh hs,high are the thresholds for the thickness of the active surface layer with thresh hs,low < thresh hs,high and thresh hs,prox is an alias for the threshold that is closer to h s .It is the objective of the layer interaction to reach a state in which h s,post is as close as possible to the middle between the thresholds and thresh hs,low ≤ h s,post ≤ thresh hs,high with h s,post = h s + (y • h s ) and y ∈ Z. (F1) However, the state of Eq. ( F1) cannot be reached if mod h s − thresh hs,prox h s < h s − thresh hs,high − thresh hs,low . (F2) This might cause instability.In simple terms, if h s was larger than the distance between the thresholds, situations might occur in which the addition or subtraction of another h s will cause h s,post to jump over both thresholds, such that it can never reach a value between the thresholds.Therefore, h s is defined as the minimum distance between the thresholds.With this minimum distance, the right-hand side of Eq. (F2) cannot become larger than zero, while the left-hand side of Eq. (F2) cannot become smaller than zero anyway.
and then kept constant.Alternatively, it can be dynamically adjusted based on a Eulerian forward approach, in which the thickness is updated at the end of each time step.The third approach (Fig. A6) is a variation of the second one.When sediment is eroded, only the volume of the active surface layer is instantaneously replaced from the subsurface, while the grain-size distribution stays the same.The sediment volume that is transported from the subsurface alluvium to the active surface layer shares the grain-size distribution of the subsurface alluvium only if a condition for the break-up of an armouring layer is fulfilled: Here θ 50 is a representative dimensionless shear stress θ for the median diameter D 50 of the active surface layer and θ c,s is a representative θ c for the active surface layer.To avoid artefacts due to a hard threshold, some fraction i sub of the sediment transported from the subsurface alluvium to the active surface layer has the grain-size distribution of the subsurface alluvium already before the break-up condition is fulfilled.The rest (1 − i sub ) has the grain-size distribution of the active surface layer: Here i sub is the relative grain-size influence from the subsurface alluvium and θ c,sub is a representative θ c for the subsurface alluvium.The value of θ c,sub can be estimated, e.g. according to Eq. ( 8), while the value of θ c,s can be estimated using the following relation of Jäggi (1992): Here D mArith s and D mArith sub are the arithmetic mean diameters of the grain-size distribution of the active surface layer (s) and of the subsurface alluvium (sub).
For non-fractional studies, the active surface layer concept can be turned off.In that case, the complete alluvium is represented by a single homogeneous layer, which directly interacts with the flow.

Comparison of sedFlow implementations with similar models
In the following subsections, various details of the implementation of sedFlow are discussed and compared to implementations in similar models.Differences between the individual models are explained in the context of their differing objectives.
F. U. M. Heimann et al.: The bedload transport model sedFlow

Fractional transport and grain-size distributions
In streams, local grain-size distributions and bedload transport influence each other in various ways.For example, flow resistance is often determined as a function of a characteristic bed surface grain diameter such as in Eqs. ( 4b) and ( 5) that use D 84 .Therefore, local grain-size distributions have a high influence on the local hydraulic radius r h , which in turn determines bed shear stress and thus bedload transport capacity.In addition, the shear stress partitioning in Eq. ( 6) is based on D 84 as well.Finally, the stream can only mobilise and transport the grains which are available in the bed and, as transport capacity is inversely related to the transported grain diameter, fine-grained bed material will yield high transport rates.These mechanisms let the bed surface grain-size distribution influence the transport capacity, which in turn alters the local grain-size distributions.At high transport capacity, the stream will erode and deplete particles with transportable diameters and will leave only coarse grains in the bed.At low transport capacity, the stream cannot transport and will, thus, even deposit grains with small diameters.This will accumulate fine-grained material in the bed.Therefore, the spatial pattern of bed surface grain-size distributions at the end of a simulation can be interpreted as a function of bed slope (coarse grains in steep sections), channel width (coarse grains in narrow sections) and channel network effects (coarse grains at confluences with steep tributaries).
In order to study dynamically evolving grain-size distributions and their effects on hydraulics and bedload transport, together with the effects of an evolving channel slope, sedFlow is optimised for the simulation of fractional bedload transport.In this context, different concepts are provided for the interaction between the active surface layer and the underlying alluvium.As grain-size distributions dynamically adjust to be consistent with the local circumstances (channel width, slope, etc.), this numeric concept might partially compensate the uncertainty in local grain-size distribution data.For a more detailed discussion of this topic see Heimann et al. (2015).To the authors' knowledge there are no in-depth studies assessing the influence of different representations of active surface layer dynamics on simulated bedload volumes.As sedFlow contains three different formulations of active surface layer dynamics in the same modelling tool, it provides the base for a future study on the effects of different active surface layer algorithms.
Within Topkapi ETH (Konz et al., 2011) fractional transport is not implemented, and within Tom Sed (Chiari et al., 2010) it is rarely used due to long calculation times.These models have a different application objective, for which dynamic grain-size distributions are of lesser relevance.
The threshold-based layer interaction (Fig. A4) implemented in sedFlow is similar to the approach of Lopez and Falcon (1999) for the evolution of the local grain-size distribution.However, Lopez and Falcon (1999) only introduced a lower threshold for the thickness of the active surface layer.This means that the subsurface alluvium remains constant even in cases of massive aggradation and that the active surface layer may reach unrealistically high thickness values, especially in cases of intense aggradation.
In SEDROUT (Hoey and Ferguson, 1994), the surface layer grain-size distribution is determined as a function of the spatial derivative of fractional transport rates and the thickness of the surface layer, which in turn is a function of the surface layer grain-size distribution.This set of equations is solved by numeric iteration within each time step.In sedFlow, this numerically extensive procedure is replaced by a constant surface layer thickness or an Eulerian forward approach, in which the layer thickness is updated at the end of each time step.These more pragmatic approaches have been selected because fast simulations are one of the main aims of sedFlow.

Adverse channel slopes
Within sedFlow adverse channel slopes (uphill slopes in the downstream direction) and their effects in terms of ponding can be considered using uniform discharge hydraulics.This approach is based on the assumption that the simulated system only consists of two extreme cases: ponding on the one hand, and parallel slopes on the other.This assumption is valid in many mountain streams and it allows for fast simulations, which have been another main objective for the development of sedFlow.However, in the intermediate case of moderate backwater effects, it may produce large errors.The implemented approach corresponds to the situation of a confined channel, in which the sudden deposition of large volumes of sediment, e.g. by debris flow inputs, may produce ponding.
Within Topkapi ETH any large volumes of deposited material, which would produce adverse channel slopes, are instantaneously distributed to downstream river reaches until all slopes are positive.This algorithm would correspond to instantaneous landslides within the channel or to debris flows with short travel distances, but such phenomena are typically not observed in mountain streams.However, instantaneous lateral sediment input is not the main focus of the Topkapi ETH model.In the context of its intended applications, the described algorithm of Topkapi ETH is an appropriate pragmatic approach to conserve mass and ensure positive bed slopes, which are used as energy slopes in the implemented kinematic wave approach.
Within Tom Sed any deposited material that would produce adverse channel slopes is fed to a virtual sediment storage, which does not contribute to elevation changes in the main channel.As long as there is sediment in this virtual storage, any erosion or deposition is applied to this storage keeping the main channel untouched.That means that the elevation of the main channel is frozen as long as there is material in the virtual storage.In some way this algorithm corresponds to a lateral displacement of the river channel due to large volumes of deposited material that are stored next to the new channel.However, in mountain rivers the amount of material that is fed from the deposits to the main channel depends on the stability of the deposit slopes.Therefore, the new channel may well lower its elevation due to erosion, even if there is still some material in storage next to the channel.However, the situation of lateral channel displacement is close to the limits of a one-dimensional simulation and the described algorithm of Tom Sed ensures positive bed slopes (used as proxy for the energy slopes) in a way which corresponds to some extent to a natural process, even though it violates conservation of mass within the main channel.

Simulation speed
sedFlow has been designed for optimal computation speeds.Besides the selection of the coding language C++, for which there are powerful compilers available, we have implemented a spatially uniform discharge (within each segment of the channel network), as well as an implicit kinematic wave flow routing approach, both aimed at providing a modelling tool for fast simulations.Both hydraulic approaches allow for coarse temporal discretisations and the implemented algorithm of Liu and Todini (2002) omits computationally demanding iterations.The ideal temporal discretisation (as fine as necessary and as coarse as possible) can be obtained from the Courant-Friedrichs-Lewy (CFL) criterion based on the speed of bedload multiplied by a user-defined safety factor.As a result, several years of bedload transport and resulting slope and grain-size distribution adjustment can be simulated with sedFlow within only few hours of calculation time on a regular 2.8 GHz central processing unit (CPU) core.
In Topkapi ETH the implicit kinematic wave flow routing is also implemented using the algorithms of Liu and Todini (2002).However, the length of the time steps used for the bedload transport simulations may differ from the ideal length, as it is defined as an integer multiple of the time steps used for the hydrologic simulations.This implementation is due to Topkapi ETH's aim to simulate different processes at different scales.
Within Tom Sed , explicit kinematic wave flow routing is implemented.Thus, time steps need to be determined using the CFL criterion based on water flow velocity.Therefore, time steps are shorter than in sedFlow or Topkapi ETH and slow down the simulations considerably.The choice of the explicit water flow routing is due to Tom Sed 's aim to simulate the effects of the shape of channel cross sections.The implicit flow routing based on the algorithms of Liu and Todini (2002) requires simple rectangular or V-shaped channels and would therefore prevent any detailed study of channel geometry.

Flexibility
To ensure flexibility of application, we selected regular spreadsheets as the file format for the data input to sedFlow.Thus, preprocessing can be done with common software applications, which are familiar to most users, so that data from any study catchment can be quickly and easily prepared for a sedFlow simulation.This contrasts with Tom Sed , which uses the extensible markup language (XML) file format for data input, as well as with Topkapi ETH, which partially requires MATLAB preprocessing.In sedFlow different equation sets can be selected and combined by the user for the main process representations (flow routing, flow resistance, initiation of bedload motion, transport capacity, etc.).The number, content and format of the output files can be defined by the user as well, in order to get the best solution for the respective study objectives.
To ensure flexibility of model development we selected an object-oriented design for the internal structure of the sedFlow code.In such a code, succeeding programmers just create new realisations for predefined code interfaces without revising the model core.It is not necessary to know (virtually) anything about the model itself.The only piece of code that the programmer will have to read is the specification of the relevant code interface, which is typically not longer than two printout pages.

Advantages and limits of the sedFlow modelling approach
The limits of the sedFlow modelling approach are defined by the implemented process representations.Hydraulics are considered in a cross-sectionally averaged way.The effects of turbulence and eddies, such as the formation of dunes and pool-riffle sequences, are not considered explicitly.Therefore, sedFlow should not be applied at spatial discretisations small enough for these processes to become relevant.Within sedFlow, discharge is assumed to be generally subcritical.It should not be applied in systems where Froude numbers larger than one occur at the scale of spatial discretisation.
As sedFlow generally assumes kinematic wave propagation, bed slopes S b should be steep enough for friction slopes S f to approximately equal bed slopes (S f ≈ S b ).Slight violations of this assumption will not cause major errors.However, if one uses the flow routing approach that allows for adverse channel slopes (uphill slopes in the downstream direction), any occurrence of moderate backwater effects (S f ≈ 0 ∧ S f ≈ S b ) will produce large errors (cf.Sect.2.1.1).
sedFlow considers sediment transport only in terms of bedload.Therefore, it should not be applied if suspension load and wash load play a major role in the studied system.
However, if the limits of applicability are respected, sedFlow provides several benefits for the user interested in the simulation of bedload transport in steep mountain streams.The model is optimised for the calculation of frac-F.U. M. Heimann et al.: The bedload transport model sedFlow tional bedload transport resulting in more robust simulations as compared to single-grain approaches.Based on an assumption that is valid for many mountain streams, sedFlow provides for the efficient simulation of adverse slopes and ponding.The selected coding language and two implemented hydraulics representations are all aimed at providing a modelling tool for fast simulations.Finally, with the user-defined structure of the model outputs, with the straightforward preprocessing of the inputs, and with the provision of different selectable options, e.g., for the flow resistance and bedload transport equations, sedFlow offers the flexibility to be easily adjusted to fit a specific application in a particular catchment.

Conclusions
The new model sedFlow complements the range of existing tools for the simulation of bedload transport in steep mountain streams.It is an appropriate tool if (i) grain-size distributions need to be dynamically adjusted in the course of a simulation, (ii) the effects of ponding, e.g., due to debris flow inputs might play a role in the study catchment or (iii) one simply needs a fast simulation of bedload transport in a mountain stream with quick and easy pre-and postprocessing.Examples of the application of sedFlow are presented in a companion article by Heimann et al. (2015).
The current version of the sedFlow code and model can be downloaded under the terms of the GNU General Public License (GPL) at the following web page: www.wsl.ch/sedFlow. .Qualitative sketch of threshold-based interaction between the active surface layer and subsurface alluvium.The water level is displayed in blue, the active surface layer with its variable thickness is light grey and bedrock is black.The subsurface alluvium consists of several strata layers with user-defined constant thickness displayed in reddish colours and one base layer with variable thickness displayed in dark grey.The horizontal dashed lines indicate the thresholds for the thickness of the active surface layer.In case of aggradation (a and c), the material input from upstream, which is displayed in green, is added to the active surface layer, which is instantaneously homogenised.Any homogenisation process is displayed by diagonal stripes.If the thickness of the active surface layer does not exceed its thresholds after aggradation or erosion (a and b), the thresholds of the active surface layer and the layers of the subsurface alluvium remain constant.If the active surface layer thickness exceeds its upper threshold after aggradation (c), the complete system of layers and thresholds is shifted upwards so that the upper boundary of the active surface layer is in the middle between its thresholds.The upper strata layers are populated with material from the homogenised active surface layer, while the material of the lower strata layers is added to the base layer, which is instantaneously homogenised.If the active surface layer thickness exceeds its lower threshold after erosion (d), the complete system of layers and thresholds is shifted downwards so that the upper boundary of the active surface layer is in the middle between its thresholds.The lower strata layers are populated with material from the base layer, while the material of the upper strata layers is added to the active surface layer, which is instantaneously homogenised. .Qualitative sketch of continuous interaction between the active surface layer and subsurface alluvium.The water level is displayed in blue, the active surface layer with its constant thickness is light grey, the subsurface alluvium consisting of one layer with variable thickness is dark grey and bedrock is black.The dashed bracket indicates the position of the active surface layer after aggradation or erosion.In case of aggradation (a), the material input from upstream, which is displayed in green, is added to the active surface layer, which is instantaneously homogenised.Any material of the homogenised active surface layer, which exceeds the constant thickness, is transferred to the subsurface alluvium, which is instantaneously homogenised as well.Any homogenisation process is displayed by diagonal stripes.In case of erosion (b), the sediment deficit of the active surface layer is replaced by material from the subsurface alluvium and the active surface layer is instantaneously homogenised.

Figure A6
. Qualitative sketch of shear-stress-based interaction between the active surface layer and subsurface alluvium.For an explanation of the symbols see the caption of Fig. A5.In this approach, the aggradation case is treated identically to the continuous update approach displayed in Fig. A5a.For erosion, three cases are differentiated, with the representative dimensionless shear stress θ 50 increasing from case (a) to case (c).If θ 50 equals or exceeds the threshold for the active surface layer θ c,s (c), the layers interact in the same way as in the continuous update approach displayed in Fig. A5b.If θ 50 does not exceed the threshold for the subsurface alluvium θ c,sub (a), the volume deficit of the active surface layer is replaced by material from the subsurface alluvium, but the grain-size distribution of the active surface layer remains constant.For the intermediate case with θ 50 greater than θ c,sub and smaller than θ c,s (b), the influence of the grain-size distribution of the subsurface alluvium on the grain-size distribution of the active surface layer is interpolated linearly between cases (a) and (c).bedload flux per unit flow width q b abr q b corrected for gravel abrasion q b lat lateral bedload influx per unit flow width q c threshold q for initiation of bedload motion

Figure 1 .
Figure 1.Simplified overview over the main process interactions within the sedFlow model.

Figure 3 .
Figure 3.Comparison of the effects of different channel representations on accumulated bedload transport estimates simulated with the Tom Sed model in three Swiss mountain rivers.SeeTable 2, text, and Stephan (2012) for details.

Figure A1 .
Figure A1.Simplified summary diagram of the bedload transport part of sedFlow, showing a selection of the most important general types of functions (bold type face) and their selectable realisations (regular type face).For the transport equations of Rickenmann (2001), versions for both fractional and uniform grain sizes are implemented in sedFlow.Among the implemented equations, only the θ -based ones of Rickenmann (2001) use θ c to determine the initiation of bedload motion.The selection of a hiding function and active surface layer dynamics is only necessary if a fractional transport equation is used.The transport equation of Wilcock and Crowe (2003) is always combined with their proposed hiding function and uses τ * rm to determine the initiation of bedload motion.

F
Figure A5.Qualitative sketch of continuous interaction between the active surface layer and subsurface alluvium.The water level is displayed in blue, the active surface layer with its constant thickness is light grey, the subsurface alluvium consisting of one layer with variable thickness is dark grey and bedrock is black.The dashed bracket indicates the position of the active surface layer after aggradation or erosion.In case of aggradation (a), the material input from upstream, which is displayed in green, is added to the active surface layer, which is instantaneously homogenised.Any material of the homogenised active surface layer, which exceeds the constant thickness, is transferred to the subsurface alluvium, which is instantaneously homogenised as well.Any homogenisation process is displayed by diagonal stripes.In case of erosion (b), the sediment deficit of the active surface layer is replaced by material from the subsurface alluvium and the active surface layer is instantaneously homogenised.

Table 1 .
Comparison of bedload transport models for steep mountain streams.Estimated calculation speeds refer to the simulation of a 20 km long study reach of a regular mountain river.

Table A1 .
Notation.The following symbols are used in this article.D mArith for the active surface layer D mArith sub D mArith for the subsurface alluvium active surface layer thickness after erosion or deposition but before the interaction between the active surface layer and the subsurface alluvium h s,post active surface layer thickness after the interaction between the active surface layer and the subsurface alluvium www.earth-surf-dynam.net/3/15/2015/Earth Surf.Dynam., 3, 15-34, 2015 32 F. U. M. Heimann et al.: The bedload transport model sedFlow
of previous time step thresh hs,low , thresh hs,high thresholds for the thickness of the active surface layer with thresh hs,low < thresh hs,high thresh hs,prox an alias for the threshold (either thresh hs,low or thresh hs,high ) that is closer to h s