Quantifying the effects of watershed subdivision scale and spatial density of weather inputs on hydrological simulations in a Norwegian Arctic watershed

The semi-distributed model SWAT (Soil And Water Assessment Tool) (Neitsch et al. 2009) was developed to predict the impacts of human activities (Gassman et al. 2007) and climate change (Dile et al. 2013) on environment and water resources in large complex watersheds. Lumped models consider the entire watershed/basin as a single system (Devi et al. 2015); on the other hand, the semi-distributed models like SWAT divide the whole watershed/basin into smaller sub-watersheds/sub-basins (Daofeng et al. 2004; Dwarakish & Ganasri 2015). It is assumed that each sub-basin is a homogeneous unit with representative parameters for the entire sub-basin (Bingner et al. 1997). Choosing the size for the sub-basins also influences the homogeneous assumption because the larger the sizes of the sub-basins, the higher variable conditions the sub-basins have (Bingner et al. 1997). When the sizes of sub-basins are reduced and the number of sub-basins are increased, it significantly influences the amounts of required input data and model parameters, the computational process (Bingner et al. 1997), as well as the calibration effort since the large number of sub-basins may require more adjusted model parameters needed to optimize the simulation results and more iterations needed for running the calibration (Rouhani et al. 2009).

Watershed delineation is considered as an important preliminary step since the accuracies of the modeling results, e.g., runoff (Rouhani et al. 2009; Gong et al. 2010; Chaplot 2014), streamflow (Norris & Haan 1993; Rouhani et al. 2009), and soil erosion and pollution (Gong et al. 2010) may be influenced by the delineation resolutions, beside the quality of input data (Chaplot et al. 2005; Ning et al. 2015; Nazari-Sharabian et al. 2020). For example, the accuracy of runoff simulation results decreases at the coarse levels of watershed discretization due to the effects of the changes in the distribution of runoff curve numbers over the entire watershed, particularly in the SWAT model (Rouhani et al. 2009). Increasing or decreasing numbers of sub-watersheds also influences the accuracy of simulation of peak flows or extreme flows (Rouhani et al. 2009). In particular, the deviation between observed and simulated peak flows/extreme flows increases when the number of sub-watershed increases. This is because higher numbers of sub-watersheds lead to higher variation of the runoff values that contribute to streamflow. In addition, the values of runoff curve numbers are automatically updated according to the variation of soil moisture condition in each sub-watershed. Therefore, when the runoff curve numbers have high fluctuation, then the values of runoff contributing to streamflow also highly fluctuate. As a result, the accuracy of simulated peak flows/extreme flows is influenced (Rouhani et al. 2009). Moreover, delay in the travel time of runoff occurring in the watersheds with large numbers of sub-watersheds may result in lower values of simulated peak flows compared to that with small numbers of sub-watersheds. The reason is because runoff from the upper sub-watersheds could reach to the outlet of the watershed only after runoff from the lower sub-watersheds has been already discharged (Rouhani et al. 2009). Furthermore, increasing the number of sub-watersheds leads to increasing channel slope and drainage density that result in higher simulation results of some water balance components (Chen et al. 2021). Also, the change in drainage density also influences the accuracy of runoff prediction (Goodrich 1992). Finally, the automated computational processes of morphological and hydrological parameters of the watershed are strongly influenced by the chosen numbers and sizes of sub-watersheds (Munoth & Goyal 2019b).

The most natural subdivision method is dividing a watershed into its natural sub-watersheds based on topography data, which are extracted from a digital elevation model (DEM). Such watershed subdivision also aims to preserve the watershed's natural boundaries, flowpaths, as well as channels for realistic flow routing (Zhang et al. 2004). With the development of GIS (geographic information system) technologies, several watershed subdivision approaches have been developed to investigate the impacts of watershed subdivision on modeling outputs (Savvidou et al. 2014), including (1) critical source area (CSA) (Thieken et al. 1999; FitzHugh & Mackay 2000; Kalin et al. 2003; Di Luzio & Arnold 2004; Arabi et al. 2006); (2) threshold drainage area (TDA) (Nour et al. 2008); (3) aggregated simulation area (Lacroix 1999); (4) representative elementary areas (Wood et al. 1988); (5) representative elementary watershed (Reggiani & Rientjes 2005); (6) hydrologic similar units (Karvonen et al. 1999); (7) functional units (Argent et al. 2006); and (8) hydrologic response units (HRUs) (Flugel 1995, 1997). In the SWAT model, watershed subdivision is basically based on the TDA, which is the minimum upstream drainage area for a channel to originate (Aouissi et al. 2013), or as a percentage of total catchment area (Di Luzio & Arnold 2004; Kumar & Merwade 2009). Each sub-basin in SWAT is further subdivided into smaller HRUs.

Many previous studies around the world applied the SWAT model to investigate the impacts of watershed subdivision on the results of hydrological simulations, including runoff (Norris & Haan 1993; Bingner et al. 1997; Jha et al. 2004; Arabi et al. 2006; Rouhani et al. 2009; Chaplot 2014; Munoth & Goyal 2019a), water balance components (Tripathi et al. 2006; Chaplot 2014; Chen et al. 2021), and streamflow (Mamillapalli et al. 1996; FitzHugh & Mackay 2000; Haverkamp et al. 2002; Jha et al. 2004; Muleta et al. 2007; Rouhani et al. 2009; Aouissi et al. 2013, 2018; Chiang & Yuan 2015; Ozdemir et al. 2017; Pignotti et al. 2017; Chen et al. 2021). Regarding the studies on runoff, there was inconsistency among previous findings. For example, in a study in the 21.3-km² Goodwin Creek Watershed, in northern Mississippi, USA, it was found that values of runoff volume generated from 10 different levels of watershed subdivision based on different values of the CSA were not significantly impacted by the chosen number and size of sub-watersheds (Bingner et al. 1997). For example, the simulation of total annual runoff volume varied less than 5% among 10 watershed subdivision schemes (Bingner et al. 1997). In contrast, little variation in the total simulated surface runoff among 12 sub-watershed delineation schemes was detected in a study in four Iowa watersheds, USA, with the areas varying from 2,000 to 18,000 km² (Jha et al. 2004). However, the variation in the total simulated surface runoff was not clearly quantitative in such a study (Jha et al. 2004). Also, changes in runoff simulation were found in a study in the 384-km² Grote Nete River catchment, in Flanders, Belgium (Rouhani et al. 2009). They pointed out that the larger number of sub-watersheds in the watershed delineation schemes resulted in higher variation of runoff that led to higher fluctuation in the values of simulated extreme flows (Rouhani et al. 2009). Nevertheless, the deviation of simulated peak flow among watershed delineation schemes was not clearly quantified (Rouhani et al. 2009). Another study in the 65,145-km² Tapi River, in India, concluded that surface runoff decreased (approximately 35%) when TDA increased (from 25 to 400 km²) (Munoth & Goyal 2019a). In contrast, a study in a 26.12-km² flat watershed, the Walnut Creek watershed, in central Iowa, USA, found that surface runoff increased (approximately 15%) when TDA increased (from 23 to 654 ha) (Chaplot 2014). However, the larger relative errors, between observed and simulated results, were also detected for the coarse watershed subdivision schemes with higher TDA values (Chaplot 2014). For example, the relative error for estimated runoff from TDA 654 ha was approximately 15%, while it was 0% from TDA 23 ha or approximately 8% from TDA 100 ha (Chaplot 2014). Therefore, it was learned from that study that using the coarse watershed subdivision schemes could produce higher values of surface runoff volume, but at the same time, the simulated results were more uncertain compared to those by using the fine schemes.

