Accuracy and efficiency have resulted in contradictions in semi-distributed hydrological models, such as Soil and Water Assessment Tool (SWAT), Topography-based Hydrological Model (TOPMODEL) and Hydrological Simulation Program-Fortran (HSPF) model (Bourdin et al., 2012). Theoretically, it is assumed that improved model accuracy could be obtained either by more precision in input data spatial resolution or finer discretisation of the watershed (Fatichi et al., 2016; Devak and Dhanya, 2017; Kan et al., 2019). Meanwhile, considerable computational time and input data preparation effort were required when implementing the model, and these did not necessarily provide a better simulation result beyond a certain level (Azizian and Shokoohi, 2015; Guan et al., 2015; Reddy and Reddy, 2015; Boithias et al., 2017; Munoth and Goyal, 2019; Park et al., 2019; Nazari-Sharabian et al., 2020; Roostaee and Deng, 2020). Consequently, the optimum watershed discretisation level as well as appropriate spatial input parameters became an on-going theoretical debate.

A SWAT model is spatially distributed and therein the watershed is partitioned into a number of sub-watersheds connected by stream networks using digital elevation model (DEM) data before further subdivision into multiple hydrologic response units (HRUs) by overlaying spatial datasets including slope, land-use, and soil maps to represent the spatial heterogeneity of the watershed (Savvidou et al., 2014; Luo et al., 2019). Model computational time is nearly proportional to the number of HRUs, since HRU is the basic calculation unit (Wang et al., 2016). In some cases, HRU numbers exceeded the computational limits on the model partly due to the large watershed scale or high level of discretisation (Chiang and Yuan, 2015). Decreasing numbers of HRUs, when a non-zero HRU threshold is applied, would contribute to reduce SWAT model simulation time and increase computational efficiency. Nonetheless, elimination of smaller HRUs disproportionately results in a poorer representation of watershed heterogeneity. Accordingly, the HRU threshold has the potential to affect the accuracy of model prediction, since each HRU corresponds to unique combinations of input parameters. Hypothetically, the lower the threshold, the less the watershed spatial information loss and the more accurate the model output.

Limited investigations have been conducted to assess the influences of HRU thresholds on SWAT hydrologic predictions and most of them focused on the single land-use threshold or land-use and soil thresholds, besides which, hydrologists have mainly investigated the effects on monthly or annual streamflow predictions, while the relationship between HRU threshold and daily streamflow remains to be clarified (Chen and Mackay, 2004; Migliaccio and Chaubey, 2008; Yacoub and Foguet, 2013; Chiang and Yuan, 2015; Her et al., 2015). Although the conversion of land-use distribution has attracted more attention from researchers, this result did not adequately explain the mechanism of hydrologic response to HRU thresholds. For instance, Chiang and Yuan (2015) and Her et al. (2015) demonstrated the variation of land-use composition due to different land-use and soil threshold settings, but the changes in soil type and slope class were ignored. Moreover, without universally accepted guideline available, threshold values were arbitrarily set by the user. Thresholds between 5% and 20% are commonly used in SWAT modelling. Hence, a more detailed analysis should be conducted to assess the effects of various HRU thresholds on model performance and identify suitable HRU thresholds.

The chief objectives of this study are to 1) investigate how various HRU thresholds affect the characteristic of watershed land-use, soil type, and slope class, 2) evaluate the effects of various HRU thresholds on SWAT daily and monthly streamflow predictions, 3) explain the main reasons for the variation of streamflow predictions and model accuracy, and 4) suggest guidelines for selecting appropriate HRU thresholds in SWAT streamflow modelling. For the obvious advantage as a physically based hydrological model that includes free availability, user-friendly interface, readily accessible inputs, and extensive application fields, SWAT has gained international acceptance as evidenced by published articles and one of the most widely used hydrological models around the world. The significance of this study is to optimize the HRU threshold sets that contribute to reducing the uncertainties for watershed streamflow simulation and water resources management as well as further acceleration of flood forecasting and restoration.

