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

Abstract

This paper describes the first version of a stand-alone 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.

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/ doi:10.5194/gmd-9-2223-2016 Author(s) 2016. CC Attribution 3.0 License.

Naoki Mizukami1, Martyn P. Clark1, Kevin Sampson1, Bart Nijssen2, Yixin Mao2, Hilary McMillan3,8, Roland J. Viger4, Steve L. Markstrom4, Lauren E. Hay4, Ross Woods5, Jeffrey R. Arnold6, and Levi D. Brekke7

1National Center for Atmospheric Research, Boulder, CO, USA

2University of Washington, Seattle, WA, USA

3National Institute of Water and Atmospheric Research, Christchurch, New Zealand

4United States Geological Survey, Denver, CO, USA

5University of Bristol, Bristol, UK

6U.S. Army of Corps of Engineers, Seattle, WA, USA

7Bureau of Reclamation, Denver, CO, USA

8San Diego State University, San Diego, CA, USA

Correspondence to: Naoki Mizukami ([email protected])

Received: 8 September 2015 Published in Geosci. Model Dev. Discuss.: 2 November 2015 Revised: 25 March 2016 Accepted: 1 June 2016 Published: 23 June 2016

Abstract. This paper describes the rst version of a stand-alone runoff routing tool, mizuRoute. The mizuRoute tool post-processes runoff outputs from any distributed hydro-logic model or land surface model to produce spatially distributed streamow 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 ner-scale river lines than the grid-based network. Streamow 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 mizuRoutes capabilities to produce spatially distributed streamow 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 stream-ow.

1 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 stream-ow 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 streamow. 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 streamow estimates at all river segments in the river net-

Published by Copernicus Publications on behalf of the European Geosciences Union.

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

2224 N. Mizukami et al.: mizuRoute version 1

work across the domain of interest at each time step, facilitating further spatial and temporal analysis of the estimated streamow. As opposed to other routing models, our goal in developing mizuRoute is to provide exibility in making routing model decisions (i.e., river network denition and routing scheme).

The paper proceeds as follows. Section 2 reviews existing river routing models. Section 3 describes how the mizuRoute tool provides exibility of routing modeling decisions and details the hillslope and river routing schemes used in mizuRoute. Section 4 provides an overview of the work-ow of mizuRoute from preprocessing hydrologic model output to simulating streamow in the river network. Section 5 demonstrates streamow simulations in river systems over the contiguous United States (CONUS). Finally, a summary and future work are discussed in Sect. 6.

2 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 Hydro-logic Engineering Center River Analysis System (HECRAS; Brunner, 2001). HEC-RAS offers various hydraulic routing schemes, ranging from simple uniform ow to one-dimensional (1-D) Saint-Venant equations for unsteady ow.HEC-RAS has been popular among civil engineers for river channel design and oodplain 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 ow equations (e.g., Miguez-Macho and Fan, 2012;Paiva et al., 2013; Clark et al., 2015), simplied 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 Musk-ingum routing (e.g., David et al., 2011). Despite their computational cost, dynamic or diffusive wave models are attractive for relatively at oodplain regions such as along the Amazon River where backwater effects on the ood wave are signicant (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 ood wave delay and attenuation, but do not simulate other streamow variables such as ow velocity and ow 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 dened 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 dened by a grid box and river segments connect neighboring grid boxes based on the ow directions.Vector-based river networks are better than coarser resolution (e.g., > 1 km) gridded river networks at preserving ne-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 ne-resolution ow 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., ODonnell et al., 1999;Fekete et al., 2001; Olivera et al., 2002; Reed, 2003; Davies and Bell, 2009; Wu et al., 2011, 2012). Earlier work (e.g., ODonnell et al., 1999; Fekete et al., 2001; Olivera et al., 2002) focused on preserving the accuracy of the ow direction at the coarser resolution and therefore on an accurate representation of the drainage area. Newer upscaling methods are designed to also preserve ne-scale ow 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).

3 Runoff routing in mizuRoute

