To what extent does river routing matter in hydrological modeling?

Spatially distributed hydrology and land surface models are typically applied in combination with river routing schemes that convert instantaneous runoff into streamflow. Nevertheless, the development of such schemes has been somehow disconnected from hydrologic model calibration research, although both seek to achieve more realistic streamflow simulations. In this paper, we seek to bridge this gap to understand the extent to which the configuration of routing schemes affects hydrologic model parameter searches in water resources applications. To this end, we configure the Variable Infiltration Capacity (VIC) model coupled with the mizuRoute routing model in the Cautín River basin (2770 km²), Chile. We use the Latin hypercube sampling (LHS) method to generate 3500 different model parameters sets, for which basin-averaged runoff estimates are obtained directly (no routing or instantaneous runoff case) and are subsequently compared against outputs from four routing schemes (unit hydrograph, Lagrangian kinematic wave, Muskingum–Cunge, and diffusive wave) applied with five different routing time steps (1, 2, 3, 4, and 6 h). The results show that incorporating routing schemes may alter streamflow simulations at sub-daily, daily, and even monthly timescales. The maximum Kling–Gupta efficiency (KGE) obtained for daily streamflow increases from 0.64 (instantaneous runoff) to 0.81 (for the best routing scheme), and such improvements do not depend on the routing time step. Moreover, the optimal parameter sets may differ depending on the routing scheme configuration, affecting the baseflow contribution to total runoff. Including routing models decreases streamflow values in flood frequency curves and may alter the probabilistic distribution of the medium- and low-flow segments of the flow duration curve considerably (compared to the case without routing). More generally, the results presented here highlight the potential impacts of river routing implementations on water resources applications that involve hydrologic models and, in particular, parameter calibration.

1 Introduction

Hydrology and land surface models are powerful tools for characterizing the terrestrial water cycle, and they provide valuable information for water resources planning under future climate scenarios (Vano et al., 2012; Mendoza et al., 2016; Melsen et al., 2018; Chegwidden et al., 2019). In applications at the catchment scale or beyond, these models are typically used in combination with river routing models that convert instantaneous runoff into realistic streamflow estimates at any location in river networks (Oki and Sud, 1998; Olivera et al., 2000; Lucas-Picher et al., 2003). Hence, streamflow estimated by the river routing model is used for several water resources applications, including flood risk assessments (Wobus et al., 2017), ecosystem health evaluations (Qiu et al., 2021), short-term streamflow forecasting (e.g., Tang et al., 2007; Emerton et al., 2016), and reservoir operations (Salas et al., 2018; Shaad, 2018).

Over the past 3 decades, many river routing models have been developed and coupled with hydrology and land surface models (Shaad, 2018). The river routing models vary in terms of the modeling reservoir, irrigation, and other human interventions in relation to river water (e.g., Hanasaki et al., 2006); the spatial resolution and type of discretization of the river network, namely grid based vs. vector based (Lehner and Grill, 2013; Mizukami et al., 2016, 2021); and, finally, the representation of flow physical processes in equations (hereafter called routing schemes). The last category spans from a simple unit hydrograph method (Lohmann et al., 1996, 1998) to storage-based routing schemes such as the Muskingum scheme (David et al., 2011), simplifications of the Saint-Venant equations like the kinematic-wave scheme (Arora and Boer, 1999; Decharme et al., 2010; Ye et al., 2013; Thober et al., 2019) or the diffusive-wave scheme (Gong et al., 2009; Yamazaki et al., 2011), local inertia equations (Bates et al., 2010; Yamazaki et al., 2013), and full dynamic wave approaches (Paiva et al., 2011).

Given the wide range of routing methods available, it is crucial to understand the benefits and limitations of each method for the specific model application (Shaad, 2018). Many studies have conducted intercomparison experiments with focus on routing schemes to evaluate their impacts on streamflow simulations. For example, Arora et al. (2001) compared a time-evolving (or variable-velocity) algorithm that uses Manning's equation against a simple storage-based routing scheme (without using a momentum equation) operating at a very different horizontal resolution. Specifically, they concluded that the variable-velocity scheme can produce higher values of peak discharge. Gong et al. (2009) demonstrated the benefits of diffusive-wave routing over a linear reservoir-routing method to get more realistic time delays in hydrograph waves in a basin located in southern China. David et al. (2011) introduced the Routing Application for Parallel Computation of Discharge (RAPID) based on the traditional Muskingum method (McCarthy, 1938), obtaining improvements in terms of the root-mean-squared error (RMSE) and the Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) when compared to a lumped runoff scheme, which accumulates upstream instantaneous runoff without any delay. Ye et al. (2013) implemented a kinematic-wave-routing scheme in the Community Land Model (CLM) version 3.5 and obtained better results compared to the original grid-based river transport model (RTM), which uses storage-based routing, in two basins in China.

More recently, Zhao et al. (2017) compared daily and monthly streamflow simulations produced with the CaMa-Flood (Yamazaki et al., 2011) model – fed with daily runoff from nine global hydrological models (GHMs) – against those obtained with the same hydrological models and their native routing schemes (which have simpler physics). They concluded that the choice of routing scheme may have large effects on the simulated streamflow and peak values. ElSaadani et al. (2018) compared streamflow simulations obtained from Variable Infiltration Capacity (VIC) runoff outputs using RAPID and the Hillslope Link Model (HLM; Mantilla, 2007), which is based on power laws that relate flow velocity, channel discharge, and upstream area at many stream gauges located in the Cedar River basin, Iowa. They noted that the choice of routing scheme has large effects on simulated hydrographs, obtaining more realistic peak times and magnitudes with the HLM model and decreasing differences in performance for larger catchments. Siqueira et al. (2018) compared a local inertia scheme against a non-hydrodynamic scheme or storage-based routing, showing that the former provided slight improvements in terms of the NSE and the Kling–Gupta efficiency (KGE; Gupta et al., 2009) over the Amazon and La Plata river basins, especially in flow timing. They highlighted that the calibration of hydrological parameters and including hydrodynamic routing are critical elements for achieving realistic streamflow simulations in South America.

Besides the complexity of the routing scheme used, the choice of routing time step may also impact streamflow calculations (Shaad, 2018). Qiu et al. (2021) characterized the effects of such a decision on hydrological variables simulated with the Soil and Water Assessment Tool (Arnold et al., 1998), which uses the variable storage coefficient routing scheme, computing flow velocity with the Manning equation. The authors used six time steps ranging from 1 min to 1 d and assessed their impacts on performance skills, including NSE and bias, finding variations in streamflow simulations that were small compared to water storage and depth.

