A network scale, intermediate complexity model for simulating channel evolution over years to decades

Excessive river erosion and sedimentation threatens critical infrastructure, degrades aquatic habitat, and impairs water quality. Tools for predicting the magnitude of erosion, sedimentation, and channel evolution processes are needed for effective mitigation and management. We present a new numerical model that simulates coupled river bed and bank erosion at the watershed scale. The model uses modified versions of Bagnold’s sediment transport equation to simulate bed erosion and aggradation, as well as a simplified Bank Stability and Toe Erosion Model (BSTEM) to simulate bank erosion processes. The model is mechanistic and intermediate complexity, accounting for the dominant channel evolution processes while limiting data requirements. We apply the model to a generic test case of channel network response following a disturbance and the results match physical understanding of channel evolution. The model was also tested on two field data sets: below Parker Dam on the lower Colorado River and the North Fork Toutle River (NFTR) which responded dramatically to the 1980 eruption of Mount St. Helens. It accurately predicts observed channel incision and bed material coarsening on the Colorado River, as well as observations for the upstream 18 km of the NFTR watershed. The model does not include algorithms for extensive lateral migration and avulsions and therefore did not perform well in the lower NFTR where the channel migrated across a wide valley bottom. Despite its parsimony, we are confident in the utility of the model for simulating channel network response to disturbance.


Introduction
Excessive river erosion and sedimentation are triggered by a variety of watershed disturbances which alter natural flow and sediment dynamics. For example, urbanization increases discharge (Hollis, 1975;Rosburg et al., 2017), channel straightening increases slope (Simon, 1989), and dam construction decreases sediment supply and modifies flow regimes (Williams and Wolman, 1984). Channel instability and sediment imbalance threatens infrastructure, degrades aquatic habitat, and impairs water quality. Landowners and environmental resource agencies often respond to these threats by attempting to stabilize channels, sometimes without success (e.g. Miller and Kochel , 2009). Stream stabilization projects fail because designers do not account for altered hydrology and sediment supply and because of the inherent uncertainty of channel response (Simon et al., 2007;Roni and Beechie, 2013;Wohl et al., 2005;Bernhardt and Palmer , 2007). It is challenging to predict how streams will adjust and what new equilibrium state -if any -they will attain.
Numerical modeling can address this issue by providing a simple and reproducible way to (1) assess channel sensitivity to disturbance and (2) predict channel adjustment. While morphodynamic modeling has advanced in recent years, most of the research has focused on large spatial and temporal scales (e.g. landscape evolution models (Lague, 2014)) or individual processes (e.g. bar formation (Nelson et al., 2015)). Models that predict channel changes at intermediate spatial and temporal scales (10s -100s km 2 watersheds; 10s -100s of years) are needed to help guide river restoration and management.
Recent research has attempted to fill this gap with regime-based models of river response (Eaton and Millar , 2017), watershed-scale accounting of sediment dynamics (Parker et al., 2015;Foufoula-Georgiou, 2014, 2015;Schmitt et al., 2016;Soar et al., 2017), and mechanistic bank erosion modeling (Langendoen et al., 2012;Stryker et al., 2017). These approaches are useful but may not account for relevant erosion processes or require significant amounts of data, making it difficult to assess uncertainty and provide results useful to managers. The aim of this study was to develop a networkscale morphodynamic model for simulating channel incision and bank erosion with limited data requirements. To achieve this goal, we use specific stream power (Bagnold , 1966), allowing us to model channel erosion and deposition without simulating detailed flow hydraulics. Additionally, the model was designed to be transparent about uncertainty, explicitly translating variability in model inputs into probabilistic predictions of channel evolution. This paper introduces this new stream power-based morphodynamic model -the River Erosion Model (REM). REM is designed primarily for modeling channel evolution in smaller watersheds (i.e. 10s -100s km 2 ) with cohesive banks, integrating a bank stability model based on Lammers et al. (2017) with novel stream power based sediment transport equations (Lammers and Bledsoe, 2018). Unfortunately, watershed-scale data on channel response are rarely available for these types of systems. We therefore test REM on a generic watershed responding to base-level fall as well as two field datasets of rivers responding to different types of disturbance. The first is a reach of the lower Colorado River which incised and coarsened after Parker Dam was constructed in 1938. The second is the North Fork Toutle River (NFTR) which has followed a complex trajectory of channel change following massive sediment deposition from the eruption of Mount St. Helens in 1980. Applying REM to these complex systems tests the basic model processes, explores uncertainty and model sensitivity, and pushes the limits of model application, determining the range of conditions for which it is most suitable.

