An improved routing algorithm for a large-scale distributed hydrological model with consideration of underlying surface impact

Climate change has received widespread attention from the scientific community and the public, as it has caused a series of serious problems, including natural disasters and extreme climate events (Kharin et al. 2018; Chen et al. 2019; Lu & Qin 2020; Ragettli et al. 2020). As tools for estimating regional and global water resources and predicting the hydrological response to climate change, large-scale hydrological models have become a hot topic recently (Müller Schmied & Döll 2017; Wang et al. 2018; Abou Rafee et al. 2019; Li et al. 2019). A large-scale routing method is an important part of a large-scale hydrological model, with which a distributed hydrological model can simulate the change in runoff at temporal and spatial scales (Beven 2001). Many routing methods based on a digital elevation model (DEM) have been developed (Wen et al. 2012; Li et al. 2014; Lu et al. 2015; Huang & Lee 2017; Ling et al. 2018; Fan et al. 2019). The main steps of runoff routing in distributed hydrological models are (1) extracting the flow network and (2) calculating the runoff routing according to the flow path (Wen et al. 2012). Generally, runoff routing based on a DEM is calculated in one of two ways (Olivera et al. 2000): cell-to-cell and source-to-sink. The cell-to-cell method calculates the water movements from each grid to its adjacent downstream grid until the outlet grid of the basin is reached. The outflow of each cell is calculated based on the inflow and the river channel's storage. The source-to-sink method, under the assumption that the water is rigid, directly calculates the water's movement from the grid where runoff is generated to the outlet of the basin.

The first of these two kinds of routing methods is based on the mass conservation equation and the channel storage function equation. It is widely used due to its simplicity (Miller et al. 1994; Sausen et al. 1994; Arora & Boer 1999; Huang & Lee 2017). In addition to its simplicity, the cell-to-cell routing method has other advantages. It has the potential to include the nonlinear losses from the surface water to the groundwater because it considers the feedback between the flow and the water storage in each cell (Naden et al. 1999). Besides this, with the cell-to-cell routing method, the water storage in any cell of interest can be queried at each time step. If the ratio of land area covered by surface water to the total land area can be calculated based on the water storage in each cell, it will contribute to land–atmosphere interaction simulation (Olivera et al. 2000). However, the choice of a proper spatial resolution is a dilemma for the cell-to-cell routing method. In distributed hydrological models, high resolution can reflect the heterogeneity of the hydrological characteristics of the basin and improve the accuracy of calibrated parameters. On the other hand, higher-resolution DEM will lead to a significant increase in the amount of required computations, especially when coupled with a meteorological model because the meteorological model has vast computational cost. However, a coarser resolution may lead to a decrease in model performance due to not considering routing within the cell (Arora et al. 2001; Gong et al. 2009). Besides this, there may be a great difference between the network that is generated at a coarse resolution and the true one, which is another disadvantage of the cell-to-cell routing method. Attempts have been made to bring a network that is generated at a coarse resolution closer to the real one. However, a gap still exists (Fekete et al. 2001; Guo et al. 2004).

Some of the source-to-sink routing methods have the same disadvantages as the cell-to-cell method, including the routing models by Ducharne et al. (2003) and Guo et al. (2004), which will provide worse performance as the cell resolution becomes coarser. The routing model by Ducharne et al. (2003) uses the mean travel time of each large routing cell, so the water in a cell will arrive at the outlet in the same day. This can be unreasonable, especially when the cell size is large. The routing model by Guo et al. (2004) uses the mean flow path length for each cell, making the routing method scale-dependent. Therefore, it would benefit a hydrological modeler to find a routing method with the capability to be scale-independent, which will preserve the spatially distributed travel time information from a finer-resolution DEM for coarser resolution cells in the large-scale hydrological model. The disadvantages of some source-to-sink routing methods have been overcome by some model developers. For instance, Wen et al. (2012) proposed a multiscale routing framework that can reduce the impact of the spatial scale by using histograms for the flow path lengths. The aggregated network-response function (NRF) routing method (Gong et al. 2009) can preserve all of the spatially distributed travel time information in high-resolution Hydro1 k DEM data for cells with different resolutions by aggregating the pixel-response function to the cell-response function (CRF). Thus, NRF has the unique advantage of being able to obtain a stable result when moving from a finer to a coarser resolution, making it a preferable routing method for large-scale hydrological models. Another advantage of source-to-sink methods is the computation efficiency (Naden et al. 1999; Olivera et al. 2000; Gong et al. 2009), which enables them to be widely applied (Olivera et al. 2000; Ducharne et al. 2003; Gong et al. 2011; Lu et al. 2015).