Although many past studies have shown that the choice of routing scheme affects streamflow simulations, efforts for improving their accuracy have been made by configuring hydrologic models and routing models independently. Hydrologists still focus on parameter calibration to improve discharge simulations, excluding river routing models or neglecting the potential impacts of river routing configurations (routing scheme and time step) and parameters if included (e.g., Beck et al., 2020; Newman et al., 2021). On the other hand, routing-model development and evaluation uses hydrologic model outputs that contain varying degrees of error, becoming especially difficult in large river basins or greater spatial domains (e.g., Mizukami et al., 2016; Zhao et al., 2017). Further, the key role of river routing parameters in reproducing observed streamflow characteristics has been previously recognized (Boyle et al., 2001; Butts et al., 2004; Sheikholeslami et al., 2021), highlighting the need for joint (i.e., hydrological and routing) parameter search strategies to characterize the benefits of routing configurations and potential compensatory effects in reproducing application-specific metrics.

In this paper, we seek to better understand the implications that the configuration of routing schemes may have when conducting hydrologic model calibration for water resources applications. To this end, we perform numerical experiments in the Cautín River basin (Araucanía, Chile) using the VIC model (Liang et al., 1994) and the vector-based routing model mizuRoute (Mizukami et al., 2016). Specifically, we disentangle the impacts of model parameters and different routing schemes (all implemented for five time steps) by combining a large sample of VIC simulations obtained from 3500 parameter sets and routing simulations with four different routing methods implemented in mizuRoute. Our end goal is to unravel how the choices of routing method and routing time step affect (i) streamflow simulated at different temporal resolutions; (ii) performance metrics; (iii) the selection of model parameters given a target calibration metric; (iv) simulated water balance and runoff partitioning (i.e., baseflow ratio); and (v) hydrological signatures used for decision making, including flood frequency curves and flow duration curves (FDCs). The results and conclusions drawn here reflect the impact that modeling decisions may have for water resources management.

2 Study domain and data

2.1 The Cautín River basin

The study domain is the Cautín River at Cajón (CatC) basin (Fig. 1), a sub-catchment of the Imperial River basin, located in the Araucanía region of Chile. The basin elevation ranges between 125 and 3104 m a.s.l., the catchment area is 2770 km², and the dominant land cover types are crop–pasture rotation (44 %) and native forest (40 %). Additionally, the basin is prone to rainfall-driven flood events during winter and, therefore, has been the subject of studies aimed at enhancing predictive capabilities (e.g., Mendoza et al., 2012).

Figure 1 (a) Location of the Cautín River at Cajón basin in Chile (CatC, 2770 km²). (b) Location of outlet and inner stream gauge stations (white circles) and contributing drainage areas (white lines). The inner stations are Muco at Muco bridge (MatPM, 651 km²), Collín at Codahue (CatCD, 259 km²), and Cautín at Rariruca (CatRR, 1305 km²). (c) Digital river network and sub-basin boundaries used in mizuRoute.

2.2 Hydrometeorological data