Stream power
Many models use the standard step method or a simple flow resistance relationship to compute flow depth, velocity, and shear stress (e.g. El Kadi Abderrezzak et al., 2008;Allen et al., 1999). In contrast, we use specific stream power to directly model channel incision and bank erosion. Specific stream power is the power available to do work in the stream, normalized by bed area (Bagnold , 1966): where ω is specific stream power [W m −2 ], Ω is total stream power [W m −1 ], γ is the specific weight of water [9,810 N m −3 ], Q is discharge [m 3 s −1 ], S is the friction slope [m m −1 ], and w is the water surface width [m]. Specific stream power is a useful variable because it is readily calculated throughout a stream network but still represents the physical processes in rivers. Because of this, it has been used to determine erosion and deposition potential (Parker et al., 2015;Vocal Ferencevic and Ashmore, 2012;Bizzi and Lerner , 2015;Soar et al., 2017), explain dominant modes of channel adjustment (Knighton, 1999;Bull , 1979), model sediment transport processes (Bagnold , 1977(Bagnold , , 1980Martin and Church, 2000;Eaton and Church, 2011), and explain historic variability and future evolution of rivers (Fryirs et al., 2012). Discharge data are typically available from gaging stations, regional regression equations, or hydrologic modeling. Channel slope and width can be obtained from high resolution digital elevation models, often created from airborne LiDAR data.

Discharge
REM is driven by a user-supplied flow record with a given time step (e.g. daily, hourly, or 15-minute). To account for overbank flooding, the model uses the Manning equation to partition flow between the channel and floodplain using the sub-area method, similar to the approach used by HEC-RAS and others (e.g. Soar et al., 2017). The channel and two floodplains are treated as separate sections (j), each with their own Manning roughness coefficient (n j ). The discharge for each section is calculated using trial values of water surface elevation. This processes is repeated until the sum of these discharges equals the known total flow: where A j is the section area [m 2 ], R j is the section hydraulic radius [m], and S is the channel slope. Only the discharge within the channel, and the corresponding flow width, are used to calculate specific stream power.

Channel incision
The model simulates incision into non-cohesive and cohesive bed material, including a mix of both bed types as described below.

Non-cohesive incision
Fundamentally, the model uses the Exner equation to simulate bed elevation changes based on a sediment mass balance: where η is the bed elevation [m], λ is the bed porosity (assumed to be 0.4), w avg is the average bottom width of adjacent cross sections, Q b is the volumetric sediment transport rate [m 3 s −1 ], and t and x are time and downstream distance, respectively. REM models sediment transport by grain size and tracks changes in bed material composition: where F k is the bed surface fraction of the k th grain size, L a is the active layer thickness [m], Q bk is the volumetric sediment transport rate of the k th grain size [m 3 s −1 ], and f lk is the interface exchange fraction which depends on whether the bed is degrading or aggrading: where f k is the bed subsurface fraction, p bk is the bedload fraction of the k th grain size, and α is a weighting parameter than ranges from 0 -1 (we assume α = 0.5). The model does not store bed stratigraphy, meaning information on buried sediment size is lost if the channel aggrades and then incises.
The active layer thickness L a is calculated as three times the surface layer D 90 . Sediment fluxes are discretized using the κ scheme with flux limiters (Hirsch, 2007): where Q be is the volumetric sediment flux out of the control volume centered on the i th cross section [m 3 s −1 ], Q bi−1 and Q bi are the volumetric sediment fluxes at the i th −1 (upstream) and i th cross sections [m 3 s −1 ], and κ is a constant that controls the discretization scheme. We use second order upwinding (κ = −1; (Hirsch, 2007)). r i is defined as: Finally, REM uses the Superbee limiter function (ψ): The model uses two new stream power based equations (Bagnold , 1980) for calculating bedload and total load sediment transport capacity (Lammers and Bledsoe, 2018): where q b is the mass sediment transport rate per unit width [kg m −1 s −1 ], Q t is the total load [ppm], q is unit discharge [m 2 s −1 ], D s is the grain size [m], ω is specific stream power [W m −2 ], and ω c is the critical specific stream power for incipient motion [W m −2 ]. This value is calculated for each grain size using a stream power based hiding function: where ω ri * is the reference dimensionless specific stream power of the i th grain size, ω r50 * is the reference dimensionless specific stream power of the median grain size, and b is an empirical exponent that varies from 0 (size independent mobilization) to 1.5 (equal threshold mobilization). Stream power is made dimensionless by: where ρ is water density [1,000 kg m −3 ], g is gravity [9.81 m s −2 ], and s is sediment specific gravity (usually 2.65). Details of this hiding function are described in Supplementary Material.

Cohesive incision
The model uses a simple excess shear stress equation to model cohesive bed erosion (Partheniades, 1965): where E is the erosion distance [m], k is the erodibility coefficient [m 3 N −1 s −1 ], ∆t is the time step [s], τ is the applied shear stress [Pa], and τ c is the critical shear stress of the bed material [Pa]. The erodibility coefficient can be supplied by the user, or it is calculated based on an empirical relationship developed by Simon et al. (2010) after work by Hanson and Simon (2001): Equation 13 calculates erosion using excess shear stress, but this model is based on a stream power approach. Since data on τ c of various soils are widely available in the literature, and there is no work that we are aware of defining critical stream power of cohesive material, we chose to use an empirical equation to calculate average bed shear stresses directly from ω (see Supplementary Material for more details): This estimated value of τ is then used to calculate cohesive erosion rates (Equation 13).

Mixed non-cohesive/cohesive incision
In streams with both non-cohesive and cohesive bed material, modeling bed elevation changes is more complicated. Sand and gravel can be deposited on top of cohesive material and transport capacity may not be representative of actual sediment movement if the stream is supply limited (e.g. no alluvium on the bed). To account for these processes, REM calculates the actual volume of sediment transported out of a cross section as the minimum of the transport capacity (Q be , Equation 6) and the sediment available for transport (sum of the incoming sediment from upstream and bank erosion and of the available noncohesive alluvium on the channel bed). The available non-cohesive sediment is calculated as: where Q bk,avail is the volume of bed sediment of the k th size class available for transport, converted to a rate [m 3 s −1 ], η cohesive is the elevation of the cohesive layer [m], ∆x is the distance to the next cross section [m], and ∆t is the time step [s]. If η cohesive = η or F k or f k = 0, there is no available bed sediment.

Bank Erosion
The model simulates two fundamental bank erosion mechanisms: fluvial erosion and mass wasting. Bank erosion is calculated at the discharge time step (e.g. daily, hourly, 15-minute, etc.), independent of the time step for bed aggradation and degradation.

Fluvial erosion
Fluvial erosion is the removal of bank soil by flowing water once the resistance threshold of the bank material has been exceeded. Similar to cohesive incision, REM models fluvial bank erosion using an excess shear stress approach (Equation 13). We use an empirical equation to calculate average wall (i.e. bank) shear stress directly from ω (see Supplementary Material for more details): where τ w is the shear stress acting on the channel bank [Pa]. A user specified fraction of the eroded bank material is added to the bed material load (i.e. sand and coarser). The remainder is exported from the watershed as washload.

Mass failure
Planar mass failure is modeled using a modified version of the Bank Stability and Toe Erosion Model (BSTEM) (Simon et al., 2000(Simon et al., , 2011. BSTEM estimates bank stability using a limit equilibrium analysis to calculate a factor of safety -the ratio of resisting to driving forces acting on the bank. The bank is predicted to be stable if the factor of safety is greater than one and unstable if it is less than one. BSTEM accounts for several processes that increase or decrease bank strength, including: (1) water pressure in soil pores (positive pressure decreasing stability and negative pressure increasing stability); (2) confining pressure of the streamflow; and (3) increased soil cohesion from plant roots. Although the simplified version of BSTEM accounts for the first two processes, we exclude vegetation effects since they have a negligible effect on BSTEM output in sensitivity analyses (Lammers et al., 2017) and increase computation time and data requirements. This gives the following factor of safety equation:  (Lammers et al., 2017)). BSTEM uses a horizontal layer method to calculate a net factor of safety for the bank, accounting for different soil layers. The simplified version follows this same approach but uses simplified bank geometry (see Section 2.7), assumes only two soil layers (main bank and toe), and assumes the failure plane intersects the bottom of the bank or top of the bank toe. For a more detailed description of BSTEM see Midgley et al. (2012); Daly et al. (2015a); Simon et al. (2000Simon et al. ( , 2011.

Coupled bank erosion modeling
Fluvial erosion and mass failure are linked processes. Fluvial erosion is typically higher at the bank toe, which steepens the bank and makes it more susceptible to failure. After a bank fails, the collapsed soil is often deposited at the base of the bank toe, temporarily protecting the bank from fluvial erosion (Thorne, 1982).
We account for these dynamic and coupled processes in two ways. First, fluvial erosion is assumed to be a maximum at the base of the toe. This node is eroded the most, with zero erosion at the top, creating a steeper toe angle. If the new toe angle exceeds 90 • (e.g. an undercut bank), the overhanging bank immediately collapses, and the bank geometry is updated accordingly. This bank steepening, coupled with bank heightening from bed erosion, increases the chance of mass failure. If the bank fails, the collapsed soil block is deposited at the bank toe -narrowing the channel -and the toe angle is reduced to conserve the mass of the failed block. If the failed block is too large to fit at the base of the toe, any extra bank material is stored in a "tank". No further fluvial erosion is allowed until the material in this "tank" is eroded (Lai et al., 2015). See Supplementary Material for more details.

Meandering
In addition to incising, meandering channels can also reduce their slope via lateral migration. REM incorporates this process by simulating meander migration from fluvial erosion, allowing the channel to increase its length, thereby decreasing its slope.
The effects of curvature on shear stress distributions can be simulated by directly modeling flow mechanics, typically using a high resolution 1-D or 2-D model (Crosato, 2007;Huang et al., 2014;Darby et al., 2002); however, REM does not directly calculate boundary shear stress distributions, meaning it cannot mechanistically account for the effects of bend geometry on bank erosion. Instead, we use an empirical equation to find the maximum shear stress on the outside of bends (Army Corps of Engineers, 1970): where τ max is the maximum bend shear stress [Pa], τ w is the wall shear stress calculated using Equation 17 [Pa], R c is the bend radius of curvature [m], and w is the channel bottom width [m].
Equation 19 is based on only five small flume datasets, and more recent analysis suggests that no single relationship adequately predicts maximum shear stress in bends (Thornton et al., 2012). Field studies, however, show that radius of curvature is a major control on channel migration rate Hickin, 1983, 1986;Hooke, 1997). We therefore used Equation 19imperfect as it may be -to account for this process.
Including meander dynamics in the model requires two user inputs for each reach. Radius of curvature and sinuosity are used to build and track changes in channel planform. We conceptualize the channel as a series of circular arc segments, where each arc is one bend. The number of bends between each cross section can be calculated from the user defined cross section spacing, radius of curvature, and sinuosity using equations describing circular arcs (see Supplementary Material for more details).

Knickpoint migration
Knickpoints or headcuts are small waterfalls or locally steep stream sections where bed erosion is especially pronounced. These vertical drops tend to migrate upstream as they erode, and can initiate substantial bank erosion (Schumm et al., 1984). We use a simple, empirical model to simulate headcut advance (Allen et al., 2018): where hc m is the headcut migration distance [m], Q cum is cumulative daily discharge [m 3 ], H hc is headcut height [m], and Ehc is an erodibility resistance parameter that is a function of soil erodibility and vegetation cover: and RCF is a root cover density factor (dimensionless, 0 -1.4). While channel beds are usually unvegetated, using RCF = 0 sometimes requires a negative K d value to accurately predict knickpoint migration rates; therefore, REM assumes RCF = 1.4 and requires users to calibrate K d to match observed migration rates (see Supplementary Material for more details). This sub-routine requires the user to input the location, elevation, height, and K d of each knickpoint. The position of each knickpoint is tracked as it migrates upstream (including into any tributaries) and bed elevations are adjusted accordingly.

Cross section geometry
REM assumes a prismatic channel, based on user-supplied bottom width, bank and toe heights and angles, and floodplain width and slope ( Figure 1). All channel geometry variables are unique for the right and left banks. Bank soil parameters (e.g. cohesion) can be distinct for the bank toe and upper bank soil but are the same for the right and left banks in a reach. For each cross section, a cohesive layer may be placed some distance below the channel bed. Aggradation and degradation only occur across the flat channel bottom.
Figure 1: Schematic of cross section (a) and network (b) geometry included in REM. w = width; H = height, α = angle, Q = water discharge, and Q s = sediment discharge.

Network structure and sediment routing
The model uses a simple reach-node network structure, where a series of channel reaches are connected by nodes ( Figure 1) (Schmitt et al., 2016;Czuba and Foufoula-Georgiou, 2015). The user specifies inputs individually by reach, and each reach may have multiple cross sections. Incoming bed material load to each cross section is the sum of sediment supplied by the upstream cross section (or cross sections at tributary junctions), sediment from local bank erosion, and any user-inputted upland sediment supply. Upland sediment and bed material load from eroded banks are assumed to be the same grain size distribution as the initial bed grain size distribution for that reach. The washload component of any bank, cohesive bed, or knickpoint erosion is immediately routed to the watershed outlet. The effects of grade controls or bank armoring can be incorporated by placing non-erodible cross sections within the channel network (i.e. cohesive soils with high τ c ). A table of required and optional model inputs is included in the Supplementary Material.

Generic model test
We simulated channel evolution in a generic watershed with six distinct reaches. The total channel length of 10.4 km corresponds to an approximate drainage area of 6.5 km 2 (Hack (1957, Eq. 3)). Initial grain size (2 mm), slope (0.003), and bank height (2 m) were constant throughout the watershed. Discharge was steady at a station but increased moving downstream. Upstream sediment supply was equal to the transport capacity of the undisturbed channel. A full table of model inputs is included in Supplementary Material. Beginning with an initially stable channel, we dropped the downstream elevation by 2.5 m, including a 1.5 m tall knickpoint, and modeled 20 years of resulting channel evolution.

Study area
Parker Dam, completed in 1938, is one of several large dams on the lower Colorado River built for water supply and power generation. Like most hydropower dams, Parker Dam altered flows and trapped sediment. The combined effects of these changes caused the Colorado River downstream from the dam to incise while the bed material coarsened (Williams and Wolman, 1984).

Data collection and modeling
Initial longitudinal profiles and grain size data for a 144 km reach downstream of Parker Dam were obtained from two U.S. Bureau of Reclamation reports (U.S. Bureau of Reclamation, 1948Reclamation, , 1950. We used a single grain size distribution for the entire reach. The pre-dam grain size data were all finer than 2 mm; however, later observations included gravel up to 32 mm, presumably unearthed as the channel incised. Since the channel coarsened over time (and REM does not account for bed stratigraphy), we adjusted the initial grain size distribution to include a small amount of coarser material. Average channel widths were calculated from 1938 aerial photographs (Norman et al., 2006) and contemporary satellite imagery (Google Earth Pro, 2017). We ran the model from 1938 -1975 using daily discharge data from USGS gage 09427520. Only bed elevation changes were modeled, bank erosion was not included. We used a cross section spacing of 2,000 m and a time step of 2,400 seconds. The total load sediment transport equation was used for all grain sizes < 4 mm and the bedload equation for all coarser grain sizes. We assumed no sediment inputs from upstream (i.e. the dam trapped all sediment). Model results were compared to measured longitudinal profiles for a 66 km subreach (from 27 -93 km downstream of Parker Dam) (Williams and Wolman, 1984). We also compared modeled D 50 to measured values from three cross sections (26, 64, and 130 km downstream of Parker Dam) (Williams and Wolman, 1984).
In additional to the single model run described above, we ran 5,000 Monte Carlo simulations varying the initial grain size distribution, channel width, floodplain geometry, roughness values, and the exponent and coefficient of the hiding function. Sobol' quasi-random numbers (using the "gsl" R package; (Hankin, 2006)) were used to generate these variables since they provide more uniform coverage than simple random numbers (Sobol' , 1976).  . This massive sediment deposit buried the channel network of the upper NFTR. Over the following months and years, channels reformed from surface runoff, pumping from Spirit Lake, and multiple lahars (volcanic debris or mudflows) .
To prevent sedimentation in the downstream Cowlitz and Columbia Rivers, two sediment retention structures where built on the NFTR. The first (N1) was built in summer 1980 and operated until it was breached in 1982. A second, more permanent sediment retention structure (the "SRS"), was built in 1987 and was essentially filled by 1998 Zheng et al., 2014). To prevent overtopping of Spirit Lake, water was released into a NFTR tributary (see TR065 and TR070 in Figure 2) at a constant rate of 5.1 m 3 s −1 from November 1982 to August 1983, causing extreme incision (up to 34 m) (Paine, 1984). For more details on the eruption and its effects, see Simon et al. (1999), Lipman and Mullineaux (1981), and Major et al. (2018).

Data collection and modeling
We modeled evolution of the upper NFTR and its tributaries from September 1983 -August 2011. We started the model 3.5 years after the eruption because there were more cross section data and this avoided several lahars and pumping from Spirit Lake which had complicated effects on channel adjustment. Following the eruption, the USGS and Army Corps of Engineers established several permanent cross sections which have been surveyed at irregular intervals since 1980. We used these data (Mosbrucker et al., 2015) for 19 cross sections on the NFTR and its tributaries to estimate initial channel and floodplain geometry (Figure 2). Each of these cross sections defined a model reach with unique inputs. Initial bed grain size distributions were estimated from field data (U.S. Army Corps of Engineers, 1988;Paine, 1984). We used the daily discharge series at the SRS constructed by Simon and Klimetz (2012) from several nearby USGS gages. These values were scaled by drainage area to give discharge in each reach. We also used bank sediment properties (τ c , k, cohesion, unit weight, and φ ) and Manning's n values estimated by Simon and Klimetz (2012). We assumed no hillslope sediment supply since upland erosion peaked soon after the eruption and remained negligible compared to in-stream sediment sources . We used a model cross section spacing of 500 m, a time step of 2,400 seconds, and the bedload sediment transport equation. Sediment specific gravity was adjusted to account for lighter volcanic material (Simon and Klimetz ,  (Mosbrucker , 2014). Modeled cross sections are differentiated into "upstream" and "downstream" which will be referenced in certain result figures. Flow is from right to left. 2012, Eq. 24). Finally, we assumed that 100% of the eroded bank material consisted of bed material load. We ran 5,000 Monte Carlo simulations to quantify uncertainty, varying initial grain size, channel width, channel roughness, hiding function parameters, and bank soil properties. Model accuracy was assessed by comparing modeled bed elevations to observations (from survey data and a 1 m DEM from 2009 (Mosbrucker , 2014)). Other parameters (e.g. D 50 and width) were not used because only sparse grain size data were available and the simplified model cross sections could not adequately represent the complex observed channel geometries.

Sensitivity analysis
For both case studies, we performed sensitivity analyses to determine which variables most influence model output. We used a density-based method that estimates parameter sensitivity based on differences between conditional and unconditional probability density functions of model output (Plischke et al., 2013). Variables with a greater effect have bigger differences in these density functions. This method has two advantages over other approaches: it requires no unique input parameter sampling design (e.g. Saltelli et al., 2010) and it requires much fewer model runs (e.g. Pianosi and Wagener , 2015). We there- fore used the output from Monte Carlo simulations to compute the sensitivity indices. Bootstrapping with 1,000 replicates was used to correct for bias and calculate uncertainty in sensitivity indices. Finally, we incorporated a dummy variable to determine the threshold for influential variables. This dummy variable is a simple set of random numbers that has no influence on the model and accounts for noise in the sensitivity analysis (Plischke et al., 2013;Khorashadi Zadeh et al., 2017). These sensitivity analyses are only applicable for each individual case study because each system has unique boundary conditions and relevant processes. Because of this, it is necessary to perform a sensitivity analysis separately for every model application to understand what variables are most influential in each case.
For the Colorado River, we quantified sensitivity for two model outputs: bed elevation and bed D 50 . For the NFTR, channel width was also included. To give a single output value for each model run, we summed the absolute value of the total change in the variable (e.g. bed elevation) for all cross sections. For the NFTR, a separate sensitivity analysis was performed for each reach. For comparison among reaches, we standardized the sensitivity indices by taking the difference between the index for each input and the "dummy" variable, divided by the dummy variable index. All analyses of model outputs were done using R version 3.4.1 (R Core Team, 2018). Figure 3 shows changes in bed elevation, channel width, and width-depth ratio for the modeled test case. The zone of disturbance migrated upstream through time, with changes in channel width lagging slightly behind changes in bed elevation. The greatest channel changes were at the far downstream end -the area with greatest disturbance. Figure 4 shows changes in stream power, bed elevation, and channel width at two locations (indicated in Figure 3(a)). For both areas, stream power was relatively constant until the knickpoint passed, after which stream power spiked before slowly decreasing. Bed elevation and width show similar trends, with abrupt changes following passage of the knickpoint. After the initial drop in channel elevation, both cross sections showed a period of aggradation followed by renewed incision. Sediment export from the watershed peaked early in the simulation and then decreased exponentially. Figure 5 shows the error in predicted bed elevation and bed D 50 . The median of the Monte Carlo simulations generally has lower error than the single model run. For the bed elevation results (Figure 5(a-c)), model error decreases over the course of the simulation, although the uncertainty increases. For the bed D 50 results ( Figure 5(e-g)), uncertainty is high for all sites but error generally decreases moving downstream. Figure 6 shows the results of the sensitivity analysis for bed elevation and bed D 50 outputs. Initial D 50 , geometric standard deviation of the grain size distribution (σ g ), and channel width have the largest influence on predicted bed elevation changes. Initial D 50 and σ g have a significant effect on the final D 50 while channel width and the hiding function parameters (ω c and b) have only a small effect. Floodplain angle has a moderate effect on both outputs. have generally low error in predicted final bed elevations normalized to the magnitude of total bed elevation change (Figure 8). Median normalized error is 43%, but is only 22% for reaches CW280 -NF130. For the remainder of the cross sections, the model did a relatively poor job of predicting changes in bed elevation.

North Fork Toutle River
There is substantial uncertainty for all sites, especially in the upper half of the watershed (e.g. > 20 m wide 90% confidence interval). The magnitude of uncertainty is generally less in the lower portion of the watershed where the magnitude of aggradation and incision was smaller.
The sensitivity results for the NFTR model runs are summarized in Figure  9. Modeled bed elevation was influenced most by bank τ c , bank cohesion, and hiding function parameters (ω c and b). Channel width and initial bed grain size (D 50 and σ g ) also had a minor effect. Modeled D 50 was influenced by similar variables, but the hiding function parameters, initial grain size, and bank cohesion had a much larger effect. For modeled channel width, bank τ c was by far the most influential but initial width and ω c also contributed to some observed model uncertainty.  Figure 3(a)). Sediment discharge at the watershed outlet is also shown (c). All variables are scaled to their starting value (horizontal lines). In (a) and (b), ω increases because slope increases slightly later in the simulation, despite the channel incising.

REM accurately predicts channel change
The generic test case and field applications show that REM can realistically and accurately simulate channel evolution -in the absence of avulsions and extensive lateral migration. First, the model test case matches physical understanding of channel evolution in response to disturbance (in this case, base level drop). The greatest channel change is observed nearest the disturbance, and rates and magnitudes of erosion decline nonlinearly with time and distance upstream (Figure 3). This is consistent with conceptual models of channel evolution (Schumm et al., 1984;Simon, 1989), and experimental (Begin et al., 1981), numerical (Simon and Darby, 1997), and field studies (Simon and Rinaldi, 2006). In general, the channel incises which destabilizes the banks, leading to rapid widening ( Figure 4). As the upstream channel begins to erode, large amounts of sediment are delivered downstream, causing aggradation. After this upstream sediment supply is cut off (i.e. upstream channel erosion has slowed or stopped), channel incision begins again. This shift between degradation and aggradation depending on sediment delivery from upstream is an important control on channel evolution, as demonstrated in both numerical modeling (Simon and Darby, 1997) and field studies (Simon and Hupp, 1992). Downstream aggradation can help stabilize these reaches and allows the channel to more rapidly attain a new stable slope (Doyle and Harbor , 2003). Disrupting this downstream sediment delivery, for example by installing grade control structures, can induce a second round of incision  downstream (Simon and Darby, 2002), similar to what the modeling showed (Figure 4).
Following a disturbance, the channel is expected to adjust rapidly, with the rate of change slowing until the channel reaches some new stable state. This results in an exponential decay in channel variables to some asymptote. These variables may include stream power (Bull , 1979;Bledsoe et al., 2002), sediment discharge (Simon, 1999;Bledsoe et al., 2002), or bed elevation (Begin et al., 1981), but all describe a reduction in the rate of energy dissipation (Simon, 1992). Modeling shows these exponential reductions in specific stream power and sediment discharge, and an exponential increase in channel width ( Figure  4). Bed elevation follows a more complex trajectory, but does decrease towards an asymptote during the second round of incision.
Modeling from the NFTR also shows this exponential decrease (or increase) in bed elevation (Figure 7), consistent with physical understanding of channel evolution. In the Colorado River modeling, the greatest incision and bed coarsening were seen closest to the dam (the disturbance), with less channel    NF405). This is the same division as the vertical line in Figure 8. Cohesion, phi, and weight results are shown for the higher value of either the bank or bank toe. Vertical dashed line is a normalized sensitivity index of zero (i.e. no influence on model output). change downstream (data not shown). Furthermore, REM accurately predicts the magnitude of channel incision in this system (bed elevation RMSE 0.7 -1.5 m for all years). Bed material coarsening is also accurately predicted, although errors are more variable (D 50 RMSE 0.1 -5 mm). In the NFTR, REM accurately predicts channel incision in the upper half of the watershed (CW280 -NF130). This portion of the channel is single thread, while the downstream portion is anastomosing or braided -features that were deliberately not incorporated into REM. Taken together, these three model tests suggest that REM can predict channel evolution across decadal time scales in single-thread systems with reasonable accuracy, matching both physical understanding of channel change and adequately predicting evolution in real-world, dynamic fluvial systems.