SWAT is a watershed or river Basin scale hydrologic model widely used for simulating the hydrological processes and the associated sediment and agricultural chemical yields. The model predicts the total amount of water flowing to the main channel from each sub-watershed: the streamflow of the watershed outlet through the channel network is then calculated (Jin et al., 2018). The governing equations for the land phase are expressed as follows (Hooghoudt, 1940; Sloan and Moore, 1984; Arnold et al., 2011):

where WYLD is total amount of water leaving the HRU and entering main channel during the time step (mm), SURQ_CNT stands for surface runoff contribution to streamflow in the main channel during time step (mm), LATQ is lateral flow contribution to streamflow (mm), GWQ stands for groundwater contribution to streamflow (mm), and TLOSS is transmission losses (mm).

where SW_ly,excess stands for the drainable volume of water stored in the saturated zone of the hillslope per unit area (mm), K_sat is the saturated hydraulic conductivity (mm/h), slp stands for average slope of the HRU (m/m), ϕ_d is the drainable porosity of the soil layer (mm/mm), and L_hill stands for the hillslope length (m).

where K_sat is the hydraulic conductivity of the aquifer (mm/day), L_gw stands for the distance from the ridge or HRU divide for the groundwater system to the main channel (m), and h_wtbl is the water table height (m). Surface runoff is estimated by SCS (Soil Conservation Service) curve number using the procedure presented below (SCS, 1972):

where SURQ_GEN is the accumulated runoff (mm), R_day is the rainfall depth for the day (mm), I_a is the initial abstractions including surface storage, interception and infiltration prior to runoff, which is commonly approximated as 0.2 S (mm), and S is the retention parameter (mm).

where CN is the curve number for the day. When the retention parameter varies with plant evapotranspiration (Neitsch et al., 2011):

where S_prev is the retention parameter for the previous day (mm), E₀ is the potential evapotranspiration for the day (mm/d), cncoef is the weighting coefficient used to calculate the retention coefficient for daily curve number calculations dependent on plant evapotranspiration, and S_max is the maximum value the retention parameter can achieve on any given day (mm). The time of concentration is calculated by summing the overland flow time and the channel flow time is defined as (Chow, 1959; Engman, 1986; Neitsch et al., 2011):

where t_conc is the concentration time (h), t_ov is the concentration time for overland flow (h), and t_ch is the concentration time for channel flow (h).

where L_slp is the HRU slope length, and n is Manning’s roughness coefficient for the HRU.

where L is the channel length from the most distant point to the sub-watershed outlet, n is Manning’s roughness coefficient for the channel, Area is the HRU area, and slp_ch is the channel slope. The amount of surface runoff released to the main channel is calculated (Neitsch et al., 2011):

where SURQ_GEN is the amount of surface runoff generated in the HRU on a given day (mm), SURQ_STOR_i−1 is the surface runoff stored or lagged from previous day (mm), and surlag is the surface runoff lag coefficient.

The upstream region of Hunhe River Basin is situated in Qingyuan County, in the northern part of Fushun City in Liaoning Province, China (Guo and He, 2013). Flood disaster frequently appeared in this region, especially in recent years. The catchment originating from the headwaters of the Hunhe River to the monitoring sites of Beikouqian has been selected as the case study watershed for testing the effects of HRU thresholds on the streamflow modelling (Fig. 1). This catchment, with an area of 1840 km², which covers almost half of Qingyuan County, is one of the two major tributaries to Dahuofang Reservoir. The topography of the area is characterised by low mountains and foothills (Ye et al, 2014). The elevation varies between 155 m and 1083 m and the average altitude is 478 m above mean sea level. The slope of the catchment ranges from 0 to 114.5%, with mean slope and median slope extracted from the DEM of 19.4% and 18.9%, respectively. Land-uses within the catchment can generally be described as mixed forests, apart from areas supporting rural residences and agricultural land along the channels. It is dominated by a continental monsoon climate typical to medium latitudes, characterised by a hot summer and a cold winter. The monthly average temperatures in January and July from 1957 to 2014 (calculated by daily mean temperature data gathered from China Meteorological Data Service Centre for the Qingyuan Weather Station) are −15.08℃ and 23.03℃ , respectively. The catchment receives an annual average precipitation of 797.93 mm (as calculated by daily precipitation data gathered from China Meteorological Data Service Centre for the Qingyuan Weather Station during the period 1957–2014), most of which occurs in months from June to August (62.6%) and shows remarkable rainfall variability.