Daily precipitation and maximum and minimum temperatures are obtained from the CR2MET v2.0 dataset (Boisier et al., 2018, available at https://www.cr2.cl/datos-productos-grillados/, last access: January 2022), which covers continental Chile with a horizontal resolution of $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$ for the 1979–2020 period. In CR2MET, precipitation data were obtained with a statistical modeling framework that uses topographic descriptors and large-scale climatic variables (water vapor and moisture fluxes) from ERA5 (Hersbach et al., 2020) as predictors and observed daily precipitation from gauge stations as predictands. For maximum and minimum daily temperatures, additional variables from MODIS land surface products were added as predictors. Daily precipitation and temperature time series are disaggregated into hourly time steps using the sub-daily distribution provided by ERA5-Land (Muñoz-Sabater et al., 2021). Relative humidity, wind speed, and shortwave radiation are derived for the same horizontal-resolution grid by spatially interpolating ERA5-Land outputs. Longwave radiation was computed with the parameterization proposed by Iziomon et al. (2003) using CR2MET air temperatures disaggregated to hourly time steps using the ERA5-Land hourly distribution.

Daily streamflow data were obtained from five stations (Fig. 1) maintained by the Chilean Water Directorate (DGA, available at the (CR)² Climate Explorer https://www.cr2.cl/datos-de-caudales/, last access: December 2021). Similarly, hourly streamflow records for the CatC basin were obtained from the official DGA website (https://dga.mop.gob.cl/servicioshidrometeorologicos, last access: December 2021).

3 Models

3.1 Hydrological model

We use the VIC model (Liang et al., 1994) to simulate state variables and fluxes at a $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$ horizontal resolution. VIC is a semi-distributed physically based hydrological model that solves energy and mass balance equations. Precipitation can be partitioned into snowfall or rainfall, and both can be stored in the canopy. The maximum amount of water intercepted by the canopy is estimated using the leaf area index (LAI; Dickinson, 1984). The soil is represented by three layers controlling the infiltration (first soil layer) and the baseflow (third soil layer). For infiltration fluxes, VIC uses the Xinanjiang formulation (Zhao et al., 1980), assuming that the infiltration capacity varies within an area (Wood et al., 1992). Excess runoff is generated in those areas where precipitation exceeds the amount of the available soil moisture storage of the first soil layer. VIC assumes that drainage is driven by gravity using the formulation proposed by Brooks and Corey (1964). In this regard, water enters the cell only from the atmosphere; i.e., VIC does not consider lateral fluxes among grid cells. Baseflow is generated in the third (deepest) soil layer using a formulation proposed by Franchini and Pacciani (1991). The snowpack is represented by two layers, where the top layer is used for energy balance computations (Andreadis et al., 2009). The reader is referred to Liang et al. (1994) for more details.

Horizontal heterogeneity is considered in each grid cell by incorporating different land cover types. Here, we use the International Geosphere–Biosphere Programme (IGBP) classification for the year 2010 from the MCD12Q1 v006 land cover product (Sulla-Menashe and Friedl, 2018) to represent all land cover types spanning at least 2 % of each grid cell area. Mean monthly LAI values for these land cover types are derived from the MOD15A2 product. Soil bulk density is estimated using the mean value from the first 2 m of soil from the SoilGrids product (Poggio et al., 2021).

3.2 River routing schemes

The mizuRoute model first performs a hillslope routing using a gamma-distribution-based unit hydrograph to delay instantaneous total runoff from the VIC model to a catchment outlet and then routes the delayed runoff for each river reach in the order defined by the river network topology. Full descriptions of hillslope routing and general routing procedures are provided in Mizukami et al. (2016). mizuRoute originally included two channel-routing schemes: (1) kinematic-wave-tracking (KWT) routing and (2) impulse response function (IRF) routing, which is similar to the model of Lohmann et al. (1996), except for the fact that mizuRoute uses a reach-to-reach routing approach instead of the source-to-sink approach. Details of both routing schemes are also provided in Mizukami et al. (2016). Here, we implement in mizuRoute two additional routing schemes commonly used for many water resources applications: diffusive-wave routing (DW – see Appendix A) and Muskingum–Cunge routing (MC – see Appendix B). All the channel-routing schemes except IRF (which uses prescribed wave celerity and diffusivity) share two parameters: Manning's n roughness coefficient and channel width (assuming rectangular channel).

4 Experimental setup

Figure 2 summarizes the approach used in this paper, which consists of the following steps:

• a.

  Sample model (VIC and mizuRoute) parameter sets and obtain, for each one, streamflow times series with five routing schemes (including instantaneous runoff or no routing as the baseline) and five temporal resolutions (Fig. 2a; see details in Sect. 4.1) at each river gauge (Fig. 1b and Table 1). We save the parameter sets that maximize each performance metric – computed with a daily time step – at each stream gauge station for each combination of routing scheme and routing time step. For the KWT scheme, we select n values of 0.01 (default option), 0.03 (i.e., the spatially constant value used by Yamazaki et al., 2011), and 0.033 (which maximizes the KGE computed with daily streamflow).

• b.

  Examine the effect of routing-model configurations (i.e., routing schemes and time steps) on simulated daily hydrographs at Cautín River at Cajón (Fig. 2b.1) and analyze the impact of excluding the river routing process on simulated streamflow at annual, monthly, daily, and sub-daily time steps (Fig. 2b.2).

• c.

  Explore the overall impacts of routing-model decisions on performance metrics (Sect. 4.2) computed with different temporal resolutions (Fig. 2c.1). Then, we examine the sensitivity of the best metric value (achievable from the simulations with all the sampled parameter sets) in relation to the river routing configuration across sub-basins (Fig. 2c.2).

• d.

  Characterize the effects of the routing configuration on simulated annual water balance (specifically, the mean annual runoff ratio) and baseflow contribution (computed from VIC output) in relation to total runoff (Fig. 2d.1) and the selected parameter values (Fig. 2d.2).

• e.

  Analyze the effects of routing configurations on flood frequency curves (see details in Sect. 4.3, Fig. 2e.1) and daily FDCs (Fig. 2e.2).

Figure 2 Overview of the analysis framework used here. (a.1) VIC model simulations are conducted at hourly time steps for 3500 parameter sets; (a.2) each runoff time series at each grid cell is aggregated to four additional time steps (2, 3, 4, and 6 h), and these new time series are routed with four schemes to produce 3500 (VIC parameters) × 5 (time steps) × 5 (Inst + four routing schemes) modeling configurations. (b) We illustrate routing effects on (b.1) simulated hydrographs during fall 2008 and (b.2) simulated streamflow at various temporal resolutions. (c) For each configuration, we compute performance metrics and examine (c.1) the impacts of routing configuration on streamflow performance computed at various time steps using the 1 % best model runs and (c.2) improvements in performance metrics across other gauge points in the Cautín River basin. (d.1) We analyze simulated mean annual water balance and baseflow contribution to total runoff and (d.2) compare the best-parameter sets for each configuration in terms of their normalized values. (e) Finally, we analyze the implications of routing configurations for (e.1) flood frequency and (e.2) flow durations curves.

Table 1 Stream gauge stations in the Cautín River basin. Annual streamflow at each station was obtained from daily records for the period April 1985–March 2020.

The steps a–d are performed using observed and simulated data for the period April 2008–March 2012, whereas step e is conducted using simulations for the period April 1985–March 2020. All the steps but c.1 (Fig. 2) use VIC and mizuRoute parameter sets that maximize performance metrics (listed in Sect. 4.1) computed with simulated and observed daily flows. The following sub-sections provide complete descriptions of the parameter-sampling strategy, streamflow simulations, performance metrics, and the computation of flood frequency and flow duration curves.

4.1 Parameter sampling and streamflow simulations

Since we aim to examine the impacts of different routing schemes on streamflow performance metrics across the parameter space, rather than seeking an optimal parameter set, we use the Latin hypercube sampling (LHS) method, which is a common strategy for sampling the parameter space and identifying behavioral sets for specific target metrics (Andréassian et al., 2014; Broderick et al., 2016; Melsen et al., 2016, 2019; Guse et al., 2017; Khatami et al., 2019). Here, we sample 3500 model parameter sets (Fig. 2a) considering the 13 VIC parameters identified by Sepúlveda et al. (2022) to be the most sensitive and the routing model parameters: one (the Manning roughness coefficient) for the KWT, MC, and DW methods, and two for the IRF method (Table 2). For each parameter set, we run VIC and mizuRoute at hourly time steps for the period April 2006–March 2012. To generate streamflow simulations at 2, 3, 4, and 6 h time steps, we aggregate hourly VIC runoff and run the mizuRoute model with all routing schemes. For example, streamflow time series at a 3 h resolution are obtained from mizuRoute simulations using 3 h VIC runoff time series, which are computed by temporally aggregating 1 h VIC outputs at each grid cell. It should be noted that this step requires an assumption of the absence of non-linear processes in time within the hydrological model. The resulting VIC runoff time series are also used to compute streamflow by spatially averaging total runoff within each (sub-)basin without using a routing scheme (hereafter referred to as instantaneous runoff, Inst, or no routing), which is a common approach used in hydrological modeling applications (e.g., Mendoza et al., 2016; Beck et al., 2020). As a result, we obtain streamflow times series at each river reach (Fig. 1c) for five routing methods (four routing schemes and Inst as the baseline) and five temporal resolutions (1, 2, 3, 4, and 6 h; Fig. 2a).

Table 2 Model parameters sampled in this study.

4.2 Streamflow performance metrics

For each model run, we evaluate the performance of streamflow simulations from VIC-mizuRoute using four metrics. The first one is the Kling–Gupta efficiency (KGE; Gupta et al., 2009; Kling et al., 2012), which quantifies performance in terms of variability, volume, and timing:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\text{KGE}\left(Q\right)& =\mathrm{1}-\sqrt{(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}+(\mathrm{1}-\mathit{\beta }{)}^{\mathrm{2}}+(\mathrm{1}-r{)}^{\mathrm{2}}}\\ & \mathit{\alpha }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{s}}}{{\mathit{\sigma }}_{\mathrm{o}}}}\mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{s}}}{{\mathit{\mu }}_{\mathrm{o}}}}\end{array},\end{array}$$

where σ is the standard deviation for simulated and observed values, μ is the mean streamflow over the n times steps, and r is the Pearson correlation coefficient between simulated and observed streamflow. The second metric is the Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970), which is computed using observed (o) and simulated (s) streamflow (Q):

$$\begin{array}{}\text{(2)}& \text{NSE}\left(Q\right)=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n{\left(Q_{\mathrm{o}}^t-Q_{\mathrm{s}}^t\right)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^n{\left(Q_{\mathrm{o}}^t-{\stackrel{\mathrm{‾}}{Q}}_{\mathrm{o}}\right)}^{\mathrm{2}}}},\end{array}$$