NRF, as a source-to-sink routing method, has been applied in basins of China, North America, and Africa (Gong et al. 2009, 2011; Li et al. 2010, 2013, 2014; Ngongondo et al. 2013). The results showed that NRF provides good daily runoff simulation performance. However, in some cases, the calibrated wave velocity is too large to be physically realistic, and NRF will fail to provide good daily runoff simulation performance if the wave velocity becomes smaller (Gong et al. 2011; Li et al. 2014). In this method, the underlying factors that affect the routing process are not comprehensively considered. In this respect, the NRF routing method needs to be improved.

This study aims to improve NRF by considering the influence of the underlying surface in order to obtain a more reasonable wave velocity and travel time and better daily runoff simulation performance. To achieve this aim, the original NRF routing method is coupled with a velocity function proposed by Sircar et al. (1991) in order to take the underlying surface's influence into consideration, in which a new routing parameter b is added to reflect how sensitive the wave velocity is to slope. The remainder of this paper is organized as follows. Following a description of the study area, we present the original and improved NRF routing methods. Then, the results of a daily runoff and peak flow simulation by the model before and after the improvement are compared. The travel time and routing parameters and the factors that influence them are analyzed and compared in different catchments to explore the relationship between catchment characteristics and model parameters. Finally, we present our conclusions.

The study area is the Yangtze River basin (24°27′–35°54′ N, 90°33′–122°19′ E). It has a drainage area of 1,800,000 km². As shown in Figure 1, the Yangtze River is the longest river in China and the third longest river in the world. It has a total length of 6,397 km.

The climate, terrain, land cover and land use, and runoff characteristics vary greatly among the catchments in the Yangtze River basin. Fifteen catchments in the Yangtze River basin were selected in this study in order to evaluate the performance of the improved routing method over a wide range of representative catchments. The characteristics of the terrain in those catchments are shown in Table 1. ‘Elevation difference’ means the height difference between the highest point and the lowest point in a catchment. Table 1 also presents the shape factor and the river length of the catchments. The basin shape factor (Morisawa 1958) is the ratio of the actual basin length to the perimeter of a circle whose area is the same as that of the basin, and measures the general runoff concentration behavior of the basin. Only the first-order stream, second-order stream, and third-order stream (Wang et al. 2019) were included in the stream length in this study. A drainage map of the river system and drainage stations can be found in Figure 1.

The Yangtze River basin has very complex terrain, as it flows through mountains, plateaus, hills, and plains. The overall distribution of elevation at the Yangtze River is high in the west and low in the east. The mean elevation of the source of the Yangtze River, which is located in the hinterland of the Qinghai–Tibet Plateau, is more than 4,000 m above sea level (a.s.l.). The middle reaches lie mostly in low mountains areas and the downstream lies mostly in a hilly plain area.

The mean annual precipitation is 1,110 mm. Due to the vast amount of territory, the complex terrain, and the typical monsoon climate, the spatial distribution and the temporal distribution of annual precipitation in the Yangtze River are very uneven (Xu et al. 2008). The main source of river runoff in the basin is the rainfall, but snowmelt and groundwater are also important sources of runoff in the upper reaches of the Yangtze River basin (Yao et al. 2014). The annual mean air temperature in the Yangtze River basin is generally high in the east, low in the west, high in the south, and low in the north. The source area of the basin experiences the lowest temperatures (Su et al. 2006).

It is of great significance to simulate the runoff in the Yangtze River basin, as it accounts for 18% of the total area of China, and has more than 400 million people living in it. The Yangtze River basin is also the most important agricultural production base in China because of its warm and humid climate, which provides good conditions for crop growth (Tian et al. 2011).

