An Alternative Approach to Overcome the Limitation of HRUs in Analyzing Hydrological Processes Based on Land Use/Cover Change

Since the concept of hydrological response units (HRUs) is used widely in hydrological modeling, the land use change scenarios analysis based on HRU may have direct influence on hydrological processes due to its simplified flow routing and HRU spatial distribution. This paper intends to overcome this issue based on a new analysis approach to explain what impacts for the impact of land use/cover change on hydrological processes (LUCCIHP), and compare whether differences exist between the conventional approach and the improved approach. Therefore, we proposed a sub-basin segmentation approach to obtain more reasonable impact assessment of LUCC scenario by re-discretizing the HRUs and prolonging the flow path in which the LUCC occurs. As a scenario study, the SWAT model is used in the Aksu River Basin, China, to simulate the response of hydrological processes to LUCC over ten years. Moreover, the impacts of LUCC on hydrological processes before and after model modification are compared and analyzed at three levels (catchment scale, sub-basin scale and HRU scale). Comparative analysis of Nash–Sutcliffe coefficient (NSE), RSR and Pbias, model simulations before and after model improvement shows that NSE increased by up to 2%, RSR decreased from 0.73 to 0.72, and Pbias decreased from 0.13 to 0.05. The major LUCCs affecting hydrological elements in this basin are related to the degradation of grassland and snow/ice and expansion of farmland and bare land. Model simulations before and after model improvement show that the average variation of flow components in typical sub-basins (surface runoff, lateral flow and groundwater flow) are changed by +11.09%, −4.51%, and −6.58%, and +10.53%, −1.55%, and −8.98% from the base period model scenario, respectively. Moreover, the spatial response of surface runoff at the HRU level reveals clear spatial differences between before and after model improvement. This alternative approach illustrates the potential bias caused by the conventional configuration and offers the possible application.