Similar to the studies in runoff, the changes in the results of streamflow simulations under different watershed delineation schemes were also inconsistent among previous studies. For example, it was found in a study in the four Iowa watersheds, USA, that streamflow components were increased only less than 7% when the number of sub-watersheds increased (from 5 to 53) for a 1,929-km² watershed, which indicated quite insensitive streamflow to the number of sub-watersheds (Jha et al. 2004). In other studies, the changes in watershed subdivision schemes had slight impacts on the results of streamflow simulation (FitzHugh & Mackay 2000; Aouissi et al. 2013; Chiang & Yuan 2015). For example, the simulated annual and monthly streamflow, in a 62-km² Pheasant Branch watershed, Dane County, Wisconsin, USA, slightly increased (approximately 12%) from the coarse schemes to the fine schemes (FitzHugh & Mackay 2000). In a study in the 418-km² Joumine watershed, northern Tunisia, the simulated annual and monthly streamflow were only few percentage of variation among watershed delineation schemes (Aouissi et al. 2013). In a study in the large-scale watershed, the Kaskaskia River watershed in Illinois, USA, with 11,350 km², the simulated average annual streamflow increased by less than 2% from the finest scheme to the coarsest scheme (Chiang & Yuan 2015). However, a study in a small-scale watershed, the 26.12-km² Walnut Creek watershed, in central Iowa, USA, found significant increase of mean streamflow (approximately 62%) when TDA increased from 23 to 654 ha (Chaplot 2014). However, that study also pointed out that simulation of streamflow by using the coarse watershed subdivision solution produced higher uncertain results than by using the finer watershed subdivision solution (Chaplot 2014). For example, the relative error between observed data and simulated results of mean streamflow was 163% for TDA 654 ha, while it was less than 6% for TDAs 23, 100, and 261 ha (Chaplot 2014). Another study in the 4,297-km² Bosque River watershed, Texas, USA, found the positive effects of changes in watershed subdivision schemes on the accuracy of streamflow prediction (Mamillapalli et al. 1996). Such study revealed that the accuracy of streamflow prediction was improved by increasing the number of sub-watersheds and/or the number of HURs (Mamillapalli et al. 1996). This finding was also similar to findings from other studies in the Weiherbach (6.3 km²) and Dietzhoelze (81.7 km²) watersheds, Germany, and the Bosque River watershed (4,297 km²) in Texas, USA (Haverkamp et al. 2002; Muleta et al. 2007). In contrast, the accuracy of peak flow prediction was concluded to be decreased (approximately 20%) when numbers of sub-watersheds increased in a study in the 384-km² Grote Nete River catchment, in Flanders, Belgium (Rouhani et al. 2009). Recently, a study in the 491,700-km² Upper Mississippi River Basin, USA, found that the fine schemes of watershed delineation (12-digit sub-basin scenario) yielded higher values of streamflow simulation (approximately 1.79–7.17%) compared to the coarser schemes (8-digit sub-basin scenario), since the fine schemes are able to capture a sophisticated level of the spatial variation of watershed features including variation of rainfall regime (Chen et al. 2021). This finding agreed with the finding of a previous study in the 152.29-km² Little Washita watershed, near Chickasha, Oklahoma, USA, since they pointed out that the number of sub-watersheds strongly impacted the simulated streamflow hydrograph (Norris & Haan 1993). In particular, the estimated peak flow increased approximately 30% when the number of sub-watersheds increased from 1 and 2 sub-watersheds up to 15 sub-watersheds (Norris & Haan 1993).

Regarding studying water balance components, a previous study in the 90.23-km² Nagwan watershed in India revealed that the size and number of sub-watersheds had significant impacts on the simulation results of evapotranspiration (ET), percolation, and soil water (SW) content, with the exception of surface runoff (Tripathi et al. 2006). For example, increasing the number of sub-watersheds (from a single watershed to the discretization of 12 and 22 sub-watersheds) resulted in increasing approximately 0.28–61.4% ET, 4.48–26.7% percolation, and 17.7–22.3% SW content (Tripathi et al. 2006). A study in a 26.12-km² flat watershed, the Walnut Creek watershed, in central Iowa, USA, found that decreasing sizes of sub-watersheds (from 654 ha down to 23 ha) or increasing numbers of sub-watersheds (from 1–4 to 96–115) resulted in increased ET (4.6%), decreased SW (5.1%), decreased percolation (2.8%), decreased surface runoff (15.1%), and decreased groundwater (2.4%) (Chaplot 2014). In this case, the trends of changes were not consistent for every water balance component, particularly in ET (Chaplot 2014). Another study in the 4,91,700-km² Upper Mississippi River Basin, USA, found that the simulation results of surface runoff, lateral flow, groundwater flow, and water yield (WYLD) increased approximately 0.98, 92, 2.73, and 2.07%, respectively, when the number of sub-watersheds and HRUs increased (Chen et al. 2021).