Model strengths and weaknesses
REM's main strength is its parsimony and utility in simulating watershed scale channel evolution processes. Watershed scale assessment is essential because channel evolution is not limited to local disturbances or dynamics. Changes in both upstream and downstream channel form and sediment delivery affect local channel response (e.g. Schumm et al., 1984;Simon, 1992;Simon and Darby, 2002). Both bed and bank erosion processes are especially important in smaller urban watersheds (Booth, 1990). Furthermore, channel hardpoints (i.e. bed and bank armoring) can significantly influence local channel evolution and adjustment in other parts of the watershed (Booth and Fischenich, 2015). REM accounts for these processes -enabling users to specify nonerodible cross sections -and may be an important tool for understanding urban channel network evolution. Other numerical models have been developed that include both bed and bank erosion, but these are typically designed for reach-scale application. For example, the CONCEPTS model (Langendoen and Simon, 2008;Langendoen and Alonso, 2008) and Darby and Thorne (1996a) model both include more detailed modeling than REM, but cannot be easily applied at the watershed scale. Alternatively, the watershed scale Soil and Water Assessment Tool (SWAT) (Allen et al., 1999;Mittelstet et al., 2016;Arnold et al., 1998) has erosion processes for cohesive channels; however, REM incorporates cohesive and non-cohesive erosion and bank failure. REM includes the most important mechanisms to realistically simulate channel evolution while still keeping data requirements to a minimum.
Another important strength of REM is its capacity to explicitly account for input variable uncertainty. It automates the use of Monte Carlo simulations, allowing users to easily quantify model uncertainty and produce probabilistic estimates of channel change. Quantifying uncertainty can be useful for decision making and assessing reliability of model outputs (e.g. Pappenberger and Beven, 2006). Model field tests illustrate this. In most cases, it appears the median of the Monte Carlo simulations predicts river behavior as well or better than the single model run (with the exception of NF130 and NF300 from the NFTR, Figure 7). This suggests that accounting for uncertainty in the inputs can actually improve model accuracy.
How much uncertainty is too much must be determined by the model user because it depends on the question(s) being asked. The model test cases show large uncertainty bounds. This may seem discouraging, but is an inescapable consequence of simulating complex and uncertain geomorphic systems (Shreve, 1975). By quantifying this uncertainty, we can at least be candid about confidence in the model's predictions. The widths of the simulated uncertainty bounds are proportional to the magnitude of modeled bed elevation ( Figures  5 and 7) and grain size ( Figure 5). This is expected -the larger the change, the greater uncertainty.
REM is only applicable for single-thread rivers. It is therefore unsurprising that it could not adequately predict channel evolution in the downstream half of the NFTR. This section of the river migrates across a wide valley bottom and -in the lower reaches -the channel braids (Zheng et al., 2017). In reality, much of the channel is 15 -20 m wide, but may be within a several hundred meter wide valley. The model cannot simulate the aggressive channel migration observed in the lower portion of the watershed and instead spreads the water out over an unrealistically wide modeled channel bottom. Figure  8 illustrates this issue, showing how error in modeled bed elevation increases substantially where the valley widens (just downstream of NF130). REM does include a meandering algorithm, but this is not entirely mechanistic and is incorporated to allow single thread meandering channels an additional mode of slope adjustment.
Other limitations are a consequence of REM's relative simplicity. The model assumes uniform flow (S o = S f ) to calculate specific stream power and relies on new empirical equations to convert stream power to shear stress for cohesive erosion modeling. This facilitates network scale analysis without detailed hydraulic modeling but may be a source of error. This also neglects local, complex flow hydraulics which can have an impact on channel changemaking it unsuitable for small scale analyses, like bridge scour. Still, REM has a strong physical basis, integrating novel stream power based sediment transport models (Lammers and Bledsoe, 2018) with a well tested bank erosion algorithm (BSTEM; (Simon et al., 2000(Simon et al., , 2011) that underwent systematic sensitivity and uncertainty analyses to identify the most parsimonious representation of essential physical processes (Lammers et al., 2017).