Hydro1 K (USGS 2001) was used in this study to delineate the basin's boundary. It is a global DEM with a resolution of 1 × 1 km and was produced by the United States Geological Survey's Earth Resources Observing System. The projection of HYDRO1 k is Lambert azimuthal equal-area projection, whose transformation equations (Weisstein) enable each grid to have an equal area.

Flow direction and flow accumulation databases based on Hydro1 k are also available, which can be used with the discharge stations to obtain the basin's boundary. Then, the estimated basin boundary can be evaluated and corrected based on the basin boundary that was delineated using the finer DEM at the resolution of 90 × 90 m.

Daily precipitation, temperature, and relative humidity data for the period 1969–2008 were taken from the daily dataset of surface climatological data of China (V 3.0) in National Meteorological Information Center (http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html). The meteorological data were interpolated using the inverse distance weighted method (Shepard 1968) to be gridded data, which were used to drive the hydrological model at resolutions of 3, and 10 arc-minutes in the study. The observed daily discharge during 1964–2008 of all catchments, which was used to calibrate and validate the model, was provided by the Bureau of Hydrology of the Changjiang Water Resources Commission.

The travel time of a flood between major hydrological stations in the Yangtze River basin was taken from Zhang et al. (2006), which was calculated based on a great deal of data and some proven methods. In this study, it was used as a benchmark to compare the travel time calculated by the original routing method with the travel time calculated by the improved routing method.

Gong et al. (2009) developed a routing method, named the aggregated NRF, which was developed with the source-to-sink concept. It builds the relationship between each grid in DEM and the outlet grid.

There are five steps in the improved NRF routing method: (1) extract the flow direction and flow net based on DEM; (2) calculate the travel time t_i between each grid and the outlet of the basin according to the flow path; (3) construct the pixel-response function (PRF) based on t_i; (4) calculate CRF based on PRF with the way of linear averaging method; and (5) calculate the discharge at the outlet of the basin based on the result of the runoff-generation model and CRF. The main advantage of NRF is that the routing accuracy is independent of cell size, because the CRF in NRF is not affected by the cell size (Gong et al. 2009).

Other details about the NRF routing methods can be seen in Gong et al. (2009). In this study, the computation cell size in each catchment was determined by its area to ensure a modest number of cells. The cell size in a larger catchment is larger. Conversely, the cell size is smaller in a smaller catchment. The cell size was set to be at catchments whose area is greater than 50,000 km². It was set to be at catchments whose area is less than 10,000 km² and otherwise.

Compared with Equation (4), Equation (6) is more general and more physically realistic. Because different basins have different underlying conditions, it cannot be known how sensitive the wave velocity is to slope. In this case, b varies in different basins. Equation (6) reflects the relationship between wave velocity and slope more flexibly and can more easily reflect the characteristics of different basins. Thus, Equation (5) is substituted with Equation (7) in the improved NRF routing method.

There are two parameters in the improved NRF routing method (referred to as ‘after’ in the following) and one parameter in the original NRF routing method (referred to as ‘before’ in the following). Table 3 shows the routing parameter value sets for calibration. They were chosen based on computer capability limitations and the physical meaning of each parameter. There are 7 routing parameter value sets before and 42 sets after, composed of six kinds of parameter b and seven kinds of parameter ⁠.

The model before and after the routing method improvement was calibrated and validated in all catchments. The warm-up period was 1 year for all catchments.

The parameters of the routing method were calibrated using the Monte Carlo algorithm (Barraquand & Latombe 1990), and the parameters of the runoff-generation model were calibrated using Covariance Matrix Adaption Evolution (CMAES) (Hansen & Ros 2010). The Monte Carlo algorithm was used to search for the best routing parameter set. First, 300 runoff-generation parameter value sets were produced using LATIN-Hypercube sampling (Mckay et al. 2000) according to the prior range of the parameters. Table 3 shows the routing parameter sets. Second, for each routing parameter set, all of the runoff-generation parameter value sets were used to drive the hydrological model, and then the NSE was calculated based on the simulated runoff and the observed runoff. Third, the best routing parameters were selected based on the NSE. Then, the Covariance Matrix Adaption Evolution algorithm was used to search for the best runoff-generation parameters under the condition that the routing parameters are fixed.