However, the fine watershed delineation does not always yield higher model performances compared to the coarse schemes (Boyle et al. 2001; Reed et al. 2004; Rouhani et al. 2009). For example, some previous studies stated that the modeling results are better when using the semi-lumped and semi-distributed model structures compared to those using the distributed models (Ogden & Julien 1994; Smith et al. 2004; Das et al. 2008). The reason is that the capacities to capture the important features of the watersheds and variation of rainfall regime of the coarse watershed delineation are better than those of the fine one. This argument is somewhat against the argument in the study in the Upper Mississippi River Basin, USA (Chen et al. 2021).

In another method of assessment, a study in two small-scale watersheds, the Dreisbach (6,23 km²) and the Smith Fry (7,30 km²) watersheds, in Maumee River Basin, Allen County, Indiana, USA, demonstrated the importance of manner of watershed subdivision on the efficiency of different best management practices (BMPs) for controlling the fate and transport of nutrients (e.g., total nitrogen and total phosphorus) and sediment within the watersheds (Arabi et al. 2006). Herein, nutrients and sediments are transported into the channels by surface runoff and lateral subsurface flow (Arabi et al. 2006). Therefore, the changes in these water balance components because of different watershed subdivision solutions could potentially affect the estimation of nutrient and sediment outputs. Besides that, the study found that watershed subdivisions caused discrepancies in watershed characteristics, e.g., drainage density, channel networks which affected nutrient and sediment yields (Arabi et al. 2006). In particular, it was expected from the study that more studies in future should be focused on the larger watersheds to verify the impacts of watershed subdivision scales on BMPs of the watersheds, since the larger watersheds may reveal different trends of changes compared to the smaller watersheds (Arabi et al. 2006). Also, it was highly expected from the study that the impacts of watershed subdivision should deserve more attention in future than those carried out in the past because of uncertainties resulting from different spatial resolutions (Arabi et al. 2006).

Beside the impacts of watershed delineation, density and spatial distribution of weather data input are also the important factors that may affect the modeling results (Chaubey et al. 1999; Bardossy & Das 2008; Aouissi et al. 2013, 2018; Chaplot 2014; Chen et al. 2021). For example, the high uncertainty in the estimated model parameters in the hydrological models results from using spatial homogeneity of rainfall and does not consider the refined variation of rainfall input (Chaubey et al. 1999). In addition, the performances of the hydrological models significantly decline when the density of integrated rain gauges is reduced (Bardossy & Das 2008). Also, the accuracy of streamflow simulation is significantly impacted by the spatial distribution of rain-gauge networks (Aouissi et al. 2013, 2018). It was found from a study in the 26.12-km² Walnut Creek watershed, in central Iowa, USA, that spatial resolution of rain-gauge networks significantly impacted water balance components in different manners (Chaplot 2014). For example, when increasing the number of rain gauges from 1 to 13, ET decreased approximately 17.7%, SW content increased approximately 41.3%, percolation increased approximately 66.67%, surface runoff decreased approximately 40.9%, and groundwater decreased approximately 42.1% (Chaplot 2014). Furthermore, it was learned from a study in the 491,700-km² Upper Mississippi River Basin, USA, that higher values of streamflow prediction were yielded from the denser weather networks compared to those from the scatter networks (Chen et al. 2021). For example, the simulated average monthly streamflow increased approximately 6.30–8.32% by using the denser climate dataset compared to using the sparser dataset in a large-scale watershed (Chen et al. 2021). However, the influences of weather networks density on hydrological simulations might vary from region to region and might depend upon the environmental characteristics or conditions of the investigated watersheds, where the differences in rainfall types (e.g., convective or advective), rainfall seasonality, the importance of snow accumulation and snowmelt processes, topographic features, and land use are identified (Bardossy & Das 2008; Chaplot 2014). The regions with complex hydrological processes might require high density of rain-gauge networks (Bardossy & Das 2008). Therefore, more studies are necessary to verify the effects of density and spatial distribution of weather data on hydrological simulations.

Obviously, the impacts of watershed subdivision and spatial resolution of weather networks on hydrological simulations were significantly investigated in numerous regions around the world, where there are differences in climate conditions, topography, land-use, and hydrological regimes. This contributed valuable knowledge for the scientific community; however, a consensus has not yet been obtained among outcomes under different environmental conditions (Mulungu & Munishi 2007; Chaplot 2014). For example, using the fine watershed delineation schemes could result in both positive and negative effects on the accuracy of hydrological simulations compared to those by using the coarse schemes. Similarly, the impacts of weather network density on hydrological simulations also varied regionally and could depend upon different environmental characteristics or conditions (Mulungu & Munishi 2007; Chaplot 2014). Therefore, the consequences of watershed subdivision and spatial density of weather networks on hydrological simulations are still a controversial issue. In addition, most of the previous studies were conducted in tropical/sub-tropical or temperate climate zones. However, studies in the Arctic region, with complex hydrological processes and sparse weather data, are still limited. Moreover, snowmelt runoff is an important component in the Arctic hydrology, since it contributes approximately 75% of the total annual flow in many Arctic watersheds (Woo 1980). Nevertheless, studies on the effects of watershed subdivision and weather network density on snowmelt runoff have not been addressed in previous studies. Also, many studies assessed the changes in streamflow, because of watershed subdivisions, mostly at the basin outlet or at certain hydro-gauging stations. However, the spatial variation of streamflow in each sub-basin, which is important for the case of flood hotspots analysis because of watershed delineation solutions, has not been well investigated. Therefore, to fill the existing knowledge gaps as well as to satisfy the expectation of previous studies, this paper conducted a study in the Arctic conditions to investigate the combined effects of watershed subdivision scale and weather network density on the results of hydrological simulations. In particular, outcomes of the present study aim to answer the following pertinent questions:

• 1.

  How much discrepancies in watershed characteristics and land-use composition could change?

• 2.

  How much model performance (in terms of statistical indicators) could be influenced?

• 3.

  How much snowmelt runoff volume and water balance components could vary spatially across the watershed?

• 4.

  How much streamflow could vary in temporal–spatial patterns across the watershed?

• 5.

  How much the sensitivity of model parameters under the Arctic conditions could be influenced?