where $Q_{\mathrm{o}}^t$ is the observed streamflow for time step t, $Q_{\mathrm{s}}^t$ is the simulated streamflow for time step t, and ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{o}}$ is the mean observed streamflow over the n time steps considered. The third metric is the NSE computed for the logarithms of the streamflow (NSE-log) to test the model's capability to simulate low flows (Krause et al., 2005). Although all these metrics range between −∞ and 1, where 1 represents a perfect simulation, their values are not comparable because they differ in terms of target streamflow characteristics and the incorporation or lack of a benchmark. For example, NSE is formulated based on a reference model simulation (i.e., mean climatology), while KGE does not have one, and NSE=0 is equivalent to $\text{KGE}=-\mathrm{0.41}$ (Knoben et al., 2019). Importantly, the above metrics are relevant not only for the hydrologic modeling community – especially for parameter calibration and evaluation (Fowler et al., 2018; Knoben et al., 2019; Clark et al., 2021) – but also for the river routing community. In fact, several examples of river routing scheme assessments can be found using the KGE (Pereira et al., 2017; Hoch et al., 2019; Qiao et al., 2019; Thober et al., 2019; Munier and Decharme, 2022), the NSE (Yamazaki et al., 2011; Ye et al., 2013; Zhao et al., 2017; ElSaadani et al., 2018; Nguyen-Quang et al., 2018; Fleischmann et al., 2019, 2020), and even the NSE-log (Paiva et al., 2013b).

Finally, we use the percent bias in the high-segment volume of the FDC (Yilmaz et al., 2008):

$$\begin{array}{}\text{(3)}& \mathit{\%}\text{BiasFHV}={\displaystyle \frac{{\sum }_{h=\mathrm{1}}^H(Q_{\mathrm{s}}^h-Q_{\mathrm{o}}^h)}{{\sum }_{h=\mathrm{1}}^H{Q}_{\mathrm{o}}^h}}\times \mathrm{100},\end{array}$$

where $h=\mathrm{1},\mathrm{2},\mathrm{\dots },H$ are the flow indices in the flow array with a probability of exceedance of less than 0.02. FDC high-segment volume (FHV) is a measure of the basin response to high-precipitation events.

The four performance metrics are calculated for the period April 2008–March 2012 (after a 2-year warm up) using all the combinations of parameter sets (3500), routing schemes (including the case without routing), and routing time steps (1, 2, 3, 4, and 6 h). Additionally, for each routing time step, the performance metrics are computed for different aggregated temporal resolutions when possible for step c.1. For example, to estimate metrics at an hourly time step, routing can only be run at a 1 h time step. Metrics computed at 3 hourly time steps use temporally averaged streamflow from 1 and 3 h mizuRoute simulations. Metrics computed at 6 hourly time steps can be computed from temporally averaged 1, 3, and 6 h mizuRoute simulations and so on. The observed streamflow for a given time step is estimated from hourly streamflow records.

4.3 Flood frequency and flow duration curves

Because high flows are relevant for engineering applications, particularly infrastructure design, we analyze the implications of routing configurations for the calculation of flood frequency curves (Fig. 2e.1). To this end, we run VIC at hourly time steps from April 1981 to March 2020 using the parameters associated with the highest KGE, NSE, and %BiasFHV values (all computed with daily flows) for each routing configuration. Then, hourly VIC total runoff is aggregated and routed at different time steps (i.e., 2, 3, 4, and 6 h), and annual maximum daily flows are obtained for the period April 1985–March 2020 (i.e., the period April 1981–March 1985 is dropped). Hence, for each routing time step, we obtain five annual time series with n=35 values (obtained from the baseline and the four routing schemes) that are used to compute maximum daily flows at return periods of 20, 50, 100, and 200 years. We use the log-normal parametric distribution – which provides the best results for the Kolmogorov–Smirnov test – for the observed time series of maximum daily flows.

Finally, we characterize the impacts of routing configurations on daily FDCs (Fig. 2e.2), which are widely used in water resources applications. Empirical FDCs are constructed from daily time series of streamflow for the period April 1985–March 2020, computing exceedance probabilities with the Weibull plotting position formula.

5 Results

5.1 Illustration of routing effects

Figure 3 illustrates the sensitivity of daily streamflow simulations to different routing-model decisions, including the routing time step of the IRF scheme (Fig. 3a), the choice of routing scheme (Fig. 3b), and the Manning's roughness coefficient of the KWT results (Fig. 3c). In each panel, simulations are displayed for the period 10 May–16 June 2008 using the parameter set (obtained from LHS) associated with the maximum daily KGE for each combination of routing scheme, routing time step, and Manning's coefficient values. In all cases, sub-daily routing simulations are aggregated to a 24 h time step. It can be noted that the choice of routing scheme and the Manning's coefficient values have a larger effect on the shape of the flood wave. Additionally, increasing routing time steps for IRF accelerates the timing of peak discharge by 1 d, though this decreases its value to 776 m³ s⁻¹ for Δt=6 h (compared to 785 m³ s⁻¹ obtained with Δt=1 h). The choice of routing scheme affects the shape of storm hydrographs, especially high flows. Finally, a delay in peak-flow simulations is obtained for larger Manning's roughness coefficients.

Figure 3 Time series with daily streamflow at Cautín at Cajón, obtained from hourly VIC runoff outputs routed with different mizuRoute configurations. (a) Application of the Impulse Response Function (IRF) with five different routing time steps, (b) effects of different routing schemes using Δt=6 h, and (c) effects of the Manning's roughness coefficient (n) when applying the kinematic-wave routing scheme with Δt=6 h (see text for details). In panel (c), the orange line (n=0.033) is associated with the parameter set that maximizes the KGE computed with daily streamflow.

Figure 4 compares streamflow obtained from mizuRoute (y axis) against instantaneous runoff (x axis) for several temporal resolutions and different routing schemes. In this case, the parameter set used to run the models is the one that maximizes the KGE among the 3500 parameter sets from the LHS. The results for hourly time steps show that the lack of routing yields much larger values (>1300 m³ s⁻¹ in some cases) compared to routed streamflow. These differences are gradually reduced when the routing time step increases to Δt=3 and 6 h, although differences can be larger than 1200 m³ s⁻¹. The impact of excluding routing lessens as the time step increases, yet it can be important even for streamflow simulations at a Δt=24 h time step. At monthly time steps, the differences between routed and instantaneous runoff are reduced considerably, though these can still be as large as 27 m³ s⁻¹ (i.e., a 10 % difference using routed runoff as the reference). Further, the differences become negligible at the annual resolution. Finally, given a specific time step, the magnitudes of differences are very similar across routing schemes, although slight differences in r² suggest that IRF and KWT affect VIC outputs more.

Figure 4 Simulated streamflow (VIC+mizuRoute) vs. instantaneous VIC runoff at Cautín at Cajón for the period April 2008–March 2012 using different time steps (rows) and routing schemes (columns): instantaneous runoff (Inst), impulse response function (IRF), kinematic-wave tracking (KWT), Muskingum–Cunge (MC), and diffusive wave (DW). Mean yearly, monthly, daily, and 12 h streamflows are computed from temporally averaged 6 h values, while 1, 3, and 6 h streamflows are obtained from mizuRoute simulations with 1, 3, and 6 h time steps, respectively. The 1:1 line is displayed in red with the coefficient of determination (R²).

5.2 Effects on performance metrics