Figure 1.  Location of the the Upper Hunhe River watershed with hydro-meteorological stations and DEM data

Spatial datasets of DEM, land-use and soil combined with attribute datasets of soil physical properties and weather data were used to prepare model input parameters. DEM data for the catchment were downloaded from the International Scientific Data Service Platform of the Chinese Academy of Science (http://www.cnic.cas.cn/zcfw/sjfw/gjkxsjjx/) at a spatial resolution of 30 m in grid format (Zhao, 2016). The land-use data were derived from ALOS (Advanced Land Observing Satellite) multispectral images (10-m resolution) for the year 2008. Eight classes of land-use pattern were interpreted according to SWAT land cover database: 1) forested land (FRST, 78.77%), 2) cultivated land (CORN, 16.54%), 3) residential land (URLD, 1.79%), 4) water body (WATR, 1.34%), 5) agricultural land (RICE, 1.20%), 6) industrial land (UIDU, 0.18%), 7) grassland (PAST, 0.11%), and 8) wetland (WETL, 0.07%). The soil map for the case study watershed was obtained from Land Resources Atlas of Liaoning Province (PECLP, 1987), and then boundaries of different soil types were manually digitised. Soil types were classified into five great groups based on genetic classification (Shi et al., 2004). Major soil types within the area were dark brown forest soil (DBFS, 47.34%) and brown forest soil (BFS, 30.30%), while bog soil (BS, 13.51%), meadow soil (MS, 6.77%), and paddy soil (PS, 2.08%) were also found in the watershed. Brown forest soil was subdivided into three sub-groups, namely cultivated meadow brown forest soil (CMBFS, 3.63%), typical brown forest soil (TBFS, 4.77%), and young brown forest soil (YBFS, 21.90%), for more accurate assignment of curve numbers over the area and more realistic representation in the model (Kumar et al., 2015). Dark brown forest soil was composed mainly of young dark brown forest soil (YDBFS). Soil physical properties were extracted from the reference (Yong and Yu, 2011) and then converted in accordance with the SWAT user soil database (Ye Yin, 2014; Zhang et al., 2014). Some 58 yr (from 1957 to 2014) of meteorological data including daily precipitation, air temperature, wind speed, sunshine duration, and relative humidity from Qingyuan Weather Station were gathered from the China Meteorological Data Service Centre (CMDC) in order to build the SWAT weather-generator database (Bieger et al., 2015). Daily river discharges from Beikouqian Station located at the outlets of this watershed were acquired from the Water Resources Department of Liaoning Province for the period 2002–2014.

In SWAT, initial HRUs were defined by overlaying land-use, soil type, and slope class maps. Each HRU is considered a hydrologically homogeneous area consisting of unique land-use-soil-slope combination within a sub-watershed. The combination and distribution of final HRUs were determined by a user-specified HRU threshold during the process of multiple HRUs definition (Winchell et al., 2010). The HRU threshold was sequentially applied to the three overarching categories, namely land-use threshold, soil threshold, and slope threshold. Any land-use that covered a percentage of the sub-watershed below the land-use threshold level was eliminated and reapportioned into the remaining land-uses regardless of the distribution of soil type and slope class on this land-use area (Winchell et al., 2010; Her et al., 2015). As an example, in Fig. 2, the hypothetical sub-watershed was divided into eight initial HRUs (Table 1): if the land-use threshold was set to 10%, the area of CORN would be eliminated and the areas of PAST and FRST would be reapportioned as follows:

Figure 2.  A hypothetical sub-watershed with three land-use categories, three soil types, and three slope classes consisting of nine land-use-soil-slope combinations (Patch number)

Accordingly, soil threshold was only examined within the remaining land-uses (Her et al., 2015). If the area percentage of soil type within a certain land-use area was below the soil threshold level, that area was eliminated and reapportioned into the other qualified soil type(s) in this land-use area. For this hypothetical sub-watershed (Table 2), if the soil threshold was set to 10%, the area of DBFS on PAST would be eliminated and the area of MS on PAST in the sub-watershed would be modified as follows:

The slope threshold is then applied to the remaining land-use-soil combinations using a similar elimination process. For this hypothetical sub-watershed (Table 3), if the slope threshold was set to 15%, area of slope class > 15° on FRST/DBFS would be eliminated and the area of slope class 5°–15° for FRST/DBFS in the sub-watershed would be modified as follows:

The combination and distribution of final HRUs for the hypothetical sub-watershed are presented in Table 4.

ArcSWAT was utilised in this application to define sub-watersheds and HRUs for setting up SWAT and running the model. A critical source area of 3580 ha was chosen to determine the stream network and the watershed was partitioned into 33 sub-watersheds. Slopes calculated from the DEM were divided into four classes: 1) < 5° (0–8.75%), 2) 5°–15° (8.75%–26.80%), 3) 15°–25° (26.80%–46.63%) and 4) > 25° (46.63%–100%). Multiple HRUs were defined within each sub-watershed by setting thresholds of 0%/0%/0% (land-use/soil/slope), which resulted in 1468 HRUs representing the watershed. Surface runoff was estimated by using SCS curve number method and the retention parameter varied temporally as a function of accumulated plant evapotranspiration.