Both routing methods were taken as the routing module of the WASMOD. To all catchments, the model was calibrated by using 7 years' data series and validated by using another 4 years' data series during the period of 1964–2008. Figure 2 shows the performance of the daily runoff simulation before and after the improvement.

Figure 2 shows that the model with the improved NRF routing method performs better than the model with the original NRF routing method for all of the study catchments to varying extents. Taking the Baihe catchment and the Waizhou catchment as examples, the Baihe catchment has an NSE of 0.87 in the calibration period and 0.74 in the validation period after the improvement, while the NSE is 0.43 in the calibration period and 0.31 in the validation period before the improvement. The Waizhou catchment has an NSE of 0.82 in the calibration period and 0.78 in the validation period after the improvement, while the NSE is 0.58 in the calibration period and 0.59 in the validation period before the improvement. The mean improvement in the NSE in both the calibration period and the validation period of the 15 catchments is 0.11, with a maximum of 0.44 in the calibration period and 0.43 in the validation period. From Figure 2, it can also be found that the improvement in NSE varies greatly in the 15 catchments of the Yangtze River basin after the improvement of the routing method.

In order to demonstrate the difference in the simulated and observed discharge hydrographs of the original and improved routing methods, daily hydrographs of a representative year 1968 for the Shimen catchment are plotted in Figure 3. It shows that with the original NRF routing method, the simulated flood peak is often later than the observed flood peak. The improved NRF can better simulate the flood peak. Tables 4 and 5 show the optimal parameters of the hydrological model with the original and the improved routing method.

In order to assess the flood simulation performance of the model with the improved NRF routing method, the averaged relative errors and averaged time lags of all flood events in each catchment were calculated using Equations (11) and (12) and are shown in Figure 4. The flood events whose peak flow was greater than three times the average runoff of catchment were selected in this study. It can be seen that the model with the improved NRF routing method provides better flood simulation performance than the model with the original NRF routing method. After the improvement, the averaged time lags of the peak flow in almost all catchments are less than one day, except for the Shigu and Waizhou catchments (Figure 4, left panel). The averaged relative error of the peak flow in all catchments decreases in nearly all of the catchments (Figure 4, right panel). The amelioration is most evident in the Waizhou catchment, where equals 0.32 with the original routing method and 0.17 with the improved routing method. The relative error of the peak flow decreases to the minimum error range of −0.1 to 0.1 after the improvement in the routing method in all of the catchments except for the Dufengkeng and Fushun catchments.

To evaluate the ability of simulating the occurrence of flood peaks, the time lags in all selected flood events in each catchment were analyzed. Figure 5 plots the frequency distributions histogram of time lags in each catchment with the original and improved NRF routing method. Frequency is the number of flood events that in a certain value is divided by the total number of flood events in this catchment. When the time lag is positive, it means that the simulated flood peak is later than the actual one. On the contrary, when the time lag is negative, it means that the simulated flood peak is earlier than the actual one. When equals zero, the hydrological model simulates the time of flood peak accurately. Figure 5 shows that with the original routing method, the frequency of positive time lag is larger than the frequency of negative time lag ⁠. That indicates that the simulated peak flow is often delayed compared with the observed peak flow with the original NRF routing method. There is a higher frequency of zero time lag of peak flow (⁠ equals zero) in most of the catchments with the improved routing method than that with the original routing method. Particularly, in the Gaochang, Shizhengjie, and Meigang catchments, the frequency of zero time lag of peak flow is greater than 0.8 in the model with the improved routing method and less than 0.4 in the model with the original routing method. Those phenomena indicate that the improved routing methods can better serve the flood simulation.

These results indicate that the model with the improved routing method performs better in both peak flow simulation and long-term runoff simulation than the model with the original routing method. Besides this, NRF's advantage of being able to preserve the travel time information in the PRF is maintained in the improved model.