The KGE, NSE, NSE-log, and %BiasFHV values obtained with the 1 % best (i.e., 35) parameter sets for each combination of routing scheme, routing time step, and metric time step are displayed in Fig. 5. To compare performance measures from different configurations, simulations were aggregated to the metric time resolution (columns). For the sake of brevity, we only show 1 and 24 h metric time steps here (the full results can be found in Fig. S1 in the Supplement). Overall, the results show a clear difference between including routing and Inst, which becomes more evident for performance metrics computed at smaller temporal resolutions. Moreover, none of the 1 % best-parameter sets for KGE and NSE with instantaneous runoff could produce better performance than the inclusion of routing. On the other hand, the choice of routing time step is comparatively less influential for a given metric time step. The maximum KGE spans 0.69–0.73 for instantaneous runoff, increasing to 0.8 or more when routing is included. Similar improvements are observed for NSE, with increments that can be larger than 0.3 in terms of NSE values (e.g., 1 h). Compared to the former metric, routing yields smaller benefits for NSE-log and less noticeable differences among routing configurations, mainly due to the minor influence of high-flow values on the metric. In all cases, a larger spread in high-flow biases is obtained with instantaneous runoff (compared to routing schemes), indicating that many VIC parameter sets do not compensate for the lack of river routing to obtain accurate high-flow simulations. Finally, the results show that the impact of representing river routing on %BiasFHV is more relevant when this metric is computed at an hourly resolution, approaching zero as the metric time step increases.

Figure 5 Impact of routing scheme and routing time step on performance metrics (rows) for the period April 2008–March 2012 at Cautín at Cajón, computed with 1 and 24 h discharge temporal resolutions (columns) and the 1 % best-parameter sets among those obtained through Latin hypercube sampling (see text for details). The results are presented for instantaneous runoff (Inst), impulse response function (IRF), kinematic-wave tracking (KWT), Muskingum–Cunge (MC) and diffusive wave (DW).

Figure 6 compares the best KGE, NSE, NSE-log, and %BiasFHV values (computed from daily flows) achievable from the large sample of model parameter sets in each basin (represented by the basin area in the x axis) given a specific combination of routing schemes and routing time steps. For completeness, the KGE components (α, β, and r) are also displayed. For KGE, NSE, and NSE-log, the maximum values increase at all streamflow gauges when the routing process is included regardless of the routing configuration. Very similar maximum KGE values are obtained with the four schemes implemented in mizuRoute, and the differences among these schemes are generally lower than 0.05 for all time steps and basins. The improvements in KGE through the inclusion of routing are explained by the enhancement of temporal correlation (r) and variability error (α). In particular, routing (especially Muskingum–Cunge and diffusive-wave schemes) helps to improve the r values of simulated instantaneous runoff by changing the timing of high peak flows. Even more so, Fig. 6 shows that, in our study domain, the correlation between streamflow simulations and observations increases with basin area when routing is incorporated.

Figure 6 Best metric value obtained with daily flows at each stream gauge station for a given performance metric (rows), routing time step (columns), and routing scheme for the period April 2008–March 2012. For completeness, the Kling–Gupta efficiency (KGE) components associated with the best KGE value are included. The results are presented for instantaneous runoff (Inst – red symbols), IRF, KWT, MC, and DW.

Figure 6 also shows considerable improvements in NSE across all catchments when routing is applied, especially in CatRR and CatC (i.e., the two largest river basins). Notably, differences between routing and Inst options are also obtained for NSE-log in all stream gauges, demonstrating the benefits of routing beyond the simulation of high flows. %BiasFHV is very close to zero for the smallest catchment, increasing with catchment size when no routing is performed; nevertheless, in every basin, at least one routing scheme can reduce high-flow biases to nearly zero.

5.3 Impact on simulated fluxes and model parameters

For the remainder of this paper, we illustrate the impacts of routing time steps by showing results for 1, 3, and 6 h (the full results are available in the Supplement). Figure 7 illustrates the impacts of the choice of routing scheme on the mean annual runoff ratio (partitioning of precipitation into runoff and evapotranspiration – x axis) and the ratio between mean annual baseflow and mean annual total runoff (y axis) for each routing time step (columns) and performance metric (all computed with daily discharge and displayed in different rows). The red symbols represent the results based on the best value for the corresponding performance metrics. To account for equifinality effects, we also include the results with the parameter sets that produce the 0.1 % best simulations for the respective performance metrics (i.e., four parameter sets). Precipitation partitioning ($Q/P$) is relatively unchanged whether routing is used or not, regardless of the routing schemes for any performance metric except %BiasFHV, for which less clear patterns in runoff ratio and flow components are obtained. Conversely, the greater impacts are seen in the baseflow ratio ($\text{bf}/Q$). For KGE and NSE, a clear separation in runoff components is obtained between instantaneous runoff and routing schemes (for any routing time step), with a much higher baseflow contribution to total runoff for Inst. When the VIC model parameters are selected based on the NSE-log metric, the differences among the routing configuration options are generally smaller compared to KGE and NSE, though the best-performing parameter set produces less baseflow.

Figure 7 Effects of performance metric (rows), routing time step (columns), and routing scheme on simulated mean annual water balance (characterized with the annual runoff ratio – x axis) and the baseflow ratio (y axis) obtained for the 0.1 % best-parameter sets at Cautín at Cajón (period April 2008–March 2012). The results are presented for instantaneous runoff (Inst), impulse response function (IRF), kinematic-wave tracking (KWT), Muskingum–Cunge (MC) and diffusive wave (DW). In each panel, the results obtained with the parameter set (among the 3500 samples) that maximizes each metric are displayed in red, and the results from a small ensemble (n=4) with the best 0.1 % VIC parameter sets are displayed in gray.

To examine the effect of river routing on the selection of the model parameters, we show, for three routing time steps and all routing schemes, the best values for KGE, NSE, NSE-log, and %BiasFHV (with the four metrics computed with daily flows) among the 3500 parameter sets from LHS (Fig. 8). The parameter values are normalized by the difference between the maximum and minimum values obtained from LHS to facilitate comparisons. Hence, a normalized value of zero indicates the lower boundary of the parameter, while a value of 1 indicates the upper limit. The results indicate that the same best-parameter set is obtained for NSE-log regardless of the selected routing scheme (when included) or the routing time step, which explains why all the routing schemes with the best-performing parameter set (red) produce the same $\text{bf}/Q$ and $Q/P$ (Fig. 7). Additionally, the absence of routing when maximizing NSE-log not only affects soil parameter values (compared to results with river routing) but also reduces the snow albedo decay parameter (α_thaw).

Figure 8 Normalized parameter values for Cautín at Cajón associated with the best performance metric (period April 2008–March 2012) obtained from the 3500 parameter sets produced with the Latin hypercube sampling given a combination of routing schemes and routing time steps. The symbols representing VIC parameters are linked with straight lines. The results are presented for daily instantaneous runoff (Inst), impulse response function (IRF), kinematic-wave tracking (KWT), Muskingum–Cunge (MC), and diffusive wave (DW).

