mizuRoute version 1 : a river network routing tool for a continental domain water resources applications

This paper describes the first version of a standalone runoff routing tool, mizuRoute. The mizuRoute tool post-processes runoff outputs from any distributed hydrologic model or land surface model to produce spatially distributed streamflow at various spatial scales from headwater basins to continental-wide river systems. The tool can utilize both traditional grid-based river network and vector-based river network data. Both types of river network include river segment lines and the associated drainage basin polygons, but the vector-based river network can represent finer-scale river lines than the grid-based network. Streamflow estimates at any desired location in the river network can be easily extracted from the output of mizuRoute. The routing process is simulated as two separate steps. First, hillslope routing is performed with a gamma-distribution-based unit-hydrograph to transport runoff from a hillslope to a catchment outlet. The second step is river channel routing, which is performed with one of two routing scheme options: (1) a kinematic wave tracking (KWT) routing procedure; and (2) an impulse response function – unit-hydrograph (IRF-UH) routing procedure. The mizuRoute tool also includes scripts (python, NetCDF operators) to pre-process spatial river network data. This paper demonstrates mizuRoute’s capabilities to produce spatially distributed streamflow simulations based on river networks from the United States Geological Survey (USGS) Geospatial Fabric (GF) data set in which over 54 000 river segments and their contributing areas are mapped across the contiguous United States (CONUS). A brief analysis of model parameter sensitivity is also provided. The mizuRoute tool can assist model-based water resources assessments including studies of the impacts of climate change on streamflow.


Introduction
The routing tool described in this paper post-processes runoff outputs from macro-scale hydrologic models or land surface models (hereafter we use "hydrologic model" to refer to both types of model) to estimate spatially distributed streamflow along the river network.The river routing tool is named mizuRoute (mizu means "water" in Japanese).The motivation for the development of mizuRoute was to enable continental domain evaluations of hydrologic simulations for water resources assessments, such as studies of the impacts of climate change on streamflow.The mizuRoute tool is suitable for processing ensembles of multi-decadal runoff outputs because the tool is stand-alone and easily applied in a parallel mode.The mizuRoute tool is also designed to output streamflow estimates at all river segments in the river net-work across the domain of interest at each time step, facilitating further spatial and temporal analysis of the estimated streamflow.As opposed to other routing models, our goal in developing mizuRoute is to provide flexibility in making routing model decisions (i.e., river network definition and routing scheme).
The paper proceeds as follows.Section 2 reviews existing river routing models.Section 3 describes how the mizuRoute tool provides flexibility of routing modeling decisions and details the hillslope and river routing schemes used in mizuRoute.Section 4 provides an overview of the workflow of mizuRoute from preprocessing hydrologic model output to simulating streamflow in the river network.Section 5 demonstrates streamflow simulations in river systems over the contiguous United States (CONUS).Finally, a summary and future work are discussed in Sect.6.