Figure 6(a) and 6(b) show the distributions of the travel time of all grids in all catchments that were calculated by the original and the improved routing methods. Figure 6(d) plots the mean grid travel time and mean slope of all of the catchments. The simulated travel time was evaluated based on two hydrological stations in each catchment: the outlet station and the specific upstream hydrological station that made the distance between the two hydrological stations as large as possible. Figure 6(c) plots the runoff travel times between the upstream hydrological station and the outlet station that were calculated by the original routing method and the improved routing method. The runoff travel times were compared with the travel time provided by Zhang et al. (2006) as a benchmark.

For most catchments, the travel time between grids that was calculated by the improved routing method was shorter than that calculated by the original routing method (Figure 6(a) and 6(b)). Besides this, the travel time between the specific upstream hydrological station and the outlet station was shorter and closer to the benchmark value for all the catchments (Figure 6(c)). The main reason why the travel time that was calculated by the improved routing method was longer than the benchmark value in some catchments is that the routing parameters and b are spatially constant over the whole catchment. However, routing refers to the process in which rainfall in a catchment evolves along a certain path into discharge at the outlet. The rainfall is first routed by a hillslope and then routed by the river before being transferred into discharge. Hillslope routing and river routing are two forms of routing (Kirkby & Beven 1993). Usually, the runoff travels faster in river routing than in hillslope routing. However, the hillslope travel time and the river travel time are not separated in both the original and improved NRF routing methods. The parameters and b reflect the state of the whole basin; however, the benchmark value reflects the travel time that is needed in the selected river section. Thus, the travel time in the river that is calculated by the improved routing method will be a little longer than the benchmark value. However, the travel time that was calculated in the Shigu catchment is much longer than the benchmark value, which is due to the Shigu catchment's underlying surface. The shape factor of the Shigu catchment is the largest of all of the studied catchments, which means that the shape of the Shigu catchment is narrow and long and the water system's development is immature. The above factors lead to much more hillslope routing and a longer travel time in the Shigu catchment. In addition, the snowmelt and groundwater are the main sources of runoff in the Shigu catchment. The snowmelt runoff moves more slowly down the hillslope than the rainfall runoff, so the optimal routing parameters that were calibrated using the observed runoff lead to a longer travel time in the Shigu catchment compared with the benchmark value, which only considers routing in the river.

The travel time that was calculated using the original routing method produces a late peak flow and a smaller amount of flow (Figure 4). An unreasonable travel time is one of the main causes of a hydrological model's poor performance in long-term runoff simulation. The improved routing method provides a more realistic travel time between the grids compared with the benchmark value.

When analyzed in conjunction with Figures 2 and 6, it can be found that the NSE and the travel time in all catchments have a close relationship. From Figure 2, it can be seen that the NSE that was calculated by the original NRF is lower in some catchments, i.e., the Waizhou, Shimen, Beibei, and Baihe catchments. As shown in Figure 6(c), the runoff travel times in these catchments, as calculated using the original routing method, were unrealistic and at least twice as long as the runoff travel times that were calculated using the improved routing model. On the other hand, for the Dufengkeng, Fushun, and Wulong catchments, we found no obvious change in the travel times that were calculated using the original and improved routing methods. The original and improved routing methods also provided a similar NSE in these catchments.

Figure 7 shows the theoretical curves of wave velocity versus slope. Figure 8(a) and 8(c) plot the calibrated curves that were calculated by the original routing method and the improved routing method (Table 3), respectively. From the analysis in the section Analysis of the routing time calculated by the improved NRF, it is known that the improved routing method is able to achieve a more reasonable travel time. The reason for this is that the improved routing method has a more flexible wave velocity–slope curve (Figure 7) that can better reflect the real world.

It can be seen from Figure 7(b) that the wave velocity of a grid is determined by the parameters b and together. If the parameter b remains constant, the wave velocity increases as increases. If the parameter remains constant, the wave velocity increases as b decreases at grids with a slope of less than 45°; however, the wave velocity increases as b increases at grids with a slope of greater than 45°. The effect that slope has on the wave velocity varies across catchments due to the differences in the hydrology of the underlying surface, which can be adjusted by the parameter b. The bigger b is, the more sensitive the wave velocity is to slope. When b is equal to zero, the slope has no effect on wave velocity, which will be a constant. A larger wave velocity results in a shorter runoff travel time in a catchment.