Introduction
In recent years, the Earth's water cycle has changed due to the impact of global changes and human activities [1], and many countries or regions worldwide are facing serious water problems and crises [2,3]. Water problems have become a key factor restricting national and regional sustainable development [4]. Current water problems are mostly related to irrational human development activities, and various studies have shown that these impacts may outweigh the effects of other global changes [1,5]. Land use change is considered to be a direct reflection of human activities, which have dramatically altered hydrological processes and spatio-temporal redistribution of regional water resources [6,7]. Moreover, current hydrological and meteorological extremes appear more frequently than in the last century due to land use/cover change (LUCC) [8,9]. Land use/cover change (LUCC) has a significant effect on watershed hydrological processes by changing canopy interception, surface roughness, the infiltration rate, soil moisture and surface evaporation. Thus, a thorough understanding of the hydrological responses of LUCC is critical to the rational planning of land use patterns and the sustainable development of watersheds.
In the 1980s, without the support of advanced computer technology, remote sensing and geographic information system technology, it was very difficult to analyze the impact of LUCC on hydrological processes (LUCCIHP). Some scholars used comparative test methods for watersheds to study the effects of vegetation, especially the impacts of forest changes on hydrological elements [10,11]. These methods involved the selection of two small watersheds with different vegetation types but that are similar in other aspects, and their runoff characteristics are compared. While experimental results can be obtained using this type of approach, it is difficult to find two watersheds that are similar [12].
Since then, new studies have used many other conventional approaches for LUCCIHP such as the watershed water balance principle [13,14], statistical approaches [15,16], and the distributed hydrological modeling method. The first two approaches do not consider the multi-scale difference of LUCCIHP and contain fewer physical process mechanisms of the water cycle. The distributed hydrological modeling method is commonly used because of its physical basis, high precision and various options. Many hydrological models have been used successfully. Table 1 presents a summary of relevant representative literature for LUCCIHP based on the watershed water balance principle, statistical approaches and the distributed hydrological modeling method, such as the large-scale Variable Infiltration Capacity (VIC) hydrological model [17], the Soil and Water Assessment Tool (SWAT) [18,19], the Geo-spatial interface for Water Erosion Prediction Project (GeoWEPP) [20], the Regional Hydro-Ecological Simulation System (RHESSys) [21], and the MIKE SHE model [22]. The SWAT model is one of the most commonly used models because of its easy setup, moderate data requirement and many other practical modules that include most water-related topics [23]. The SWAT model divides the sub-basins into multiple HRUs, which have unique land use, soil and slope characteristics [18]. The conventional approach of considering LUCCIHP in the SWAT model is to directly replace the land use/cover map with another for different periods [7,24,25]. The SWAT model since version 2009 with LUP (land update) module can automatically update the land use/cover data of different periods, and keep the same number of sub-basins and HRUs [23,26,27]. The SWAT_LUP updated land use distribution by updating the HRU_FR (Proportion of HRU area in sub-basin) variable during the model run. However, there is a significant drawback in the models based on the concept of hydrological response unit (HRU) that the spatial distribution of LUCC is considered in watershed discretization but not in runoff routing. In addition, the threshold specified in the HRU definition has a small effect on LUCCIHP and does not consider different spatial distributions between land use types [28]. Therefore, a severe over-or underestimation of LUCCIHP might occur when applying this type of model structure.
This paper has three main aims: (a) to reveal the variation of land use/cover type transfer-in and transfer-out; (b) to dissect the importance of coupling the spatial information of LUCC with watershed discretization; and (c) to explain the impacts of LUCC on hydrological processes, and compare whether differences exist between the conventional approach and the improved approach for LUCCIHP.
For Objective (a), defining a suitable land use/cover classification system is the basis for LUCC analysis. There are also differences in the extent of land surface coverage information reflected by different land use/cover classification systems [29]. At present, five units have been released by the mainstream global land use/cover classification system: FAO (Food and Agriculture Organization)/ UNEP (United Nations Environment Programme) [30], USGS (United States Geological Survey) [31], CORINE (Co-ORdinated INformation on the Environment) [32], IGBP (International Geosphere-Biosphere Programme) [29] and IGSNRR (Institute of Geographic Sciences and Natural Resources Research) [33]. The IGBP classification system includes 17 types, of which 11 are natural vegetation types, three are land use types and three are non-vegetation growth land types [29]; this system can also consider synthetic land cover types and land use types. In response to the situation in the study area, we further refined the vegetation types at different coverage levels to specify the variation of hydrological processes. Similar to our research, many scholars also defined the land use/cover classification system according to the actual situation of the study area [7,13,[16][17][18]. Too few classes may generalize the detailed information of land use/cover changes and further reduce its effect on hydrological processes. In addition, we analyzed the land use/cover change from two aspects of time and space, which can reflect changes in the amount of time and can reflect the distribution of the changes in space. Some scholars only analyzed the main types of land use/cover change at the catchment scale and qualitatively analyzed LUCCIHP based on the catchment scale [13,16]. Other scholars assumed some land use/cover change scenarios to analyze the impact mechanism of land use/cover change [18,22,25]. Others temporal and spatial analyses of land use/cover changes were performed before carrying out LUCCIHP [7,17]. However, most only analyzed the time variation at the river catchment scale, including the conversion area, transformation ratio and degree of change of each land use/cover type during the study period. Very few researchers have focused on the temporal and spatial changes in land use/cover at the basin scale, and there is a lack of quantitative analysis at the sub-basin scale [18]. To analyze LUCCIHP accurately, this study not only analyzes the temporal and spatial changes of land use/cover at the catchment scale, but also places the changed regions in the corresponding sub-basins.
For Objective (b), as everyone knows, the hydrological modeling method has a more solid physical basis than traditional methods such as the watershed water balance principle and statistical methods [34]. These conventional methods can only analyze the impact of land use/cover change on runoff at the catchment scale [13]. In hydrological models, such as VIC, MIKE SHE, etc., the smallest unit of computation is a grid, which has better accuracy in hydrological process spatial simulation.
However, since grid computing increases computational complexity and runs for a long time, it is suitable for smaller watersheds [35]. For large watersheds, hydrological models such as SWAT, GeoWEPP, and RHESSys that use HRU as the minimum calculation unit have greater advantages, and these models can well simulate the hydrological processes and reduce the computational load and save model running time [20,21,25]. The SWAT model has more land use/cover related parameters and options, so it has more capacities for analyzing LUCCIHP.
There are three scale units in the SWAT model: the watershed, sub-basin and HRU. A study has shown that smaller catchment areas and more sub-basins have higher simulation accuracy [36]. However, another study showed that, if the catchment area is relatively small, the simulation accuracy will not increase, but will stabilize [37]. However, more sub-basins cannot replace HRUs. In a sub-basin, land use/cover data, soil data and slope data are superimposed to form a plurality of patches. There is only one type of land use/cover, soil type and slope type in each patch, i.e., the HRU. Therefore, as the smallest unit in the SWAT model, the HRU has a unique land use/cover type, and many related parameters.
In addition to meeting the water balance, grasping the flood peak and guaranteeing base flow, corresponding parameters should be given different values according to different land use/cover types due to their impact on hydrological processes. The calculation of SWAT surface runoff is based on the SCS curve method, in which a curve number (CN) is introduced, which is a comprehensive parameter reflecting the characteristics of the basin before rainfall. Land use/cover types correspond to different CN values, which also have differences in runoff generated under different rainfall conditions [25]. In addition, different types of land use/cover corresponding to the amount of evapotranspiration are also quite different, e.g., the plant uptake compensation factor (EPCO) parameter values are specified for different land use/cover types.
The following steps were performed in this study: (i) the temporal and spatial variation characteristics of watershed land use/cover for the Aksu River Basin, China, were analyzed; (ii) DEM, land use, soil, long term (2000-2007) streamflow and climate data were used to build the SWAT model, and perform calibration and validation; (iii) the spatial distribution of LUCC was analyzed to identify sub-basins with significant changes, and the sub-basins were segmented; (iv) the runoff variation of the watershed was comparatively analyzed at the catchment scale; (v) the runoff components were analyzed at the sub-basin scale; and (vi) the spatial difference in surface runoff was analyzed at the HRU scale based on the conventional approach and the improved approach. Following these steps, an alternative approach was designed to compare whether the differences exists between the conventional approach and the improved approach for LUCCIHP using the modified sub-basin segmentation to reveal the importance of considering the spatial distribution of LUCC before watershed discretization.

Study Area
The Aksu River Basin is located west of the Xinjiang Uygur Autonomous Region, east of the Republic of Kyrgyzstan, and southeast of the Republic of Kazakhstan, enclosed between latitudes 40 • 16 -42 • 28 N and longitudes 75 • 4 -80 • 18 E with an area of 4.29 × 10 4 km 2 ( Figure 1). The Aksu River is the largest water system on the southern slope of the Tianshan Mountains, with two main headwater tributaries, the Kumalak River and Toxkan River. The Toxkan River originates from the Artbash Mountains, with a length of 457 km and a catchment area of 2.91× 10 4 km 2 . The Kumalak River is 293 km long and originates from the Hantengri Mountains, with a catchment area of 1.38 × 10 4 km 2 . Each tributary has a hydrological gauging station: Shaliguilank (SLGLK) on the Kumalak River and Xiehela (XHL) on the Toxkan River. The rivers are mainly fed by snow and glacier melt water from the mountainous regions. The Aksu River Basin has a complex regional climate because of the rugged terrain with elevation between 1118 and 7126 m. The high mountainous areas (Elevation (E) > 3500 m) have low temperature, moderate precipitation, and glaciers and seasonal snow. In the middle elevation area (2000 m < E < 3500 m), there is a relatively clear boundary between warm and cold, and the highest precipitation occurs here. The lower mountain areas (E < 2000 m) have low precipitation and a large day-night temperature difference. Since the Aksu River Basin is embedded in the Tianshan Mountains and is next to the Taklimakan Desert, the arid environment formed under a continental climate. The downstream region is characterized by low precipitation, high potential evaporation, large day-night temperature difference, and a high number of sunshine hours.
This river is the most important tributary of the Tarim River [38] and plays a vital role in food security and the ecological health of the entire Tarim River Basin. Recent research examined the temporal variation of climate variables [39], LUCC [40], and runoff variation [41,42], indicating that LUCC has an impact on runoff variation in the Aksu River Basin.

Materials
Digital elevation map (DEM), soil map, land use/cover map and basic climate datasets are required to build the SWAT model. DEM data at 90 m resolution were obtained via the Shuttle Radar Topography Mission (SRTM, http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp) which can provide watershed slope and elevation information. The soil map was derived from the Harmonized World Soil Database version 1.1 (HWSD, http://westdc.westgis.ac.cn/data/611f7d50-b419-4d14-b4dd-4a944b141175) at 1 km spatial resolution from the Food and Agriculture Organization (FAO) and International Institute for Applied Systems Analysis (IIASA). The land use/cover maps for 2000 and 2010 were generated based on Landsat MSS/TM products, which were produced by the Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences.
The basic climate data (daily precipitation, daily maximum and minimum temperatures) were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn/). We used daily discharge time series data which were measured in m 3 /s from 2000 to 2007 at the XHL and SLGLK hydrological gauging stations. The data were acquired from the Xinjiang Tarim River Basin Management Bureau for the purpose of SWAT model calibration and validation.

Land Use/Cover Change
Land use/cover change (LUCC) includes changes in both land use type and land cover type over time, as well as in space. LUCC analysis can provide the basic information for LUCCIHP. We use ARCGIS 10.2 software to analyze the spatial and temporal characteristics of LUCC by using functions of Clip, Identify, and Zonal statistics as table. Then, the transition matrix of land use/cover change in topical sub-basins were created [43] including the following standards and processes: (i) The two periods of land use/cover maps were standardized based on the unified projection and geographic coordinate system, which were foundation of LUCC analysis and model building. (ii) According to the land use/cover classification system of the International Geosphere-Biosphere Programme (IGBP) [44], we classified the two periods of land use/cover maps into thirteen categories: evergreen coniferous forest, deciduous broadleaf forest, shrubland, meadows, grassland, sparse grass, water, permanent snow and ice, upland, paddy field, bare land, built-up and wetland. This step reflects the conversion of various land types, including conversions within the same vegetation type, such as the conversion of sparse grasslands and grasslands. (iii) The number and spatial changes of watershed land use/cover change, including location, area and type in sub-basins, were analyzed.
A LUCC map for the period between 2000 and 2010 was produced ( Figure 2) to reveal the characteristics of change, such as the extensive expanding in upland and bare land and the shrink of sparse grass. The images suggest that LUCC had strongly disturbed the ecosystem and natural balance in the local area. The extent of the LUCC distribution is presented in Figure 2b  To focus on the hydrological effects of grassland degradation, the lands of meadows, grassland and sparse grassland are collectively referred to as grasslands. The spatial locations of the main LUCC areas between 2000 and 2010 are shown in Figure 2a. The LUCC type with the greatest change was the conversion of grassland to bare land, accounting for 26.28% of the total transfer area. Subsequently, the conversion of grassland to upland counts for 7.79% of LUCC.
Our results indicated that the largest increases in land use/cover from 2000 to 2010 were in bare land and dryland, which originated mainly from snow/ice cover retreat, grassland degeneration, and reclamation of grassland or shrubland. Since the 1990s, global warming has also affected the snow/ice cover area through mechanisms such as decreasing the proportion of solid precipitation, enhancing snowmelt [45,46], and increasing bare land. In addition, cultivated land in oases has expanded rapidly since 1950, and made this basin one of the most important production areas for grains, cotton, and fruits in Xinjiang [47]. The newly added cultivated land was mainly converted from grassland [40]. Therefore, permanent snow, ice and grassland have become the main sources of transfer-out to other land types.

Base Period Model Building and Calibration
The SWAT model [48] is a semi-distributed hydrological model that acts as a continuous simulator of water, sediment, pesticide, and nutrient transport at multiple scales. The SWAT model is based on physical algorithms to simulate surface and subsurface processes, such as precipitation interception, snow fall and snow melt, plant uptake, evaporation, infiltration, percolation into aquifers, lateral flows, ground flow, and surface runoff processes. In the SWAT model, a watershed can be segmented into several sub-basins regarding the digital elevation and threshold value of the catchment area. Soil maps, land use/cover maps and slope maps are overlaid to generate the HRUs, which are the smallest calculated units in the sub-basin model with unique properties of land use, soil and slope [18]. Therefore, HRUs are lumped together and there is no interaction among HRUs in one sub-basin. Only at the sub-basin level are the spatial relationships among sub-basins specified [48]. Therefore, if the interaction of one land use area with another is important, rather than defining the relevant land use areas as HRUs, they should be defined as sub-basins. To reduce the noise information of land use/cover change, the HRU area ratio was set to 1%, which means that land use/cover patches up to 1% areas cannot be ignored.
We use the DEM data, soil data and 2000 land use/cover data to build the basic model (for the traditional method analysis): the watershed was divided into 35 sub-watersheds within the catchment area of 5000 hectares. The largest sub-basin area was 3149.81 km 2 , and the smallest sub-basin area was 0.10 km 2 . If the HRU division is too large, many land use/cover patches that cover less area than the specified proportion will be ignored, that is, only HRUs that meet the criteria are retained in the sub-basin. A study has shown that, when the proportion of the HRU area is small based on a certain threshold, there will be no significant distinction of land use/cover information [28]. To reduce the loss of land use/cover change information in the model, the land use threshold was set to 1%, the soil threshold to 5% and the slope threshold to 5% in the HRU definition step, and the watershed was divided into 986 HRUs. Therefore, it is suggested that the catchment area and HRU area be set proportionally according to the actual situation of the basin size and land use/cover change. After entering the maximum temperature, daily minimum temperature and daily precipitation data from 2000 to 2007, we used the same period of daily discharge data at the SLGLK and XHL hydrological stations to calibrate and validate the Aksu River Basin SWAT model. Among them, the first two years  Figure S1).
To ensure the reliability of model simulation and obtain a set of reasonable parameters, we adopted the automated calibration method of sequential uncertainty fitting (SUFI-2) [49] based on SWAT-CUP software to calibrate daily runoff data. The Latin hypercube sampling method is integrated into the SUFI-2 procedure, which can generate the optimum solution from a multidimensional distribution and ensure that it is not a local minimum [50]. This is referred to as the 95% prediction uncertainty, or 95PPU. These 95PPUs are the model outputs in a stochastic calibration approach. Thus, an envelope of good solutions is expressed by the 95PPU and a set of certain parameter ranges. The p-value is an indicator to assess the parameter sensitivity by comparing the mean value between simulated values and measurements [51]. This gives relative sensitivities based on linear approximations and, hence, only provides partial information about the sensitivity of the objective function to model parameters. For further details, refer to the SWAT-CUP manual (http://swat.tamu.edu/media/ 114860/usermanual_swatcup.pdf). The smaller the average p value is, the more sensitive the parameter is. The snowfall temperature (SFTMP, • C), snow melt base temperature (SMTMP, • C), melt factor for snow on 21 June (SMFMX, mm H 2 O/ • C-day), melt factor for snow on 21 December (SMFMN, mm H 2 O/ • C-day), precipitation lapse rate (PLAPS, mm H 2 O/km), temperature lapse rate (TPLAPS, • C/km), and baseflow alpha factor (ALPHA_BF, 1/days) are more sensitive than other parameters in the Aksu River Basin SWAT model [38] (Table 2). In addition, we assigned different values to the parameter curve number (CN) for each type of land use/cover ( Table 2). The Nash-Sutcliffe coefficient (NSE) [52], ratio of the root mean square error to the standard deviation of measured data (RSR) and percent bias (PBias) were chosen to evaluate the "goodness-of-fit" of modeled results [53,54]. The formulas of these indices are presented in Equations (1)- (3): where Q i,obs and Q i,sim are the average value of the observed and simulated values, respectively; Q i,obs and Q i,sim are the observed and simulated values at time step i, respectively; and n is the number of time steps. Normally, when NSE > 0.5%, RSR < 0.7% and PBias < 25%, the results of the modeling discharge are supposed to be satisfactory [54].

Modified Approach for LUCCIHP
As mentioned in Sections 1 and 2.2, conventional analysis of LUCCIHP has a potential bias if the changed land use pattern has a jigsaw and interlocking spatial pattern. In such cases, the original parallel configuration of HRUs might not represent reality and may cause bias [55]. Because the conventional SWAT analysis method lacks pertinent considerations in model building and parameter calibration, a simplified analysis of LUCCIHP cannot effectively grasp the true hydrological effects of LUCC. In particular, with respect to parameter configuration, if all land use/cover types are given the same parameter values, the results of the simulated hydrological response will not be obvious after a change in land use/cover type [25]. When land use/cover change occurs in a sub-basin, the type of some HRUs may change. Since the hydrological processes between HRUs are calculated in parallel, the flow elements of all HRUs are aggregated to the sub-basin outlet [55]. In the conventional approach, the relative spatial location of changed HRUs is not considered in response to land use/cover change types in the sub-basin.
A modified sub-basin segmentation approach was tested to assess this issue, quantify potential bias, and provide a possible approach. We compared the impact results of LUCC using the conventional approach and a modified approach. The land use/cover maps for 2000 and 2010 were used for assessment.
Building and calibrating models as described in Section 3.2.1 will help to improve the accuracy of LUCCIHP analysis. The improved method first provides the model base for analyzing LUCCIHP, and can clearly reflect the response of hydrological processes when land use/cover change occurs. By analyzing the spatial location of LUCC, we can identify sub-basins with significant changes that refer to the percentage of the total area that underwent land use/cover change. Sub-basins with large land use/cover change areas can be considered for segmentation. Segmenting a sub-basin is similar to discretizing a watershed into more sub-basins, but there are essential differences when considering the key nodes where the land use/cover changed dramatically. Setting the catchment area threshold to a very small size can result in many small sub-basins. We used a similar idea to integrate the spatial location of land use/cover change and to determine the location of nodes in the sub-basin. According to this new idea, the HRU where land use/cover change occurred was segmented into different sub-basins to fully consider LUCCIHP. A detailed description of the method is provided below.

a. The conventional analysis approach
The conventional analysis approach (Figure 3a) comprises three basic steps: (1) Set up the SWAT model with general routing with calibration and validation processes based on the land use/cover map before LUCC, that is build the SWAT model using DEM, soil, and the base period of land use/cover data and dividing HRUs, and then utilize existing climate data and runoff data to calibrate and validate model. (2) Backup the base period model and replace the new periods of land use/cover map after LUCC in the original model, and then carry out the standard procedure as previous step with previous parameter and variable sets of base period model. (3) Analyze the hydrological processes by the results from different land use/cover scenarios. When land use/cover types in watershed changes, the number of sub-basins will not change. However, due to the proportion change of the original land use/cover, the HRUs in sub-basin will change accordingly in shape, size, location and number [26]. Some hydrological elements in HRU will be calculated independently. Therefore, the LUCC in sub-basin altered to some extent the total water of each of hydrological elements.

b. Improved sub-basin segmentation approach
According to the SWAT model structure, the hydrological process elements in each HRU in certain sub-basins (i.e., surface runoff, lateral flow, ground flow, and evapotranspiration (ET)) are computed in parallel at the HRU level and simultaneously routed to the outlet of this sub-basin [18]. When replacing the new land use/cover map, the previous HRUs in sub-basins were reconstructed and redistributed, including in shape, size, number and location, because of the change of proportion of every land use/cover type, but new HRUs are not routed with consideration of the spatial position relationship among HRUs. Therefore, the improved sub-basin segmentation approach adds one step before the watershed delineator. The whole analysis process of LUCCIHP is described in Figure 3b. The improved analysis approach adds one more key step rather than original analysis approach. In this step, the spatial position relation of the land use unit is pre-processed to locate the main areas of LUCC. That is, we firstly analyzed the land use/cover change based on ARCGIS software in shape, position and area; then selected some larger change proportion of land use/cover patches and identified which sub-basins it belongs to; and then added notes between the changed land use/cover patches and unchanged patches under watershed delineation. Based on the key zones of change and the existing river network, more nodes (outlets of sub-basin) were added and specified to increase the complexity of stream routing (Figure 4). Then, we again divided HRUs according DEM, soil, and the new period of land use/cover data and run the model with previous parameter and input variable sets of base period model. Finally, we analyzed the hydrological processes under different land use/cover scenarios. The principle of sub-basin segmentation is to use the relative spatial position of LUCC to divide the sub-basin into further multi-sub-basins based on the relationship between upstream and downstream. For example, although the land use/cover has changed (grassland HRU changed to farmland HRU), the hydrological elements of all HRUs (e.g., surface runoff, lateral flow, and ground flow) will directly collect at the outlet of sub-basin 1 without considering the spatial position relationship. After LUCC occurred, the outlet of sub-basin 1 does not change, but the hydrological process elements of the forest HRU and farmland HRU preferentially move toward the outlet of sub-basin 1-1, then pass sub-basin 1-2 to the original outlet of sub-basin 1 ("Improved" in Figure 4). Therefore, this approach can reallocate the hydrological elements in a sub-basin.

Impact of Land Use/Cover Change on Hydrological Processes
The impact of LUCC, such as cultivated land or bare land expansion, on the hydrological response has already been studied in previous research [56]. The impact of land use/cover change on hydrological processes is mainly reflected in the change in water cycle elements and water quality [7,24]. Land use/cover changes affect the surface and near-surface hydrological cycle elements, including runoff, evaporation, interception, soil water, etc. [57]. Furthermore, land use/cover change influences on hydrological processes at multiple scales, including the catchment scale, regional scale, sub-basin scale, patch scale, HRU scale, grid scale, etc. According to Objective (c), combined with the model characteristics, we chose to analyze LUCCIHP at the catchment scale, sub-basin scale and HRU scale, and to analyze the difference between the conventional approach and the improved approach.
Based on the catchment scale, we analyzed the runoff variation of the outlet of the whole watershed, which can reflect the comprehensive effect of LUCC. In addition, we analyzed the runoff composition change resulting from surface runoff, lateral flow and groundwater flow at the sub-basin scale, which can reflect how runoff composition is affected by typical land use/cover change type. The first two scale analyses could reflect the comprehensive effect of LUCC on the amount of water. Based on the HRU scale, LUCCIHP can be reflected from the spatial perspective. The spatial distribution response of surface runoff was analyzed, and the response of evapotranspiration and recharge entering aquifers during time step (total amount of water entering shallow and deep aquifers during time step, mm H 2 O) distribution to LUCC can be found in the Supplementary Materials ( Figures S2 and S3). The multi-scale and multi-factor analyses can prove that the sub-watershed segmentation approach has a significant effect on LUCCIHP from different scale angles.

SWAT Model Building and Improvement
Using the spatial overlay analysis method to overlay the land use/cover change map and watershed discretization map, we analyzed the intensity of land use/cover change for each sub-basin and identified some sub-basins with the most significant land use/cover change. The land use/cover area ratio of 4, 17, 19, 20, and 22 sub-basins exhibited the greatest total change areas. Based on this result, before the watershed discretization, we divided the original 35 sub-basins into 40 sub-basins with five sub-basin segmentation nodes and obtained 1015 HRUs. The parameterization scheme of the base period model (LU2010) was revised to that of the new period model (LU2010-im) to ensure the two approaches were comparable. Through comparative analysis of NSE, RSR and Pbias, model simulations before and after model improvement showed that NSE was increased by up to 2% (from 0.46 to 0.47), while RSR decreased from 0.73 to 0.72, and Pbias decreased from 0.13 to 0.05.
The largest land use/cover change area in SUB19 was grassland to arable land. After segmentation of SUB19, the change area was located upstream, and the downstream area did not change. In contrast, the largest land use/cover change area in SUB4 was bare land, distributed in the downstream region. After segmentation of SUB4, the upstream sub-basin did not change, whereas the downstream sub-basin did. SUB22 was similar to SUB19, and SUB17 and SUB20 were similar to SUB4 ( Figure 5 and Tables S1-S5).   (Figure 6a). Based on the comparison of simulated annual average daily streamflow results between the base period model (LU2000, which is driven by land use/cover data for 2000) and the changed period models (LU2010, which is driven by land use/cover data for 2010 (not improved); and LU2010-im, which is driven by land use/cover data for 2010 (improved)), the latter (with less grassland cover and more bare land) had higher streamflow in spring. The impact of LUCC was considerable in the hydrographs ( Figure 6). The LU2010-im with the same change trend (reduced grassland and expanded bare land) showed a different response: the peak flow (400 m 3 /s) in spring based on LU2010-im was higher than that of the base period model (350 m 3 /s) and lower than that based on LU2010 (480 m 3 /s), whereas the peak flow in summer based on LU2010-im was similar to that of LU2010 and much lower than that of the base period model.

Variation of Flow Components at the Sub-Basin Level
To analyze the effects of LUCC on streamflow on the sub-basin level, we considered the variation of surface runoff, lateral flow, and groundwater flow in detail. We selected sub-basins with clear LUCC (sub-basins 4 and 22 in the upstream; sub-basin 17 in the mid-stream; and sub-basins 19 and 20 in the downstream) as typical sub-basins to analyze LUCCIHP. The proportion of the flow components of the LU2000 model can be summarized as follows: groundwater flow dominated in the mountainous areas; lateral flow was the major flow component in the downstream; and the three flow components were equally important in the mid-stream (Figure 7). The proportion of flow components was distinct due to the various combination of the key factors, such as land use and land cover, slope and aspect, soil properties, and rainfall intensity.
Compared with the LU2000 model, the LU2010 and LU2010-im models generated greater surface runoff, less lateral flow and less groundwater flow (Figure 7). Larger changes in sub-basin 22 in the downstream region of the Aksu River were observed, with +19.73% surface runoff, +1.32% lateral flow and −21.04% groundwater flow. These phenomena were clearly observed in the mid-stream sub-basin, e.g., sub-basin 17, but there was only a small change in the upstream mountainous areas (e.g., sub-basin 4). The opposite change trend was observed in sub-basin 19, with −0.62% surface runoff, +0.01% lateral flow and +0.61% groundwater flow, but this was not severe. Despite generally similar change trends between LU2010 and LU2010-im, the extent of change varies greatly. The increased surface runoff based on LU2010-im was less than that based on LU2010 (e.g., sub-basins 17 and 20). However, a greater level of increase occurred in other sub-basins of LU2010-im (e.g., sub-basins 4, 19, and 22). Consequently, compared with LU2010, there was an increasing trend of lateral flow in most sub-basins of LU2010-im except sub-basin 4. There was no substantial difference in groundwater flow.

Spatial Distribution of Surface Runoff at the HRU Level
The surface runoff at the sub-basin level was strongly influenced by LUCC, as shown in the previous section. Since land use/cover is the most direct representation of the underlying surface, its changes immediately influence the surface runoff processes and further affect streamflow components. The spatial distribution of surface runoff at the HRU level (mean value of 2002-2007) for the LU2000, LU2010, and LU2010-im models in typical sub-basins is illustrated in Figure 8. Spatial distribution of surface runoff altered markedly between LU2000 and LU2010, and between LU2010 and LU2010-im, with strong geological effects. Surface runoff in high-elevation sub-basins increased in some HRUs and decreased in other HRUs (e.g., sub-basins 4 and 22). In contrast, surface runoff in low-elevation sub-basins was considerably lower. The spatial distribution of surface runoff at the HRU level for the LU2010 model acted differently to that of the LU2010-im model after sub-basin segmentation. The surface runoff volume for the LU2010-im model was higher than that of LU2010 in some sub-basins (e.g., sub-basins 4, 17, and 22), which are mainly situated in high-and mid-elevation areas. In contrast, the surface runoff at the HRU level for the LU2010-im model in the lower sub-basins (e.g., sub-basin 19) showed overall decreasing trends. Furthermore, surface runoff at the HRU level for LU2010-im was between that of LU2000 and LU2010, indicating that the LU2010 model without the sub-basin segmentation approach overstated the LUCCIHP.
Spring snowmelt runoff in the Aksu River Basin is the major trigger of spring floods. Therefore, analysis of the spatial distribution change of surface runoff at the HRU level in April is favorable for understanding the effects of LUCC on surface runoff and the water resources in spring. The increasing trend of surface runoff volume in April for LU2010 and LU2010-im was clear in both the conventional approach and the improved approach (Figure 9). The surface runoff volume based on the LU2010 model was greater than that of the LU2000 model in high-and mid-altitude areas (e.g., sub-basins 4, 17, and 22). However, since the surface runoff was not the dominant hydrological component in lower sub-basins, the surface runoff in lower sub-basins (e.g., sub-basins 19 and 20) showed little change in April. Although the LU2010-im model had a similar change trend to LU2010, the distribution of surface runoff showed distinct differences. Following application of the sub-basin segmentation approach, the spatial variation of surface runoff between LU2000 and LU2010-im was less than that between LU2000 and LU2010. Furthermore, the surface runoff distributed in the upper region in the LU2010-im model was higher than that based on the LU2010 model. In lower regions, surface runoff based on LU2010-im was less than that based on LU2010. In mid-altitude areas, the surface runoff based on the LU2010-im model was lower than that based on the LU2010 model.

Discussion
The results show the improved model slightly improves the simulation results (NSE < 0.5, RSR > 0.7). This may be because land use/cover change in this basin is not very dramatic. According to our study, the change ratio of land use/cover in this basin is less than 2%. The regions with severe changes are mainly distributed in some sub-basins [40]. Moreover, this study does not focus on improving the accuracy of model simulation, but improving the accuracy of LUCCIHP scenario analysis. Therefore, the values of NSE and RSR are not satisfactory.
Our results also show a relationship between grassland cover degradation and streamflow increase, which may reduce the surface roughness as well as strengthen the peak flow. In spring in particular, the overflow from snowmelt is the main source of streamflow in this basin [41]. Vegetation degradation promotes peak flow in spring.
Our model results reveal that the streamflow at the basin outlet was increased by the reduction in grassland area and the increase in bare land and farmland. A similar hydrological response has been verified in previous studies. Haregeweyn et al. (2015) showed that decreases in vegetated land and increases in cultivated land in an Ethiopian highland catchment led to an increase in the annual runoff yield by 101 mm [58].
There is no doubt that LUCC will alter hydrological processes; however, different simulation approaches may produce different analysis results. In our study, the impact degree of LUCC on hydrological elements varied as a result of the consideration of the spatial position of LUCC between 2000 and 2010. Although the streamflow increased at the basin outlet in spring under LUCC, our results show that the peak flow based on the LU2010-im model was slightly lower than that based on the LU2010 model. This may be because conversion of grassland to bare land and grassland degeneration in some sub-basins of the mid-stream increase surface runoff, but grassland and farmland remain the dominant land use/cover types in these sub-basins (e.g., from 43.11%, 43.47% and 13.42% to 56.57%, 37.89% and 5.53% in surface runoff, lateral flow and groundwater at SUB17, respectively). Insertion of the LUCC nodes before watershed discretization allows streamflow to pass through the sub-basin to the outlet. The different routing scheme of channel flow and surface runoff might alter the arrival time and amount of peak flow.
Our model results also reveal remarkable spatial changes in surface runoff at the HRU level. With conversion of grassland to farmland or bare land, there were clear increasing trends of surface runoff on the spatial distribution in most sub-basins. This phenomenon was more apparent in April due to the water holding capacity of plants and roots of grassland in spring; however, this ability is lost after conversion to bare land or farmland. In contrast, after sub-basin segmentation, the sub-basins with grassland type HRUs in the LU2000 model were converted to sub-basins with bare land or cultivated land in the LU2010-im model. These new sub-basins can clearly indicate the spatial relationship between upstream and downstream areas, and overcome the drawback of the parallel relationship of HRUs in the original sub-basin.
Some related studies have indirectly identified the above phenomenon using other methods [18,25,59]. With respect to setting a different catchment size threshold value, Romanowicz et al. (2005) divided the watershed into different numbers of sub-basin scenarios and determined that the greater the number of sub-basins, the higher the NSE value is. However, optimal modeling results were not obtained when the catchment size threshold value (CSTV) was minimized. Similar to this study, some researchers analyzed the effect of LUCC on surface runoff at the catchment scale but did not perform sub-basin segmentation [7,18,25]. The results of Baker and Miller indicated that, when the proportion of land cover varied, the surface runoff did not change correspondingly [25], which may be partly because the change of land cover in the sub-basin area was not properly described. However, the impact of LUCC on the hydrological processes may have been better presented in this study after sub-basin segmentation.

Conclusions
The improved sub-basin segmentation approach has potentially advanced the simulation accuracy of how hydrological processes respond to LUCC in large-scale watersheds. In the current study, we first divide the sub-basin into further multi-sub-basins by analyzing the relative spatial position of LUCC and then improve the influence degree of LUCC on hydrological processes (NSE increased from 0.46 to 0.47, RSR decreased from 0.73 to 0.72, and Pbias decreased from 0.13 to 0.05) at three different scales: watershed, sub-basin, and HRU. Our results indicate that grassland degeneration and expansion of bare land/farmland resulted in higher surface runoff (from 43.10% to 56.57% in SUB17) and peak flows in spring (from 350 m 3 /s to 480 m 3 /s), and less lateral flow (from 43.47% to 37.89% in SUB17) and groundwater (from 13.42% to 5.53% in SUB17) than in the base period model, and a weaker change trend in the LU2010-im model (400 m 3 /s, 53.47%, 42.20% and 4.33%). At the HRU level, although an increasing trend of surface runoff was observed from 2000 to 2010 in the spatial distribution in most sub-basins, the magnitude of change differed between the LU2010 model and the LU2010-im model. The improved sub-basin segmentation approach is easy to perform and effectively overcomes the limitation of HRUs to analyze LUCCIHP by purposefully adding some outlets based on LUCC to improve the simulation accuracy. This method can serve as a reasonable approach for watershed spatial management, sustainable land use planning and ecological environmental protection.