Målselv watershed in northern Norway, distributing from 68°21′N to 69°17′N, was selected as the study area (Figure 1). This is a large-scale watershed with an area of approximately 5,913 km². Målselv has the features of mountainous terrain with the ground surface elevation ranges from 0 to 1,718 m. This area is located in the cold climate zone with the average annual air temperature varying from −5 to 6 °C. Rainfall regime is also highly variant across the watershed. The long-term average annual rainfall in this area fluctuates from ∼500 to 1,500 mm.

The routing phase describes several processes occurring in the stream including movement of water, sediments, flow mass in the channel, and transformation of chemicals in the stream and streambed.

The QSWAT interface, version 1.9, was used in this study. This is a coupling product of the hydrological model SWAT version 2012 and the open-source QGIS version 2.6.1.

In order to run the SWAT model, several temporal-spatial input data including time series of climate data such as precipitation, maximum and minimum air temperature, wind speed, relative humidity, solar radiation, and spatial (grid) data such as land use, soil, and topography (e.g., DEM) were required. In addition, time series of river discharge were needed for model calibration and validation. These data were collected from several sources. Details of data types, data resolution, and sources of data are summarized in Table 1.

The technique of watershed subdivisions in the SWAT model is based on the values of the TDA. Four different TDA schemes, using TDA values of 200, 2,000, 5,000, and 10,000 ha, were developed in this study.

Watershed subdivision and HRU creation were performed for each option of TDA values. Multiple HRUs were generated for each sub-basin, from the inputs of land use, soil, and slope classes, based on an HRU threshold (Figure 2). This threshold considers the percentage of the representative land use/soil/slope for each sub-basin. The HRU thresholds from 5 to 15% were widely used in many studies (Sexton et al. 2010; Srinivasan et al. 2010; Han et al. 2012; EPA 2013). In this study, the designed HRU thresholds for land use, soil, and slope were 5% for each. According to this threshold, only types of land use, soil, and slope, which are higher than 5% of the sub-basin area, were considered. In addition, the terrain slope was classified into five classes such as 0–5, 5–10, 10–25, 25–30, and >30%.

For each level of watershed subdivision, the model ran on a monthly time step from 1995 to 2012, including a 3-year warming up period (1995–1997). The 10-year period, from 1998 to 2007, was used for model calibration, and the remaining 5 years, from 2008 to 2012, were used for model validation.

The Sequential Uncertainty Fitting Version 2 (SUFI-2) algorithm in the SWAT Calibration Uncertainties Program (SWAT_CUP) (Abbaspour et al. 2007) was used for model calibration, model validation, parameters sensitivity, and uncertainty analyses. Each SWAT model corresponding to each watershed subdivision scheme was calibrated separately. A total of 18 model parameters, which were recommended as the most sensitive parameters for streamflow calibration (Abbaspour et al. 2007, 2015), were used in calibration and validation processes. Those model parameters are classified into five different subgroup processes of hydrological cycle, including (1) surface runoff (e.g., CN2.mgt and CANMX.hru); (2) lateral flow (e.g., ESCO.hru, SOL_AWC.sol, SOL_BD.sol, and SOL_K.sol); (3) snowmelt (e.g., SMTMP.bsn, TIMP.bsn, SMFMN.bsn, SMFMX.bsn, and SFTMP.bsn); (4) channel water routing (e.g., CH_K2.rte and CH_N2.rte); and (5) groundwater (e.g., ALPHA_BF.gw, GW_REVAP.gw, GWQMN.gw, REVAPMN.gw, and GW_DELAY.gw). The finest scheme used a total of 2,500 simulations to detect the optimal model parameters, while each of the other coarser schemes used a total of 2,000 simulations. In addition, cross-validation was approached to test whether or not the calibrated model parameters achieved from the finest watershed subdivision scheme could also perform well in other coarser schemes. All the possible results, which are found during the calibration process, are distributed in the so-called 95PPU band. The two statistical indicators, such as P-factor and R-factor, were used to measure the uncertainty of the calibration results. Herein, the values of P-factor range from 0 to 1, of which a threshold of 0.7 or 0.75 is suggested for river discharge calibration. The optimal values of R-factor, which presents the thickness of the 95PPU band, should be close to zero. For river discharge calibration, the value of R-factor is suggested to be smaller than 1.5. When the thickness of the 95PPU band is large, it means that the possibility of the model to capture most of the observed data is high; however, the model uncertainty is also high.

To identify whether or not a parameter is significant in sensitivity analysis, a t-test was used. This method uses two indicators, namely p-value and t-stat, to measure and rank the sensitive level of each calibrated model parameter. The hypothesis of the t-test method is that the larger the absolute values of t-stat, and the smaller the p-values, the more sensitive the parameters are determined. In addition, a parameter is considered significant in sensitivity analysis if the p-value calculated for that parameter is smaller than a value of 0.05. Finally, all calibrated model parameters are ranked for their sensitivity levels according to the magnitudes of t-stat and p-value.

The threshold values of the statistical coefficient R², NSE, and RSR for monthly simulation are shown in Table 2 (Santhi et al. 2001; Van Liew et al. 2003; Moriasi et al. 2007; Premanand et al. 2018).

To investigate the effects of watershed subdivisions on the hydrological simulations, the present study focused on the evaluation of the simulation results of water balance components and streamflow. Regarding water balance components, the annual mean values of total areal rainfall (PCP), actual ET, surface runoff (SUR_Q), snowmelt runoff, lateral runoff (LAT_Q), groundwater recharge amount (PERCO), groundwater contribution to streamflow (GW_Q), and WYLD (YIELD = SUR_Q + LAT_Q + GW_Q – Transmission losses) contributing to streamflow were calculated. In addition, to compare the spatial variation of such water balance components across the entire watershed among different TDA schemes, the GIS maps were produced. Furthermore, the long-term monthly average streamflow at five different hydro-gauging stations was analyzed. Also, the GIS maps of the spatial variation of the long-term annual mean streamflow were produced. The ArcGIS software version 10.6.1 was used in this study for generating the GIS maps and for spatial analysis. The results are discussed in the following section.