Simulations were performed from 2002 to 2014 including the first three-year of warming-up period. The model was calibrated for the period 2005–2009 and validated from 2010 to 2014 for streamflow at Beikouqian Station. Model performance was assessed by comparison of simulated and observed streamflow using standard statistical parameters: the correlation of coefficient (R²) and Nash-Sutcliffe coefficient (NS):

where Q_m and Q_s stand for measured and simulated discharges (m³/s), respectively, and i is the ith measured or simulated data.

A total of 20 scenarios were investigated by running the model with various combinations of HRU thresholds in this study. HRUs were delineated by applying five threshold levels for land-use (Scenario No. 1–5), five threshold levels for soil (Scenario No. 6–10), five threshold levels for slope (Scenario No. 11–15), and five threshold levels for all the three categories (Scenario No. 16–20), where thresholds of 0%/0%/0% were used as the baseline scenario (Table 5). All scenarios were given the same initial parameters, input data and calibrated parameters of the baseline scenario. Model output uncertainty due to various HRU threshold levels was compared using R² and NS (2005–2014), as well as the relative error (RE) defined as:

where Q_SCNR is the predicted stream discharges (taken from the output.rch file) with different HRU threshold levels (m³/s), and Q_base represents stream discharges for the baseline scenario (m³/s).

Model performance in terms of streamflow was firstly assessed for the baseline scenario. Although accurate predictions were not the goal of this study, the calibration and validation served to establish a realistic base line for making relative comparisons of the model output (Chen and Mackay, 2004). Fig. 3 shows the observed and simulated discharges at the outlet of the study watershed in the simulation period (2005–2014). R² and NS values were 0.60 and 0.58 for daily streamflow, while these values for monthly streamflow reached 0.92 and 0.90, respectively. Overall, the predicted stream discharges were reasonably matched with the observations although the simulated discharges were slightly over-predicted during the simulation period.

Figure 3.  Hydrographs of daily (a) and monthly (b) stream discharge observations and SWAT simulations for the baseline scenario during the simulation period at Beikouqian Station (the y-axis shows daily average values)

As depicted in Fig. 4, the number of HRUs decreased sharply with increasing HRU thresholds. Thresholds of 5%/5%/5% resulted in only 31.40% of the total number of HRUs in the baseline scenario, while thresholds of 25%/25%/25% (the maximum threshold combination used), resulted in the number of HRUs accounting for only 6.06% of those in the baseline scenario. The number of HRUs also decreased when slope, soil, and land-use thresholds were increased, and the number of HRUs was more sensitive to the land-use threshold than the slope and soil thresholds. For example, thresholds of 5%/0%/0% resulted in 701 HRUs for the whole watershed, even less than the number acquired under the maximum threshold levels of slope and soil (0%/0%/25% and 0%/25%/0%).

Figure 4.  Relationship between various thresholds and the number of hydrologic response units