The runoff routing in mizuRoute provides more exibility in continental domain routing applications. The mizuRoute tool enables model exibility in two ways: rst, mizuRoute can be used to simulate streamow for both grid- and vector-based river networks. Given either type of river network data, mizuRoute offers an option to route ow along all the river segments in the river network data or route runoff at an outlet segment specied by a user. With the latter option, routing computations are performed only upstream of the spec-ied outlet, which reduces the computational cost. Second, the modular structure of the mizuRoute tool offers the exibility to congure 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 exibility 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-

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

N. Mizukami et al.: mizuRoute version 1 2225

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).

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 hydro-graph. 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 ow 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.

3.1 Hillslope routing

Hillslope routing accounts for the time of concentration (Tc) 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 two-parameter 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

(t : a, ) =

1[Gamma1] (a) a ta1e

(s : a, ) [notdef] R (t s)ds, (2)

where q is delayed runoff or discharge [L3 T1] at time step t [T], R is HRU total runoff depth [L3 T1] from the hydro-logic 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 efcient manner, which is important for continental-scale domains.

3.2 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 ood wave propagation through a river channel. The 1-D conservation equations for continuity (Eq. 3) and momentum (Eq. 4) are

@q @x +

@A

@t = 0, (3) @v

@t + v

@v

@x + g

t

, (1)

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 ashiness (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:

q (t) =

tmax

[integraldisplay]

0

@y

@x g(S0 Sf) = 0, (4)

where q is discharge [L3 T1] at time step t [T] and location

x [L] in a river network, A is cross sectional ow area [L2],

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

Figure 2. Visualization of waves. The top panel (a) plots a discharge of the wave [m3s1] against its location (distance [km] from the beginning of the 1st segment) at the beginning of ve consecutive daily time steps. A vertical line indicates the river segment boundary. The bottom panel (b) plots a discharge of the wave [m3 s1]

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.

waves that reside in the 16 segments (16 river segments are shown in the inserted map) at the beginning of ve 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 predened 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

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

2226 N. Mizukami et al.: mizuRoute version 1

v is velocity [LT1], y is depth of ow [L], S0 is channel slope [], Sf is friction slope [], and g is gravitational constant [LT2]. The continuity equation (Eq. 3) assumes that no lateral ow is added to a channel segment. The following sub-sections describe the two routing schemes.

3.2.1 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 ow 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 C [LT1] is a function

of channel width w[L], Mannings coefcient n [], channel slope S0 [], and discharge q [L3 T1]. Further details are provided in Appendix A. Among the four variables, the channel slope S0 is provided by the river network data and discharge is computed with hillslope routing for the head-water basin, or/and updated via routing from the upstream segment. The other two variables, Mannings coefcient n and river width w, are much more difcult to measure or estimate. The river width is determined with the following widthdrainage area relationship (Booker, 2010):

w = Wa [notdef] Abups, (5) where Wa is a width factor [], Aups is the total upstream basin area [L2], and b is an empirical exponent equal to 0.5.The width factor Wa and the Mannings coefcient 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 rst time step, only the wave from local HRUs exists. Figure 2a visualizes the discharge of

N. Mizukami et al.: mizuRoute version 1 2227

Table 1. Routing model parameters.

Parameters Routing methods Descriptions Values used in Sect. 5

a Hillslope Shape factor [] 2.5 []

Hillslope Timescale factor [T] 86 400 [s] n KWT Manning coefcient [] 0.01 []

W KWT River width scale factor [] 0.001 [] C IRF-UH Wave velocity [LT1] 1.5 [ms1]

D IRF-UH Diffusivity [L2 T1] 800 [m2 s1]

Figure 3. Overview of streamow 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 agged as exited. The exit time then becomes the time the wave entered the downstream segment. Otherwise, the wave is agged 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 ow depth, and thus an increase in the discharge at a xed 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:

Cmerge =

[Delta1]q [Delta1]A =

[Delta1]q

w[Delta1]y , (6)

where Cmerge is the merged wave celerity [LT1], [Delta1]q and [Delta1]A are differences in discharge [L3 T1] and cross sectional ow area [L2 T1], respectively, between slower and faster waves. Note that Eq. (6) is the mathematical denition of the wave celerity. Since we assume a rectangular channel, whose width is constant for each segment, the merged celerity Cmerge is a function of ow depth y, which is computed with Eq. (A4).

4. Finally, the time-step-averaged discharge (streamow) 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.

3.2.2 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 ows from gridded land surface models such as the Variable Inltration 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 dened. 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 briey 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):

