Understanding the impacts of catchment characteristics on the shape of the storage capacity curve and its influence on flood flows

Flooding is one of the leading natural hazards (Hirabayashi et al. 2013; Blaikie et al. 2014; Winsemius et al. 2016), causing significant casualties and property damage globally every year (DHI 2012; Blaikie et al. 2014). Flood peak flow is one of the most important variables in flood modelling and forecasting. The impacts of rainfall amount, intensity and the antecedent soil moisture on flood generation have been well studied (Yu et al. 1999, 2000; Blöschl et al. 2015; Berghuijs et al. 2016). However, the impact of catchment features on flood peak flow is neither systematically studied nor well explored, particularly in large-sample datasets. This lack of understanding hinders improved understanding of catchment storm-flood processes, and furthermore, prediction in ungauged basins (Sivapalan et al. 2003).

Additionally, understanding the mechanisms of storm-flood processes and their link to catchment features is essential for decision-making on flood prevention and mitigation (Palmer et al. 2015), as well as land-use planning (Andréassian 2004; Sriwongsitanon & Taesombat 2011). Understanding relationships between flood flows and catchment features is also beneficial for prediction of future hydrological changes (Xu & Singh 2004; Xu et al. 2005; Milly et al. 2008; Montanari et al. 2013), and the impact of land use changes (Stonestrom et al. 2009). Therefore, there is an urgent need to improve our understanding of linkages between floods and catchment characteristics.

Partitioning rainfall into runoff and infiltration is the central question of all hydrological models (Singh & Woolhiser 2002; Beven 2012). It also plays an important role in land surface modelling (Sellers et al. 1996; Dai et al. 2003) and biogeochemical modelling studies (Austin et al. 2004) due to the fact that energy and nutrients are transported with the movement of water (Mo & Zhang 2013). The dominant perceptual models (Beven 2012) of rainfall-runoff process mechanisms are infiltration excess overland flow (Hortonian overland flow (HOF)) (Horton 1933), and saturation excess flow (SEF) (Hewlett 1961; Dunne & Black 1970). HOF, in conjunction with the unit hydrograph, dominated engineering hydrology for several decades in the early 20th century (Beven 2004a). However, findings from more and more field experiments (Sklash & Farvolden 1979; Beven 2004b; Burt & McDonnell 2015) questioned the role of HOF, even in arid catchments (Liu et al. 2012; Li et al. 2014). In contrast, the SEF theory, also known as the variable contributing area (VCA) model, assumes that only certain saturated areas of the catchment, contribute to runoff yield (Hewlett 1961; Dunne & Black 1970; Hewlett & Troendle 1975). The saturated area varies with different soil moisture conditions, and is strongly impacted by the spatial heterogeneity of the storage capacity at a catchment scale. This theory has been implemented in various conceptual hydrological models, i.e. Xinanjiang (Zhao et al. 1980), HBV (Bergström & Forsman 1973; Bergström & Lindström 2015), TOPMODEL (Beven & Kirkby 1979), GR4 J (Perrin et al. 2001). All these models have demonstrated good performance in streamflow simulation in many catchments (Singh & Woolhiser 2002). The Xinanjiang, HBV and GR4 J models apply a conceptual storage-contributing area relationship to estimate the event runoff yield. The TOPMODEL assumes that topography determines the relation of storage-contributing area, and Dynamic Topmodel improves this idea by including landscape information (Beven & Freer 2001).

The amount of rainfall and the antecedent soil moisture are two essential variables to determining the amount of runoff in SEF-type conceptual models. However, it should be noted that different catchments will likely generate different amounts of runoff even with the same amount of rainfall and precedent soil moisture, due to the different spatial distribution of unsaturated storage, which is controlled by land surface features. Many efforts have been made to study the impacts of extreme precipitation events and antecedent soil moisture on flood peak flows (Blöschl et al. 2015; Berghuijs et al. 2016), however, very few studies have systematically explored the influence of various catchment land surface features (except for land cover, see Andréassian 2004; Blöschl et al. 2007) on flood flows. In conceptual hydrological models, such as Xinanjiang, HBV, and GR4 J, catchment features are expressed by two essential parameters: the root zone storage capacity (S_uMax) and the shape of the storage capacity curve (often termed β). The former was found to be controlled by climate (Gao et al. 2014a), and determines the long-term runoff coefficient and evaporation. The latter reflects the spatial distribution of the storage capacity (Zhao 1992; Sivapalan & Woods 1995), and determines the peak flow given the amount of rainfall and antecedent soil moisture. As a result, the shape of the storage capacity curve is useful in estimating the amount of runoff yield and more importantly, the flood peak flow. Further details about the shape of the root zone storage capacity curve β are given below in Methodology: Storage capacity curve and flood.