For KGE, the KWT, MC, and DW schemes yield the same VIC parameter sets, which are different from the parameter set obtained with instantaneous runoff. Conversely, all the routing schemes yield different VIC parameter sets for NSE, except for MC and DW. For KGE and NSE, the choice of routing time step does not affect the best-parameter set given a routing scheme, with the exception of the combination of KGE–IRF. It should be noted that, for both NSE and KGE, excluding routing produces higher W_s (fraction of maximum soil moisture where non-linear baseflow occurs), higher d_max (maximum total soil thickness), and lower b (infiltration parameter) regardless of the routing time step, which may contribute to the larger baseflow fraction seen in Fig. 7. Finally, different routing configurations (routing methods and routing time step) result in unique best-parameter sets if %BiasFHV is used as the performance metric, which is in contrast to KGE, NSE, and NSE-log.

5.4 Implications for flood frequency and flow duration curves

Figure 9 shows the flood frequency curves from the annual time series of maximum daily flows using model parameters obtained with KGE, NSE, and %BiasFHV as target metrics. Note that the curve for daily instantaneous runoff is the same for each metric regardless of the time step. As expected, differences in flood estimates between instantaneous runoff and any other routing scheme are considerable, surpassing 400 m³ s⁻¹ in some cases (see results for T=200 years with KGE and %BiasFHV). Additionally, the dispersion among routing schemes increases with larger Δt for KGE and, in particular, for %BiasFHV, which can be explained by differences obtained for model parameter values (Fig. 8). Even if the same parameter sets are obtained for a target metric and a suite of routing schemes, regardless of the choice of Δt (e.g., KWT, MC, and DW for KGE 00 see Fig. 8), variations in routing time step affect the time series of daily flows and, therefore, flood quantiles obtained with a specific routing scheme, as well as inter-method differences for a given return period. For example, when the target metric is KGE, Q(T=100 years) values obtained with KWT and DW for Δt=1 h are 1202 and 1189 m³ s⁻¹, respectively, while Δt=6 h yields 1196 and 1171 m³ s⁻¹ using the same schemes. Interestingly, differences in frequency curves arising from the choice of routing scheme with NSE decrease with larger Δt, although the best-parameter set does not depend on the routing scheme (see Fig. 8).

Figure 9 Frequency curves for annual maximum daily flows (y axis) in Cautín at Cajón, derived from numerical simulations conducted with different routing schemes, routing time steps (columns), and performance metrics (rows). All frequency curves are computed from annual time series of n=35 annual maximum daily flows (April 1985–March 2020) using a log-normal density function. The results are presented for instantaneous runoff (Inst), impulse response function (IRF), kinematic-wave tracking (KWT), Muskingum–Cunge (MC), and diffusive wave (DW).

Figure 10 shows daily FDCs obtained with different routing schemes; routing time steps (columns); and model parameters that maximize KGE, NSE, or %BiasFHV (rows). It can be noted that, for KGE, the disagreement arising from the choice of routing scheme is generally small for medium and low flows as opposed to the disagreement obtained with NSE and %BiasFHV. For KGE and NSE, no appreciable differences in FDCs are observed among routing time steps, whereas the opposite is observed for %BiasFHV. For example, FDCs with very different mid-segment slopes (which is a signature for flashiness of runoff) and low-flow volumes (i.e., segment for P_exc>70 %) are obtained for MC with Δt=3 and Δt=4 h (see Fig. S4 in the Supplement). This can be explained by the fact that the choice of routing time step does not impact the parameter sets obtained with NSE and KGE (as it happens with %BiasFHV).

Figure 10 Mean daily flow duration curves for the period April 1985–March 2020 in Cautín at Cajón derived from different routing schemes, routing time steps (columns) and performance metrics (rows). The results are presented for daily instantaneous runoff (Inst), Impulse Response Function (IRF), kinematic-wave tracking (KWT), Muskingum-Cunge (MC) and diffusive-wave (DW) routing schemes.

6 Discussion

6.1 Implications for hydrological modeling

In this paper, we use the LHS approach to evaluate the impact of routing on streamflow performance metrics across the parameter space. Our results suggest that, regardless of the routing scheme, including the river routing process improves the overall streamflow performance (Figs. 5 and 6). In other words, the lack of river routing in the modeling chain may not be fully compensated for through the calibration of hydrologic model parameters. Nevertheless, such a conclusion may depend on the hydrological regime of the catchment and the distributed spatial configuration of the river routing implementation. The Cautín River basin has a rainfall-dominated runoff regime, with high-flow peaks associated with heavy rainfall events during the winter season and a slight influence of snowmelt during the spring season. The catchment response time and peak discharge depend on the runoff-routing process in the river network; hence, its explicit inclusion in hydrological modeling may yield better results, especially for performance metrics influenced by high flows (Clark et al., 2021).

The results presented here show that the implementation of river routing is also relevant for medium and low flows. For example, including river routing provided higher values for NSE-log (Figs. 5 and 6) – improving the simulation of low discharges – and modified the shape of the middle- and low-flow segments in the FDC (Fig. 10), which are characteristic signatures of flashiness in runoff response and long-term baseflow, respectively (Yilmaz et al., 2008). The effects of river routing are also reflected in the partitioning of total runoff between baseflow and surface runoff. In fact, the results presented here show that the parameter search process compensates for the lack of routing by modifying other fluxes and state variables (Khatami et al., 2019) to increase streamflow-oriented performance metrics. For example, when the target metrics are KGE, NSE, or NSE-log, excluding routing (represented by squares) forces VIC to compensate for the absence of this process by delaying the runoff response with a larger contribution of baseflow to total runoff compared to any routing schemes. In such cases, the contribution of baseflow to total runoff increases by >30 % when river routing is excluded, which is achieved by modifying soil parameters – including W_s, one of the most sensitive to baseflow in this type of hydroclimate (Sepúlveda et al., 2022) – to delay the streamflow response. This result suggests that including routing processes may impact the outcomes from drought-oriented studies since baseflow is the primary flux sustaining streamflow during water scarcity periods (Karki et al., 2021).

Conversely, we found smaller variations in the partitioning of precipitation between evapotranspiration and runoff in the absence of river routing (Fig. 7), especially for KGE, NSE, and NSE-log. Hence, when models are configured to maximize these metrics to conduct hydroclimatic or water balance analyses at the annual timescale, the incorporation of routing processes is relatively less important.