As the slope threshold increased, the slope class distribution of the watershed was substantially changed, while the composition of land-use and soil type was kept constant. Greater portions of minor slope classes (< 5° and > 25°) were regrouped into slope class of 5°–15°, which was the dominant slope class of the watershed (Table 6). For the slope threshold of 25%, the percentage of slope class 5°–15° was exaggerated by 17.21% compared to the baseline scenario, while the areas of slope class < 5° and > 25° in the watersheds were decreased by 34.60% and 99.65%, respectively. The areal proportion of slope class 15°–25° increased until using thresholds of 0%/0%/10%, then decreased.

The distribution of soil type and slope class of the watershed was continuously changed when increasing the soil thresholds, while the land-use composition was kept constant. The direction of soil distribution changes varied among different soil types (Table 7). The areal proportion of YDBFS (the dominant soil type) expanded with an increase in soil threshold, while an opposite trend was observed in the case of BS. The area of CMBFS decreased from 66.88 km² for the baseline scenario to 63.71 km² for the 5% soil threshold but increased to 66.28 km² for the soil threshold of 10%, then decreased with an increase in soil threshold. The areas of YBFS and MS varied erratically, while the other minor soil types, namely TBFS and PS, showed an insignificant change and the areal proportions thereof fluctuated as the soil threshold increased. Table 7 also shows that the slope class composition showed a slight change as the soil threshold was altered.

The areas of land-use categories, soil types and slope classes from different land-use thresholds, depicted in Table 8, which showed the discrepancies between them and the baseline scenario. When the land-use thresholds were increased, minor land-uses were decreased significantly, while FRST and CORN (the two major land-use categories in the watershed) even became the only land-use categories therein (25%/0%/0%). Table 8 also shows that the soil type configuration was altered less than the soil threshold. The major soil types (YDBFS and YBFS) tended to increase when the land-use threshold increased, whereas the areas of BS, MS, and PS decreased slightly. The area of other soil types (CMBFS and TBFS) was generally well preserved. In addition, Table 8 shows that the slope class configuration was more sensitive to the land-use threshold than the soil threshold. The area of slope class < 5° decreased from 300.99 km² to 188.48 km², where the area of slope classes 5°–15°, 15°–25°, and > 25° increased to 43.16 km², 65.56 km², and 3.79 km², respectively.

Land-use, soil, and slope class distributions from different HRU threshold settings are listed in Table 9. The land-use composition was the same as the land-use threshold, while the distributions of soil type and slope class were significantly changed. Areal proportions of minor soil type and slope class decreased significantly and were integrated into the preserved dominant classes with an increase of HRU threshold. For instance, thresholds of 25%/25%/25% resulted in only 38.46% of BS remaining (covering only 5.20% of the watershed). Similarly, the proportion in slope class < 5° decreased from 16.36% to 3.28% and that in slope class > 25° decreased to zero as HRU thresholds were increased from 0% to 25%.

Evaluation statistics showed an insignificant change for all slope threshold scenarios, which reflected there was no clear relationship between slope threshold and streamflow. According to the model predictions, average daily stream flows for the five slope threshold scenarios were similar to the baseline scenario, with RE between −0.33% and 0.13% (Fig. 5(a)). Except for the thresholds of 0%/0%/10%, R² and NS values were observed to be 0.60 and 0.58 for daily streamflow and 0.92 and 0.90 for monthly streamflow, respectively. Model evaluation statistics indicated a close agreement between the simulated and observed stream discharges for all slope threshold scenarios.

Figure 5.  Evaluation statistics for different hydrologic response unit threshold scenarios

As shown in Fig. 5a, increasing the soil threshold from 0 to 25% caused a persistent over-prediction of average daily streamflow, and the maximum RE was 3.08% for the 25% soil threshold. Figs. 5b–5e show decreasing values of R² and NS as the soil threshold was increased from 0 to 25%: this indicated that the baseline scenario was the best simulation.