Model sensitivity
Sensitivity is a function of (1) how much an input influences model output and (2) how much the input varies. Sensitivity analyses can therefore reveal information about model structure and suggest which variables should be most accurately quantified to obtain the most reliable results. REM sensitivity analyses largely confirm the validity of the model as important parameters are known to be linked to important channel evolution processes and are consistent with results reported in the literature.
Bed elevation is most controlled by D 50 , σ g , width, and floodplain angle (Colorado River, Figure 6) plus hiding function parameters and bank τ c and cohesion (NFTR, Figure 9). The size and erodibility of the bed material directly influences the extent of incision. Bank erodibility has a secondary effect by either allowing the channel to widen and reducing incision, or limiting widening and forcing the channel to incise more (Simon, 1992). Other numerical models have shown that bed D 50 has a significant effect on modeled channel profiles (El Kadi Abderrezzak et al., 2008;El Kadi Abderrezzak and Paquier , 2009); however, Darby and Thorne (1996b) found that D 50 had a minimal effect compared to discharge.
Predicted bed D 50 was most influenced by initial grain size distribution in the Colorado River case study ( Figure 6), but hiding function parameters were equally or more influential for the NFTR (Figure 9). Others have also shown that hiding function parameters (in their case, critical shear stress and the hiding factor) control modeled grain sizes (Ruark et al., 2011;Hoey and Ferguson, 1994). The NFTR results also show that bank τ c and cohesion had an influence on modeled D 50 . Sediment from bank erosion has the same grain size distribution of the initial bed sediment. As the bed coarsens, bank erosion therefore becomes a source of finer grains.
Channel width was controlled most by bank τ c (Figure 9). This suggests that fluvial erosion, not mass failure, was the dominant bank erosion process. Darby and Thorne (1996b) also found that τ c had a much greater influence on channel widening than bank cohesion. The three variables controlling bank failure (cohesion, φ , and weight) all had similar relative importance, unlike other sensitivity analyses of bank erosion models that found cohesion was the dominant control on bank stability (Lammers et al., 2017;Van de Wiel and Darby, 2007;Parker et al., 2008;Samadi et al., 2009). These studies also did not show that τ c was important, possibly because they did not model cumulative bank erosion and therefore did not incorporate the threshold effect of τ c determining when erosion occurs.
Despite its relative simplicity, REM is dependent on field data which may be difficult to collect at a network scale; however, the sensitivity results provide guidance on which variables should be most accurately quantified to yield the best model results. This is especially important for bank τ c which has a strong influence on the model, is subject to considerable uncertainty, and is difficult to measure in the field (Wynn et al., 2008;Konsoer et al., 2016;Daly et al., 2015b). Bank τ c may need to be estimated through model calibration to provide more reliable model inputs than field measurements.