Existing river routing models
The water resources and Earth system modeling communities have developed a wide spectrum of river routing schemes of varying complexity (Clark et al., 2015).For example, the US Army Corps of Engineers (USACE) has developed a stand-alone river modeling system called Hydrologic Engineering Center -River Analysis System (HEC-RAS; Brunner, 2001).HEC-RAS offers various hydraulic routing schemes, ranging from simple uniform flow to onedimensional (1-D) Saint-Venant equations for unsteady flow.HEC-RAS has been popular among civil engineers for river channel design and floodplain analysis where surveyed river geometry and physical channel properties are available.At the continental to global scale, unit-hydrograph approaches have been used (e.g., Nijssen et al., 1997;Lohmann et al., 1998;Goteti et al., 2008;Zaitchik et al., 2010;Xia et al., 2012), though more recent, large-scale river models use fully dynamic flow equations (e.g., Miguez-Macho and Fan, 2012;Paiva et al., 2013;Clark et al., 2015), simplified Saint-Venant equations such as the kinematic wave or diffusive wave equation (e.g., Arora and Boer, 1999;Lucas-Picher et al., 2003;Koren et al., 2004;Yamazaki et al., 2011;Li et al., 2013;Yamazaki et al., 2013;Gochis, 2015;Yucel et al., 2015) or non-dynamical hydrologic routing methods such as Muskingum routing (e.g., David et al., 2011).Despite their computational cost, dynamic or diffusive wave models are attractive for relatively flat floodplain regions such as along the Amazon River where backwater effects on the flood wave are significant (Paiva et al., 2011;Yamazaki et al., 2011;Miguez-Macho and Fan, 2012).At the other end of the spectrum, simpler, non-dynamic routing schemes, such as the unit-hydrograph approach, estimate the flood wave delay and attenuation, but do not simulate other streamflow variables such as flow velocity and flow depth.
One of the key issues for large-scale river routing, besides the choice of the routing scheme, is the degree of abstrac-tion in the representation of the river network (Fig. 1).A vector-based representation of the river network refers to a collection of hydrologic response units (HRUs or spatially discretized areas defined in the model) that are delineated based on topography or catchment boundary.River segments in the vector-based river network, represented by lines, meander through HRUs and connect upstream with downstream HRUs.On the other hand, in the grid-based river network, the HRU is defined by a grid box and river segments connect neighboring grid boxes based on the flow directions.Vector-based river networks are better than coarser resolution (e.g., > 1 km) gridded river networks at preserving fine-scale features of the river system such as tortuosity and drainage area, therefore representing more accurate sub-catchment areas and river segment lengths.
For large-scale applications, many studies have developed and evaluated methods to upscale fine-resolution flow direction grids (∼ 1 km or less) to a coarser resolution (∼ 10 km or more) to match hydrologic model resolution and/or reduce the cost of routing computations (e.g., O'Donnell et al., 1999;Fekete et al., 2001;Olivera et al., 2002;Reed, 2003;Davies and Bell, 2009;Wu et al., 2011Wu et al., , 2012)).Earlier work (e.g., O'Donnell et al., 1999;Fekete et al., 2001;Olivera et al., 2002) focused on preserving the accuracy of the flow direction at the coarser resolution and therefore on an accurate representation of the drainage area.Newer upscaling methods are designed to also preserve fine-scale flow path length (e.g., Yamazaki et al., 2009;Wu et al., 2011).More recent river routing models have also begun to employ vector-based river networks (Goteti et al., 2008;David et al., 2011;Paiva et al., 2011;Lehner and Grill, 2013;Paiva et al., 2013;Yamazaki et al., 2013).

Runoff routing in mizuRoute
The runoff routing in mizuRoute provides more flexibility in continental domain routing applications.The mizuRoute tool enables model flexibility in two ways: first, mizuRoute can be used to simulate streamflow for both grid-and vectorbased river networks.Given either type of river network data, mizuRoute offers an option to route flow along all the river segments in the river network data or route runoff at an outlet segment specified by a user.With the latter option, routing computations are performed only upstream of the specified outlet, which reduces the computational cost.Second, the modular structure of the mizuRoute tool offers the flexibility to configure multiple routing schemes.The current version of mizuRoute includes two different river routing schemes: (1) kinematic wave tracking (KWT) routing and (2) impulse response function -unit-hydrograph (IRF-UH) routing, mimicking the Lohmann et al. (1996) model.This flexibility offers new capabilities not present in existing routing models.One capability is to provide an opportunity to explore routing model uncertainties originating from the repre-  (∼ 12 km) gridded river network and vector river network from United States Geological Survey (USGS) Geospatial Fabric for the upper part of Snake River basin (Viger, 2014).sentation of the river system and routing scheme differences (equations and parameters) separately.
The mizuRoute tool uses a two-step process to route basin runoff.First, basin runoff is routed from each hillslope to the river channel using a gamma-distribution-based unit hydrograph.This allows for the representation of ephemeral channels or channels too small to be included in the river network.Second, using one of the two channel routing schemes, the delayed flow from each HRU is routed to downstream river segments along the river network.The routing time step is the same as that of the runoff output from the hydrologic model, typically an hourly or daily time step.The following sub-sections provide descriptions of the two routing steps.