@q

@t = D

@2q

@x2 C

@q

@x , (7)

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

2228 N. Mizukami et al.: mizuRoute version 1

Table 2. River network information required in mizuRoute.

Variables Vector data type Descriptions

seg_id River segment line ID of segmenttosegment River segment line ID of immediate downstream segment Length River segment line length of segment [m]

Slope River segment line Slope of segment [m m1]

hru_id HRU polygon ID of HRUSeg_hru_segment HRU polygon ID of segment to which the HRU discharges hru_area HRU polygon Area of HRU [m2]

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

where parameters C and D are wave celerity [LT1] and diffusivity [L2T1], 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

q =

t

and the outlet segment (Eq. 9). The unit-hydrograph convolution with delayed ow (i.e., hillslope-routed ow) is computed for each upstream segment and then all the routed ows from the upstream segments are summed to obtain the streamow 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.

4 mizuRoute workow

The overall workow 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

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

[integraldisplay]

0

U (t s)h(x,s)ds, (8)

where

h(x,t) =

2tpDt exp[parenleftBigg]

(Ct x)2

4Dt

x

[parenrightBigg]

(9)

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

N. Mizukami et al.: mizuRoute version 1 2229

Figure 5. Spatially distributed daily streamow on 15 July 1986 in the GF river network simulated with mizuRoute. Gray lines indicate ow less than 100 m3 s1. The streamow time series shown in Fig. 6 are extracted at three USGS gauges AC.

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 dened based on information about the immediate downstream segment, the river routing schemes in mizuRoute require identication of all the upstream river segments. For this purpose, we have developed a program that identies all the upstream segments for each segment in the river network data based solely on

network topology. This identication 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 efciency of the main routing tool, especially when the routing is performed for multiple hydro-logic model outputs for a large river system. In addition to the identication of all upstream segments, the topology program identies upstream HRUs, upstream areas (cumulative area of all the upstream HRUs), total upstream distance from each segment to all the upstream segments, etc.

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

2230 N. Mizukami et al.: mizuRoute version 1

Figure 6. Daily mean streamow (DMQ) at the three selected gauges in the GF river network. See Fig. 5 for the locations of the three gauges.

5 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

Web End =http://wwwbrr.cr.usgs.gov/projects/SW_ http://wwwbrr.cr.usgs.gov/projects/SW_MoWS/GeospatialFabric.html

Web End =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

Web End =http://gdo-dcp.ucllnl.org/ http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.html

Web End =downscaled_cmip_projections/dcpInterface.html ). In that project, the VIC model was forced by the spatially downscaled temperature and precipitation outputs at 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).

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 (Mannings coefcient, 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 ow 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 streamow, not to attain the most accurate simulation.

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

N. Mizukami et al.: mizuRoute version 1 2231

Figure 7. Monthly mean of CMIP5 projected streamow 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 20102039, Middle P2 20402069 and bottom P3 2070 2099). The dash line denotes streamow estimated from runoff output from VIC forced by M02 historical data while gray lines indicate projected streamow based on future runoff outputs from VIC forced by 28 CMIP5 RCP8.5 data.

In addition, sensitivity of the streamow 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 ow. 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 ow is not discussed because it depends largely on the performance of the hydro-logic model that produces the distributed runoff elds 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 denitions. The performance of the IRF-UH approach in routing ows compared to observed ows has been discussed by Lohmann et al. (1996).

5.1 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 ne-scale river segments and corresponding HRUs from the rst 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 km2 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 grid-

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

2232 N. Mizukami et al.: mizuRoute version 1

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 Wa with three xed Manning coefcients n (from left to right: n = 0.005, 0.02, and 0.05). The bottom panel

shows sensitivity to Manning coefcient n with three xed width factor w (from left to right: Wa = 0.0005, 0.0050, and 0.0100).

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.

5.2 Spatially distributed streamow in the river network