Land-use thresholds exerted a limited influence on average daily streamflow in this watershed. As the land-use threshold level was increased, average daily streamflow tended to increase, leading to a slight increase in RE (Fig. 5a). Figs. 5b–5e show that model evaluation statistics decreased as the land-use threshold was increased from 0 to 25%. Daily R² and NS values decreased by 4.85% and 5.89%, respectively, when the land-use threshold level was increased to 25% whereas monthly model evaluation statistics decreased to a lesser extent, and did so monotonically, with R² exceeding 0.9 and NS of 0.887 to 0.901.

As the HRU threshold increased from 0 to 5%, 10%, 15%, 20%, and 25%, the average daily streamflow was over-predicted by 0.45%, 1.25%, 2.78%, 3.83%, and 4.51% compared to the baseline scenario, respectively (Fig. 5a). R² and NS followed the opposite trend to RE with increasing HRU threshold. The decrement of daily R² induced by HRU threshold was 1.00%, 2.84%, 5.85%, 11.04%, and 12.88%, respectively, whilst daily NS was reduced by 1.39%, 3.81%, 7.28%, 13.87%, and 14.73%, respectively. Monthly R² and NS values were also decreased due to the increasing HRU thresholds but such increases were deemed negligible compared with those in the daily values.

Given the mountainous nature of the watershed, the daily streamflow was expected to change significantly with increasing slope thresholds, however, the five scenarios predicted similar values at the outlet of the study watershed. Daily surface runoff generated in HRUs was not directly related to changes in slope based on Equs. (4)–(6). These changes only affected overland flow time and channel flow time of concentration that may cause slight changes in the daily contribution of surface runoff to streamflow in the main channel (Equs. (7)–(10)). The most important changes in daily streamflow when changing the slope of the HRUs arise mainly due to daily lateral flow contribution to streamflow. Following Equ. (2), decreasing the slope would cause a decrease in slp and an increase in hillslope length that combine to cause a decrease in daily lateral flow to the main channel. In this study watershed, steep slope class and flat ground were continuously converted into intermediate slope class as the slope threshold increased from 0 to 25%, and areal proportion of slope increment was approximately equal to the slope decrement, resulting in a relatively stable lateral flow. Consequently, the average daily streamflow of the watershed underwent smoothed fluctuations with changes in slope threshold and followed the same trend with the average slope of the watershed.

In contrast to slope threshold, average daily streamflow was persistently over-predicted across the five soil threshold scenarios and the baseline scenario. Indeed, this increment is mostly attributed to the differences of saturated hydraulic conductivity between different soil types. Mean values of saturated hydraulic conductivity for strata in the watershed were 162.41, 24.79, 9.65, 72.73, 64.45, 13.66, and 12.58 for each soil type, corresponding to YDBFS, CMBFS, TBFS, YBFS, BS, MS, and PS, respectively. Consequently, area-weighted means of saturated hydraulic conductivity for the watershed increased consistently from 104.06 to 104.75, 105.42, 107.91, 111.04, and 111.46 caused by the change of soil type distribution when changing the soil thresholds. As described in Equs. (2) and (3), LATQ and GWQ increase as the area-weighted average values of watershed saturated hydraulic conductivity were increased, since other input parameters in the equations varied slightly or not at all. Even though increasing curve numbers in response to changes in soil type composition would cause a slight decrease in surface runoff according to Equs. (4) and (5), it is not enough to counteract the increases in LATQ and GWQ.

As the land-use threshold level increased, average daily streamflow followed the same trend as the soil threshold but to a lesser extent. On the one hand, increasing the land-use threshold altered the composition of watershed land-use as well as its soils and thus affected the surface runoff through the input parameters of curve numbers. Since FRST has the lowest curve number, surface runoff was reduced when other classes of land-use were redistributed to FRST based on Equs. (4) and (5). Likewise, YDBFS and YBFS have higher infiltration rates than BS, MS, and PS resulting in lower hydrologic group curve numbers and surface runoff in comparison to the baseline scenario when areal proportions of YDBFS and YBFS increased. On the other hand, decreases in the area of BS, MS and PS and increases in the area of YDBFS and YBFS lead to increased LATQ and GWQ which were altered much more than surface runoff. More importantly, slope class < 5° was consistently regrouped into steeper slope classes, increasing the average slope of the HRUs, further resulting in greater lateral flow and increasing the overall trend in average daily streamflow as land-use thresholds increased.