Hillslope routing
Hillslope routing accounts for the time of concentration (T c ) of a HRU to estimate temporally delayed runoff (or discharge) at the outlet of the HRU from runoff computed by a hydrologic model.
For hillslope routing, mizuRoute uses a simple twoparameter Gamma distribution as a unit-hydrograph to route instantaneous runoff from a hydrologic model to an outlet of HRU.The Gamma distribution is expressed as where t is time [T ], a is a shape parameter [-] (a > 0), and θ is a timescale parameter [T ].Both the shape and timescale parameters affect the peak time (mode of the distribution: (a − 1)θ ) and flashiness (variance of the distribution: aθ 2 ) of the unit-hydrograph and depend on the physical HRU characteristics.Convolution of the gamma distribution with the runoff depth series is used to compute the fraction of runoff at the current time, which is discharged to its corresponding river segment at each future time as follows: where q is delayed runoff or discharge [L 3 T −1 ] at time step t [T], R is HRU total runoff depth [L 3 T −1 ] from the hydrologic model, and tmax is the maximum time length for the gamma distribution [T ].This two-parameter gamma distribution has been widely used in unit-hydrograph-based models for water resources engineering applications (e.g., Bhunya et al., 2007;Nadarajah, 2007).Kumar et al. (2007) presented methods to estimate the two parameters in the gamma distribution based on geomorphological information.The gamma distribution offers a parsimonious way to describe a wide range of hillslope-to-channel responses in a computationally efficient manner, which is important for continental-scale domains.

River channel routing
Two different river channel routing schemes are implemented in mizuRoute: (1) KWT routing and (2) IRF-UH routing.Both schemes are based on the one-dimensional (1-D) Saint-Venant equations that describe flood wave propagation through a river channel.The 1-D conservation equations for continuity (Eq. 3) and momentum (Eq.4) are where q is discharge [L 3 T −1 ] at time step t [T] and location , and g is gravitational constant [LT −2 ].The continuity equation (Eq. 3) assumes that no lateral flow is added to a channel segment.The following sub-sections describe the two routing schemes.

Kinematic wave tracking (KWT)
In contrast with several other kinematic routing models that solve a kinematic wave equation with numerical schemes (e.g., Arora and Boer, 1999;Lucas-Picher et al., 2003;Koren et al., 2004), the KWT method computes a wave speed or a celerity for the runoff (or discharge) that enters an individual stream segment from the corresponding HRU at each time step using kinematic approximation (Goring, 1994;Clark et al., 2008).The runoff, represented as a particle, is propagated through the river network based on a travel time (the celerity divided by the segment length).Note that the wave celerity differs from the flow velocity, as the wave typically moves faster than water mass (McDonnell and Beven, 2014).
In the kinematic wave approximation with the assumption that the channel is rectangular and hydraulically wide (channel width y), the wave celerity and discharge q [L 3 T −1 ].Further details are provided in Appendix A. Among the four variables, the channel slope S 0 is provided by the river network data and discharge is computed with hillslope routing for the headwater basin, or/and updated via routing from the upstream segment.The other two variables, Manning's coefficient n and river width w, are much more difficult to measure or estimate.The river width is determined with the following width-drainage area relationship (Booker, 2010): where W a is a width factor [-], A ups is the total upstream basin area [L 2 ], and b is an empirical exponent equal to 0.5.The width factor W a and the Manning's coefficient n are treated as model parameters as shown in Table 1.
The KWT routing starts with ordering all the segments in the processing sequence from upstream to downstream.The KWT routing is performed at each segment in the processing order at each time step.The procedures of the KWT routing method are as follows: 1. Obtain the information on the waves that reside in the segment at a given time step: This includes the waves routed from the upstream segments, the wave that remains in this current segment form the previous time step, and the wave generated from the runoff from local HRUs during the current time step.Three state variables of the waves are kept in the memory: (i) discharge; (ii) the time at which the wave enters the segment; and (iii) the time at which the wave is expected to exit the segment.At the first time step, only the wave from local HRUs exists.Figure 2a visualizes the discharge of waves that reside in the 16 segments (16 river segments are shown in the inserted map) at the beginning of five time steps against the wave locations in the segments.
2. Remove waves in order to reduce the memory usage and the processing time for the wave routing (the next step).
The number of the waves in the segment is limited to a predefined number (20 by default).In Fig. 2, the threshold for the number of waves in each segment is set to 100.To determine which waves can be removed, the difference between the discharge of the wave and the linearly interpolated discharge between its two neighboring waves is computed for all waves, and the wave that produces the least difference (from the interpolated discharge) is removed so that loss of wave mass is minimized.This process is repeated until the number of waves is below the threshold.
3. Route waves through a given river segment.In the routing routine, the celerity of each wave in the segment is computed with Eq. (A6) and the time at which each Figure 3. Overview of streamflow simulation with mizuRoute.The green cylinder and blue box denote data and computational process, respectively.
wave is expected to exit the river segment is updated.
If the exit time occurs before the end of the time step, the wave is propagated to the downstream segment and flagged as "exited".The exit time then becomes the time the wave entered the downstream segment.Otherwise, the wave is flagged as "not-exited", and remains in the current segment.Figure 2b shows the discharge of the waves against the exit times of the corresponding waves at segments 4 and 13.As a reference, the end of each time step is shown as a vertical line.
The routing routine checks for (and corrects) the special case of a kinematic shock.A kinematic shock is a sudden rise in the flow depth, and thus an increase in the discharge at a fixed location, and occurs when a faster-moving wave successively overtakes multiple slower-waves to build a steep wave front.It occurs in models due to the kinematic approximation; in reality, diffusion would act to reduce the steepness of the wave front.Two neighboring waves are evaluated to check if a slower wave is overtaken by a faster wave before the waves exit the river segment.If this occurs, those two waves are merged into one, and the celerity of the merged wave is updated with the following equation: where C merge is the merged wave celerity [LT −1 ], q and A are differences in discharge [L 3 T −1 ] and cross sectional flow area [L 2 T −1 ], respectively, between slower and faster waves.Note that Eq. ( 6) is the mathematical definition of the wave celerity.Since we assume a rectangular channel, whose width is constant for each segment, the merged celerity C merge is a function of flow depth y, which is computed with Eq. (A4).
4. Finally, the time-step-averaged discharge (streamflow) is computed by temporal integration of the discharge of all the waves that exit the segment during the time step.Temporal integral of wave discharge is visualized in Fig. 2b) as the area enclosed by the discharge curve formed by all the waves that exit during the time step.

