A Channel Network Model for Sediment Dynamics Over Watershed Management Time Scales

A mountain watershed network model is presented for use in decadal to centurial estimation of source‐to‐sink sediment dynamics. The model requires limited input parameters and can be effectively applied over spatial scales relevant to management of reservoirs, lakes, streams, and watersheds (1–100 km2). The model operates over a connected stream network of Strahler‐ordered segments. The model is driven by streamflow from a physically based hydrology model and hillslope sediment supply from a stochastic mass wasting algorithm. For each daily time step, segment‐scale sediment mass balance is computed using bedload and suspended load transport equations. Sediment transport is partitioned between grain size fractions for bedload as gravel and sand, and for suspended load as sand and mud. Bedload and suspended load can deposit and re‐entrain at each segment. We demonstrated the model in the Elwha River Basin, upstream of the former Glines Canyon dam, over the dam's historic 84‐year lifespan. The model predicted the lifetime reservoir sedimentation volume within the uncertainty range of the measured volume (13.7–18.5 million m3) for 25 of 28 model instances. Gravel, sand, and mud fraction volumes were predicted within measurement uncertainty ranges for 18 model instances. The network model improved the prediction of sediment yields compared to at‐a‐station sediment transport capacity relations. The network model also provided spatially and temporally distributed information that allowed for inquiry and understanding of the physical system beyond the sediment yields at the outlet. This work advances cross‐disciplinary and application‐oriented watershed sediment yield modeling approaches.