The methods of watershed delineation have significant impacts on the levels of terrain complexity as well as the topographic attributes. The number of sub-basins declined 75, 90, and 96% in the coarse TDAs 2,000 ha (115 sub-basins), 5,000 ha (48 sub-basins), and 10,000 ha (18 sub-basins), respectively, compared to the finest TDA 200 ha (459 sub-basins) (Table 3; Figure 3). When numbers of sub-basins increase and the sub-basin sizes decrease, the accurate level of the representative land uses for the watershed will be high. For example, the TDAs 200 ha (with 459 sub-basins) and 2,000 ha (with 115 sub-basins) presented a total of 11 main land-use groups for the entire watershed, while other coarser schemes such as TDAs 5,000 and 10,000 ha lost two and four land-use groups, respectively. In addition, the areas of each land-use group varied as the number of sub-basins declined, of which some land-use groups decreased in their areas while others increased, but the magnitudes of decreasing were greater than those of increasing (Table 4). For example, land-use groups of barren or sparsely vegetated, mixed forest, grassland, shrubland, bare ground tundra, water, and wooded wetland had declining trends, whereas the remainder had slight increasing trends. In particular, under TDA scheme 2,000 ha, the area of wooded wetland decreased approximately 70% compared to the finest scheme 200 ha, and this was the highest ratio among the decreased land-use groups. Regarding TDA scheme 5,000 ha, the highest percentage of decreasing area (approximately 75%) was in the group of barren or sparsely vegetated. For the coarsest scheme 10,000 ha, mixed forest had the highest percentage of declined area with approximately 55%. Noticeably, two land-use groups, evergreen needleleaf forest and wooded wetland, disappeared under TDA scheme 5,000 ha, while TDA scheme 10,000 ha lost four land-use groups, namely barren or sparsely vegetated, evergreen needleleaf forest, bare ground tundra, and wooded wetland. In contrast, other land-use groups including deciduous broadleaf forest, evergreen needleleaf forest (excluding TDA 5,000 ha and TDA 10,000 ha), mixed tundra, and wooded tundra slightly decreased in their areas compared to the finest scheme 200 ha. For land-use groups with increased areas, the highest percentage of increasing was only 6.5% and this was for the group of wooded tundra under TDA 10,000 ha. Obviously, land-use groups with small areas significantly decreased, even disappeared in the coarse schemes. This is because the threshold for the HRU definition in the present study was 5% for land-use/soil/slope. Therefore, the land uses with their areas smaller than 5% of the sub-basin areas were not defined, since they were regrouped into the major land-use groups. The decrease in the areas of the minor land-use groups in the coarse watershed subdivision schemes was also validated in some previous studies in the USA (Bingner et al. 1997; Chiang & Yuan 2015; Chen et al. 2021). Table 3 provides the summary of the discrepancies in watershed characteristics resulting from the changes in TDA schemes.

In addition to the changes in the presence of land uses over the watershed resulting from the changes (decreasing) in the number of sub-watersheds, other topographic attributes were also changed. For example, the increases of sub-basins’ sizes resulted in the changes in the average elevation and average overland slope of the sub-basins, which may affect the surface runoff process (Table 3). The finest TDA 200 ha had the highest average sub-basin slope compared to other coarse schemes. Such an increase in the overland slope could result from a better representation of spatial variation of surface elevation by discretization to smaller sub-watershed. In addition, the finest TDA scheme generated denser stream networks. For example, the TDA 200 ha produced approximately 1,922 km total stream length with 5 levels of stream order, while other coarser schemes 2,000, 5,000 and 10,000 ha produced approximately 1,015 km stream length (−47%) with 4 levels of stream order, 688 km stream length (−64%) with 4 levels of stream order, and 488 km stream length (−75%) with 3 levels of stream order, respectively (Figure 4). The decrease in stream length could affect some important in-stream processes. However, generating more sub-channels may not be realistic, since the channels from a very detailed level of sub-watershed description may only represent the low-lying areas in nature but they may not be existing/real channels. Moreover, the average stream slope remarkably declined (approximately 68–83%) from the finest scheme to the coarsest scheme. Drainage density dropped from 0.33 km km⁻² (TDA 200 ha) to 0.17 km km⁻² (TDA 2,000 ha), 0.12 km km⁻² (TDA 5,000 ha), and 0.08 km km⁻² (TDA 10,000 ha). The decrease in drainage density of the coarse TDAs may affect the accuracy of runoff simulation. However, one of the advantages of using the coarse watershed subdivisions is requiring less inputs, computation time (e.g., for running the model, calibration, and validation), as well as computer resources (e.g., reducing storage space). For example, in this study, in order to run 500 simulations for one iteration in the calibration process, it took around 56 h for the TDA 200 ha, while it was only 19.1, 7, and 3.2 h for TDA 2,000, TDA 5,000, and TDA 10,000 ha, respectively. Therefore, using a lower number of sub-watersheds could benefit some users in the cases of limitations in time and the available resources. However, for the watersheds with high variation of land uses, it may require a detailed analysis of sub-watersheds for sophisticated description of the important features of the watersheds. The decreases in drainage density, total stream length, stream order level, average sub-watershed slope, and average stream slope in the coarse schemes compared to those in the fine schemes in the present study also agreed with findings from previous studies, but the level of declines depended upon the drainage area of the study area as well as the designed TDA values (Bingner et al. 1997; Jha et al. 2004; Chiang & Yuan 2015; Munoth & Goyal 2019a; Chen et al. 2021).

During the calibration period, at Lille Rostavatn and Målselvfossen hydro-gauging stations, model performances were relatively stable under the changes of number of sub-watersheds. Although the number of integrated weather grids decreased due to the decrease of numbers of sub-watersheds, the selected weather grid points may be the correct representatives for the watershed as well as for the sub-watersheds surrounding these two hydro-gauging stations. The negligible impact of watershed subdivision on model performance was also validated in previous studies (Aouissi et al. 2013, 2018). In contrast, model performances at three remaining stations fluctuated under different watershed subdivisions. For example, at Høgskarhus station, model performance increased from TDAs 200 to 2,000 ha, then slightly decreased when numbers of sub-watersheds decreased. At Skogly station, model performance declined gradually when numbers of sub-watersheds decreased. At Lundberg, model performance highly fluctuated. For instance, it was stable from TDAs 200 to 2,000 ha, then decreased with TDA 5,000 ha, and afterward increased with TDA 10,000 ha. The decreases in model performances in the coarse schemes in the present study agreed with conclusions from the previous studies (Mamillapalli et al. 1996; Haverkamp et al. 2002; Tegegne et al. 2019). Obviously, model performances were heterogeneous among hydro-gauging stations under the changes of TDA schemes (Figure 5(a)). The reasons could be the complexity of hydrological processes as well as topographic characteristics in the Arctic.