Each routing parameter value set refers to a certain curve between slope and calibrated wave velocity. The original NRF only produces three kinds of curves for the 15 catchments (Figure 8(a)). The curves of 13 catchments coincide because they have the same optimal value of ⁠. However, the improved NRF produces a different calibrated curve for each catchment (Figure 8(c)). Tables 6 and 7 present the optimal values of that were calibrated by the original routing method with a strict range of 4–10 m/s (cali1), as in Table 3, and with a wider prior range of 4–31 m/s (cali2), respectively. Due to computer capability limitations, the increment of parameter above 10 m/s is 3 m/s. Figure 8(e) plots the calibrated curves between slope and calibrated wave velocity in cali2. In combination with the optimal value of shown in Tables 6 and 7, it can be seen that the optimal value of in cali1 is the same as 10 m/s, reaching the upper limit of the prior range for nearly all of the catchments. When the prior range for calibration is fixed to be 4–10 m/s, the larger the value of ⁠, the shorter the travel time, which is closer to the reality in most catchments. This results in a larger NSE value. Thus, the calibrated value of reaches the upper limit of the prior range for calibration in 13 catchments, resulting in only three kinds of curves between slope and calibrated wave velocity in 15 catchments. Moreover, in many basins, the optimal value of in cali2 is larger than 10 m/s and has a maximum of 31 m/s; thus, the grid wave velocity is too large to be realistic (Figure 8(e)). Figure 8(b), 8(d) and 8(f) show the distributions of the wave velocity of all grids in all catchments before the improvement with the strict prior range of parameter for calibration; after the improvement with the strict prior range of parameter for calibration; before the improvement with the wider prior range of parameter for calibration. It also shows that the grid wave velocity is too large to be realistic before the improvement with the wider prior range of parameter for calibration. The model can provide good runoff simulation performance because the excessively large value of compensates for the problem of the routing time being too long due to the unreasonable distribution of the grids' wave velocity. This indicates that the original routing method cannot adapt to the variation in the underlying surface of the different basins. This is consistent with the results from Gong et al. (2011), who found that the optimal value of was physically unrealistic as it was determined to be greater than 26 m/s with the original Hydro1 k-driven NRF routing method.

The original NRF routing method uses a constant (0.5) to reflect the influence that the slope has on the wave velocity; however, it fails to adapt to the variation in the underlying surface of the different basins. The improved routing method makes it possible to fit the actual physical relationship between the slope and wave velocity in a more flexible way. The reasonable distribution of wave velocity leads directly to a more accurate simulation of the travel time, which in turn leads to a better simulation of the runoff.

The improved NRF routing method belongs to hydrological methods rather than hydraulics methods. There are some other hydrological routing methods calculating the travel time of each grid and then routing the runoff to discharge. But the parameter to calibrate and influence factors that have been taken into consideration are different in different methods.

For the routing method in Ducharne et al. (2003), only the slope between two grids will influence the velocity. For the routing method in Guo et al. (2004), not only the slope, but also the hydraulic radius, the roughness of the cell and the scale factor about the cell size will influence the wave velocity. For the routing method in Bunster et al. (2019), it proposed a travel time formulation that accounts for the dynamics of the upstream contributions, compared to the previous methods. In that method, the dynamics of the upstream contributions and slope will influence the wave velocity. Du et al. (2009) proposed a time variant routing method, the spatially distributed direct hydrograph travel time method (SDDH), in which time variant runoff, grid length, the roughness of the cell, and slope are used to calculate the wave velocity.

Some data such as hydraulic radius may not be available when applying the large-scale hydrological model. Thus, the routing method requiring too much input information may not preferable in a large-scale hydrological model. In the improved NRF routing method, all the necessary information can be taken from the DEM. Compared with the routing methods that only take the slope into consideration, the improved NRF can be adapted to more basins. Overall, the improved NRF routing method is simple and effective for a large-scale hydrological model.