Future improvements and applications
There are a number of modifications that could improve model predictions.
Coupling REM with an upland erosion model would provide more realistic estimates of sediment inputs and channel response (e.g. Stryker et al., 2017). Furthermore, floodplains can be significant sediment sinks (Kronvang et al., 2007;Fryirs and Brierley, 2001); although, floodplain sedimentation likely has a larger effect on fine sediment delivery (e.g. Walling et al., 1998) than the bed material load that controls channel incision and aggradation. Adding these processes may improve model predictions, but this extra complexity also increases data requirements and uncertainty. It is important to balance the need to incorporate relevant processes while retaining the simplicity that makes REM applicable at the watershed scale.
REM has a number of potential applications, both in river management and research. For example, channel erosion can be a significant -but difficult to quantify -source of fine sediment and phosphorus pollution in watersheds (Fox et al., 2016). REM could be used to estimate loading of these pollutants at watershed scales. Urban stormwater management (or mismanagement) is a leading cause of channel degradation (Walsh et al., 2016(Walsh et al., , 2005. While certain stormwater design standards can help mitigate channel degradation (e.g. Tillinghast et al., 2011), REM may allow a more comprehensive analysis of channel stability when coupled with a stormwater management model. REM also has a number of research applications. The search for an "optimal" or "equilibrium" channel form has intrigued scientists for decades (e.g. Langbein and Leopold , 1964;Yang et al., 1981;Millar , 2005;Huang et al., 2014). Tools like REM can be used to explore this concept in more detail, looking beyond the "optimal" channel cross section and examining interactions between parts of a network and their influence on watershed scale channel evolution.

Conclusions
We present a new model for simulating channel evolution at the watershed scale. This model is based on specific stream power and does not require detailed hydraulic modeling. Results from a generic test case of channel response to base level lowering match physical understanding of channel evolution. The model also accurately predicts channel incision and bed coarsening for a reach of the lower Colorado River below Parker Dam. In the North Fork Toutle River, the model accurately predicted channel incision and widening in the upper portion of the watershed where the channel remained single thread. Model predictions were poor in the lower watershed where the river migrated significantly across the valley floor -a behavior that REM is not designed to simulate. Results from these case studies suggest the model can provide useful predictions of watershed-scale channel erosion, while recognizing it is limited to single thread channels. Importantly, the model can also account for uncertainty in input variables -allowing for a probabilistic assessment of channel change. More model testing is required to fully understand its capabilities and limitations. For example, REM's ability to simulate cohesive incision, knickpoint migration, or meandering was not tested because of a lack of sufficient field data. Further testing is also warranted on the smaller watersheds (i.e. 10 -100 km 2 ) for which REM was designed.
Understanding how and how much rivers may change under future climate and land use variability is an essential question for sustainable river management. Other tools have been developed to estimate watershed sediment dynamics (Czuba et al., 2017;Schmitt et al., 2016;Czuba and Foufoula-Georgiou, 2014) and erosion and deposition potential (Soar et al., 2017;Parker et al., 2015). In smaller, urbanizing watersheds, however, channel changes are driven by both bed and bank erosion processes (Booth, 1990) and strongly influenced by channel armoring and other channel "improvements" (Booth and Fischenich, 2015). By accounting for these processes, REM can provide insight into urban stream evolution. Additionally, the model can be used to test different mitigation strategies; for example, by simulating how the river erodes under different stormwater and/or stream restoration scenarios to support cost effective and successful solutions to address excessive channel erosion.

Acknowledgments
This work was partially funded by the National Science Foundation, Integrative Graduate Education and Research Traineeship (IGERT) [Grant No. DBE-0966346] 'I-WATER: Integrated Water, Atmosphere, Ecosystems Education and Research Program' at Colorado State University. Additional funding was provided by the United States Environmental Protection Agency (USEPA) [grant RD835570]. Its contents are solely the responsibility of the grantee and do not necessarily represent the official view of the USEPA. Further, USEPA does not endorse the purchase of any commercial products or services mentioned in the publication. We are grateful to Peter Nelson for providing guidance on model development and Mazdak Arabi and Sara Rathburn for constructive comments on the manuscript. Additional information is available in "Supplementary Material". REM code and all model inputs and outputs from this paper are available at www.github.com/rodlammers/REM. Darby, S. E., and C. R. Thorne (1996a), Numerical simulation of widening and bed deformation of straight sand-bed rivers. I: Model development, Journal of Hydraulic Engineering, 122 (4), 184-193.