The rst example demonstrates mizuRoutes capability to produce spatially distributed streamow estimates over the continental domain. Figure 5 shows daily mean streamow 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 streamow.

The mizuRoute tool outputs the time series of the stream-ow estimates at all the river segments in the river network in the NetCDF output le, and modeled streamow for the point of interest (e.g., streamow 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 streamow 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 ow simulations with the two river routing schemes are very similar, but the day-to-day differences in estimated streamow due to the different routing choices become visible.

The next demonstration of mizuRoutes capability is to produce an ensemble of projected streamow estimates from the runoff simulations using CMIP5 data. Figure 7 shows the monthly mean of 28 projected streamow 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.

5.3 Sensitivity of streamow 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 streamow simulations. In this paper, qualitative analysis was performed using VIC simulated runoff with M02 data and using different river routing

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

N. Mizukami et al.: mizuRoute version 1 2233

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 xed 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 xed diffusivity D values (from left to right: D = 200, 1000, and 3000 ms

2).

graph to diffusivity D, the degree of hydrograph sensitivity to celerity C is consistent across different diffusivity values (bottom panel of Fig. 9).

5.4 Comparison between grid-based and vector-based river network

This section illustrates the effect of river network denitions (grid- or vector-based network) on simulated streamow 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 streamow 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 streamow 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 sub-basin (Fig. 10b). This leads to discrepancies of basin areas

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

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 Wa in Eq. (5) (top panels) and the Manning coefcient n (bottom panels) for the KWT scheme. As expected, wider channels (larger Wa value) delay the hydrograph, because the larger ow area results in slower velocities. This effect is enhanced with larger Manning coefcient n, because more friction slows the water ow. A similar effect is seen in the sensitivity experiments for Mannings 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 IRFUH 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 inuence on peak timing of IRF although it changes the degree of ashiness (Eq. 9). Due to the low sensitivity of the hydro-

2234 N. Mizukami et al.: mizuRoute version 1

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. Comparison of IRF-UH routed ow 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).

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

N. Mizukami et al.: mizuRoute version 1 2235

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 ow amounts are similar at Lees Ferry for both networks, they differ for the two sub-basins. For example, the simulated streamow 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.

6 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 streamow on a vector-based river network using the USGS GF river network over the CONUS. The stream-ow 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 streamow estimations from multiple hydro-logic models. As an example of a practical application of mizuRoute, an ensemble of streamow 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 streamow 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 difcult 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 ood 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 uid. 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 ow 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.

7 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

Web End =http://dx.doi.org/ http://dx.doi.org/10.5281/zenodo.56043