The impacts of routing-scheme choice exhibit less clear patterns if the model chain is calibrated with integrated time series metrics such as KGE, though differences remain in the performance metrics (Fig. 5), high-flow analyses (Fig. 9), and FDCs (Fig. 10). Here, we obtain very similar results with MC and DW, including runoff partitioning, best-parameter sets (including Manning's roughness coefficient n), and flood frequency curves if the target metrics are KGE and NSE. Both schemes use the same routing parameters and essentially simulate wave diffusion. MC mimics physical diffusive phenomena via numerical diffusion in the explicit numerical solution, while the DW routing explicitly incorporates the diffusion process in the diffusion equation. Despite the fact that very good results were achieved with DW, limited impacts are expected for the Cautín River basin because the slopes of the river reaches therein range from 0.0004 to 0.274 m m⁻¹. In fact, the largest benefits of DW are expected for flatter river systems (slope<0.001m m⁻¹; e.g., Kazezyılmaz-Alhan et al., 2007), where flood wave diffusion processes can dominate. It can be argued that a more physically realistic routing scheme will better simulate the hydrograph. However, IRF routing, which uses the simplest algorithm among the schemes used in this study, may reproduce the results from the other routing schemes after the calibration.

Although the routing time step does not yield important effects on performance metrics in our experimental setup, it can affect the choice of VIC parameters (e.g., see results for KGE and %BiasFHV – Fig. 8), which is in line with previous hydrologic modeling research. For example, Kavetski et al. (2011) found that the temporal data resolution may alter parameter values in conceptual hydrological models. More recently, Melsen et al. (2016) found that the parameter values may greatly vary if calibration metrics are computed at hourly, daily, or monthly time steps. Accordingly, variations in the VIC model time step – which is fixed to Δt=1 h here – may also alter the selection of parameters and performance measures (see Sect. 6.2).

6.2 Limitations and future work

Here, we only focused on the choice of routing scheme and routing time step, though there are many other decisions that could be explored in the implementation of river routing models. For example, we did not examine the effects of surface storage elements like reservoirs, wetlands, and flood plains on river flow dynamics. Additionally, we did not estimate the spatial variability of the Manning's roughness coefficient (n) across the Cautín River basin. Due to data restrictions, many past studies used spatially constant values of n (Arora and Boer, 1999; Lucas-Picher et al., 2003; Yamazaki et al., 2011; Siqueira et al., 2018) or adopted indirect approaches. For instance, Decharme et al. (2010) estimated n as a linear function of the river width W; Miguez-Macho and Fan (2012) used satellite land cover to assign the Manning's roughness coefficient; and Verzano et al. (2012) estimated n variability in space based on topography, the location of the urban population, and river sinuosity. These or other techniques could be applied in combination with field data to estimate spatial n fields that can be subsequently calibrated through spatial regularization strategies (e.g., Mendoza et al., 2012).

An important assumption of our experiment is the lack of non-linear processes in time within the hydrological model in order to aggregate hourly runoff to coarser time steps. Such a decision was required to isolate the impact of the hydrologic model configuration from river routing decisions and to achieve a clean experimental design, though we recognize that the choice of hydrological model time step may also alter performance metrics (e.g., Bruneau et al., 1995; Wang et al., 2009).

In this study, we did not include the full dynamic wave scheme, which might yield improvements compared to the routing schemes tested here, especially for very large flood events downstream of the bases or for the flatter parts of basins. Paiva et al. (2013b) validated a full hydrodynamic model in stream gauges within the Amazon River basin, finding that discharge and water levels were simulated accurately, outperforming the Muskingum–Cunge approach. The same model was evaluated against satellite observations, showing good performance in terms of water levels and inundation extents (Paiva et al., 2013a). Hence, future assessments of routing schemes may include more detailed comparisons against remotely sensed data, including such additional hydraulic variables over ungauged catchments with different hydrological regimes (e.g., snowmelt-driven and mixed regimes) and physiographic characteristics (e.g., contributing area, average slope, land cover types). Further, it would be interesting to examine the interplay between structural uncertainty (e.g., channel geometries, river network or drainage density, and floodplain) and parametric uncertainty in river routing models, a topic that has been widely explored in the hydrologic modeling literature (e.g., Ajami et al., 2007; Günther et al., 2019; Newman et al., 2021).

In this work, we contrast the performance of different model configurations using metrics that describe specific aspects of the system's response. However, the interpretation of differences in metric values is not straightforward (Clark et al., 2021). For example, an improvement of 0.2 in KGE may be explained by a better simulation of streamflow timing and variability at the cost of larger volume bias (e.g., see results for Cautín at Rariruca, 1305 km² – Fig. 6), which stresses the need to understand the information content in the metrics used for model diagnostics. Finally, we only considered a single-objective (e.g., NSE, KGE) parameter search based on a Monte Carlo sampling scheme. Future studies could characterize the impacts of river routing schemes exploiting single-objective optimization algorithms (e.g., Duan et al., 1992; Tolson and Shoemaker, 2007) or address multi-objective problems using Pareto principles (Yapo et al., 1998; Pokhrel et al., 2012; Shafii and Tolson, 2015). Although different behavioral parameter sets and, therefore, different internal fluxes could be obtained, we hypothesize that similar conclusions could be drawn regarding the benefits of river routing representation to achieve realistic streamflow simulations. Nevertheless, further research is needed to understand implications for catchments with different hydroclimatic regimes and physiographic characteristics (e.g., Murillo et al., 2022; Muñoz-Castro et al., 2023).

7 Conclusions

Despite the general consensus in the hydrology and Earth system modeling communities about the relevance of river routing schemes for realistic streamflow simulations, there is little knowledge of the extent to which this process is relevant. Additionally, hydrologic model calibration research has been done in a way that neglects the impacts of river routing model configurations, and routing model development has been conducted in a way that ignores the effects of hydrologic model parameters. In this paper, we try to reduce these gaps by performing modeling experiments at the Cautín River basin (Chile), coupling the VIC model with four different routing schemes implemented in the mizuRoute model to produce streamflow simulations at various time steps with an ensemble of 3500 parameter sets.

Our main conclusions are as follows:

1. Runoff routing alters streamflow simulations considerably at sub-daily and daily time steps, with slight (negligible) impacts at the monthly (annual) time step.

2. Including a river routing model may provide better hydrologic model calibration results compared to the case without routing.

3. The timing of streamflow simulations may improve for larger contributing areas if runoff routing is performed.

4. For popular performance metrics (i.e., KGE, NSE, and NSE-log), including routing processes may yield different parameter sets compared to the case without routing, with notable impacts on the baseflow contribution to total runoff. Additionally, different routing schemes may yield different hydrologic parameter sets.

5. Including routing models decreases annual maximum daily flow values in frequency curves and, depending on the target streamflow metric, the disagreement in flood quantile estimates among schemes may increase for larger routing time steps.

6. When the calibration metric is NSE(Q_24h) or %BiasFHV, including routing models may affect the probabilistic distribution of medium and low daily flows considerably.

Appendix A: Diffusive-wave routing

The flood wave propagation through a river channel is described with the one-dimensional Saint-Venant equations:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle \frac{\partial Q}{\partial x}}+{\displaystyle \frac{\partial A}{\partial t}}=\mathrm{0}\text{(A2)}& {\displaystyle }& {\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)+gA{\displaystyle \frac{\partial Z}{\partial x}}-gA(S_{\mathrm{o}}-S_{\mathrm{f}})=\mathrm{0},\end{array}$$