According to Fig. 5, HRU thresholds have a greater influence than the threshold of any single category on average daily streamflow for the same threshold level. This can be attributed to the combined effects of model input parameters variation, which is directly related to the changes in watershed land-use, soil type, and slope class composition as discussed previously. The key factors affecting stream discharge when changing the HRU thresholds were saturated hydraulic conductivity, average slope of the HRU, and curve number, which were different from the factors affected by drainage area thresholds (or critical source area, CSA). Lin et al. (2020) and Lin et al. (2021) pointed out that average sub-watershed areas, reach lengths, channel density, average channel slope and width varied substantially due to drainage area thresholds and consequently resulted in a considerable variability in estimated streamflow and sediment yields.

Generally, average streamflow was not greatly affected by HRU threshold (measured by RE), which agrees with the results of Her et al. (2015) and Wang et al. (2016). In addition, various HRU thresholds had little effect on monthly R² and NS values. Model performance for daily discharges was more sensitive than monthly results in terms of evaluation statistics, which indicated that daily simulation results exhibited greater fluctuations than the observed values compared to monthly simulation results. It should be noted that the components of streamflow (i.e., surface runoff, lateral flow, and groundwater) were more responsive to changes in HRU threshold as found by Chiang and Yuan (2015) implying that HRU thresholds may have a significant effect on sediment and nutrient-load predictions.

The above results also indicate that various HRU thresholds significantly affected model efficiency since the computational time is related to the number of HRUs. Besides, the number of HRUs was more responsive to HRU threshold than the threshold of any single category for the same threshold level. For daily streamflow predictions, thresholds of 5%/5%/5% resulted in numbers of HRUs between the thresholds of 20%/0%/0% and 25%/0%/0% which are much less than the maximum thresholds of slope and soil (0%/0%/25% and 0%/25%/0%) but returned better model performance than when using thresholds of 20%/0%/0%: this was considered to have been the optimum HRU threshold for the watershed studied in the present work (Fig. 6). For monthly streamflow predictions, thresholds of 25%/25%/25% would achieve greater computational efficiency.

Figure 6.  Relationship between evaluation statistics and hydrologic response units numbers for different hydrologic response unit threshold scenarios

It was difficult to predict the direction and magnitude of land-use, soil, and slope changes under the influence of HRU threshold in other watersheds, which was directly related to the streamflow prediction, because the elimination and reapportion of minor HRU combinations were calculated separately across the sub-watershed scale and each sub-watershed would have stochastic distributions of land cover, soil, and slope. Besides, increased model output uncertainty would be introduced depending on other factors such as spatial data classification criteria, property disparity between different spatial data categories, parameter adjustment method, and meteorological variability with additional precipitation data from rain-gauge stations being used. Although the aforementioned factors would be likely to change the discrepancy between the simulated result and the observed value, the effects of HRU thresholds on streamflow modelling were generally limited and are not compatible with the effect on sediment and nutrient-load predictions (Chiang and Yuan, 2015; Her et al., 2015).

This study shows that applying a non-zero HRU threshold can allow management of the number of HRUs and address the issue caused by exceeding model computation limits. Since stream discharge was determined by the input parameters extracted from spatial data category, an alternative approach was to combine spatial data categories with similar properties. According to the basic characteristics of HRU threshold effects on SWAT streamflow modelling (as discussed previously), recommendations for the selection of a suitable HRU configuration depend on project objectives. It is found that the most commonly used or even larger thresholds are appropriate for monthly stream discharge prediction which balances the need for computational efficiency with model output accuracy because monthly streamflow was not found to be greatly affected by HRU threshold. As for daily stream discharge prediction, a single category of HRU threshold no more than 25% or HRU threshold level no more than 15% was the better option for achieving satisfactory model performance and the thresholds of 5%/5%/5% was the optimum HRU threshold level of the watershed to allow full consideration of model accuracy and efficiency. Additionally, HRU thresholds of 0%/0%/0% (or the lowest thresholds possible) should be used when the modelling is aimed at assessing hydrological responses to changes in land-use or evaluating the effects of management practices to preserve each unique landscape feature.