The general tendencies of the changes in model performances in the validation period were repeated to those in the calibration period at Målselvfossen, Lundberg, and Høgskarhus, except at Skogly and Lille Rostavatn (Figure 5(b)). Details of the statistical coefficients for calibration and validation at all five hydro-gauging stations under different TDA schemes are presented in Tables 5–9.

Rainfall is one of the main inputs of water balance components. It was observed from the present study that the annual mean values of areal rainfall increased with the decreases in the number of sub-watersheds (Figure 6(a)). However, the magnitudes of deviations in annual mean values of areal rainfall were not significant among the TDA schemes. For example, the gap in the annual mean values of areal rainfall was only 24 mm between TDAs 200 and 10,000 ha. It could be interpreted that the coarse scheme had higher rainfall input than the finer scheme because fewer weather grid points in the coarse scheme generated more uniform rainfall across the watershed. In particular, the average and minimum values of annual rainfall from 14 weather grids in the TDA 10,000 ha were 1,207 and 826 mm, respectively, which were higher than those from 21 weather grids in the TDA 200 ha, with 1,185 and 750 mm, respectively. Therefore, the integrated rainfall amount from the coarse scheme was higher than that from the fine scheme, although the number of integrated weather grids from the coarse scheme was less than that from the fine scheme. Also, the denser weather grid points of the fine TDA scheme produced lower rainfall amount, since rainfall had high variation among weather grids. This could be true for the mountainous watershed, since rainfall is usually high variation. Figure 7(a)–7(d) illustrates the spatial variation of areal rainfall resulting from different resolutions of watershed discretizations. Obviously, the higher number of sub-watersheds produced finer variation of areal rainfall across the watershed. The finest scheme was able to display some locations in the watershed with low rainfall amount. For example, based on Figure 7(a), the minimum value of annual mean rainfall in some sub-basins in the upstream, from the finest TDA scheme 200 ha, was 758 mm which was lower at 87 mm compared to that from other coarser schemes.

Annual mean potential ET (PET) (Figure 6(b)) and actual ET (Figure 6(c)) slightly increased from TDAs 200 to 5,000 ha, then slightly decreased to TDA 10,000 ha. This could be because the coarse scheme 10,000 ha simplified the land-use and cropland variations that resulted in lower ET amount. For example, the coarse scheme 10,000 ha lost 4 land-use groups compared to 11 land-use groups in the fine scheme 200 ha (Table 4). In contrast, annual mean lateral flow (LATQ) dropped from TDAs 200 to 5,000 ha, after which it went up to TDA 10,000 ha (Figure 6(e)). The annual mean surface runoff (SURFQ) (Figure 6(d)) and WYLD (Figure 6(f)) had similar trends. In general, SURFQ and WYLD had an upward trend from TDAs 200 to 10,000 ha, but the trend dropped down to TDA 5,000 ha. Therefore, this could reveal that the TDA 5,000 ha could be a threshold which sketched the line of the discrepancies between watershed subdivision schemes. The water balance components could increase/decrease if the number of sub-watersheds was lower/higher than this threshold.

Although the annual mean values of PET and ET decreased, and surface runoff and WYLD increased when TDA increased, the magnitudes of these changes were not significant. The trends of changes in such water balance components from the present study agreed with findings from the previous studies in the 21.3-km² Goodwin Creek Watershed, in northern Mississippi, USA (Bingner et al. 1997), the four Iowa watersheds (2,000–18,000 km²), USA (Jha et al. 2004), and the 26.12-km² Walnut Creek watershed, in central Iowa, USA (Chaplot 2014). The slight increase in surface runoff and WYLD in the present study could suggest that rainfall amount had a slight increasing trend, and PET and ET had a slight decreasing trend when the number of sub-watersheds decreased. Therefore, additional water was more than water loss. However, the findings from the present study also contradicted the findings in the study at the 65,145-km² Tapi River, India, since they concluded that surface runoff decreased when TDA increased (Munoth & Goyal 2019a). They stated that reduction in the number of streams was the reason for decreasing runoff in the coarse schemes. However, loss of streams could affect sediment yield more than runoff volume (Bingner et al. 1997), since stream loss would result in loss of deposition process in streams (Jha et al. 2004).

According to the GIS maps of the spatial variation of ET (Figure 7(e)–7(h)), SURQ (Figure 7(i)–7(l)), and WYLD (Figure 7(m)–7(p)), the higher number of sub-watersheds produced finer variation of water balance components. This was in agreement with the study in the 384-km² Grote Nete River catchment, in Flanders, Belgium (Rouhani et al. 2009). In addition, hotspots (e.g., the places with extremely high/low values of water balance components in comparison to their surroundings) were presented more clearly in the finer schemes than those in the coarser ones.

Snowmelt is the main source of spring runoff in the Arctic. Determining the vulnerable locations due to high snowmelt runoff volume is highly important for risk management, e.g., flash flood, erosion, and landslide. In general, four different TDA schemes produced similar locations with high or low snowmelt volume across the watershed (Figure 8).

However, the maximum and minimum amounts of snowmelt runoff as well as the areas of sub-basins dominating by snowmelt are somewhat different among four TDAs. The reason could be that the fine TDA schemes have higher land-use composition, which could result in decreasing snow albedo and accelerate the snowmelt process (Szczypta et al. 2015). Therefore, the higher values of maximum snowmelt volume were found in the fine schemes. Moreover, the fitted model parameters were somewhat different among four TDA schemes during the calibration process (Supplementary Tables S1– S5). This could result in the differences in calculated snowmelt runoff across the watershed.