Impulse response function -unit hydrograph (IRF-UH)
The IRF-UH method mimics the river routing model of Lohmann et al. (1996), which has been used to route flows from gridded land surface models such as the Variable Infiltration Capacity model (VIC; Liang et al., 1994).The only difference between the current tool and the Lohmann routing tool is the way in which the river network is defined.The Lohmann routing model is designed as a grid-based model as shown in Fig. 1 to ease the coupling with grid-based land surface models.In mizuRoute, the same IRF-UH method can be used either on a vector-or grid-based river network.The descriptions of IRF-UH are given briefly as follows.
The mathematical developments of IRF-UH are based on 1-D diffusive wave equation derived from the 1-D Saint-Venant equations (Eqs.3 and 4): www.where parameters C and D are wave celerity [LT −1 ] and diffusivity [L 2 T −1 ], respectively.The complete derivation from Eqs. ( 3) and (4) to Eq. ( 7) is given in Appendix B. Equation ( 7) can be solved using convolution integrals where and U (t − s) is a unit depth of runoff generated at time t − s.This solution is a mathematical representation of the IRF used in unit-hydrograph theory.Wave celerity C and diffusivity D are treated as input parameters for this tool (Table 1), and ideally they can be estimated from observations of discharge and channel geometries at gauge locations.Given a river segment or outlet segment, a set of unique unit-hydrographs is constructed for all the upstream segments based on the distance between the upstream segment and the outlet segment (Eq.9).The unit-hydrograph convolution with delayed flow (i.e., hillslope-routed flow) is computed for each upstream segment and then all the routed flows from the upstream segments are summed to obtain the streamflow at the outlet segment.As opposed to the KWT routing, the IRF-UH routing does not require the segment sequence for the routing computation.In other words, the routing can be performed in any order of the segments within a river network and for a given segment the unit-hydrograph convolution can be also performed in any order of its upstream segments.

mizuRoute workflow
The overall workflow of mizuRoute is illustrated in Fig. 3.There are two main, separate data preprocessing steps that are executed prior to the routing computation.First, if the hydrologic model simulations are performed with spatial discretization that differs from the HRU used in the river network data, it is necessary to map the runoff output from the hydrologic models to the river network HRUs.This process is done by taking the area-weighted runoff of the intersecting hydrologic model HRUs.We developed the python scripts to identify the intersected hydrologic model HRUs for each river network HRU and their fractional areas to the river network HRU area to assist with this process.
The second data pre-processing step is augmentation of the river network data set.Typical topological information in this data set is the immediate downstream segment for each segment.While a river network can be fully defined based on information about the immediate downstream segment, the river routing schemes in mizuRoute require identification of all the upstream river segments.For this purpose, we have developed a program that identifies all the upstream segments for each segment in the river network data based solely on network topology.This identification of upstream segments only has to be done once for each unique river network data set.Therefore, the program can be used as a preprocessor, which improves the efficiency of the main routing tool, especially when the routing is performed for multiple hydrologic model outputs for a large river system.In addition to the identification of all upstream segments, the topology program identifies upstream HRUs, upstream areas (cumulative area of all the upstream HRUs), total upstream distance from each segment to all the upstream segments, etc.