where Q is discharge (L³T⁻¹) at time step t (T) and location x (L) in a river reach, A is cross-sectional flow area (L²), Z is flow depth (L), S₀ is channel slope (–), S_f is friction slope (–), and g is gravitational constant (LT⁻²). The continuity equation (Eq. A1) assumes that no lateral flow is added to a channel segment. A friction slope is expressed using channel conveyance K_c:

$$\begin{array}{}\text{(A3)}& S_{\mathrm{f}}={\displaystyle \frac{Q\left|Q\right|}{K_{\mathrm{c}}}}.\end{array}$$

In large-domain river routing, one-dimensional full Saint-Venant equations, or fully dynamic wave equations, are typically simplified by neglecting some force terms in the momentum equation (Eq. A2). The kinematic-wave approximation is obtained by neglecting acceleration and pressure gradient terms, assuming that river bed slope and energy slope are equal. This assumption is the basis of the kinematic-wave-tracking algorithm (Mizukami et al., 2016). If a rectangular channel with a channel width w is used, the diffusive-wave equation can be obtained by neglecting acceleration terms (first and second terms in Eq. A2) and by combining Eqs. (A1) and (A2) (Sturm, 2021):

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial Q}{\partial t}}=D{\displaystyle \frac{{\partial }^{\mathrm{2}}Q}{\partial x^{\mathrm{2}}}}-C{\displaystyle \frac{\partial Q}{\partial x}}\end{array}$$

Where:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }C={\displaystyle \frac{\mathrm{1}}{K_{\mathrm{c}}}}{\displaystyle \frac{\mathrm{d}{K}_{\mathrm{c}}}{\mathrm{d}A}}={\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}A}}\\ {\displaystyle }& {\displaystyle }D={\displaystyle \frac{K_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}qw}}={\displaystyle \frac{Q}{\mathrm{2}w{S}_{\mathrm{o}}}},\end{array}$$

where K_c is conveyance; and parameters C and D are wave celerity (LT⁻¹) and diffusivity (L² T⁻¹), respectively.

To solve the diffusive-wave equation for discharge Q, Eq. (A4) is discretized using weighted averaged finite-difference approximations across two time steps in space (i.e., second-order central difference in the first term in Eq. (A4) and first-order central difference for the second term in Eq. (A4)). The resulting discretized diffusive-wave equation is

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}(\mathit{\alpha }{C}_a-\mathrm{2}\mathit{\beta }{C}_d)& \cdot Q_{j+\mathrm{1}}^{t+\mathrm{1}}+\left(\mathrm{2}+\mathrm{4}\mathit{\beta }{C}_d\right)\cdot Q_j^{t+\mathrm{1}}\\ & -\left(\mathit{\alpha }{C}_a+\mathrm{2}\mathit{\beta }{C}_d\right)\cdot Q_{j-\mathrm{1}}^{t+\mathrm{1}}\\ =& -\left[\left(\mathrm{1}-\mathit{\alpha }\right)C_d-\mathrm{2}\left(\mathrm{1}-\mathit{\beta }\right)C_d\right]\cdot Q_{j+\mathrm{1}}^t\\ & +\left[\mathrm{2}-\mathrm{4}\left(\mathrm{1}-\mathit{\beta }\right)C_d\right]\cdot Q_j^t\\ & +\left[\left(\mathrm{1}-\mathit{\alpha }\right)C_a+\mathrm{2}\left(\mathrm{1}-\mathit{\beta }\right)C_d\right]\cdot Q_{j-\mathrm{1}}^t\\ C_a={\displaystyle \frac{C\cdot \mathrm{\Delta }t}{\mathrm{\Delta }x}}& C_d={\displaystyle \frac{D\cdot \mathrm{\Delta }t}{{\left(\mathrm{\Delta }x\right)}^{\mathrm{2}}}}\end{array},\end{array}$$

where α is the weight factor for the first-order space difference approximation of the second term in Eq. (A4), and β is a weight factor for the second-order space difference approximation of the first term in Eq. (A4). If both weights are set to 1, the finite difference becomes a fully implicit scheme, while setting both weights to zero results in a fully explicit scheme.

If internal nodes are defined within each reach (here we used five), Eq. (A5) becomes a system of linear equations that can be expressed in tridiagonal matrix form and solved with the Thomas' algorithm. In this paper, we use a fully implicit finite-difference approximation (i.e., $\mathit{\alpha }=\mathit{\beta }=\mathrm{1}$). The solution of the implicit method requires downstream and upstream boundary conditions, being the latter inflow from upstream reaches. We use the Neumann boundary condition, which specifies the gradient of discharge between the current and downstream reaches. Note that in diffusive-wave routing, celerity (C) and diffusivity (D) are updated at every time step based on the discharges (Q) and flow area (A) as opposed to IRF routing, in which celerity and diffusivity are provided as model parameters.

Appendix B: Muskingum–Cunge

In the Muskingum–Cunge (MC) routing approach, the desired streamflow value is computed as the weighted (C₁, C₂, and C₃) average of known discharge values at upstream and downstream positions at current and previous time steps:

$$\begin{array}{c}\text{(B1)}& {\displaystyle }{Q}_{j+\mathrm{1}}^{t+\mathrm{1}}=C_{\mathrm{1}}\cdot Q_j^t+C_{\mathrm{2}}\cdot Q_{j+\mathrm{1}}^t+C_{\mathrm{3}}\cdot Q_j^{t+\mathrm{1}},{\displaystyle }{C}_{\mathrm{1}}={\displaystyle \frac{\mathrm{2}KX+\mathrm{\Delta }t}{\mathrm{2}K(\mathrm{1}-X)+\mathrm{\Delta }x}},C_{\mathrm{2}}={\displaystyle \frac{\mathrm{2}K(\mathrm{1}-X)-\mathrm{\Delta }t}{\mathrm{2}K(\mathrm{1}-X)+\mathrm{\Delta }x}},\\ {\displaystyle }{C}_{\mathrm{3}}={\displaystyle \frac{-\mathrm{2}KX+\mathrm{\Delta }t}{\mathrm{2}K(\mathrm{1}-X)+\mathrm{\Delta }x}}.\end{array}$$

The parameters K and X are defined as

$$K={\displaystyle \frac{\mathrm{\Delta }x}{C}},X=\mathrm{0.5}-{\displaystyle \frac{Q}{\mathrm{2}{S}_{\mathrm{o}}C\mathrm{\Delta }x}}.$$

Here, both parameters are computed with discharge Q updated at every time step based on the average of the inflow at the current time step and the inflow and outflow at the previous time step. Note that celerity is also a function of discharge. Since Muskingum–Cunge is an explicit method, the routing time step can affect the numerical stability of the solution. To stabilize the solution, the sub-routing time step is determined at every simulation step so that the Courant condition ($C\cdot \mathrm{d}T/\mathrm{d}x$, where C is the wave celerity (L T⁻¹), dT is the routing time step (T), and dx is the channel length (L)) is less than unity.