Furthermore, it is obvious from Figure 8 that the fine TDAs 200 ha (Figure 8(a)) and 2,000 ha (Figure 8(b)) generated finer and higher ranges of annual mean snowmelt volume across the watershed. Also, TDAs 200 and 2,000 ha could point out some hotspot locations of snowmelt volume within the watershed (marked with dark red colors). The relatively high annual mean value of snowmelt runoff volume and its large impacted areas were detected in the central section of the watershed. However, magnitude and impacted areas were inconsistent among four TDAs. For example, the annual mean snowmelt of 450–500 mm accounted for the largest area with 2,462 km², calculated from TDA 200 ha, while the annual mean snowmelt of 500–550 mm accounted for the largest impacted area with 3,338 km², achieved from TDA 2,000 ha. TDA 5,000 ha had the largest impacted area (2,896 km²) regarding the annual mean snowmelt of 400–450 mm. Similar to TDA 200 ha, TDA 10,000 ha detected the largest impacted area (1,927 km²) regarding the annual mean snowmelt of 450–500 mm. Based on the spatial distribution of simulated snowmelt runoff, it is recommended from the present study that more inspection should be focused on the central parts of the watershed as well as the locations of snowmelt hotspots for better risk management due to high snowmelt volume. Additionally, the central sections of the watershed and locations of snowmelt hotspots are the high mountain areas; therefore, the risks for flash flood or landslide could be high.

For each TDA scheme, the simulation result of monthly streamflow during a 10-year period, from 1998 to 2007, was averaged for each hydro-gauging station, and all the results were plotted as shown in Figure 9. According to Figure 9, all TDA schemes had capacities to simulate the observed tendency of monthly mean streamflow at all five hydro-gauging stations. However, the accuracy of simulation of monthly mean peak discharge was heterogeneous among TDA schemes as well as among five hydro-gauging stations. For example, the finest scheme TDA 200 ha performed quite well in simulating peak flow at Lundberg and Målselvfossen, while the coarsest scheme 10,000 ha was able to capture the peak flow at Høgskarhus and also at Målselvfossen. Lille Rostavatn was only the station where all TDA schemes yielded similarly the simulations of peak flows.

Therefore, the present study found that the accuracy of streamflow simulation did not totally depend on the levels of watershed discretization, but may also be controlled by other factors such as the geographic location and/or topographic characteristics of the sub-basins surrounding the investigated hydro-gauging stations. Unlike the homogeneous effects of watershed subdivisions on streamflow simulation in the previous studies, heterogeneous effects were found in the present study. For example, a previous study concluded that the accuracy of streamflow prediction was only increased (Mamillapalli et al. 1996) or decreased (Rouhani et al. 2009) when the number of sub-watersheds increased. Thus, it could be revealed that because of the complexity of hydrological cycles in the Arctic conditions as well as the topographic characteristics of the watershed, it could result in the heterogeneous effects of watershed subdivisions on streamflow simulations at different locations within the watershed. Furthermore, the heterogeneous effects of watershed subdivisions on the simulation results of streamflow hydrograph from the present study could help to explain the impacts of topographic variation compared to the homogeneous effects of watershed subdivisions in a flat watershed, e.g., the 152.29-km² Little Washita watershed, USA (Norris & Haan 1993). For example, that study found that the simulated peak flow linearly increased with the increase of number of sub-watersheds (Norris & Haan 1993).

Figure 10 illustrates the spatial variation of annual mean streamflows resulting from different TDA schemes. Under the finest TDA 200 ha, five levels of stream order were generated and displayed the high variation of the annual mean streamflow across the entire watershed (Figure 10(a)). Levels of stream order gradually declined in the coarser schemes TDA 2,000 ha (Figure 10(b)), TDA 5,000 ha (Figure 10(c)), and TDA 10,000 ha (Figure 10(d)) because of the decrease in number of sub-basins and the increase in sub-basins’ sizes. As a consequence, flow data in the coarser schemes were only detected in the major streams. In contrast, the finest TDA scheme was able to detect some minor streams and their flow values within the watershed; however, such minor streams may not be realistic and may not exist on ground.

Additionally, discrepancies in stream order levels from the fine schemes to the coarse schemes resulted in differences in the spatial distribution of streamflow values. For example, range of values and spatial distribution of streamflow data Q1 (for stream order level 1), Q2 (for stream order level 2), Q3 (for stream order level 3), Q4 (for stream order level 4, except TDA 10,000 ha), and Q5 (for stream order level 5, and only defining in TDA 200 ha) were inconsistent from the finest scheme to the coarsest scheme. This might affect the management of water resources by zones subdivisions. The coarser schemes simplified the stream networks and may lose some important in-stream processes. Loss of minor streams and their flow data may also affect local flood risk analysis (if this task would be planned to carry out in the future studies). For example, total stream length (e.g., minor streams) and drainage density generated from the coarsest TDA 10,000 ha reduced approximately 75% compared to those from the finest TDA 200 ha.

In the SWAT model, the nearest neighbor search (NNS) method is used to pick up the weather grid point representative for each sub-basin. According to this technique, one weather grid point closest to the centroid of a sub-basin is selected. The areal rainfall was then calculated for each sub-basin based on the rainfall data from the selected weather grid point. It was observed in this study that the number of weather grid points dropped from 21 grids in TDA 200 ha to 20 grids (−5%), 18 grids (−14.3%), and 14 grids (−33.33%) in TDAs 2,000, 5,000, and 10,000 ha, respectively (Figure 4). Fewer weather grid points produced more uniform distribution of areal rainfall across the watershed (Figure 7(b)–7(d)). As a result, the spatial variation of other water balance components was also influenced (Figure 7(e)–7(p)), since rainfall is the main input to generate water resources in the watershed.

Global sensitivity analysis in the SUFI-2 algorithm in the SWAT-CUP program detected the most sensitive model parameters among 18 model parameters used for calibration. It was found in the present study that the number of sensitive model paramters and their ranks changed significantly under different TDA schemes. Figure 11 provides the sensitivity changes of each calibrated parameter, which is classified into four different hydrological subgroup processes, while Figure 12 provides the sensitivity rank for all 18 calibrated parameters, based on the values of t-stat (magnitude of sensitivity) and p-value (significance of sensitivity), under each TDA scheme. Herein, the sensitivity rank, from highest to lowest, is the bottom-up order.

Regarding the surface runoff subgroup process (Figure 11(b)), runoff curve number CN2 showed its higher impacts (e.g., higher sensitivity rank) on the formation of runoff than the remaining parameter when the number of sub-watersheds changed (decreased). This is because CN2 depends strongly on land-use and soil characteristics of the watershed (Tegegne et al. 2019). Therefore, the discrepancies in land-use/soil composition by changing TDAs (e.g., the coarsest TDA 10,000 ha lost 4 land-use groups (−36.36%) compared to the finest TDA 200 ha) resulted in high fluctuation of CN2 values. The findings from the present study contradicted findings from the previous study in the Yongdam watershed in South Korea and in the Gilgelabay watershed in Ethiopia, since the sensitivity of CN2 in those studies increased when the number of sub-watersheds increased (Tegegne et al. 2019).