N. Mizukami et al.: mizuRoute version 1
Figure 6.Daily mean streamflow (DMQ) at the three selected gauges in the GF river network.See Fig. 5 for the locations of the three gauges.

CONUS-wide mizuRoute simulations
The purpose of this section is to demonstrate the capabilities of mizuRoute to route multi-decadal runoff outputs from hydrologic model simulations over the continental domain.We use the United States Geological Survey (USGS) Geospatial Fabric (GF) vector-based river network (Viger, 2014; http://wwwbrr.cr.usgs.gov/projects/SW_MoWS/GeospatialFabric.html) over the CONUS.We routed the daily runoff simulations archived by Reclamation (2014) as part of their project "Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections" (http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.html).In that project, the VIC model was forced by the spatially downscaled temperature and precipitation outputs at The river routing scheme uses both KWT and IRF-UH.The routing parameters for each scheme (see Table 1) need to be predetermined.The channel parameters included in the KWT routing method (Manning's coefficient, n, and river width, w) can be determined by a survey of river channel geometry and river bed condition if the spatial scale of the model domain is very small, but this is usually infeasible for large spatial domains such as the entire CONUS used here.For the IRF-UH method, the determination of celerity and diffusivity with Eq. (B8) requires information on flow and channel geometry, so for simplicity we follow Lohmann et al. (1996) and treat celerity and diffusivity as parameters.For both schemes, parameter estimation methods need to be developed to determine appropriate values for large-scale applications.For this simulation, the parameter values are set somewhat arbitrarily to reasonable values, with the objective to demonstrate the capabilities of mizuRoute to produce spatially distributed streamflow, not to attain the most accurate simulation.In addition, sensitivity of the streamflow estimates to the river routing parameters is examined at selected locations.Different routing model choices (routing scheme and parameters) will differently affect the attenuation of runoff (i.e., the magnitude of peak and rate of rising and recession limbs) and the timing of the peak flow.We also discuss the effect of different river networks (grid-based and vector-based networks) on the results of the runoff routing.
Note that the accuracy of the routed flow is not discussed because it depends largely on the performance of the hydrologic model that produces the distributed runoff fields and hydrologic model outputs are input to the routing model.The hydrologic simulations can have large errors, which makes a direct comparison with observations less meaningful.For this reason, we focus on an inter-comparison between the two channel routing schemes or two river network definitions.The performance of the IRF-UH approach in routing flows compared to observed flows has been discussed by Lohmann et al. (1996).

The geospatial fabric network topology
The GF data set was developed primarily to facilitate CONUS-wide hydrologic modeling with the USGS Precipitation Runoff Modeling System (PRMS; Leavesley and Stan-nard, 1995).To reduce the computational burden of the hydrologic simulations, the GF data set is generated by aggregating fine-scale river segments and corresponding HRUs from the first version of National Hydrography Dataset Plus (NHDPlus v1; HorizonSystemsCorporation, 2010), while still representing small catchments (equivalent in area to 12 digit Hydrologic Unit Code ∼ 100 km 2 or smaller basin).The GF data set includes line and polygon geometries representing river segments and their HRUs, respectively, along with their attribute information including the connectivity between segments (topological information) and their physical attributes such as channel length and area of the HRU.Table 2 lists the river network vector information necessary for mizuRoute.The GF data set (both geometry and attribute information) is stored in Environmental System Research Institute (ESRI) Geodatabase Feature Classes and the topological and physical data (Table 2) in the attribute table is converted to NetCDF format to start with the augmentation of river network topology (Fig. 3).The GF data set include 54 929 river segments and 106 973 HRUs (including the right and left bank of each segment).Figure 4 displays the distribution of river segments in the GF vector data with colors coded by the total upstream HRU area of each river segment.To use the GF vector-based river network, the 1/8  ded runoff outputs from VIC forced by CMIP-5 data were mapped to each GF HRU by taking the areal-weighted average of the intersecting area between grid boxes and the GF HRUs.Although this paper illustrates runoff routing using GF, the mizuRoute tool can work with any other river network data as long as it includes information about the correspondence between HRUs and river segments as well as segment-to-segment topology.

Spatially distributed streamflow in the river network
The first example demonstrates mizuRoute's capability to produce spatially distributed streamflow estimates over the continental domain.Figure 5 shows daily mean streamflow distribution estimated based on VIC-M02 runoff with KWT and IFR-UH routing methods for 15 June 1986 as an example.As shown in Fig. 5, both routing schemes produce qualitatively the same spatial pattern of the daily streamflow.The mizuRoute tool outputs the time series of the streamflow estimates at all the river segments in the river network in the NetCDF output file, and modeled streamflow for the point of interest (e.g., streamflow gauge location) can be extracted from the NetCDF based on the ID of the river segment (i.e., seg_id) where the point of interest is located.Figure 6 shows daily streamflow from 1 January 1995 to 31 December 1999 extracted at three locations from the NetCDF output: (a) Snake River below Ice Harbor Dam, (b) Colorado River at Lees Ferry, and (c) Apalachicola River Blountstown.Temporal patterns of flow simulations with the two river routing schemes are very similar, but the day-to-day differences in estimated streamflow due to the different routing choices become visible.
The next demonstration of mizuRoute's capability is to produce an ensemble of projected streamflow estimates from the runoff simulations using CMIP5 data.Figure 7 shows the monthly mean of 28 projected streamflow estimates (using CMIP5 RCP 8.5 scenario) extracted at the three locations over three periods: (P1) from 2010 to 2039, (P2) from 2040 to 2069, and (P3) from 2070 to 2099.In this example, the results from the KWT scheme are shown in Fig. 7.

Sensitivity of streamflow estimates to river routing parameters
Analysis of the sensitivity of simulated hydrographs to channel routing parameters (Table 1) is performed to examine the effect of parameter values on the streamflow simulations.In this paper, qualitative analysis was performed using VIC simulated runoff with M02 data and using different river routing parameter values (two parameters for each scheme).We carried out the parameter sensitivity analysis at the three locations in Fig. 5, but found the characteristics of the parameter sensitivity are the same at all three.Therefore, we present the results for the Colorado River at Lees Ferry, where a single, distinct snowmelt runoff peak illustrates the impact of the routing parameter values on the peak timing.Figure 8 shows the effect of the width factor W a in Eq. ( 5) (top panels) and the Manning coefficient n (bottom panels) for the KWT scheme.As expected, wider channels (larger W a value) delay the hydrograph, because the larger flow area results in slower velocities.This effect is enhanced with larger Manning coefficient n, because more friction slows the water flow.A similar effect is seen in the sensitivity experiments for Manning's n (bottom panel of Fig. 8).
Figure 9 shows the sensitivity of a simulated hydrograph from 1 October 1990 to 30 September 1991 to the two IRF-UH parameters at Colorado River at Lees Ferry (top panel for sensitivity to celerity C and bottom panel for sensitivity to diffusivity D).Interestingly, the effect of diffusivity D is small while celerity C affects the timing of the hydrograph peak.This is because celerity C directly changes peak timing without attenuation of IRF, while diffusivity D has little influence on peak timing of IRF although it changes the degree of flashiness (Eq.9).Due to the low sensitivity of the hydro-graph to diffusivity D, the degree of hydrograph sensitivity to celerity C is consistent across different diffusivity values (bottom panel of Fig. 9).

Comparison between grid-based and vector-based river network
This section illustrates the effect of river network definitions (grid-or vector-based network) on simulated streamflow using the upper Colorado River basin (outlet: Colorado River at Lees Ferry) and two sub-basins (outlets: Colorado River near Cameo and East River near Almont; See Fig. 10).The daily simulated streamflow from March to August 1999 is shown in Fig. 11.The IRF-UH routing scheme in mizuRoute was used to route the VIC-M02 runoff through both GF river network and 1/8 • grid-based river network.The simulated streamflow time series at the two sub-basins were extracted from the routing results over the entire upper Colorado River basin.Model elements can have only one downstream outlet, Colorado River at Lees Ferry for this simulation.As a result, fractional areas of model grid cells on internal basin boundaries cannot be accounted for.In other words, internal basin boundaries for sub-basins follow grid box edges and a grid cell is either inside or outside a subbasin (Fig. 10b).This leads to discrepancies of basin areas  for sub-basins and total runoff volume that is routed to the gauge as indicated in Fig. 11.Even though the basin areas and therefore flow amounts are similar at Lees Ferry for both networks, they differ for the two sub-basins.For example, the simulated streamflow at Colorado River near Cameo is larger for the vector-based network than for the grid-based network (middle panel in Fig. 11) because of a mismatch in drainage area (Fig. 10).The vector-based river network preserves a more accurate drainage shape or area for sub-basins than the 1/8 • grid-based network.

Summary and discussion
This paper presents mizuRoute (version 1.0), a river network routing tool that post-processes runoff outputs from any hydrologic or land surface model.We demonstrated the capability of mizuRoute to produce multi-decadal, spatially distributed streamflow on a vector-based river network using the USGS GF river network over the CONUS.The streamflow time series are easily extracted at any locations in the network, facilitating hydrologic modeling evaluation, and other hydrologic assessments.The tool is independent of the hydrologic simulations, making it possible to produce ensembles of streamflow estimations from multiple hydrologic models.As an example of a practical application of mizuRoute, an ensemble of streamflow projections was produced at USGS gauge points on the river systems across the CONUS from 97 runoff simulations from Downscaled CMIP5 Climate and Hydrology Projections (Reclamation, 2014).Section 5.2 shows some of the streamflow simulations based on the runoff generated with VIC forced by CMIP5 data.
Based on the simulations presented in the Sect.5.3, the routing parameters can affect the simulated hydrograph especially for the KWT method.Though more detailed investigations of those effects need to be performed to fully understand the routing model behaviors, the parameter sensitivity is substantial.More sophisticated methods to estimate routing model parameters need to be developed.River physical parameters are difficult to obtain in a consistent way at the continental scale, but recent developments of the retrieval algorithms for river physical properties (channel width, slope etc.) with remote sensing data are promising (e.g., Pavelsky and Smith, 2008;Fisher et al., 2013;Allen and Pavelsky, 2015), and we expect to see advances in capabilities to estimate the hydraulic geometry of rivers over the coming years (Clark et al., 2015).
One limitation of mizuRoute is that the channel routing schemes -KWT and IRF-UH -are both 1-D approaches that do not explicitly track physical parcels of water.The 1-D approach does not allow for explicit modeling of inundation extent, which can occur during flood events.Also, the wave particles that are used in the KWT approach travel at the speed of the wave (celerity) rather than the mean velocity of the fluid.Therefore, direct use of KWT for water quality modeling such as stream temperature is not recommended.Extension of mizuRoute to simulate stream temperature and water quality can be done in one of two ways: adaptation of the existing routing methods or inclusion of an additional routing scheme that is more directly suitable for tracking water masses and their constituents.
Toward future enhancements of mizuRoute performance, both routing schemes lend themselves well for parallelization.Computing speed can be improved by implementing parallel processing directive (e.g., open MP) for routing routines.While kinematic wave routing has to be done sequentially from upstream to downstream, the processing can be parallelized through appropriate choices of the domain decomposition.For example, sub-basins that contribute to flow along a mainsteam segment can be processed in parallel because the basins are independent.On a CONUS-wide river network, individual river basins (e.g., the Colorado River and Mississippi River basins) can be processed simultaneously.For IRF-UH routing, the routing computation is performed for individual river segments independently (see Sect. 3.2.2);therefore, the parallelization for river segment loops can be made possible.Lastly, routing of an ensemble of runoff outputs such as the CMIP5 projected runoff is easily parallelized.

Code availability
The source codes for the river network topology program and the hillslope and river routing along with test data are available along with the user manual on GitHub (http://dx.doi.org/10.5281/zenodo.56043).Those codes are developed in For-tran90 and require installation of a NetCDF 4 library (http: //www.unidata.ucar.edu/downloads/netcdf/index.jsp).In addition, there are several pre-processing python scripts to map runoff outputs from hydrologic models to other types of HRUs.These pre-processing scripts are also available in GitHub.Those Python scripts process ESRI Shapefiles and NetCDF data and require GDAL, SHAPELY, NetCDF4 packages.∂y ∂s = S 0 − S f .(B1) Now, Manning's equation can be expressed in terms of channel conveyance, K c (carrying capacity of river channel), where Substituting S f from Eq. (B2) into Eq.(B1) and differentiating with respect to time, the momentum equation (Eq.A1) becomes where parameters C and D are wave celerity [LT −1 ] and diffusivity [L 2 T −1 ], respectively.Here, we assume the flow is uniform (i.e., S f = S 0 ).

Figure 1 .
Figure 1.Comparison of 1/8• (∼ 12 km) gridded river network and vector river network from United States Geological Survey (USGS) Geospatial Fabric for the upper part of Snake River basin(Viger, 2014).

Figure 2 .
Figure 2. Visualization of waves.The top panel (a) plots a discharge of the wave [m 3 s −1 ] against its location (distance [km] from the beginning of the 1st segment) at the beginning of five consecutive daily time steps.A vertical line indicates the river segment boundary.The bottom panel (b) plots a discharge of the wave [m 3 s −1 ] against its exit time (day from 1 October) for the 4th and 13th segments.A vertical line indicates the boundary of routing time step.The inserted map shows the 16 river segments used for the plots.

Figure 4 .
Figure 4. GF river network color coded by upstream drainage areas.Gray lines indicate the total upstream drainage areas less than 12 000 km 2 .

Figure 5 .
Figure 5. Spatially distributed daily streamflow on 15 July 1986 in the GF river network simulated with mizuRoute.Gray lines indicate flow less than 100 m 3 s −1 .The streamflow time series shown in Fig. 6 are extracted at three USGS gauges A-C.
1/8 • (∼ 12 km) resolution from 97 global climate model outputs from 1950 through 2099.The details of the Coupled Model Intercomparison Project Phase 5 (CMIP5) are described by Taylor et al. (2011).Additionally, historical runoff simulations were produced at 1/8 • resolution by the VIC model forced by meteorological forcings from Maurer et al. (2002) from 1950 through 1999 (Maurer meteorological data and the simulated runoff with Maurer data are referred to as M02 and VIC-M02 runoff, respectively).

Figure 7 .
Figure 7. Monthly mean of CMIP5 projected streamflow at three locations indicated in Fig. 5. (left column-Snake River below Ice Harbor Dam, middle column -Colorado River at Lees Ferry, and right column -Apalachicola River near Blountstown).The river routing scheme is KWT.Monthly mean values are computed over three future periods (Top -P1 2010-2039, Middle -P2 2040-2069 and bottom -P3 2070-2099).The dash line denotes streamflow estimated from runoff output from VIC forced by M02 historical data while gray lines indicate projected streamflow based on future runoff outputs from VIC forced by 28 CMIP5 RCP8.5 data.

Figure 8 .
Figure 8. Sensitivity of simulated runoff at Colorado River at Lees Ferry (Location B in Fig. 5) to the two KWT parameters.The top panel shows sensitivity to width factor W a with three fixed Manning coefficients n (from left to right: n = 0.005, 0.02, and 0.05).The bottom panel shows sensitivity to Manning coefficient n with three fixed width factor w (from left to right: W a = 0.0005, 0.0050, and 0.0100).

Figure 9 .
Figure 9. Sensitivity of simulated runoff at Colorado River at Lees Ferry (location B in Fig. 5) to IRF-UH parameters.The top panels show sensitivity to diffusivity D with three fixed celerity C values (from left to right: C = 1.0, 2.0, and 3.0 ms −1 ).The bottom panels show sensitivity to celerity C with three fixed diffusivity D values (from left to right: D = 200, 1000, and 3000 ms −2 ).

Figure 10
Figure 10.1/8 • (∼ 12 km) grid-based river network for the upper Colorado River basin (a) and two sub-basins inside the upper Colorado-Colorado River near Cameo and East River at Almont (b).HRUs for each GF river segment are not shown for clarity in (b).

Figure 11 .
Figure 11.Comparison of IRF-UH routed flow with two river networks (1/8 • grid-based river network and GF vector-based river network) at three locations in the upper Colorado River basin (See Fig. 10 for river network and basin boundaries).
continuity equation (Eq. 3) can be re-rewritten by differentiating both sides of the equation with respect to distance x as, conveyance, K c , is a function of flow depth, y, or flow area, A, the differentiation part of the second term of Eq. (B5) can be written as

Table 2 .
River network information required in mizuRoute.