Web End =10.5281/zenodo.56043 ). Those codes are developed in Fortran90 and require installation of a NetCDF 4 library (http://www.unidata.ucar.edu/downloads/netcdf/index.jsp

Web End =http: http://www.unidata.ucar.edu/downloads/netcdf/index.jsp

Web End =//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 Shapeles and NetCDF data and require GDAL, SHAPELY, NetCDF4 packages.

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

2236 N. Mizukami et al.: mizuRoute version 1

Appendix A: Derivation of wave celerity equation used in KWT

The kinematic wave approximation to the full Saint-Venant equations (Eqs. 3 and 4) uses the continuity equation combined with a simplied momentum equation. The simplied momentum equation is based on the assumption that the friction slope is equal to the channel slope and that ow is steady and uniform. Under this assumption, Eq. (4) is reduced to S0 = Sf. In other words, the gravitational force that moves

water downstream is balanced with the frictional force acting on the riverbed. With this assumption, the discharge q can be expressed using a uniform ow formula such as Mannings equation:

q = A

basis of the IRF-UH method. The development of the IRFUH method starts with the derivation of the diffusive wave equation from the 1-D Saint Venant equations (Eqs. 3 and 4) by neglecting inertia terms (the second term in Eq. (4) and assuming steady ow (eliminating the rst term in Eq. 4). The momentum equation (Eq. 4) can therefore be reduced to

@y

@s = S0 Sf. (B1)

Now, Mannings equation can be expressed in terms of channel conveyance, Kc (carrying capacity of river channel),

q = Kc [notdef]

pSf, (B2)

where Kc = k/n [notdef] A [notdef] R h. Substituting Sf from Eq. (B2) into

Eq. (B1) and differentiating with respect to time, the momentum equation (Eq. A1) becomes

2q K2c

@q

@t

k

nR hS1/20, (A1)

where k is a scalar, whose value is 1 for SI units and 1.49 for Imperial units, n is the Manning coefcient, Rh is hydraulic radius [L], which is dened as the cross sectional ow area A [L2] divided by the wetted perimeter P [L], and is a constant coefcient ( = 2/3).

We assume the channel shape is rectangular and the geometry is constant throughout one river segment, with width w, A = wy, and P = w + 2y. Assuming the channel is wide

compared to ow depth (i.e., w y), the hydraulic radius

Rh is expressed as

Rh =

A

P =

2q2

K3c

@Kc

@t =

@2y

@x@t . (B3)

Also, the continuity equation (Eq. 3) can be re-rewritten by differentiating both sides of the equation with respect to distance x as,

@2q

@x2 + w

@2y

wy w + 2y

= y. (A2)

By substituting Eq. (A2) into Eq. (A1), the Manning equation is re-written as

q = w

k

ny +1S1/20. (A3)

For each stream segment within which the channel width w is constant, the wave celerity C is given by

C =

dq dA =

@x@t = 0 (B4)

Combining Eqs. (B3) and (B4) results in

2q K2c

@q

@t

2q2

K3c

@Kc

@t =

1 w

@2q

@x2 . (B5)

Because the channel conveyance, Kc, is a function of ow depth, y, or ow area, A, the differentiation part of the second term of Eq. (B5) can be written as

@Kc

@t =

dKc dA

@A

@t =

dKc dA

dq d(wy)

=

dq

wdy . (A4)

By substituting Mannings equation (Eq. A3) into Eq. (A4), the wave celerity C can be given by

C = ( + 1) [notdef]

k n

@q

@x . (B6)

Finally, inserting Eq. (B6) into Eq. (B5), results in the one-dimensional diffusive wave equation

@q

@t = D

@2q

@x2 C

@q

@x , (B7)

where

C =

pS0 [notdef] y (A5)

or expressed as a function of discharge q as

C = ( + 1) [notdef] (w)

+1 [notdef]

q Kc

dKc dA =

dq dA

, (B8)

where parameters C and D are wave celerity [LT1] and diffusivity [L2T1], respectively. Here, we assume the ow is uniform (i.e., Sf = S0).

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/

k n

pS0

[parenrightbigg]

1 +1

[notdef] q

+1 . (A6)

D =

K2c

2qw =

q 2wS0

Appendix B: Derivation of 1-D diffusive equation

We describe in detail the derivation of diffusive wave equations from Saint-Venant equations (Strum, 2001) that are the

N. Mizukami et al.: mizuRoute version 1 2237

Acknowledgements. This work was nancially supported by theU.S Army Corps of Engineers Climate Preparedness and Resilience Program.

Edited by: D. Roche

References

Allen, G. H. and Pavelsky, T. M.: Patterns of river width and surface area revealed by the satellite-derived North American River Width data set, Geophys. Res. Lett., 42, 2014GL062764, doi:http://dx.doi.org/10.1002/2014gl062764

Web End =10.1002/2014gl062764 http://dx.doi.org/10.1002/2014gl062764

Web End = , 2015.

Arora, V. K. and Boer, G. J.: A variable velocity ow routing algorithm for GCMs, J. Geophys. Res.-Atmos., 104, 3096530979, 1999.

Bhunya, P. K., Berndtsson, R., Ojha, C. S. P., and Mishra, S. K.:

Suitability of Gamma, Chi-square, Weibull, and Beta distributions as synthetic unit hydrographs, J. Hydrol., 334, 2838, 2007. Booker, D. J.: Predicting wetted width in any river at any discharge,

Earth Surf. Process. Landf., 35, 828841, doi:10.1002/esp.1981, 2010.

Brunner, G. W.: HEC-RAS, River Analysis Ssytem UsersManual, US Army Corps of Engineers, Hydrologic Engineering Center, 320, 2001.

Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G., Schmidt, J., and Uddstrom, M. J.: Hydrological data assimilation with the ensemble Kalman lter: Use of stream-ow observations to update states in a distributed hydrological model, Adv. Water Res., 31, 13091324, 2008.

Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster,D., Gochis, D. J., Hooper, R. P., Kumar, M., Leung, L. R., Mackay, D. S., Maxwell, R. M., Shen, C., Swenson, S. C., and Zeng, X.: Improving the representation of hydrologic processes in Earth System Models, Water Resour. Res., 51, 59295956, doi:http://dx.doi.org/10.1002/2015WR017096

Web End =10.1002/2015WR017096 http://dx.doi.org/10.1002/2015WR017096

Web End = , 2015.

David, C. H., Maidment, D. R., Niu, G.-Y., Yang, Z.-L., Habets,F., and Eijkhout, V.: River Network Routing on the NHDPlus Dataset, J. Hydrometeorol., 12, 913934, 2011.

Davies, H. N. and Bell, V. A.: Assessment of methods for extracting low-resolution river networks from high-resolution digital data, Hydrol. Sci. J., 54, 1728, 2009.

Fekete, B. M., Vrsmarty, C. J., and Lammers, R. B.: Scaling gridded river networks for macroscale hydrology: Development, analysis, and control of error, Water Resour. Res., 37, 1955 1967, 2001.

Fisher, G. B., Bookhagen, B., and Amos, C. B.: Channel planform geometry and slopes from freely available high-spatial resolution imagery and DEM fusion: Implications for channel width scalings, erosion proxies, and uvial signatures in tectonically active landscapes, Geomorphology, 194, 4656, 2013.

Gochis, D. J., Yu, W., and Yates, D. N.: The WRF-Hydro model technical description and users guide, version 3.0. NCAR Technical Document, NCAR, 120 pp., Boulder CO, 2015.

Goring, D. G.: Kinematic shocks and monoclinal waves in the Waimakariri, a steep, braided, gravel-bed river, Proceedings of the International Symposium on waves: Physical and numerical modelling, University of British Columbia, Vancouver, Canada, 336345, 1994.

Goteti, G., Famiglietti, J. S., and Asante, K.: A Catchment-Based Hydrologic and Routing Modeling System with explicit river channels, J. Geophys. Res.-Atmos., 113, D14116, doi:http://dx.doi.org/10.1029/2007JD009691

Web End =10.1029/2007JD009691 http://dx.doi.org/10.1029/2007JD009691

Web End = , 2008.

Horizon Systems Corporation: NHDPlus Version 1 (NHDPlusV1) User Guide, available at: http://www.horizon-systems.com/NHDPlus/NHDPlusV1_documentation.php

Web End =http://www.horizon-systems.com/

http://www.horizon-systems.com/NHDPlus/NHDPlusV1_documentation.php

Web End =NHDPlus/NHDPlusV1_documentation.php (last access: 19 June 2016), 2010.

Koren, V., Reed, S., Smith, M., Zhang, Z., and Seo, D.-J.: Hydrology laboratory research modeling system (HL-RMS) of the US national weather service, J. Hydrol., 291, 297318, 2004.Kumar, R., Chatterjee, C., Singh, R. D., Lohani, A. K., and Kumar,S.: Runoff estimation for an ungauged catchment using geomorphological instantaneous unit hydrograph (GIUH) models, Hydrol. Process., 21, 18291840, 2007.

Leavesley, G. H. and Stannard, L. G.: The precipitation-runoff modeling system-PRMS. Computer Models of Watershed Hydrology, edited by: Singh, V. P., Water Resouces Publications, 281310, 1995.

Lehner, B. and Grill, G.: Global river hydrography and network routing: baseline data and new approaches to study the worlds large river systems, Hydrol. Process., 27, 21712186, 2013.

Li, H., Wigmosta, M. S., Wu, H., Huang, M., Ke, Y., Coleman, A.M., and Leung, L. R.: A Physically Based Runoff Routing Model for Land Surface and Earth System Models, J. Hydrometeorol., 14, 808828, 2013.

Liang, X., Lettenmaier, D. P., Wood, E. F., and Burges, S. J.: A simple hydrologically based model of land surface water and energy uxes for general circulation models, J. Geophys. Res., 99, 1441514428, 1994.

Lohmann, D., Nolte-Holube, R., and Raschke, E.: A large-scale horizontal routing model to be coupled to land surface parametrization schemes, Tellus A, 48, 708721, doi:http://dx.doi.org/10.3402/tellusa.v48i5.12200

Web End =10.3402/tellusa.v48i5.12200 http://dx.doi.org/10.3402/tellusa.v48i5.12200

Web End = , 1996.

Lohmann, D., Raschke, E., Nijssen, B., and Lettenmaier, D. P.: Regional scale hydrology: I. Formulation of the VIC-2L model coupled to a routing model, Hydrol. Sci. J., 43, 131141, 1998.Lucas-Picher, P., Arora, V. K., Caya, D., and Laprise, R.: Implementation of a large-scale variable velocity river ow routing algorithm in the Canadian Regional Climate Model (CRCM), Atmos.-Ocean, 41, 139153, 2003.

Maurer, E. P., Wood, A. W., Adam, J. C., Lettenmaier, D. P., and Nijssen, B.: A Long-Term Hydrologically Based Dataset of Land Surface Fluxes and States for the Conterminous United States, J.Climate, 15, 32373251, 2002.

McDonnell, J. J. and Beven, K.: DebatesThe future of hydrological sciences: A (common) path forward? A call to action aimed at understanding velocities, celerities and residence time distributions of the headwater hydrograph, Water Resour. Res., 50, 53425350, 2014.

Miguez-Macho, G. and Fan, Y.: The role of groundwater in the Amazon water cycle: 1. Inuence on seasonal streamow, ooding and wetlands, J. Geophys. Res.-Atmos., 117, D15113, doi:http://dx.doi.org/10.1029/2012JD017539

Web End =10.1029/2012JD017539 http://dx.doi.org/10.1029/2012JD017539

Web End = , 2012.

Nadarajah, S.: Probability models for unit hydrograph derivation, J.Hydrol., 344, 185189, 2007.

Nijssen, B., Lettenmaier, D. P., Liang, X., Wetzel, S. W., and Wood,E. F.: Streamow simulation for continental-scale river basins, Water Resour. Res., 33, 711724, 1997.

www.geosci-model-dev.net/9/2223/2016/ Geosci. Model Dev., 9, 22232238, 2016

2238 N. Mizukami et al.: mizuRoute version 1

ODonnell, G., Nijssen, B., and Lettenmaier, D. P.: A simple algorithm for generating streamow networks for grid-based, macroscale hydrological models, Hydrol. Process., 13, 1269 1275, 1999.

Olivera, F., Lear, M. S., Famiglietti, J. S., and Asante, K.: Extracting low-resolution river networks from high-resolution digital elevation models, Water Resour. Res., 38, 1231, doi:http://dx.doi.org/10.1029/2001WR000726

Web End =10.1029/2001WR000726 http://dx.doi.org/10.1029/2001WR000726

Web End = , 2002.

Paiva, R. C. D., Collischonn, W., and Tucci, C. E. M.: Large scale hydrologic and hydrodynamic modeling using limited data and a GIS based approach, J. Hydrol., 406, 170181, 2011.

Paiva, R. C. D., Buarque, D. C., Collischonn, W., Bonnet, M.-P.,

Frappart, F., Calmant, S., and Bulhes Mendes, C. A.: Large-scale hydrologic and hydrodynamic modeling of the Amazon River basin, Water Resour. Res., 49, 12261243, 2013.Pavelsky, T. M. and Smith, L. C.: RivWidth: A Software Tool for the Calculation of River Widths From Remotely Sensed Imagery, Geosci. Remote Sens. Lett. IEEE, 5, 7073, 2008.Reclamation: Downscaled CMIP3 and CMIP5 Hydrology Projections Release of Hydrology Projections, Comparison with Preceding Information and Summary of User Needs, Technical Memorandum No. 86-68210-2011-01, 110 pp., available at: http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/techmemo/BCSD5HydrologyMemo.pdf

Web End =http://gdo-dcp.ucllnl.org/downscaled_cmip_ http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/techmemo/BCSD5HydrologyMemo.pdf

Web End =projections/techmemo/BCSD5HydrologyMemo.pdf (last access:19 June 2016), 2014.

Reed, S. M.: Deriving ow directions for coarse-resolution (1 4 km) gridded hydrologic modeling, Water Resour. Res., 39, 1238, doi:http://dx.doi.org/10.1029/2003WR001989

Web End =10.1029/2003WR001989 http://dx.doi.org/10.1029/2003WR001989

Web End = , 2003.

Strum, T. W.: Open Channel Hydraulics, 1st Edn., McGraw-Hill,493 pp., 2001.

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485498, 2011.

Viger, R. J.: Preliminary spatial parameters for PRMS based on the Geospatial Fabric, NLCD2001 and SSURGO, US Geological Survey, doi:http://dx.doi.org/10.5066/F7WM1BF7

Web End =10.5066/F7WM1BF7 http://dx.doi.org/10.5066/F7WM1BF7

Web End = , 2014.

Wu, H., Kimball, J. S., Mantua, N., and Stanford, J.: Automated up-scaling of river networks for macroscale hydrological modeling, Water Resour. Res., 47, W03517, doi:http://dx.doi.org/10.1029/2009WR008871

Web End =10.1029/2009WR008871 http://dx.doi.org/10.1029/2009WR008871

Web End = , 2011.

Wu, H., Kimball, J. S., Li, H., Huang, M., Leung, L. R., and Adler, R. F.: A new global river network database for macroscale hydrologic modeling, Water Resour. Res., 48, W09701, doi:http://dx.doi.org/10.1029/2012WR012313

Web End =10.1029/2012WR012313 http://dx.doi.org/10.1029/2012WR012313

Web End = , 2012.

Xia, Y., Mitchell, K., Ek, M., Cosgrove, B., Shefeld, J., Luo,L., Alonge, C., Wei, H., Meng, J., Livneh, B., Duan, Q., and Lohmann, D.: Continental-scale water and energy ux analysis and validation for North American Land Data Assimilation System project phase 2 (NLDAS-2): 2. Validation of model-simulated streamow, J. Geophys. Res.-Atmos., 117, D03110, doi:http://dx.doi.org/10.1029/2011JD016051

Web End =10.1029/2011JD016051 http://dx.doi.org/10.1029/2011JD016051

Web End = , 2012.

Yamazaki, D., Oki, T., and Kanae, S.: Deriving a global river network map and its sub-grid topographic characteristics from a ne-resolution ow direction map, Hydrol. Earth Syst. Sci., 13, 22412251, doi:http://dx.doi.org/10.5194/hess-13-2241-2009

Web End =10.5194/hess-13-2241-2009 http://dx.doi.org/10.5194/hess-13-2241-2009

Web End = , 2009.

Yamazaki, D., de Almeida, G. A. M., and Bates, P. D.: Improving computational efciency in global river models by implementing the local inertial ow equation and a vector-based river network map, Water Resour. Res., 49, 72217235, 2013.

Yamazaki, D., Kanae, S., Kim, H., and Oki, T.: A physically based description of oodplain inundation dynamics in a global river routing model, Water Resour. Res., 47, W04501, doi:http://dx.doi.org/10.1029/2010WR009726

Web End =10.1029/2010WR009726 http://dx.doi.org/10.1029/2010WR009726

Web End = , 2011.

Yucel, I., Onen, A., Yilmaz, K. K., and Gochis, D. J.: Calibration and evaluation of a ood forecasting system: Utility of numerical weather prediction model, data assimilation and satellite-based rainfall, J. Hydrol., 523, 4966, 2015.

Zaitchik, B. F., Rodell, M., and Olivera, F.: Evaluation of the Global Land Data Assimilation System using global river discharge data and a source-to-sink routing scheme, Water Resour. Res., 46, W06507, doi:http://dx.doi.org/10.1029/2009WR007811

Web End =10.1029/2009WR007811 http://dx.doi.org/10.1029/2009WR007811

Web End = , 2010.

Geosci. Model Dev., 9, 22232238, 2016 www.geosci-model-dev.net/9/2223/2016/