Parameter regionalization is an approach to understanding the spatial distribution of parameters (Merz & Blöschl 2004; Kling & Gupta 2009) and the correlations between optimized parameter values and catchment physical features (Samaniego et al. 2010). Undertaking parameter regionalization is not only beneficial for understanding the linkages between parameters and catchment features, but also helpful in predicting hydrological processes in ungauged basins (Hrachowitz et al. 2013). However important questions remain. For example, the β parameter was hypothesized to be determined by topography and geology (Zhao et al. 1980), but Sivapalan & Woods (1995) found that soil depth distribution also plays a determinative role in β estimation. Therefore, further studies are required to explore the spatial distribution pattern of β, and its linkage to the physical features of catchments.

The representativeness and comprehensiveness of sample data used for hypothesis testing is an important point to consider in hydrological studies considering the uniqueness of different catchments. Since findings from a limited number of catchments may be only applicable at local and regional scales, wider generalization to other catchments may be problematic. Large sample data sets are therefore essential for hypothesis testing and the drawing of robust conclusions for diverse environments (Gupta et al. 2014). Additionally, the knowledge gained from regional studies needs to be confirmed and validated by sufficiently large sample datasets due to the complexity of catchment characteristics. Moreover, over the past decades, great efforts have been made to collate and share hydrological data sets, leading to increasing availability of relevant data sets for a large number of catchments at continental and global scales, thereby facilitating large-sample empirical studies. These hydrological datasets include, for example, MOPEX (Model Parameter Estimation Experiment) for the contiguous United States (CONUS) (www.nws.noaa.gov/oh/mopex/), the datasets of the UK (https://nrfa.ceh.ac.uk/), France (www.cnrm.meteo.fr/isbadoc/projects/rhoneagg/), and Australia (Peel et al. 2000). In addition, global datasets include GEWEX (www.gewex.org/) and GLACE (gmao.gsfc.nasa.gov/research/GLACE-2/). The MOPEX data (Duan et al. 2006) contain good quality daily data of key variables for 438 catchments in the USA, and will be used in this study.

In this study, the Xinanjiang storage capacity curve (Zhao et al. 1980; Liang et al. 1994), widely used in hydrological and land surface modelling, was coupled with the FLEX model framework (Fenicia et al. 2011; Gao et al. 2014b; Gharari et al. 2014). This approach was applied to simulate the hydrological processes of the large-sample MOPEX catchments, and to calibrate the β parameter which determines the flood peak flow given certain soil moisture conditions. Subsequently, correlations between β and typical catchment characteristics were explored. These included topography, soil hydraulic characteristics, climate, vegetation cover, urban area, geology (soil depth to bedrock), basin area and area of impermeable surfaces. We aimed at answering the following three questions in this study: (1) Is there any regionalized pattern of the parameter of the shape of storage capacity curve (β)? (2) Which catchment characteristics affect the β value and the catchment peak flows, and to what extent? (3) What is the physical meaning behind the linkages between the β value and catchment characteristics?

The MOPEX data were collated for hydrological model parameter estimation experiments (Duan et al. 2006), and contain relevant data for 438 catchments in the USA, including, daily precipitation, daily maximum and minimum air temperature, and daily streamflow from 1948 to 2003. Thirty-four catchments with short streamflow observations (less than 3 years of continuous streamflow data) were excluded, while the remaining 404 catchments were used in this study. Among these catchments, rainfall-runoff flooding is the dominant flood/runoff generation mechanism, but snowmelt-runoff flooding also occurs (Berghuijs et al. 2016). Temperature data were used to estimate snowmelt and its potential evaporation. Precipitation and calculated potential evaporation are the two major input forcings for the rainfall-runoff model. Daily streamflow was used to calibrate the free parameters in the hydrological model. Further details of the modelling process are given later in section Methodology: Storage capacity curve and flood.

Various types of landscape information have been collated and further details are provided below and in Table 1.

The digital elevation model (DEM) of the CONUS at 90 m resolution (Figure 1) was downloaded from the Earth Explorer of United States Geological Survey (USGS, http://earthexplorer.usgs.gov/). Slope (in degrees) (Figure 1) was derived from the DEM using ArcMap toolbox. For each catchment, the average elevation and slope values were calculated and used to represent the catchment elevation and slope for further analysis.

Depth to bedrock was regarded as the integrated geologic indicator for hydrological processes, since bedrock as the lower boundary of subsurface stormflow substantially impacts the saturated storage capacity and then runoff generation processes and flooding (Tromp-van Meerveld & McDonnell 2006). The depth to bedrock data (Schwarz & Alexander 1995) (Figure 1) were accessed from STATSGO (State Soil Geographic, http://www.soilinfo.psu.edu/index.cgi?soil_data&conus&data_cov&dtb). The average depth to bedrock for each catchment was calculated by ArcGIS toolbox and used for further analysis.

Soil texture is synthetically represented by the K factor, since the K factor is a lumped soil erodibility factor which represents the soil profile reaction to soil detachment caused by raindrops and surface flow (Renard et al. 2011). Generally, the soils (high in clay and sand) have low K values, and soils with high silt content have larger K values. The K factor for each catchment was calculated from soil survey information available from USGS (Wolock 1997).

Land use data were obtained from National Land Cover Database (NLCD, http://www.mrlc.gov/nlcd.php). Due to the complexity of land use information, we selected the proportion of forest cover and the urban area proportion (Figure 1) as two integrating indicators to distinguish the driving force either by nature or human impacts. Forest plays an essential role in hydrological processes (Jasechko et al. 2013), especially in the runoff generation processes (Brooks et al. 2010). The proportion of urban area is a proxy land use change related to human activities (Grimm et al. 2008), with consequent impacts on catchment hydrology (Lee & Heaney 2003).

Area data were obtained from the MOPEX dataset. Stream density (km/km²) is the total length of all the streams and rivers in a drainage basin divided by the total area of the drainage basin. Stream density was obtained from Horizon Systems Corporation (Table 1).

Daily potential evaporation was estimated using the Hargraves equation with maximum and minimum daily temperature from MOPEX as input (Hargreaves & Samani 1985). The aridity index (long term potential evaporation divided by long term precipitation) was selected to represent the climate for a given catchment (Budyko 1974; Fu 1981), and long-term mean annual precipitation and estimated potential evaporation from MOPEX datasets was used to compute the aridity index for individual catchments.

We used the FLEX modelling framework to simulate the hydrological processes of the MOPEX catchments in this study. The FLEX framework permits flexibility in model structures for different modelling purposes (Fenicia et al. 2008; Gao et al. 2014b; Gharari et al. 2014). For example, the snowmelt and interception processes have been included in this study for relevant MOPEX catchments. The FLEX framework has been used successfully to simulate streamflow in diverse catchments, with satisfactory performance (Fenicia et al. 2008; Gao et al. 2014b; Gharari et al. 2014; Li et al. 2014; Zhao et al. 2016).

The FLEX model contains five conceptualized reservoirs, including the snow reservoir (S_w), interception reservoir (S_i), unsaturated reservoir (S_u), and the two recession reservoirs (S_f and S_s). Snowmelt was calculated by a widely used parsimonious temperature-index method (Hock 2003; Gao et al. 2012, 2017). Before infiltration occurs, precipitation/snowmelt is first intercepted by the vegetation canopy. In this study, the interception was estimated by a threshold parameter (I_max), below which all precipitation was intercepted and evaporated (De Groen & Savenije 2006). The effective rainfall/snowmelt is split by the Xinanjiang storage capacity curve (further details are given in the next section), into generated runoff (R_u) and infiltration. The generated runoff (R_u) is further split into two fluxes, the flux to the fast response reservoir (R_f) and the flux to the slow response reservoir (R_s), by a splitter (D). The delayed time from storm to the flood peak is estimated by a convolution delay function, with a time delay of T_lagF. Subsequently, the fluxes into two different response reservoirs (S_f and S_s) is controlled by two linear equations between discharge and storage, representing the subsurface storm flow and the groundwater flow, respectively. The two discharges (Q_f and Q_s) were merged to generate the simulated streamflow (Q_m). The model parameters are shown in Table 2, while the equations are given in Table 3. A more detailed description of the model structure can be found in Gao et al. (2014b). The FLEX model has a total of 10 free parameters to be calibrated that are shown in red in Figure 2.

In Figure 2, precipitation (P) is split into snowfall (P_s) and rainfall (P_f) based on threshold temperature (T_t). Snowmelt is calculated based on the temperature-index method, using degree-day factor (F_DD). Rainfall is intercepted by an interception reservoir (S_i) with a capacity parameter I_max (mm). P_e (mm d⁻¹) is the effective rainfall to soil surface. The unsaturated reservoir (S_u) is used to split the P_e (mm d⁻¹) into generated runoff R_u (mm d⁻¹) and infiltration by two free parameters S_uMax (mm) and β (−). C_e (−) controls the relationship between soil moisture and actual evaporation. R_u recharges two linear reservoirs (S_f and S_s) by a parameter (D) and delay time parameter (T_lagF). The response processes of subsurface storm flow Q_f (mm d⁻¹) and groundwater runoff Q_s (mm d⁻¹) are modelled with two recession parameters, K_f (d) and K_s (d).

Different shapes of the curve will generate different amounts of runoff under the same rainfall and antecedent soil moisture conditions (Zhao et al. 1980). The prior range of β is recommended to be between 0.1 and 0.4 by the original model developer (Zhao et al. 1980). However, Liang et al. (1994) found that the optimized β can be 0.001 in certain catchments, while Wood et al. (1992) found that the optimized β could be larger than 1 in some cases. Sivapalan & Woods (1995) recommended setting the prior range as 0.01 to 5.

In order to illustrate the impact of varying β on the storage capacity curve and runoff generation, we conducted a scenario analysis (Figure 3) using the following assumptions: S_uMax = 100 mm, antecedent soil moisture (S_u) = 50 mm, instantaneous precipitation is set as 20 mm, and all initial conditions and other relevant parameters are kept the same. If β = 0, the storage capacity curve becomes a simple rectangular bucket, and no runoff will be generated until the rectangular bucket is full. Therefore, there is no runoff yield in this case because the rectangular bucket is not saturated, even after 20 mm rainfall. If β is larger than 0 and less than 1, the storage capacity curve will be a convex curve, which indicates that most runoff will be infiltrated and runoff generation is limited when soil moisture is low. In this case, flat peak flow is generated, as shown in the right panel of Figure 3. If β is larger than 1, the storage capacity curve will become a concave curve, which means the fraction of runoff to rainfall is much larger than the case of β with smaller values. In this condition, the shape of flood hydrograph will be higher, thinner and sharper (Figure 3). In summary, Figure 3 demonstrates that the β parameter significantly affects the runoff response to rainfall, including the flood volume and the flood peak flow. Therefore, we can use the β value, as a lumped parameter, to quantify the storm-flood response. More specifically, the catchments with higher β values are more likely to generate floods with larger flood volumes and sharper and thinner flood hydrographs, and vice versa.

A multi-objective parameter optimized algorithm (MOSCEM-UA) (Vrugt et al. 2003) was applied for calibration purposes. The parameter sets on the Pareto-frontier were assumed to be behavioral parameter sets and can equally represent model performance under certain criteria. Therefore, we used all the parameter sets on Pareto frontier in the calibration. The number of complexes in MOSCEM-UA was set as 10 reflecting the number of parameters, and the number of initial samples was set to 200 and a total number of 50,000 model iterations for all the catchment runs. Since the model structure has been validated in the USA (Berghuijs et al. 2014; Gao et al. 2014b) and internationally (Fenicia et al. 2008; Gharari et al. 2014; Gao et al. 2016; Zhao et al. 2016), only these parameters were calibrated and no further validation was implemented.

Linear correlation analysis was conducted to understand the association between the calibrated β parameter and catchment features. Associations were quantified by the correlation coefficient (r), and significance was evaluated by the significance test (indicated as p), to diagnose whether r is different from zeros.

Figure 4 shows a boxplot of two objective functions (I_KGE and I_KGL) to summarize the model performance. The boxplot shows that the median value (50% percentile) of the I_KGE of the 404 MOPEX catchments is 0.72. The maximum and minimum values of the I_KGE whisker are 0.92 and 0.49, respectively. The 75% and 25% percentiles of the I_KGE boxplot are 0.78 and 0.63, respectively. For the boxplot of the I_KGL, the median value is 0.66, with 0.89 and 0.43 for the maximum and minimum of the whisker, respectively, and 0.71 and 0.59 for the 75% and 25% percentiles, respectively. The results show that the FLEX model can perform well to simulate the hydrological processes of most MOPEX catchments (Duan et al. 2006; Ye et al. 2014).

Figure 5 shows the spatial distribution patterns of the two objective functions. In relation to the spatial distribution of I_KGE, catchments in the east and west coast of the CONUS outperform simulations in the interior regions, most probably due to the impacts of the arid interior climate. This is in line with previous hydrological studies that more arid zones are characterized by poorer model performances (Wheater et al. 2007; Berghuijs et al. 2014). The model performance in relation to the I_KGL function is less clear. The I_KGL function is more sensitive to baseflow simulation, where human activities, including dams, irrigation and water allocation (Wang & Hejazi 2011), have more influence than natural processes.

The β value of each catchment was optimized and calibrated for each of the 404 MOPEX catchments and for each, the averaged β value on all Pareto-frontier parameter sets was regarded as the β of that catchment. In Figure 6, the β value of each catchment is shown together with the slope map. The β values exhibit a strong spatial distribution pattern in the 404 MOPEX catchments, and generally, the β parameter has a strong positive correlation with slope. Figure 6 clearly shows that the catchments with relatively small β values (0.01–0.5) have flat slopes, and these catchments are mostly located in Great Plain the interior region of the CONUS. This means the flatter catchments tend to generate relatively dampened flood peaks, as indicated by smaller β values. In contrast, the catchments located in the Rocky Mountains and the Appalachian Mountains have lager β values (>0.5). This indicates that catchments with steeper slopes have relatively high, sharp and thin flood peaks, which is in line with previous studies (Ruiz-Villanueva et al. 2010; Bodoque et al. 2015). The regionalized pattern of the β parameter can help us to understand and predict hydrological processes and flood peak flow in ungauged basins.

Figure 7 shows the correlations between the calibrated β values and the catchment landscape features considered for all 404 catchments. Specifically, the β value increases with increases in the catchment averaged slope, elevation, impermeable area, and percentage of forest. These positive correlations can be interpreted from the perspective of catchment rainfall-runoff processes. Particularly, steeper slopes are prone to generate sharper (higher and thinner) flood peak flow represented as larger beta values, since steeper catchments have larger gravitational potential energy causing the quicker release of generated runoff and thus resulting in more concentrated flooding in a short time. Similar to the slope, elevation also has a positive correlation with β values, which is probably due to the fact that most areas with steep slopes are located in regions of higher elevation. Additionally, larger urban area corresponds with thinner and higher peak flow and larger β values, since the urbanized areas are normally characterized by impermeable surfaces (Grimm et al. 2008), reduced soil storage capacity and altered rainfall infiltration and surface runoff generation (Lee & Heaney 2003).

It is worth noting the positive correlation between forest proportion and beta values, which means a larger proportion of forest cover is associated with higher flood peak flow. The impact of forest cover on flood peak flow is still controversial in the hydrology community (Robinson et al. 2003; Andréassian 2004), because some studies have found that forest can reduce peak flow for small and medium floods (Sriwongsitanon & Taesombat 2011). Therefore, a negative correlation between β and forest cover was expected. We suggest three factors may be contributing to this positive correlation: (1) the spatial concurrence of mountainous region and forest cover, linked to a greater likelihood of anthropogenic activities such as deforestation and reclamation in flat areas. This may refer that the correlation of β and forest cover is due to the coherent of forest and mountainous region rather than forest cover itself; (2) the positive correlation between forest cover and precipitation, as Duan & Bastiaanssen (2013) found by satellite data; and (3) the complexity of simultaneously considering time and space. From a temporal perspective, deforestation likely increases flood peak flow, and reforestation probably decreases peak flow. In spatial terms, there are still few studies that address the association between flood peak flow and area proportion of forest cover with large sample dataset (Weingartner et al. 2003), and no clear relation has been found between forest proportion and flood discharge. Further research is therefore needed to explore the relation between forest cover and flood peak flow.

In contrast, the calibrated β parameter has a negative correlation with the soil erodibility factor (K factor), aridity index, the depth to bedrock and catchment area. Specifically, soils with a low K factor are very likely those with a higher proportion of clay and less detachment, which are prone to reduce infiltration and generate higher and sharper flood peak flows. The aridity index, representing climate, determines the correlation between precipitation and potential evaporation, water supply and evaporation demand. The negative correlation between aridity index and β indicates that in arid catchments soil is generally in an unsaturated condition, which allows more rainfall to infiltrate and thus reduces the peak flow and β values. The catchments with a shallower soil layer (the depth to bedrock) have less storage capacity, which makes them prone to flood generation. Catchment area also correlated negatively with β values, indicating that catchments with smaller area are apt to generate sharper flood peaks, probably due to the smaller capacity of channels in small catchments to regulate flood peaks. Catchments with larger areas generate flatter hydrographs than small catchments, because of the larger regulating capacity of the river channel network and the averaging effect of floods from different tributaries. This could also be related to there being a greater chance that rainfall is nonuniformly high across smaller catchments, and in larger catchments there is a greater chance for other processes to buffer the flood flows. The stream density had no significant correlation with β and the range of responses suggests that stream density can either increase flood peaks or reduce flood peaks depending on specific local geomorphology and channel bathymetry.

In summary, among all considered catchment features, slope has the largest impact on β and flood peaks, with a correlation coefficient (r) of 0.58 (p < 0.001) (Figure 7). The percentage of forest cover has the second largest impact on the β parameter, with r of 0.56 (p < 0.001). The K factor, elevation, impermeable area, aridity index, the depth to bedrock and catchment area also have significant impacts on β (p < 0.001), but with lower correlation coefficients (from 0.2 to 0.48) (Figure 7). Vegetation, soil characteristics, urban area, climate, geology and catchment area also impact flood peaks, but with weaker correlations. A clear spatial distribution pattern of the β parameter was also found. These conclusions are generally in line with previous parameter regionalization studies in Austria (Merz & Blöschl 2004) and the southeast of England (Wagener & Wheater 2006).

This research highlights the impact of topography on flood peaks. Globally, mountain flooding is an increasing risk to human populations, as settlement in these areas increases (Ruiz-Villanueva et al. 2010). Between 2000 and 2010, floods in mountain regions resulted in 1079 casualties every year in China, accounting for 65–92% of total flooding casualties. Mountain basins are often poorly gauged, and given the potential impacts of flash flooding in these areas, this study will be beneficial to predicting floods in mountainous regions. The findings of this study also have practical implications. For example, to reduce the peak flow for flood mitigation, we should take measures to reduce relative slope (e.g. building terraced fields), reduce the proportion of impermeable areas (using green roofs, permeable roads, green spaces and other ecosystems to absorb rainwater in urban areas), or increase the soil infiltration capacity (e.g. reforestation), For example, our results suggest that implementing the ‘sponge city’ in China could reduce flood impacts (Xia et al. 2017).

However, it should be noted that there are several limitations in our current study that should be addressed in further studies. Only nine catchment attributes were considered in this study, and further attributes should be explored in future research. In addition, for each of the considered features we used a single value (catchment-average) to represent each catchment in this study, and neither the variation nor the spatial distribution of corresponding features (e.g. distribution pattern of topography, vegetation and soil texture within catchments) were taken into account. How these spatial distribution patterns could be represented should be investigated further.

In this study, we used the MOPEX dataset encompassing 404 catchments in the USA to explore the correlations between the shape parameter (β) of the storage capacity curve and various typical catchment features. We found that the β parameter, has clear regional patterns, influenced by topography, vegetation, soil, urban area, geology, and catchment area after controlling for the flood flows given rainfall and antecedent soil moisture. Generally, the catchments with steeper slopes, less erodible soil texture, larger impermeable areas, shallower soil depth, and smaller catchment areas will have larger β values, and thus are prone to sharper flood peaks. Among these features, topography has the largest impact on peak flow. This study is beneficial to understanding parameter regionalization related to flood peak flow and the impact of catchment features thereon.

This research is supported by the National Natural Science Foundation of China (Grant No. 91547202, 51479216), the State Key Laboratory of Cryospheric Sciences, Cold and Arid Regions Environment and Engineering Research Institute, Chinese Academy Sciences (SKLCS-OP-2016-04), and the Key Laboratory of Mountain Hazards and Earth Surface Process, Chinese Academy of Sciences (KLMHESP-17-02), and the National Key Research and Development Program of China (No. 2016YFC0402600). We thank Wouter Berghuijs for sharing the catchment characteristic information of the MOPEX catchments, and Dingbao Wang for generously sharing the shapefile of the MOPEX map. Thanks also go to Eric Moody for his review and grammatical corrections for the final manuscript. We would like to thank two anonymous reviewers and editor Nevil Wyndham Quinn for their valuable and constructive comments and suggestions for improving this manuscript.