Parameters of soil properties, which governed the lateral flow process, fluctuated under changes of TDAs (Figure 11(c)). However, these soil parameters were not ranked in high sensitivity levels. Based on soil map, there are not many soil types in the Målselv watershed, and sandy loam dominates the major area. Therefore, the fluctuation of soil parameters resulted from changes in land-use composition rather than from changes in soil textures/types. The consequences of land-use change on soil properties such as soil moisture, infiltration capacity, and water storage were well validated in previous studies (Moges et al. 2013).

Several parameters in the snowmelt subgroup process had higher sensitivity rank orders compared to parameters of other subgroup processes (Figure 11(d)). The changes in watershed subdivision scales also resulted in the fluctuation of these parameters. The reason could be the changes in land-use structure in the watershed. For example, high/low density of land cover significantly affects the amount of snowfall accumulated on the ground, snow interception, sublimation process of snow, snow surface radiation balance (by the shading of trees), and snow albedo (Szczypta et al. 2015).

Changes in TDAs also resulted in the fluctuation of the sensitivity ranks of parameters in channel water routing (Figure 11(e)). However, the sensitivity ranks of the Manning's value (n) for the main channel (CH_N2) were more stable than those of the effective hydraulic conductivity in the main channel (CH_K2) in the channel subgroup process. The variation of parameters in the channel water routing process could result from the discrepancies in terrain complexity and the topographic attributes such as average overland sub-basins slope, streams length, average streams slope, and streams width (Table 3).

Most of the parameters in the groundwater subgroup process had low sensitivity rank and highly fluctuated (except GW_DELAY) when the number of sub-watersheds decreased (Figure 11(f)). The reasons for the fluctuation of parameters in the groundwater subgroup process could be the consequences of the changes in land-use structures and the associated soil properties, which affected the processes of infiltration, groundwater recharge, water transmission in soil, as well as evaporation from soil layers.

Watershed delineation is an important preliminary step for setting up the hydrological models. The application of the TDA technique for watershed subdivisions has a significant impact on hydrological simulation results. This study used the TDA technique to discretize the Målselv watershed, an Arctic watershed in the North of Norway. Four different TDA schemes, from the finest to the coarsest one, including 200, 2,000, 5,000, and 10,000 ha, were designed. The aim was to evaluate the impacts of different TDA schemes and spatial density of weather grid integration on hydrological simulations. The simulation results of water balance components, snowmelt runoff volume, and streamflow were evaluated. The main conclusions were drawn from the present study as follows:

• 1.

  The complexity of terrain and topographic attributes changed significantly with increasing TDA values such as decreasing total stream length, average stream slope, and drainage density. Stream order levels also declined as TDA increased.

• 2.

  The spatial density of weather grid point integration decreased from −5 to −33.33% in the coarse schemes 2,000, 5,000, and 10,000 ha compared to the finest scheme 200 ha. Fewer numbers of weather grid points produced more uniform distribution of areal rainfall across the watershed, particularly TDA 10,000 ha.

• 3.

  Watershed subdivision did not strongly affect the model performances during the calibration period at Lille Rostavatn and Målselvfossen hydro-gauging stations compared to the remaining stations in upstream.

• 4.

  Watershed subdivisions affected the spatial variation of areal rainfall across the watershed, but the total rainfall amount for the whole watershed was only slightly changed.

• 5.

  The changes in annual mean values of water balance components from the finest TDA 200 ha to the coarsest TDA 10,000 ha: rainfall (increased), PET and ET (decreased), surface runoff (increased), lateral flow (decreased), and WYLD (increased). However, the magnitudes of changes were not significant. The coarsest TDA 10,000 ha simplified land-use composition (e.g., loss of 4 land-use groups compared to TDA 200 ha with 11 groups) resulted in decreasing PET and ET.

• 6.

  The fine TDAs produced finer variation of snowmelt runoff volume and higher values of maximum snowmelt runoff volume across the watershed. Most TDAs detected similarly the most vulnerable areas by high snowmelt runoff volume, which are located in the central sections of the watershed and dominated by highly mountainous terrain. Therefore, high inspection should be focused on the central parts for better risk management, especially the risks of flash flood, erosion, and landslide.

• 7.

  All four TDA schemes had similar capacities to replicate the observed tendency of monthly mean streamflow hydrograph at all five hydro-gauging stations within the watershed. However, monthly mean peak flow had slight discrepancies among TDA schemes and among five hydro-gauging stations. Lille Rostavatn station was the only one where all four TDA schemes produced similar simulations of monthly mean peak flows.

• 8.

  The finest TDA scheme 200 ha could generate a high spatial variation of streamflow and could reveal the flow values of the minor streams, while other coarser schemes only showed streamflow in the major streams. Therefore, using the coarse schemes could influence the accuracy/detail of local flood risk analysis or finding flood hotspots in case such studies would be planned in future.

• 9.

  Watershed subdivisions strongly affected the calibration process regarding the changes in numbers and disordering the sensitivity ranks among the 18 calibrated model parameters.

• 10.

  One of the advantages of using the coarse TDA schemes was highlighted in the calibration process since the calibration time for approximately 5,913 km² watershed, and for simulation with monthly time step, decreased (approximately 5.7–34%) when the number of sub-watersheds decreased (from 459 to 18). This would benefit the users who have limitations in time and the available resources for running the model.

In brief, according to findings from the present study, it would be recommended that choosing the suitable threshold of TDA for watershed subdivision could depend on several factors including the purposes of studies, the expected level of accuracy in the output variables, the limited/allowed time, or the available resources that using the fine schemes or the coarse schemes would be an appropriate choice. Moreover, findings from the present study could help to enrich knowledge in designing the TDA for watershed delineation in hydrological models, particularly in the Arctic environment.

The authors acknowledge the Department of Technology and Safety, University of Tromsø – The Arctic University of Norway – for their financial support in this study. We recognize the in-kind contribution from Foundation CSDI WaterTech.

All relevant data are included in the paper or its Supplementary Information.