Introduction
In September 2014, the largest dam removal in history was completed with the removal of two dams along the Elwha River, located in Washington State's Olympic National Park (Figure 1). This historic record was attributed, in part, to the combined volume of sediment impounded behind the two dams,~21 million m 3 Randle et al., 2015). The massive volume of sediment to be released upon dam removal posed major challenges to the removal strategy and uncertainties in the future ecosystem recovery. The scrutiny of reservoir sediment impacts in the Elwha watershed (e.g., Downs et al., 2009;Randle et al., 2015;Warrick et al., 2015) illuminates how critical sediment management is in all phases of a dam's life cycle-design, construction, operations and maintenance, and removal. In each stage, miscalculating reservoir sediment yields can have dramatic consequences on a system's storage capacity, hydropower production, operations, maintenance, costs, safety, and ecosystem impacts. However, most methods used in dam planning and operations do not adequately represent the complex sediment dynamics of river and reservoir systems (U.S. Society on Dams, 2015). Existing methods are especially deficient in steep mountainous regions, where a large portion of global hydropower dams are existing or planned (Lehner et al., 2011a;Lehner et al., 2011b;Zarfl et al., 2015).
Methods commonly used for watershed management and the planning and design of water infrastructure do not adequately integrate the spatial and temporal variations in watershed sediment dynamics over typical infrastructure lifespans (~40 to 100 years). For example, the common approach of empirical at-a-station sediment rating curves uses streamflow and sediment data collected at a single location and for a limited temporal frequency and extent. As a result, peak flow magnitudes that dominate sediment yields and sediment resulting from land management and disturbances may not be adequately represented. The prediction of hillslope sediment supply is often from simple empirically based methods such as the Universal Soil Loss Equation (Wischmeier & Smith, 1978) and its variants (e.g., Revised Universal Soil Loss Equation; Renard et al., 1991). However, these approaches were developed for agricultural fields and are poorly suited for mountainous watersheds (Shen et al., 2009;U.S. Society on Dams, 2015). They also have limited sediment transport mechanisms (e.g., sheet or rill erosion) and only provide mean annual sediment yields. More sophisticated mechanistic models are used such as Water Erosion Prediction Project (Flanagan & Nearing, 1995), designed for rural agricultural watersheds, or HEC-6 (Hydrologic Engineering Center, 1993), designed for scour and deposition in rivers and reservoirs. However, these models do not account for climate variability and extremes, lack the holistic representation of source-to-sink geomorphic processes, and often require estimation of numerous unknown parameters (U.S. Society on Dams, 2015).
Watershed network sediment routing models developed for academic research purposes, on the other hand, are typically developed to investigate domain-specific questions (e.g., fluvial sediment transport, hillslope geomorphology). However, they often lack adequate representation of complexity in other process domains (e.g., distributed hydrology) as well as space and time scale considerations that are relevant to watershed planning and management. The seminal work of Dunne (1997a, 1997b) demonstrated a stochastic sediment input routine to drive sediment routing and storage in channel networks using spatially homogeneous rainfall. Studies by Czuba et al. (2017), Czuba and Foufoula-Georgiou (2014) and Schmitt et al. (2016Schmitt et al. ( , 2018aSchmitt et al. ( , 2018b have further advanced understanding on the network-scale connectivity of sediment sources and sinks. Various other sediment routing models have had sophisticated representations of mixed-size bedload transport such as the network model of Czuba (2018), the network model of Gasparini et al. (2004), The Unified Gravel-Sand model (Cui, 2007), and Morphodynamics and Sediment Tracers in 1-D model (Lauer et al., 2016). However, each of these modeling frameworks lacks integration of both bedload and suspended load as well as advanced representation of either hillslope input, distributed hydrology, or long-term watershed sediment yields.

A Channel Network Model for Sediment Dynamics
Predicting sediment balances in mountainous watersheds requires mechanistic understanding of the intermittency and magnitudes of sediment flux into the channel network and propagation of sediment as it interacts with the river channel. There is disconnect between methods and models used to predict sediment yields for water management and engineering applications and those used for researching fluvial fluxes, forms, and patterns and their linkage to climate and environmental change. Considering these limitations, we developed a parsimonious channel network sediment model to address two key research questions: 1. How can domain-specific research approaches to modeling network sediment dynamics be combined to predict watershed sediment yields over temporal and spatial scales relevant to planning and management? 2. What improvements does a network-scale sediment routing model provide compared to at-a-station based transport estimates in predicting reservoir sediment yields? Our new model connects processes of hillslope sediment supply with fluvial sediment deposition, accumulation, and transport routed across watershed network channel segments. Two novel attributes of our model are that (1) it simulates both bedload and suspended load transport, partitioned between gravel, sand, and mud fractions, and (2) the streamflow forcing is from a physically based distributed hydrology model. The network sediment transport is driven by intermittent mass wasting events represented by a stochastic arrival process, run independently for watershed area above each network segment. The network sediment model operates on a daily time step and outputs annual sediment yields for each segment.
The network sediment model's novel attributes as well as its temporal and spatial scales are relevant to dam and watershed management needs. Predictions of distributed fluvial sediment transport capacity can help, for example, to identify stream reaches where river management (e.g., erosion mitigation) and frequent in situ monitoring may be needed. Predictions of annual sediment fluxes can help with short-term watershed planning and operations. For example, understanding the annual movement of bedload can inform whether channel bed gravel is sufficiently replenished for optimal fish habitat. Modeling of both bedload and suspended load is valuable since they both impact river ecology, geomorphology, and infrastructure, however, often in different ways. Streamflow forcing from an advanced hydrologic model is critical for to capturing the spatial and temporal complexity of watershed sediment dynamics. This is particularly valuable in watersheds with complex topography and climatology, where sediment dynamics are not adequately represented by simplified hydrologic relationships (e.g., streamflow as a function of contributing area). This paper introduces our new network sediment routing model and its implementation in the Elwha River watershed upstream of the Glines Canyon dam. An ensemble of 28 model experiments was conducted over the 84-year dam lifespan using historic meteorological forcing data. Historic hydrologic, sediment transport, and reservoir sedimentation data were also used to parameterize, calibrate, and validate the model. Three model parameters were systematically varied in the model experiments. We also compared the network model performance in predicting lifetime reservoir sedimentation to the traditional approach of at-a-station sediment rating curves. We tested both empirical and physically based approaches to developing these sediment rating curves and streamflow duration curves.

Paper Outline
In the text that follows, section 2 describes the network sediment model processes as well as the methods and rationale used in developing them. Details are included on the model framework, stochastic hillslope sediment supply, suspended load and bedload transport equations, hydrologic forcing, and model experiments. Section 2 also describes the development and application of the sediment rating curve approaches. Section 3 describes the Elwha River watershed background, historic hydrologic and sediment data, hydrologic model input data, empirical sediment rating curves, and measured reservoir sedimentation estimates. Section 4 provides the results and discussion for the hydrologic forcing model, sediment rating curves, and network sediment model. Section 5 provides the paper summary and conclusions. The supporting information provides further details on the development and of the network model and rating curves. A complete list of notations used throughout the paper is at the end of this text. The sediment model code is written in Python programming language. It is publicly available online along with the data used in this study (at github. com/cbev/network_sediment_model).

Methods
An overview of the network sediment model processes follows (section 2.1) and details of model processes are in subsequent sections (sections 2.1.1 to 2.2.3). This is followed by discussion of the bedload transport formula performance (section 2.3) and the hydrologic forcing model implementation (section 2.4). Then, the testing of the network sediment model with 28 model instances is described (section 2.5). In the network model testing, three parameters were systematically varied. These parameters were the hillslope sediment production rate (rate d ) and fractions of sand (λ ss ) and mud (λ ms ) on the channel bed surface that are available for suspended load transport. These parameters are described in their respective sections (sections 2.1.1 and 2.1.2). Finally, we present our approach for using sediment rating curves to predict lifetime reservoir sediment yields in order to compare them to network model outputs (section 2.6).

Network Model Framework
Our sediment model ( Figure 2) connects processes of hillslope sediment supply with fluvial sediment transport and deposition across a Strahler-ordered channel network. Each segment has geomorphic and fluvial state attributes treated as either static or dynamic. Static attributes of each segment are stored in the geomorphic and fluvial state matrix and include segment Strahler order number o j ; segment length l j ; average segment slope S j ; average segment width W j computed from hydraulic geometry (section 2.2.1); downstream segment number ds j ; total contributing area CA T,j ; local contributing area CA L, j ; representative channel bed grain sizes D mean, j , D gb, j , D s,j , D 50, j , and D 90,j (section 2.2.2); and channel bed grain roughness n g,b (H 50 ) j (section 2.2.3). Dynamic attributes of each segment include the volumes of sediment incoming, outgoing, stored on the channel bed, and bedload transport capacity (Table 1 (a)). These sediment volumes are classified according to their grain size and transport mode i (Table 1 (b)). The channel bed volumes are used to compute the dynamic bed surface fractional proportions of gravel F gb ; sand F s ; and sand and silt (i.e., nongravel size sediments) F ssm . Dynamic segment attributes are updated in each daily time step in conjunction with the segment-scale sediment mass balance.
For each daily time step¸a sequence of model processes is conducted at each segment. The processes start at segments with a Strahler order of one and continue downstream to the highest order. The initial model process step is the sediment supply component, which is a simple algorithm that emulates stochastic mass wasting (section 2.1.1). If mass wasting occurs for a given time step, a sediment volume sourced from the segment's local contributing area is transferred into the channel. The partitioning of the deposited volume between gravel, sand, and mud is based on model parameters (α gb , α s , and α ms ). Fractions of sand and mud in the deposit go directly into suspension based on a model parameter (λ mw ). The entire fraction of gravel and remaining fraction of sand in the deposit are added onto the channel bed. In addition, sediment volumes transported out of segments directly upstream (if the segment order is greater than one) are added into the suspended sediment column or onto the channel bed.
Next, suspended load (section 2.1.2) and bedload (section 2.1.3) transport are estimated for the volumes of sediment in suspension and on the channel bed. The transport processes are driven by streamflow outputs from a distributed hydrology model (section 2.4). Hydraulic geometry relationships are also used to compute the transport formula variables of streamflow velocity U j,t , streamflow depth H j,t , and grain shear stress τ g,j,t . Suspended load transport, partitioned between sand and mud volumes, is computed based on models by Patil et al. (2012) and Abad et al. (2008). Bedload transport capacity, partitioned between gravel and sand volumes, is computed from a modified, two-fraction version of the Wilcock and Crowe (2003) formula.

Stochastic Hillslope Sediment Supply
In this study, the primary purpose of the stochastic hillslope sediment supply component is to be a perpetual sediment driver for testing the network sediment routing. The two main premises of the sediment supply forcing component are: (1) the Olympic Mountains are in erosional steady state (Montgomery & Brandon, 2002) with a constant denudation rate and (2) most fluvial sediment is supplied by threshold-dependent and episodic mass wasting events such as landslides, gullying, and debris flows (Buffington & Woodsmith, 2003;Warrick et al., 2011). In the algorithm, the probability that each segment j will have a mass wasting event is from an exponential distribution (e.g., Stark & Passalacqua, 2014): In 1, t L,j is the time lag between mass wasting events for each segment and t L is the mean time lag for the entire network. After each mass wasting event, t L,j is reset for the respective segment. When a mass wasting event occurs, the sediment volume deposited in the segment is the amount that has generated over the preceding lag time t L,j. The generation rate of the sediment volume deposit depends on the local contributing area of each segment CA L,j and a constant soil production rate d . Hence, the total mass wasting deposit volume V d,j is This is an approximate model to provide intermittent sediment supply to channel segments controlled by soil production. In any given year there could be more than one mass wasting event or no event. This level of detail for sediment input is deemed sufficient for this study as our focus is on the network sediment routing model. The total mass wasting volume is divided into proportions of gravel, sand, and mud (i = gb, s, ms) based on multipliers of α gb , α s , and α ms which collectively sum to 1 (i.e., α gb V d,j = V d,gb,j ; α s V d,j = V d, s,j ; α ms V d,j = V d,ms,j ; V d,gb,j + V d,s,j + V d,ms,j = V d,j ) ( Figure 2). From the proportions of sand and mud in the deposit volume, a fraction is directly added into the suspended load of the segment based on the multiplier of λ mw (λ mw V d,s,j , λ mw V d,ms,j ). These sand and mud volumes from the mass wasting deposit are added to the incoming fluvial suspended load volumes from the segment(s) immediately upstream. These suspended load volumes can potentially be transported further downstream based on suspended load transport computations (section 2.1.2). The mass wasting deposit fractions of sand and mud that were not directly added to the suspended load ([1-λ mw ]V d,s,j , [1-λ mw ]V d,ms,j ) along with the entire mass wasting deposit volume of gravel are added to the channel bed (V b,s,j , V b,ms,j , V b,gb,j ). The channel bed also contains gravel, sand, and mud volumes from the previous time step as well as incoming fluvial bedload volumes from the segment(s) immediately upstream. Hence, the changes in suspended load volume (ΔV o,j ) and channel bed sediment volume (ΔV b,j ) following a hillslope mass wasting event are The parameters rate d , t L , α gb , α s , α ms , and λ mw are user-specified and network model results are sensitive to their values. However, there is limited published research to derive their values from; therefore, values were selected largely based on scientific reasoning and model performance. We varied rate d based on values published for the study area soil denudation rates (0.18-0.22 mm/year; sections 2.5 and 3.1). We set t L to be constant at 20 years based on an average landslide return interval for upland catchments in the North Cascades , which are similar to the Olympic Mountains. Additional parameters that we kept constant were λ mw at 0.10 and α gb , α s , α ms each at 0.33. We chose these values since the parameters are empirical and the model performed strongly for these values in the study area.

Suspended Load Transport
For each daily time step, volumes of suspended load (sand and mud) can enter a segment from hillslope mass wasting (Δ[V o,ss,j +V o,ms,j ], 3) and the segments immediately upstream V o,ss,upstream j and V o,ms, upstream j ( Figure 2). The depth-averaged concentration of sand and mud in the water column C ssm, j (kg/m 3 ) is computed using the suspended load volumes, modeled streamflow volume V w,j , water density ρ, and ratio of sediment to water density sg (2.65), as Then, based on C ssm and other local conditions, in-channel settling and erosion of suspended sand and mud may occur prior to suspended sediment transport out of the segment. We made the simplifying assumption that suspended load transport is in local equilibrium and therefore the local settling and erosion rates are from equilibrium transport (Kondolf & Piégay, 2003). Settling and erosion formulae are based on the network model framework of Patil et al. (2012) and equations in Abad et al. (2008). However, in the suspended load erosion rates of these formulae, we incorporate three additional multipliers of λ ss, λ sm , and F ssm . These multipliers implicitly account for shielding of sand and mud on the channel bed by gravel and were necessary for preventing anomalously high suspended sediment transport rates. The multipliers λ ss and λ ms account for the fractions of sand and mud, respectively, on the channel bed that are erodible (i.e., in the active surface layer). The multiplier λ ss also accounts for the fraction of sand on the bed surface that has a small enough diameter to be transported as suspended load. The multipliers λ ss and λ ms are user-specified and network model results are sensitive to their values. However, there is limited published research to derive their values from. In the model testing for this analysis (section 2.5), we varied values of λ ss and λ ms selected based on scientific reasoning and model performance.
The multiplier F ssm accounts for shielding of sand and mud by gravel on the bed surface. It is updated prior to the computation of suspended load as.
2.1.2.1. Sand Following Abad et al. (2008), we approximate the net rate of change in suspended sand as E ss − S ss where E ss is the potential suspended sand erosion rate and S ss is the potential suspended sand settling rate (both in kg/m 2 /s). The approximation of S ss is Here, v ss is the ss settling velocity (m/s), Z R is the Rouse number, and INT (Z R )is the integral of the suspended load concentration profile in the water column. For v ss , we use the following expression (Camenen, 2007) In 8, v is the kinematic viscosity of water (10 −6 m 2 /s), and A, B, and m are the calibration coefficients for natural sand material (Julien, 1995), which are 24, 1.5, and 1, respectively. D ss* is the dimensionless particle diameter defined as Here, g is the acceleration of gravity. The definition of Z R is v ss k·u * , where k is the dimensionless von Kármán constant (0.41) and u * is shear velocity u * = (τ g /ρ) 0.5 . We use the following numerical approxima- Here, the coefficients C 0 , C 1 , C 2 , C 3 , C 4 , C 5 , and C 6 are 1.1038, 2.6626, 5.6497, 0.3822, −0.6174, 0.1315, and −0.0091, respectively. These coefficient values are for a boundary layer thickness to streamflow depth ratio (z b /H) of 0.05 (see Abad and Garcia, 2006, for details).
The approximation of E ss is The first equation in 11 follows Garcia and Parker (1991) with two multipliers we have added (λ ss and F ssm ). The empirical model testing parameter λ ss was varied between 0.010 and 0.030, which were the end members of the range that performed strongly in the study area. The second equation in 11 is the constraint for available sand on the bed surface that can be transported as suspended load. In 11, A E is 1.3•10 −7 and Z u is Here, Re p is the Reynolds number for sediment, defined as Re p ¼ v ss D ss v . The coefficients α 1 and α 2 depend on Re p as (α 1 , α 2 ) = (1, 0.6) for Re p > 2.36 and (α 1 , α 2 ) = (0.586, 1.23) for Re p < 2.36.

Mud
The settling of mud (S ms ) and erosion of mud (E ms ) (both in kg/m 2 /s) are modeled following the approach described in Patil et al. (2012). Each depends on the grain shear stress, t g , relative to the critical shear stress of mud, t c,ms . We calculate t c,ms as (Mitchener & Torfs, 1996) In 13, ρ b,f is the bulk density of mud (or fine sediment) on the bed. Settling of mud occurs when t g < t c,ms , and is computed using the falling velocity of mud v ms as In 15, the constant C vm depends on C ssm as C vm = 0.15 for C ssm ≤ 0.15, C vm = 2.11 for C ssm > 2.11 and C vm = C ssm for 0.15 < C ssm ≤ 2.11 (Hwang, 1989).
Erosion of mud occurs when t g > t c,ms and is computed as In 16, e 0 = 0.2 · 10 −3 kg/m 2 /s. The empirical model testing parameter λ ms was varied between 0.025 and 0.075, which were the end members of the range that performed strongly in the study area. The second equation of 16 is the constraint for available mud on the bed surface that can be transported as suspended load.
Hence, the outgoing suspended sediment load volume changes (m 3 ) are the sum of the erosion and settling rates of sand and mud The subsequent channel bed volume is updated as

Bedload Transport
Bedload transport capacity is computed for gravel and sand using a modified version of the Wilcock and Crowe (2003) bedload transport formula, which is further described below. The gravel and sand bedload volumes transported out of each segment are the minimum of the bedload transport capacity and channel bed volume. The gravel and sand bed volumes consist of the sediment remaining from the previous time step, deposited from hillslope mass wasting (Δ[V b,gb,j +V b,s,j ], 4), and transported from the segments immediately upstream (V o,gb,upstream j , V o,sb, upstream j ) ( Figure 2). The bed surface grain size fractions of gravel, sand, and mud are computed using these volumes as well as the suspended sand and mud that are on the bed from hillslope mass wasting (Δ[V b,ms,j ], 4) and preceding suspended sediment transport computations (Δ[V b, s,j +V b,ms,j ], 18). We make the simplifying assumption that the bed surface grain size fractions of gravel (F gb ) and sand (F s ) are the same as the gravel and sand volume proportions on the bed (i.e., the bed is well mixed) and thus are computed as While the proportions of gravel and sand on the channel bed can change, we make the simplification that their median grain sizes are assumed constant for each segment. These simplifications for the bed surface grain size distribution neglect active layer dynamics in which bedload transport and grain size sorting is restricted to the bed surface layer (e.g., Czuba, 2018;Gasparini et al., 2004). (2003) is a commonly used surface-based bedload transport formula for mixtures of sand and gravel. The formula computes the transport rate for each grain size fraction i on the channel bed surface as

Wilcock and Crowe
For each grain size on the channel bed surface, q bi is the volumetric bedload transport rate per unit channel width (m 2 /s); W i * is the dimensionless bedload transport parameter; F i is the fraction of the grain size on the bed surface; τ ri is the reference shear stress, or the value of t g at which W i * is equal to a reference value of 0.002 (Parker & Klingeman, 1982;Wilcock, 1997). Calculation of t ri using Wilcock and Crowe (2003) is Here, D i is the median grain diameter for each i; τ rmean is the reference shear stress (t ri ) for the mean grain size on the bed surface D mean ; and τ * r mean is the dimensionless τ rmean . Finally, the bedload transport capacity volume (m 3 ) for each fraction V c,i,j is calculated as The Wilcock and Crowe (2003) equations ((21a)-(21c) and (22a)-(22d)) can be applied over an extensive, multifraction grain size distribution to highly resolve multiple transport rates. However, given the uncertainty in the grain size distribution of the channel bed and sediment supply across the study area river network, we used a simplified two-fraction approach. A two-fraction approach, which only differentiates between sand (i = sb) and gravel (i = gb), has been used effectively in literature (e.g., Gasparini et al., 2004;Wilcock & Kenworthy, 2002). We compared Wilcock and Crowe (2003) in its multifraction form (with i for 16 grain size bins) and the two-fraction form (i=sb, gb) using data from the study area, and the difference in results was trivial. Hence, in our model framework, the Wilcock and Crowe (2003) equations are applied in a two-fraction form that we label hereinafter as W&C-2F.
The volumes of bedload transported as gravel and sand, respectively, are assumed to be the bedload transport capacity 23 unless constrained by the volume available on the bed. This is computed as The subsequent channel bed volume is updated as

Network Model Initialization
Channel width, streamflow depth, bed surface grain size diameters, bed grain roughness, and grain shear stress of each segment are parameterized based on observations collected at a single location in the study area watershed. This location, Lake Mills (LM) gage, and its observations are described in section 3.2, however, are referred to in this section.

Hydraulic Geometry
At-a-station hydraulic geometry relationships for LM gage of depth H LM and velocity U LM of streamflow Q LM were developed using the standard approach of power regression. The functions contain empirically derived coefficients (a H , a U ) and exponents (b H , b U ) ( Table 2) and have the following form: The LM gage was in a wide bedrock canyon and observations of channel width W LM had little variation with streamflow. Therefore, we set the LM gage channel width to be constant at the median value of measurements, 41.6 m.
Following Czuba (2018), we scaled streamflow depth H and channel width W for each segment relative to those at LM gage (H LM and W LM ) along with the total contributing areas of LM gage CA T,LM and respective segment CA T,j . We used typical exponents for downstream hydraulic geometry scaling (Leopold & Maddock, 1953;Park, 1977).
Per 28, channel width varies across the channel network; however, since width is locally constant, it does not vary with streamflow. Per 29, channel depth also varies across the channel network; and, per 26, depth varies with changes in streamflow. Rectangular cross sections were assumed throughout the channel network and model duration. Thus, streamflow velocity for each segment was computed for each daily time step by dividing streamflow with the product of depth and width. The hydraulic radius of segments were approximated as the depth.

Grain Sizes
The approaches for modeling shear stress and bedload transport require the characterization of representative grain sizes on the channel bed surface including the median D 50 , 90th percentile D 90 , mean D mean , median gravel D gb , and median sand D S . To assign fixed grain sizes of D 50 , D 90 , D mean , and D gb for each channel  (2017), Snyder et al. (2013), and Barry et al. (2004). We assumed a constant bankfull Shields stress τ * bf for each segment j and solve for grain size at each segment based on its predicted τ * bf . We began with the equation for bankfull Shields stress, In 30, we substituted bankfull depth H bf,j with 29; substituted D 50 with D x,j , which can be D 50,j , D 90,j , D mean,j , or D gb,j ; and substituted τ * bf with τ * bf ;x which can be τ * bf ;50 ; τ * bf ;90 , τ * bf ;mean , or τ * bf ;gb (corresponding to the respective D x,j ). Then, we rearranged 30 with these substitutions as In 31, it is assumed that the bracketed term as well as τ * bf ;50 ; τ * bf ;90 , τ * bf ;mean , and τ * bf ;gb were each respectively constant across all channel segments (Pfeiffer & Finnegan, 2017). Thus, we computed the four τ * bf ;x values at the LM gage using known grain size data from the site (D 50,LM , D 90, LM , D mean, LM , and D gb, LM , section 3.2). Then, we rearranged 31 to solve for grain sizes at each segment in the network (D 50,j , D 90,j , D mean,j , D gb,j ) as Values in the denominator of 32 were constant across all segments in the network. Values in the numerator of 32 are functions of topography that were assumed as static for each respective segment. Hence, D 50,j , D 90,j , D mean,j , and D gb,j can conveniently be computed for each segment.
We assumed D S was constant for all segments and therefore equal to the value computed at the LM gage (section 3.2). This is justified since sand grain collisions are viscously dampened and therefore not subject to abrasion during downstream transport (Jerolmack & Brzinski, 2010).

Shear Stress and Surface Roughness
The model sediment transport equations we used depend on grain shear stress τ g , which is the fraction of boundary shear stress τ that acts only on sediment grains (rather than woody debris, etc.). To isolate grain shear stress from boundary shear stress, we considered the channel roughness n. Channel roughness similarly has a fraction of grain roughness n g that is due to bed grains only and thereby relevant to sediment transport. Channel roughness is typically expressed from Manning's equation as n = H 2/3 S 1/2 U −1 , which can be rearranged as H = (nU) 3/2 S -3/2 . Following Wilcock et al. (2009), this expression for H can be substituted into the equation for boundary shear stress, τ = ρgHS. Then, to compute τ g , n can then be replaced with n g . This resulted in For our model, we developed an approach for computing grain roughness based on published formula and empirical findings. There are many empirical models for directly computing grain roughness, however they were developed across variety of site conditions and therefore each may or may not be relevant to a given application. From the studies of Rickenmann and Recking (2011) and Wang and Dawdy (2014), we obtained nine equations for computing grain roughness in gravel-bedded streams that are either widely used or represent a low or high degree of complexity (Bathurst, 1985;Bray, 1979;Ferguson, 2007;Griffiths, 1981;Hey, 1979;Keulegan, 1938;Leopold & Wolman, 1957;Limerinos, 1970;Smart & Jäggi, 1983;Strickler, 1923; section S1.2). We computed grain roughness using LM gage observations in each of the nine equations and found that the computed grain roughness was nearly constant with streamflow depth in each equation. The computed grain roughness values from the nine equations fell in a narrow range of 0.022-0.027 (median of~0.024). This median value is most closely represented by the equation of Bray (1979) n g,b , which is When we use equation 34 and compare predicted bedload by the W&C equation given fixed bed substrate size as observed at the Lake Mills gage site with bedload measurements (section 2.3 below) the W&C under (over) estimate bedload in low (high) flows. Total roughness computed from Manning's equation and LM gage observations (of depth, slope, and velocity) showed significantly more variability with streamflow depth compared to grain roughness. The range of total roughness was 0.04 to 0.15 and a depth weighted mean n was~0.079. The computed variability of n g,b is a lot smaller (range: 0.022-0.027). We argue that this mismatch in bedload could be a result of disparity between the variation of estimated total n from observations and the calculated n g,b from 35.
An alternative approach to estimating n g,b is proposed that scales the variability of Manning's total roughness with depth to the magnitude of grain roughness computed from 34 (section S1.2). We label this scaling relationship as n g , E (E for the Elwha watershed), and compute it as In 35, n g,b (H 50 ) j is from 34 evaluated at the median flow depth, H 50,j (which is computed using 29 with H LM of the median streamflow discharge Q 50 from 26). Thus, n g,b (H 50 ) j is a static attribute for each segment. n * is Manning's total roughness n normalized by n. Computing n=n is only plausible where there are a broad range of observations available (in our case, at LM gage). So, we conducted a power regression for n * as a function of flow depth using LM gage observations, emulating an at-a-station hydraulic geometry relationship for normalized roughness, as We assume that the coefficient (c ng ) and exponent (d ng ) are constant throughout the network (Table 2).
Hence, in the model, we used the static attribute of n g,b (H 50 ) j , computed from 26, 29, and 34, along with n *,j computed from 36 to calculate n g,E,j with 35. Then, n g,E,j is substituted for n g,j in 33 to compute τ g,j .

Lake Mills Gage Bedload Formula Performance
Our parameterizations of the W&C-2F formula for computing bedload discharge showed strong performance when compared to observations at LM gage. The performance of the W&C-2F formula with grain roughness parameterization by Bray (1979) (n g,b , 34) and our empirically developed formula (n g,E, 35) is demonstrated in Figure 3. Replacing n g,b with n g,E improves the bedload discharge predictions r 2 from 0.44 to 0.51. This improvement is largely due to improved predictions at higher bedload discharge values (>~10 −5 m 2 /s), where the n g,b parameterization overestimated bedload discharge. Since median to high sediment discharges dominate long-term yields, we used the n g,E (rather than n g,b ) grain roughness parameterization in our network sediment model. We tested both grain roughness parameterizations (n g,E and n g,b ) in our sediment rating curve analyses (section 2.6).

Hydrologic Forcing
Our sediment network model uses streamflow inputs defined at the upstream node of each segment of the channel network, for each daily time step. The streamflow inputs are from the Distributed Hydrology Soils Vegetation Model (DHSVM) (Wigmosta et al., 1994). DHSVM is an open-source distributed, physically based hydrology model that combines spatially explicit characteristics of a watershed with meteorological forcing data. DHSVM has been extensively applied in the Pacific Northwest (e.g., Bandaragoda et al., 2019;Cuo et al., 2011;Frans, 2015;Mauger et al., 2016;Murphy, 2016;Sun et al., 2018;Wigmosta et al., 1994).
To derive the study area soil depth grid and channel network topology, we used open source research software specifically developed for setting up DHSVM (PNNL, 2018). Soil depths were assigned based on local elevation and combined curvature variables within a specified a soil depth range. Segment attributes were assigned by DHSVM following channel network extraction. Historical meteorological forcing data was disaggregated from a daily time step to 3-hr time step and incoming longwave and shortwave radiation and relative humidity were estimated. This was done using the Mountain Microclimate Simulation Model algorithms (Thornton & Running, 1999) as implemented in Bohn et al. (2013).
DHSVM was run in the study area at 30 m grid cell spatial resolution and 3-hr time step from years 1915-2011 (97 years). Streamflow was output at each segment at the 3-hr time step, and then aggregated to daily streamflow predictions. We calibrated DHSVM with a one-at-a-time method to fit the modeled streamflow outputs to observations at LM gage for water years 2006 and 2007. We limited the calibration to a computationally inexpensive approach based on parameters known to be sensitive for DHSVM (Kelleher et al., 2015). We focused on adjusting the soil texture and depth parameters (final parameters in Table S2), as they had the strongest influence (compared to other parameters) on the modeled hydrographs. However, the input soil texture data in the study area is relatively low resolution (three generic soil classes), which could limit model performance where physically based soil parameterization is significant. We evaluated the DHSVM performance during the remaining (noncalibration) periods of streamflow observations periods (water years 1994-1998, 2004-2005 and 2004-2011).

Network Model Instances
We tested our model performance and calibrated parameters by comparing model outputs to observations of bedload and suspended load at LM gage. In this paper we present and ensemble of 28 model instances (unique parameter sets) from varying three model parameters within limited ranges. We varied the hillslope sediment production rate (rate d ) as 0.18, 0.20, and 0.22 mm/year. We varied the fractions of sand (λ ss ) and mud (λ ms ) on the channel bed surface that are available for transport as 0.010, 0.015, 0.025, and 0.030 for λ ss , and 0.025, 0.050, and 0.075 for λ ms . These are the plausible ranges of parameter values we see realistic for the study site conditions. Thus, we limited our simulations to show the cumulative sediment accumulation for the study period to this given parameter combination.
We conducted a single model spin-up run to establish the initial quasi-equilibrium conditions of the channel bed for all segments (Cui et al., 2006). This spin-up run and each of the 28 model instances were run over 97 years ) at a daily time step. We evaluated model outputs for the 84 years  that the dam in the study area existed.
The reservoir is not explicitly incorporated into the network model. However, there are model network segments that overlay the former reservoir surface and dam, which we used to estimate the final total reservoir sedimentation. We summed the channel bed sediment volumes (V b,i,j ) for segments overlaying reservoir surface (j ∈ reservoir surface) at the final time step (t final ) along with summing the sediment volumes transported out of segments (V o,i,j ) crossing over the GC dam (j ∈ GC dam) throughout the model evaluation period. Suspended sediment was multiplied by the empirical reservoir trap efficiency rte. Sediment volumes were converted from sediment in transport (with density ρ s ) to deposited (bulk) sediment in LM using the bulk densities of fine (ρ b,f ) and coarse (ρ b,c ) sediments. In equation form, lifetime reservoir sedimentation is as follows.

Lake Mills Gage Sediment Rating Curve Testing
We developed lifetime reservoir sedimentation predictions using the traditional approach of at-a-station sediment rating curves and streamflow duration curves in order to compare to network model outputs.
We tested two forms of sediment rating curve equations, regression, and physically based, driven by two forms of streamflow duration curves, measured and modeled, at the LM gage. We compared the lifetime reservoir sedimentation predicted by the rating curves to our network model outputs. The two streamflow duration curves for the LM gage were constructed using: (i) the extent of complete water year observations available at LM gage and (ii) model outputs from DHSVM over the 84-year dam lifespan. The regression sediment rating curve equations were bedload and suspended load equations developed from streamflow and sediment observations by Curran et al. (2009) (section 3.4). The two physically based sediment rating curve equations are W&C-2F (equations 21a-21c and 22a-22d) for bedload. The two equations use with grain shear stress computed from the two different roughness parameterizations tested in this study (n g,b and n g,E ; section 2.2.3). For these two bedload rating curves, we also used parameterizations for the LM gage streamflow depth, streamflow velocity, channel width, and channel bed grain sizes presented in previous sections. The channel bed gravel and sand fractions were assumed constant at values computed using the Lake Mills channel bed grain size distributions (section 3.2). Since predam removal measurements of the LM gage local channel slope were not available, we used predam removal Lidar data from Entrix, Inc. (2009) and criteria from Paola and Mohrig (1996) to estimate the slope as 0.005 (section S1.1).
We integrated the four daily sediment rating curves (one suspended load and three bedload) with both daily streamflow duration curves (measured and modeled) to compute the mean annual sediment loads. Then, we multiplied the annual loads by the reservoir lifespan to get the lifetime reservoir sedimentations. We used observations of LM gage suspended load grain size fractions to develop a linear equation, which approximated the proportions of sand (coarse sediment) and silt (fine sediment) in the suspended load as a function of streamflow (section S1.3). As with the network model outputs (section 2.5), we adjusted suspended load yields using the reservoir trapping efficiency and converted the transported sediment volumes to deposited sediment volumes using the bulk densities of fine and coarse sediments.

Elwha Watershed Regional Setting
The Elwha River watershed model study domain (Figure 1) spans the 636 km 2 watershed upstream of the former Glines Canyon (GC) dam, a 64 m tall concrete arch dam located 21.6 km upstream of the Elwha River mouth. The GC dam construction was completed in 1927 and removal began in 2011, a project lifetime of~84 years. The GC dam's reservoir, Lake Mills (LM), had the capacity for 50 million m 3 of storage at the time of construction . The study area is within the Olympic National Park, which is protected forested land with minimal human impacts aside from the former dam and reservoir, hiking trails, campsites, and historic homesteads.
The Elwha watershed has a maritime climate with dry summers and wet winters. The watershed elevation ranges from sea level to 2,250 m and mean annual precipitation ranges from 1,400 mm at low elevations to 5,600 mm at high elevations. Precipitation is predominantly rain at elevations below~1,200 m, but in upper 10.1029/2019MS001852 Journal of Advances in Modeling Earth Systems elevations is largely dominated by snow in the winters (Duda et al., 2011). There are small glaciers (total area of <2.1 km 2 ) in the headwaters of the Elwha watershed that have been rapidly retreating since 1980 (Pelto, 2018). The 5th, 50th, and 95th percentile daily flows at LM gage are 8.9, 29, and 92 m 3 /s (USGS, 2016).
Geologically, the study area lies within the Olympic Subduction Complex, which is primarily composed of metasedimentary phyllites and schists . The Olympic Mountains are considered to be in tectonic steady-state, although rock uplift rates vary spatially across the mountain range (Montgomery & Brandon, 2002). In the Elwha watershed, estimated denudation rates over geologic time range from 0.3 mm/year at the coast to~0.6 mm/year in the mountainous uplands (Batt et al., 2001;Brandon et al., 1998). Within the study area, the estimated modern erosion rate is~0.18 mm/year (Brandon et al., 1998). At the headwaters, uplift rates and the dominant erosional processes that include debris flows, rockfalls, and shallow and deep landslides strongly influence the steep slopes and channel morphology . The rapid erosion processes occur mostly during storm events and supply the majority of sediment to headwater valleys (Buffington & Woodsmith, 2003). The Elwha River is broadly classified as a gravel-bedded river and alternates between narrow bedrock canyons and wider alluvial reaches for much of its length. This alternating boundary composition leads to differential sediment transport behavior along the river .

Streamflow, Sediment, and Channel Geometry Data
The streamflow, sediment, and channel geometry measurements used in this study were collected by the U. S. Geological Survey (USGS) at USGS gaging station number 12044900, Elwha River above Lake Mills (hereinafter referred to as LM gage) (Figure 1). The LM gage was located in a bedrock canyon approximately 400 m upstream of LM and 4.3 km upstream of GC dam (Curran et al., 2009). The LM gage drainage area of 513 km 2 is~81% of the GC dam drainage area. The LM gage is the only long-term gaging station that existed upstream of the GC dam, with data collected from March 1994 to November 2015. Data were collected at LM gage by the USGS as part of their nationwide water monitoring program (USGS, 2016) and for special studies done in anticipation of dam removal (Childers et al., 2000;Curran et al., 2009). The temporal extent and number of samples for streamflow, channel geometry, bedload, channel bed/gravel bar surface distribution, and suspended sediment observations that were used in our study are summarized in Table 3.
The observed grain size distributions of the channel bed and one gravel bar were used for computing the model parameters for LM gage channel bed grain diameters (D 50,LM , D 90, LM , D mean, LM , D gb, LM , D S; section 2.2.2) and fractions (F S, LM, F gb, LM ). The bedload discharge grain size distributions were also used to evaluate the channel bed parameters. Each of the 42 bedload discharge observations consisted of multiple samples across the channel transect, totaling to 156 samples. Grain size distributions were measured for all 156 samples and their distribution curves are shown in Figure 4. Similarly, the channel bed surface observation consisted of five samples across the channel transect and their distribution curves are shown in Figure 4. The grain size distribution curves for the two gravel bars are also included in Figure 4. In computing the channel bed grain size parameters, we included all five channel bed surface observations and the one downstream gravel bar observation. We excluded the one upstream gravel bar sample as an outlier since its distribution is relatively coarse compared to the channel surface, downstream gravel bar, and bedload grain size distributions. We interpolated the representative grain sizes and fractions for each respective sample, and then took the average of the six samples as the model parameter. The final model parameters are D S of 0.93 mm, D 50, LM of 9.7 mm, D mean, LM of 12.5 mm, D gb, LM of 13.2 mm, D 90, LM of 27.5 mm, F S, LM of 0.37, and F gb, LM of 0.63. We used suspended load grain size distribution data at the LM gage to estimate the median suspended sediment grain size to be 0.25 mm. The corresponding streamflow, suspended sediment load, and suspended sediment grain size distribution data were also used to develop the empirical relationship for predicting the proportion of sand and mud in suspended sediment load as a function of streamflow (section S1.3).

Hydrologic Forcing Model Input Data
The DHSVM input data were a 30 m digital elevation model of the Puget Sound region of Washington from Finlayson et al. (2000); soil texture data from NRCS STATSGO2 (Soil Survey Staff, 2016); and vegetation data from the USGS National Land Cover Database 2011 edition (USGS, 2014). We determined the soil depth range from STATSGO2 data and manual calibration to be 0.75-3 m. We also specified a minimum channel initiation source area of 1 km 2 . We obtained this by developing a slope-area plot with the watershed digital elevation model and estimating the transition from hillslope to fluvial regime to be at a drainage area of 1 km 2 ( Figure S4.; for background see Montgomery & Dietrich, 1994;McNamara et al., 2006;and Montgomery & Foufoula-Georgiou, 1993). The stream network of the study area includes 293 channel network segments ranging from lengths of 0.042 to 6.68 km.

Journal of Advances in Modeling Earth Systems
For the meteorological forcing of DHSVM, we evaluated various gridded historical meteorological data sets (Daly et al., 1994;Kalnay et al., 1996;Livneh et al., 2013;Livneh et al., 2015;Salathé et al., 2014;Skamarock et al., 2008) using open source automated retrieval, preprocessing, and visualization software (Phuong et al., 2019). We compared the data sets to observations of temperature and precipitation in the Elwha watershed, and checked that annual precipitation from the gridded data sets in the LM gage tributary was greater than annual streamflow from observations for the LM gage. The Salathé et al. (2014) data set best matched observations however it only spanned years 1950-2010, which does not cover the lifespan of the GC dam. Only the Livneh et al. (2013) data set covered the GC dam lifespan, encompassing data for years 1915-2011. Hence, we developed a hybrid product of the Salathé et al. (2014) and Livneh et al. (2013) data sets (section S3) at 1/16°spatial resolution (~6km grid size) and daily temporal resolution.

Empirical Sediment Rating Curves and Reservoir Sedimentation Volume Estimates
In our testing of sediment rating curves for estimating reservoir sedimentation, we used sediment rating curve equations from Curran et al. (2009) Curran et al. (2009) to estimate the mean reservoir trap efficiency of LM to be 0.86. We used this reservoir trap efficiency in our study. However, we consider it a limitation of our approach since the value is based on only 2 years of data. We expect that the trap efficiency would vary with changing streamflow and reservoir storage ratios.
The LM sedimentation volume at the time of the GC dam removal was estimated by Randle et al. (2015) to be 16.1 ± 2.4 · 10 6 m 3 . The LM sedimentation composition was also estimated by Randle et al. (2015) to be 44% fine sediment (silt, clay <0.062 mm) with a bulk density of 1.13 g/cm 3 ; and 56% coarse sediment (sand, gravel, cobble >0.062 mm) with a bulk density of 1.71 g/cm 3 . These bulk densities were used in our computations for reservoir sedimentation 27 and suspended load transport of mud 13.

Hydrologic Forcing Model
The hydrologic forcing model performance, which was assessed using the Nash-Sutcliffe model efficiency (NSE; Nash & Sutcliffe, 1970) and coefficient of determination (r 2 ), was acceptable at the monthly and annual time scales but weak at the daily time scale. The NSE is a common metric used to assess the predictive power of hydrologic models. Its value can range from <0 to 1. Generally, a hydrologic model with NSE > 0.5 for a monthly time step is considered acceptable (Moriasi et al., 2007). The NSE and r 2 over the DHSVM hydrologic model calibration period (water years 2006-2007) for daily flow were 0.53 and 0.57 and for monthly flow were 0.67 and 0.71. The NSE and r 2 over the evaluation period (water years 1994-1998, 2004-2005, and 2004-2011) for daily flow were 0.30 and 0.43 and for monthly flow were 0.50 and 0.53. The NSE and r 2 over the evaluation period for annual flow were 0.84 and 0.85.
The model performed poorly in predicting extremely high and low streamflows, but performed stronger for moderate, higher-frequency flows. The model also overpredicted snow accumulation in most water years. Thus, modeled flow was generally lower than observations in cold months and higher than observations in warm months (i.e., snowmelt period). For flows of 2-40% exceedance probability (EP), the predictions were up to 26% greater than observations. For flows of 40-90% EP the predictions were up to 36% lower than observations. For flows over 90% EP (i.e., extremely low flows) the model predictions were poor, however model performance at extremely low flows has minimal impact on sediment transport predictions. Peak flow magnitudes (~<2% EP) were underrepresented by up to 54%, which thereby misrepresented peak sediment yields.
The limitations in the model performance are primarily attributed to the meteorological forcing data (section S3). The gridded meteorological data sets we used are generally more appropriate for capturing the average regional climatology rather than isolated extreme events. The data sets are gridded at a relatively low resolution (~6km) compared to the scale of complexities in the terrain and hydroclimatology of the Elwha watershed. They largely misrepresented the regional temperature and precipitation lapse rates (Minder et al., 2010). They also did not adequately capture regional atmospheric river events (Neiman et al., 2008;Currier et al., 2017) that are expected to be responsible for large, short pulses of streamflow and sediment discharge. Future improvements in gridded hydrometeorology data sets are expected to significantly improve hydrology model performance, especially for >2% EP flood events and snow accumulation.

Sediment Load Estimates From Rating Curves at Lake Mills Gage
Lifetime reservoir sediment yields predicted from at-a-station sediment rating curves driven by measured streamflow were primarily influenced by suspended load yields and peak streamflow discharges ( Figure 5 and Table 4, top section). The rating curve total reservoir sediment yields exceeded the measured range largely because the rating curves do not take sediment supply limitations into account. At ≤2% exceedance probability (EP), the suspended load rating curve exceeded the bedload discharge rating curves ( Figure 5 inset). Similarly, rating curve suspended load yields dominated the total reservoir sedimentation predictions (suspended load yield was 73-84% of total yield) and caused the rating curve predictions to largely exceed measurements (Table 4, top section). Most of the suspended load rating curve yield was mud, which exceeded the measured reservoir mud fraction estimate (upper end of range) by 51%. The sand fraction of the suspended load rating curve yields was toward the upper end of the measured reservoir sand volume estimate (measured reservoir sand estimate was bedload and suspended load combined). Also at ≤2% EP, bedload yield predictions from the equation of W&C-2F with n g,b were highest followed by predictions from the regression equation and then W&C-2F with n g,E. (Figure 5 inset). Similarly, rating curve total reservoir sediment yields followed the order of the bedload rating curves at peak flows (Table 4, top section). The rating curve predictions of gravel from W&C-2F with n g,b exceeded the measured reservoir gravel fraction estimated range while the W&C-2F with n g,E prediction was below the estimated range.
Lifetime reservoir sediment yields predicted from the at-a-station sediment rating curves driven by modeled streamflow ( Figure 6 and Table 4, bottom section) were 50-59% lower than the yields from curves driven by measured streamflow. The rating curve total reservoir sediment yield was within the measured range when bedload was computed using W&C-2F with n g,b , however, lower than the measured range for the other two cases (regression and W&C-2F with n g,E ). The lower rating curve sediment yields are mainly because suspended load rating curve peak values from modeled streamflow are substantially less than those from measured streamflow. As previously discussed (section 4.1) the hydrologic forcing model substantially underpredicted peak flows. The yield from the suspended load rating curve using modeled streamflow had mud and sand fraction volumes below the measured reservoir mud and sand volume estimates (Table 4, lower portion). The yields from the bedload rating curve using modeled streamflow marginally changed compared to those driven by measured streamflow, and rating curve gravel yields were again outside the ranges estimated from measurements.
Overall, the limitations and uncertainties in using the streamflow rating curve and at-a-station sediment rating curves were consistent with their known limitations (e.g., Vaughan et al., 2017). This verified the need for our network model approach incorporating distributed sediment supply and transport dynamics.

Network Sediment Model 4.3.1. Predicted Volumes of Lifetime Reservoir Sedimentation and Model Parameterization
Overall the results showed that the network model performance in predicting the measured total reservoir sedimentation and its estimated grain size volume fractions (particularly gravel and mud) was superior compared to the at-a-station sediment rating curves. Along with these improved sediment yields, the network model provided spatially and temporally distributed information allowing for inquiry and understanding of the physical system beyond the sediment yields at the outlet. Of the 28 network model instances, 25 were within the uncertainty range of the measured total reservoir sedimentation (13.7 · 10 6 -18.5 · 10 6 m 3 ; Randle et al., 2015) (Table 5 and Figure 7). The remaining three instances were within 12% of the lower end of the measured uncertainty range. The model-predicted mud volumes all fell within the uncertainty range of the measured reservoir mud fraction estimate (5.6 · 10 6 -8.4 · 10 6 m 3 ). The model-predicted gravel volumes mostly fell within uncertainty range of the measured reservoir gravel fraction estimate (2.4 · 10 6 -3.6 · 10 6 m 3 ), except two which were within 1% of the range. The model-predicted sand volumes mostly fell within the uncertainty range of the measured reservoir sand fraction estimate (4.7 · 10 6 -7.3 · 10 6 m 3 ). However, all of the model-predicted sand volumes fell below the estimated sand fraction (7.0 · 10 6 m 3 ) and eight were below the estimated uncertainty range. Seven of the low sand predictions were within 9% of the uncertainty range and one was 19% lower.
As expected, the model-predicted total reservoir sedimentation volumes generally correlated with the soil production rates in the hillslope mass wasting algorithm (rate d ; Table 5). The three model predictions below the measured total volume all had a soil production rate of 0.18 mm/year and had the lowest modeled sand Figure 5. Daily measured streamflow duration curve (black curve; right side y axis) and corresponding at-a-station sediment rating curves (red, green, blue, and cyan curves; left side y axis) for Lake Mills gage. In Randle et al. (2015), the "coarse fraction" consisting of boulders, cobble, gravel, and sand were reported together as 9.0 ± 1.9 × 10 6 m 3 . It was also reported that 1/3 of the coarse fraction was boulders, gravel, and cobbles, and 2/3 was sand. In this table, we report the result of applying these fractions to the estimated coarse fraction volume as well as the error range (±1.9 × 10 6 m 3 ).
b Suspended load predictions are from the regression-based approach but are added the bedload predictions in order to compared to the total predicted lifetime sediment load. Figure 6. Daily modeled streamflow duration curve (black line; right side y axis) and corresponding at-a-station sediment rating curves (blue and red curves; left side y axis) for Lake Mills gage. Also shown are network model sediment outputs (green and cyan dots, left side y axis) plotted against corresponding modeled streamflow exceedance probability for daily outputs at the network segment containing Lake Mills gage.   3.0 (2.4-3.6) 6.0 (4.7-7.3) 7.0 (5. fraction volumes. Although in our limited sensitivity analysis we did not vary the mass wasting deposit grain size composition (α gb , α s , α ms ), average frequency (t L ), and sand and mud fractions directly added into suspension (λ mw ), their constant values were important to model performance. Further research is needed to provide guidance on the mass wasting deposit parameter value estimation and/or improved model representation of the relevant physical processes.
The network model performance was also sensitive to parameterizations that dictated suspended load: λ ss , λ ms , and F ssm in the suspended load equations and λ mw in the mass wasting algorithm. The order of magnitudes for λ ss , λ ms , and λ mw prevented suspended sediment transport from being anomalously high. Without the four parameters (i.e., setting them all equal to one), the model essentially transported all suspended load as soon as it was supplied to the watershed. This caused nearly all sand to be transported as suspended load and thus the bedload consisted of nearly all gravel. This, furthermore, led to excessive gravel accumulation on shallow-sloped channels and caused the model to underpredict the gravel deposited in the reservoir. Although the volumes of sand and mud deposited in the reservoir fell within the estimated uncertainty ranges without incorporation of these four parameters (λ ss , λ ms , λ mw , and F ssm ), the overall model processes and results were physically unrealistic. An increase in λ ss correlated with a decrease in the gravel fraction of the predicted reservoir sedimentation volume (Table 5). There was not a clear pattern in the model results from varying λ ms , suggesting that this parameter is relatively insensitive in the range we used and could be tested more broadly in future work.

Interannual Reservoir Sedimentation Variability
In addition to the prediction of lifetime reservoir sedimentation volume, reservoir operations and maintenance can also benefit from the understanding of interannual sediment loads and their correlation with streamflow. This is helpful, for example, in predicting ongoing reservoir storage loss and developing reservoir sediment flushing routines. Here, we discuss interannual sediment and streamflow yields from Model Instance 23 as an example of model performance. We use Model Instance 23 for these and subsequent results sections (sections 4.3.3 and 4.3.2) since instance 23 was one of the strongest performing predictions of the  Table 5. modeled reservoir sedimentation volumes. For Model Instance 23, the range of variability in interannual total sediment loads was~2.1 times the mean, and in the grain size fraction volumes was 1.7-3.6 times the respective means ( Figure 8). Patterns of annual magnitude of streamflow and total sedimentation volume were poorly correlated (r 2 =0.19). Correlation between gravel volumes and the other fractions were poor (r 2 =0.27 for sand, r 2 =0.11 for mud); however, correlation between sand and mud fractions were strong (r 2 =0.84). The cases with poor correlation were reasonable considering the network model's distributed and stochastic nature of sediment supply as well as the time lag in sediment transport of different grain size fractions throughout the watershed.

Sediment Rating Curve Comparison at Lake Mills Gage
Suspended load from the network model was driven by suspended sediment supply rather than the hydrologic forcing. Compared to the suspended sediment rating curve driven by modeled streamflow (Figure 6), the network model suspended sediment predictions at LM gage are dramatically less sensitive to streamflow. It is notable that at peak streamflow (~<2% EP), suspended load from the network model did not steeply increase as the suspended load rating curve from the regression-based equation. This is because the network model accounts for sediment supply limitations.
As expected, the network model bedload predictions were sensitive to the hydrologic forcing ( Figure 6) due to the dependence of the W&C-2F equations on shear stress. The scatter in the network model bedload predictions relative to the bedload rating curve is largely due to the difference in gravel and sand fractions on the channel bed. The bed fractions evolve in the network model however not in the bedload rating curve. We expect that more advanced representation of channel geomorphic conditions in network model (e.g., dynamic slope and grain sizes, advanced active layer dynamics) would further reflect the dynamicity of bedload.

Network Sediment Dynamics
In watershed management, spatial and temporal understanding of network-wide sediment dynamics also valuable in addition to long-term yields. For example, identifying locations and time periods of high sediment supply and/or transport can help in planning more effective monitoring and mitigation measures.
To serve this interest, our network sediment model allows for intercomparison of sediment dynamics for different segments, years, grain size fractions, and sediment volumes (Figure 9, example of Model Instance 23). Our ability to analyze network-wide model behavior also allowed us to verify that the simulated sediment transport and storage dynamics were physically realistic and reflected the geomorphology of the study area as it has been described in previous studies (e.g., Warrick et al., 2011).
The stochastic nature of the model's mass wasting algorithm is reflected in the mass wasting volume depth and frequency at different channel segments (Figure 9, row 1). The differences in suspended load and bedload behavior (e.g., supply and transport limitations) at a given location and over time and space are reflected in the bedload transport capacity (Figure 9, row two), volumes transported out ( Figure 9, row three), and volumes stored on the bed (Figure 9, row four). Mass wasting deposit volumes (i.e., sediment supply) are conserved in the volume transported out and/or stored on the bed at each respective segment. Sediment volumes transported out of segments show that the relative magnitude of suspended load sand and mud increased relative to bedload gravel and sand moving downstream (Figures 9a-9c). Thus, the model captured how suspended load transported relatively rapidly with streamflow and compounded in higher order streams. Bedload, on the other hand, transported relatively slowly and was more often stored on the channel bed, particularly in low slope reaches of the channel network. Sediment volumes stored on the channel bed showed that the bed-primarily its gravel fraction-was generally more dynamic moving downstream, and bedrock streams (i.e., no sediment stored on the channel bed) were more common in upper reaches. There was little channel bed storage of mud which suggested that mud was largely supply limited. By contrast, bed storage of gravel and sand varied in time and space, showing that these larger grain size fractions varied between being supply and transport limited.

Primary Model Limitations and Future Work Recommendations
Because of the partial transport limitations of gravel and sand (section 4.3.4), we expect that their model-predicted volumes would have been closer to their measured sedimentation estimates (i.e., higher) if the forcing hydrology model had improved representation of peak streamflow events (section 4.1). The multipliers we incorporated into the suspended sediment modeling approach to represent gravel shielding and grain size partitioning of sand on the bed (λ ss, λ ms, and F ssm ) as well as the entrainment from mass wasting processes (λ mw ) were significant simplifications of the physical processes. However, these factors were critical to producing suspended load predictions at the LM gage that were not anomalously high compared to observations (section 4.3.1). The need for these multipliers highlights future work opportunities to improve upon our relatively simplistic approach of combining bedload and suspended load dynamics at the network scale as well as the partitioning of sand transport between bedload and suspended load. This future work would need to encompass improved process representation of their combined deposition, storage, erosion, and transport dynamics. Additional processes not included in our model that could be incorporated in future work would be accounting for the lateral exchange between sediment in the channels and floodplains (e.g., Lauer et al., 2016) as well as the channel bed active surface layer (e.g., Czuba, 2018;Gasparini et al., 2004).
The stochastic sediment supply algorithm was valuable for emulating the steady-state erosion and episodic mass wasting events of the study area. However, our model representation of network sediment supply was a first-order approximation, which can be improved with process-based hillslope mechanisms. These improvements include spatially explicit parameterizations for hillslope mass wasting (e.g., Benda & Dunne, 1997a;Czuba et al., 2017;Strauch et al., 2018) and additional sediment supply processes such as erosion. From a management perspective, it would also be beneficial to demonstrate the locations of dominant sediment supply as well as the connectivity of the hillslope sediment supply deposits to downstream sinks (e.g., Czuba, 2018;Schmitt et al., 2016Schmitt et al., , 2018aSchmitt et al., , 2018b.

Conclusions
Our watershed network sediment modeling approach and findings address an unmet need in applied engineering and management of predicting sediment yields in mountainous watersheds. Our model was developed by integrating empirical data sets with modeling approaches from hydrology, geomorphology, and sediment transport research fields. Compared to at-a-station sediment rating curves, our network model showed strong improvement in predicting short-and long-term reservoir sediment yields as well as understanding spatial and temporal dynamics of watershed sediment processes.
The key novelties of our network sediment modeling approach were (1) combining bedload and suspended load transport modeling and (2) using streamflow forcing from a physically based distributed hydrology model. We also incorporated a simplified stochastic, semidistributed hillslope mass wasting algorithm, which can be replaced by process-based hillslope input mechanisms. In taking these unprecedented steps to network sediment modeling, we illuminated both synergies and gaps in modeling physical processes that are critical to watershed sediment yields over engineering and management time scales. We demonstrated the particular needs for improved understanding and modeling of (1) landslide sediment deposit volumefrequency, grain size fraction distributions, and the immediate fate of each grain size fraction (e.g., deposited on channel bed, transported in suspension); (2) combined bedload and suspended load transport dynamics; (3) improved partitioning of sand between bedload and suspended load; and (4) hydrologic modeling of peak flow events. These improvements would likely help to address the underprediction of sand in the network model. Improvements would also reduce the uncertainty due to the simplifying model parameters that we introduced in modeling suspended load (λ ss , λ sm , λ mw , and F ssm ). We intend to motivate and inform future work in these specific research areas as well as interdisciplinary network sediment modeling approaches. These advances can further serve applied engineering and management needs for better understanding of watershed sediment yield predictions under increasing environmental change.