The wave velocity of a grid is affected not only by the grid's slope but also by many other factors, including the basin shape factor, rainfall intensity, vegetation cover, branching ratio, basin area, climate, and form of drainage. The factors that influence the travel time are so numerous that they may compensate for each other to some extent. It is difficult to calculate the wave velocity of each grid precisely. So, these influencing factors were generalized to be reflected by the parameter ⁠, the parameter b, and slope in this study. The parameter b varies across basins due to differences in basin characteristics, and it represents the underlying surface.

The underlying surface influences the wave velocity and travel time in many ways. There is no strict and deterministic functional relationship between basin characteristic and mean wave velocity in a large-scale hydrological model. Multivariate regression analysis (Alexopoulos 2010; Zhang et al. 2014) is a mathematical analysis method that uses the least squares method to model the relationship between a dependent variable and multiple independent variables. In this study, we used it to examine the relationship between the mean wave velocity of grids (m/s) calculated with the improved routing method and several influencing factors. The result equation will not be used to calculate the wave velocity directly. But this method can identify basin characteristics that have a greater impact on the mean wave velocity. Some variables that are related to wave velocity are chosen (Table 1): mean discharge at the outlet station (m³/s); area A (km²); shape factor F; mean slope S (°); river length (m); and range of elevation (m). The area, shape factor, and river length are basin geometric factors which affect the path and sequence of runoff and then influence the mean wave velocity. The slope and range of elevation influence the wave velocity, because potential energy will be converted into kinetic energy. In the process of runoff routes, runoff will continually merge and diverge, and the speed between water will affect each other. Thus, the mean discharge at the outlet station is also taken as an independent variable.

There are other factors that are not included in Equation (13), and for R² it is not equal to 1. However, we have yet to find more accurate and detailed rules. The mean wave velocity that was predicted by Equation (13) was compared with the mean wave velocity that was calculated with the improved routing method as shown in Figure 9, which shows an overall good result.

In the future study, more data on the factors that influence wave velocity, such as soil, vegetation cover, and the branching ratio, should be inputted into the method in order to examine the relationship between routing parameters and basin characteristics. If appropriate, the improved routing model has the potential to be applied to data-sparse areas by a parameter transformation based on the distribution in a subzone of climate, soil, and vegetation, or even to calculate the travel time without parameters, which would enable this model to be used to simulate the hydrological effect of climate and land-use change.

In this paper, a flow routing algorithm was developed by combining an aggregation NRF with a velocity function. The improved routing NRF method uses a parameter b that can reflect how sensitive the wave velocity is to a slope. The routing methods before and after the improvement were applied to 15 catchments in the Yangtze River using a daily WASMOD-M model based on Hydro1 k. The main conclusions can be summarized as follows:

• (1)

  The original NRF routing method was found to provide unsatisfactory runoff simulation performance in most of the studied catchments, with an unreasonable calibrated travel time and wave velocity of the grids, for it cannot adapt to the different underlying conditions in the different catchments.

• (2)

  After the improvement in the NRF routing method, the relationship between wave velocity and slope was found to be more flexible. It can simulate the complexity of the wave velocity distribution more effectively in all types of basins, resulting in a reasonable runoff travel time.

• (3)

  The runoff travel time that was calculated by the model with the improved NRF routing method was found to be shorter and more reasonable, which leads to better runoff simulation performance.

• (4)

  The model with the improved NRF routing method was found to yield better results with respect to long-term daily runoff time series and peak flow amounts than the model with the original NRF routing method. The improved model was found to perform better in catchments with different characteristics. The improved model also retains the best advantage of NRF; i.e., it records all of the travel time information from the PRF.

This study was partly supported by the National Natural Science Fund of China (grant no. 51539009) and the National Key Research and Development Program (grant no. 2017YFA0603702). The authors are grateful to the National Meteorological Information Center for providing the meteorology datasets and the Bureau of Hydrology of the Changjiang Water Resources Commission for providing the datasets for the Yangtze River basin.