Assessing Variation in Water Balance Components in Mountainous Inland River Basin Experiencing Climate Change

Quantification of the changes of water balance components is significant for water resource assessment and management. This paper employed the Soil and Water Assessment Tool (SWAT) model to estimate the water balance in a mountainous watershed in northwest China at different spatial scales over the past half century. The results showed that both Nash-Sutcliffe efficiency (NSE) and determination coefficient (R2) were over 0.90 for the calibration and validation periods. The water balance components presented rising trends at the watershed scale, and the total runoff increased by 30.5% during 1964 to 2013 period. Rising surface runoff and rising groundwater flow contributed 42.7% and 57.3% of the total rising runoff, respectively. The runoff coefficient was sensitive to increasing precipitation and was not significant to the increase of temperature. The alpine meadow was the main landscape which occupied 51.1% of the watershed and contributed 55.5% of the total runoff. Grass land, forest land, bare land, and glacier covered 14.2%, 18.8%, 15.4%, and 0.5% of the watershed and contributed 8.5%, 16.9%, 15.9%, and 3.2% of the total runoff, respectively. The elevation zone from 3500 to 4500 m occupied 66.5% of the watershed area, and contributed the majority of the total runoff (70.7%). The runoff coefficients in the elevation zone from 1637 to 2800 m, 2800 to 3500 m, 3500 to 4000 m, 4000 to 4500 m, and 4500 to 5062 m were 0.20, 0.27, 0.32, 0.43, and 0.78, respectively, which tend to be larger along with the elevation increase. The quantities and change trends of the water balance components at the watershed scale were calculated by the results of the sub-watersheds. Furthermore, we characterized the spatial distribution of quantities and changes in trends of water balance components at the sub-watershed scale analysis. This study provides some references for water resource management and planning in inland river basins.


Introduction
Water resources are affected by the alteration of a number of water balance components which have different responses to climate change [1,2].Effective management of water resources in cold and alpine regions requires detailed attention to be paid to the spatial distribution of the water balance components [3,4].Water balance components include: precipitation, evapotranspiration, total runoff, surface runoff, groundwater flow and the soil water content.The quantification of the various water balance components of hydrological processes in a watershed remains a challenging topic [5][6][7][8], because several unknown components, e.g., evapotranspiration and soil water content, included in the water balance are difficult to measure [9].However, the calculation of water balance components is significant for water resources assessment and management especially for water-scarce regions with regard to assessment of the impact of climate change.The water balance components can be realistically simulated using a continuous watershed hydrological model [10][11][12].Meanwhile water balance models could provide an accurate estimate of runoff, relative changes in soil moisture, evapotranspiration, etc. [13,14].
In recent years, the water balance components have been estimated at different scales based on water balance models.For example, Gregory et al. [15] simulated the water balance components at the global scale; and Herrmann et al. [16] determined the spatial difference of water balance components at the regional scale.Generally, researchers mainly focused on the water balance at the basin or the watershed scale [3,17,18].Several studies have discussed the water balance across land cover type, vegetation pattern, and ecosystem scale.For example, Jian et al. [19] quantified the changes of each water balance component for different land cover types in a small catchment.Thompson et al. [20] proposed a synthesis framework to compute the water balance from vegetation patterns and catchment scales.Knowles et al. [21] estimated the relative contributions of the two ecosystems to the water balance.These reports facilitate our current understanding of the changes of the spatial distribution of water balance components at specific scales.However, there is a major question regarding how the relationships in water balance components vary across different spatial scales in response to climate change?Thus, a multi-scale analysis focused on changes in the spatial distribution of water balance components is necessary to fully understand the interaction between the hydrological processes and environment.
The mountainous regions in inland river basins are very sensitive to climate change [22][23][24].We choose the upper reaches of Heihe River Basin (HRB) as the study area.The HRB is the second largest inland river basin in China.There are more than 1.3 million people and 266,000 ha of irrigated agricultural land in the middle and lower reaches of the river basin depending on the runoff coming from the upper watershed, which is crucial to maintain the oasis and agricultural ecosystems in the middle and lower reaches of the river.Because of the increasing water resources and eco-environmental problems, the HRB has attracted many researchers' attention in China.In the HRB, there are several studies simulating hydrological processes [25][26][27][28][29][30].However, there have been few studies conducted to estimate the climate change impact on water balance components in the mountainous watershed at different spatial scales.This study aims to understand the distribution characteristics of hydrological processes using the SWAT model and to analyze the water balance components and runoff recharge behavior at different spatial scales (watershed, landscape, elevation zones, and sub-watershed) under the background of climate change in the upper reaches of HRB, which can serve as a reference for regional water resource assessment and planning in the HRB.

Study Area
The HRB is the second largest inland river basin in China, located in the central part of the Hexi Corridor, which used to be an important section of the ancient Silk Road.The geographical coordinates range from 98 • to 101 • 30 E and 38 • to 42 • N, which is a representative inland river basin in the semi-arid and arid regions of northwest China and covers an area of approximately 130,000 km 2 .The Yingluoxia watershed is located at the main stream of HRB which is a mountainous watershed having an area of approximately 10,018 km 2 , mostly comprising mountainous terrain (Figure 1).It is the main region of runoff generation in HRB.Most of the water resources in the middle oasis and lower desert ecosystem are therefore recharged by runoff from the upper reaches.The mountainous runoff of the Yingluoxia watershed is monitored by the Yingluoxia (YLX) hydrological station.The climate of the study area is warm and wet in summer and cold and dry in winter with large spatial and temporal variability.The average annual precipitation tends to increase from west to east from 200 mm to 700 mm.The precipitation mainly occurs in wet seasons usually from June to September accounting for nearly 70% of the annual precipitation, only 3.5% of the annual precipitation occurs from December to February, and less than 10% occurs from November to March.The annual air temperature of the watershed varied from −5 to 4 • C during the period from 1961 to 2013.Altitudinal landscape zones exist in the watershed, including a desert steppe zone from 1637 to 2800 m, shrubbery forest grassland zone from 2800 to 3500 m, sub-alpine shrubbery meadow zone from 3500 to 4000 m, alpine cold-and-desert meadow zone from 4000 to 4500 m, and alpine permafrost-snow-ice zone from 4500 m to 5062 m [31].The annual air temperature of the watershed varied from −5 to 4 °C during the period from 1961 to 2013.Altitudinal landscape zones exist in the watershed, including a desert steppe zone from 1637 to 2800 m, shrubbery forest grassland zone from 2800 to 3500 m, sub-alpine shrubbery meadow zone from 3500 to 4000 m, alpine cold-and-desert meadow zone from 4000 to 4500 m, and alpine permafrost-snow-ice zone from 4500 m to 5062 m [31].

Soil and Water Assessment Tool (SWAT) Model
The SWAT model was developed by the United States Department of Agriculture.It is a semidistributed model which is based on physical processes [32].The model has some advantages in predicting climate change effects on water-related and hydrological processes over long periods, and it is a continuous-time, river basin scale model.The SWAT model can run with different time scales, including annual, monthly, and daily time scales [33].In SWAT, a basin is divided into many subbasins, which are then subdivided into many hydrologic response units (HRUs) in development that are composed of homogeneous land cover, soil type, and terrain features.The entire water balance is simulated based on each HRU, including: precipitation interception, precipitation partitioning, snowmelt water, redistribution of soil water content, evapotranspiration, lateral subsurface flow from the soil profile, and return flow from shallow aquifers.
The hydrological processes simulated by SWAT are on the basis of the following water balance equation:

Soil and Water Assessment Tool (SWAT) Model
The SWAT model was developed by the United States Department of Agriculture.It is a semi-distributed model which is based on physical processes [32].The model has some advantages in predicting climate change effects on water-related and hydrological processes over long periods, and it is a continuous-time, river basin scale model.The SWAT model can run with different time scales, including annual, monthly, and daily time scales [33].In SWAT, a basin is divided into many sub-basins, which are then subdivided into many hydrologic response units (HRUs) in development that are composed of homogeneous land cover, soil type, and terrain features.The entire water balance is simulated based on each HRU, including: precipitation interception, precipitation partitioning, snowmelt water, redistribution of soil water content, evapotranspiration, lateral subsurface flow from the soil profile, and return flow from shallow aquifers.
The hydrological processes simulated by SWAT are on the basis of the following water balance equation: where SW t is the final soil water content (mm), SW 0 is the initial soil water content on day i (mm), t is the time (days), P day is the amount of precipitation on day i (mm), Q s is the amount of surface runoff on day i (mm), E a is the amount of ET on day i (mm), w seep is the amount of water entering the vadose zone from the soil profile on day i (mm), and Q g is the amount of return flow on day i (mm).

Snow Melt Algorithm
SWAT distinguishes precipitation as rain or snow by the mean daily temperature.The threshold temperature used to categorize precipitation as rain or snow is defined by the user.If the mean daily air temperature is lower than the threshold temperature, then the precipitation will be modeled as snow and the water equivalent will be added to the snow pack.If the mean daily air temperature is higher than the threshold temperature, then the precipitation will be modeled in the form of liquid rain.The snow pack will increase with additional snowfall or will decrease with snowmelt or sublimation, and the water stored in the snow pack will be reported as a snow water equivalent.
In this case study, the temperature index with elevation band approach was applied to simulate snowmelt [34,35].The snowmelt in SWAT is calculated as a linear function of the difference between the average snow pack-maximum air temperature and the base or threshold temperature for snowmelt.The snowmelt equation on a given day is: where SNO mlt is the amount of snowmelt (mm), b mlt is the melt factor for the day (mm/day-• C), SNO cov is the fraction of the HRU area covered by snow, T snow is the snow pack temperature ( • C), T mx is the maximum air temperature ( • C), and T mlt,sno is the base temperature above which snowmelt is allowed ( • C).The melt factor is allowed a seasonal variation with maximum and minimum values occurring on summer and winter solstices: where b mlt is the melt factor for the day (mm/day-• C), b mlt6 is the melt factor for 21 June (mm/day-• C), b mlt12 is the melt factor for 21 December (mm/day-• C), and d n is the day of the year.

Glacier Melt Algorithm
A distributed temperature-index approach was used for glacier melt modeling.It calculates the melting of glacial ice based on a linear relationship between ablation and positive air temperature.This approach takes into account the effect of potential solar radiation, thus including the variability in topographic shading, aspect, and slope [34].The local melt rate is calculated as: where F M is the melt factor, r ice is the ice radiation factor.I pot denotes the potential solar direct radiation, and T is the mean daily air temperature, whereby temperature below zero indicates that no melting occurs.The potential solar direct radiation, I pot , is a function of top-of-atmosphere solar radiation, solar geometry, topographic characteristics use, and atmospheric transmissivity: where S = 1360 W•m −2 is the solar constant, (d m /d) 2 is the eccentricity correction factor of the Earth's orbit, Ψ = 0.75 is the mean atmospheric clear-sky transmissivity, P is the atmospheric pressure, P 0 is the mean atmospheric pressure at sea level, Z is the zenith angle (which is approximated as a function of latitude, solar declination, and hour angle), θ is the slope angle, φ sun is the solar azimuth angle, and φ slope is the slope azimuth angle.Potential direct solar radiation is set to zero between sunset and sunrise and for the surface that is shaded by surrounding topography.

Data Availability
The ASTER GDEM was used as the digital elevation model (DEM) data with a 30 m resolution.The meteorological and hydrological observation data were all downloaded from the portal of Environmental and Ecological Science Data Center for West China (http://westdc.westgis.ac.cn).Two periods of vegetation map were obtained (version 1.0, 2001; version 2.0, 2014), with the vegetation types including: alpine meadow, grassland, forest land, bare land, and glacier, as shown in Figure 2. The soil map (scale 1:1,000,000) was obtained from the China Soil Scientific Database (CSSD).Related main soil properties including soil particle-size composition, soil texture, soil depth, and organic carbon content for the different layers of each soil type were obtained directly from the CSSD.Bulk density, available water capacity of the soil layer, and saturated hydraulic conductivity were calculated by a SPAW (Soil-Plant-Atmosphere-Water) model [36].The modified logistic growth model with two parameters was used to convert the soil particle-size from international standard to the US standard [37].The other soil properties that were difficult to be acquired were set to the default value in the SWAT soil database.
where S = 1360 W•m −2 is the solar constant, (dm/d) 2 is the eccentricity correction factor of the Earth's orbit,  = 0.75 is the mean atmospheric clear-sky transmissivity, P is the atmospheric pressure, P0 is the mean atmospheric pressure at sea level, Z is the zenith angle (which is approximated as a function of latitude, solar declination, and hour angle), θ is the slope angle, φsun is the solar azimuth angle, and φslope is the slope azimuth angle.Potential direct solar radiation is set to zero between sunset and sunrise and for the surface that is shaded by surrounding topography.

Data Availability
The ASTER GDEM was used as the digital elevation model (DEM) data with a 30 m resolution.The meteorological and hydrological observation data were all downloaded from the portal of Environmental and Ecological Science Data Center for West China (http://westdc.westgis.ac.cn).Two periods of vegetation map were obtained (version 1.0, 2001; version 2.0, 2014), with the vegetation types including: alpine meadow, grassland, forest land, bare land, and glacier, as shown in Figure 2. The soil map (scale 1:1,000,000) was obtained from the China Soil Scientific Database (CSSD).Related main soil properties including soil particle-size composition, soil texture, soil depth, and organic carbon content for the different layers of each soil type were obtained directly from the CSSD.Bulk density, available water capacity of the soil layer, and saturated hydraulic conductivity were calculated by a SPAW (Soil-Plant-Atmosphere-Water) model [36].The modified logistic growth model with two parameters was used to convert the soil particle-size from international standard to the US standard [37].The other soil properties that were difficult to be acquired were set to the default value in the SWAT soil database.Daily weather data series: precipitation, maximum and minimum air temperature, average wind speed, relative humidity, and sunshine duration were obtained from the records of the six meteorological stations as shown in Figure 1 and Table 1.Monthly river flow data were obtained from the YLX hydrological station for the observation period from 1961 to 2013.In this study, the model was calibrated based on monthly time scale using the data from January 1961 to December 1988, and the model validated by using the data from January 1986 to December 2013 of the YLX hydrological stations.The data during January 1961 to December 1963 and January 1986 to December 1988 were set as 'warm-up' periods respectively for calibration and validation.Daily weather data series: precipitation, maximum and minimum air temperature, average wind speed, relative humidity, and sunshine duration were obtained from the records of the six meteorological stations as shown in Figure 1 and Table 1.Monthly river flow data were obtained from the YLX hydrological station for the observation period from 1961 to 2013.In this study, the model was calibrated based on monthly time scale using the data from January 1961 to December 1988, and the model validated by using the data from January 1986 to December 2013 of the YLX hydrological stations.The data during January 1961 to December 1963 and January 1986 to December 1988 were set as 'warm-up' periods respectively for calibration and validation.

Model Setup
After preparation of forcing data, four major steps must be done to setup the model: (i) watershed discretization and sub-watershed characteristics derivation; (ii) define the HRU; (iii) run the model and analyze the parameter sensitivity; and (iv) calibrate and validate the SWAT model, including uncertainty analysis.
During the watershed delineation, we divided the YLX watershed into 43 sub-watersheds using ArcSWAT tool.The sub-watersheds water balance components were used to determine the basin-scale data in this study.In the HRU definition step, four classes of slope with large ranges 0% to 5% (7.6% of the total area), 5% to 25% (36.6% of the total area), 25% to 45% (24.5% of the total area), >45% (31.3% of the total area) were addressed because the study area is an alpine terrain with numerous rugged topographies.As a result, the whole watershed was delineated into 2650 HRUs.The entire water balance was simulated based on each HRU, and then the quantities of the water balance components were added to the outlet of the sub-watershed.The sub-watershed scale water balance components were routed to the outlet of the whole watershed (watershed scale).The third step was to run the model using the prepared weather data inputs and HRU information defined in the previous steps.The simulation was run first for the calibration period from January 1961 to December 1988.The fourth step in the modeling process was calibration and validation.The process of calibration and validation was necessary and will be detailed in the next section.

Model Calibration and Uncertainty Analysis
The Sequential Uncertainty Fitting algorithm (SUFI-2) was applied to perform the calibration and uncertainty analysis [38].All sources of uncertainty such as the model structure, driving variables (e.g., rainfall), parameters, and measured data were included in the ranges of parameters, which were calibrated to bracket most of the measured data in 95% prediction uncertainty (95PPU).95PPU was calculated at the 2.5% and 97.5% levels of the cumulative distribution of output variables.P-factor and R-factor were used to quantify the goodness of calibration/uncertainty performance.A P-factor is the percentage of measured data bracketed by the 95PPU band, and an R-factor is the average width of the band divided by the standard deviation of the corresponding measured variable.Theoretically, the value of the P-factor ranges between 0 and 100%, while that of R-factor ranges between 0 and infinity.A P-factor of 1 and R-factor of 0 is a simulation that exactly corresponds to the measured data.
Based on the sensitivity analysis and our previous study [26], 16 parameters were chosen, whose final fitted parameter values were shown in Table 2.

Model Assessment
The Nash-Sutcliffe Efficiency (NSE), the determination coefficient (R 2 ), and the percent bias (PBIAS) are frequently used measures in hydrological modeling studies [39], which are calculated as: where S i is the simulated streamflow, O i is the observed streamflow at time step i, and S avg and O avg are the average simulated and observed streamflow values in time period 1, 2, . . ., n, respectively.

Model Performance
The model performance was shown in Table 3 and Figure 3.We considered that SWAT simulated the observed monthly river discharge series fairly through NSE at 0.91 and 0.92 of the YLX hydrological station during calibration and validation periods, respectively.The high NSE values demonstrated that the results of model simulation met 'very good performance' according to the performance classification of Moriasi et al. [39].Moreover, the values of R 2 in the two periods showed that the simulated results had a high correlation with the observed data, as shown in Figure 4.The values of PBIAS indicated that the model slightly underestimated when compared to the observed data.However, these underestimations were still within a reasonable scope.
Water 2016, 8, 472 8 of 16 hydrological station during calibration and validation periods, respectively.The high NSE values demonstrated that the results of model simulation met 'very good performance' according to the performance classification of Moriasi et al. [39].Moreover, the values of R 2 in the two periods showed that the simulated results had a high correlation with the observed data, as shown in Figure 4.The values of PBIAS indicated that the model slightly underestimated when compared to the observed data.However, these underestimations were still within a reasonable scope.hydrological station during calibration and validation periods, respectively.The high NSE values demonstrated that the results of model simulation met 'very good performance' according to the performance classification of Moriasi et al. [39].Moreover, the values of R 2 in the two periods showed that the simulated results had a high correlation with the observed data, as shown in Figure 4.The values of PBIAS indicated that the model slightly underestimated when compared to the observed data.However, these underestimations were still within a reasonable scope.Moreover, the 95PPU model results are shown in Figure 3, wherein the uncertainty of the modeling is shown as the shaded region that includes most types of uncertainties in the model prediction.The result showed that the P-factor was 88% and R-factor was 0.76 in the calibration period and also showed the two performance criteria were quite satisfactory.This indicated that SWAT is indeed a robust tool in the upper stream of HRB with high performance accuracy.

A Changing Trend in Water Balance Components at the Watershed Scale
In this paper, the water balance components we mentioned included the precipitation (P), the evapotranspiration (ET), the total runoff (Q), the surface runoff (Qs), the groundwater flow (Qg), and the soil water content (SW), wherein Q = Qs + Qg.The precipitation and air temperature data of the six meteorological stations were also used in the variation trends analysis.The annual variations of precipitation and air temperature were averaged by the data of TL, YNG, and QLX meteorological stations which were located in the deep of Qilian Mountain during 1964-2013 as shown in Figure 5a; the precipitation showed an upward trend, increasing at an average rate of 1.05 mm/a (significance level p < 0.05), and the annual average air temperature increased significantly by 0.034 • C/a (significance level p < 0.01).It indicated that climate change has impacted upon this region with more moisture and warming during last 50 years [40,41].According to the relationship between the annual runoff coefficient (annual runoff/annual precipitation) with precipitation and air temperature (Figure 5b), we found that there was a relatively higher correlation between runoff coefficient and precipitation, compared to that of air temperature.This indicated that precipitation is the major factor that influences the increase of streamflow in the YLX watershed.The impact of the air temperature on streamflow was non-significant.Moreover, the 95PPU model results are shown in Figure 3, wherein the uncertainty of the modeling is shown as the shaded region that includes most types of uncertainties in the model prediction.The result showed that the P-factor was 88% and R-factor was 0.76 in the calibration period and also showed the two performance criteria were quite satisfactory.This indicated that SWAT is indeed a robust tool in the upper stream of HRB with high performance accuracy.

A Changing Trend in Water Balance Components at the Watershed Scale
In this paper, the water balance components we mentioned included the precipitation (P), the evapotranspiration (ET), the total runoff (Q), the surface runoff (Qs), the groundwater flow (Qg), and the soil water content (SW), wherein Q = Qs + Qg.The precipitation and air temperature data of the six meteorological stations were also used in the variation trends analysis.The annual variations of precipitation and air temperature were averaged by the data of TL, YNG, and QLX meteorological stations which were located in the deep of Qilian Mountain during 1964-2013 as shown in Figure 5a; the precipitation showed an upward trend, increasing at an average rate of 1.05 mm/a (significance level p < 0.05), and the annual average air temperature increased significantly by 0.034 °C/a (significance level p < 0.01).It indicated that climate change has impacted upon this region with more moisture and warming during last 50 years [40,41].According to the relationship between the annual runoff coefficient (annual runoff/annual precipitation) with precipitation and air temperature (Figure 5b), we found that there was a relatively higher correlation between runoff coefficient and precipitation, compared to that of air temperature.This indicated that precipitation is the major factor that influences the increase of streamflow in the YLX watershed.The impact of the air temperature on streamflow was non-significant.Over the past 50 years, the annual average precipitation and total runoff of the watershed have increased with climate change, with an increasing trend of 5.8 mm/10a and 8.2 mm/10a, respectively.According to the calculation results, the ET of the watershed accounted for 70.9% of the precipitation.The surface runoff and groundwater flow changes remained consistent with the total runoff, showing an increasing trend with the rates of 3.5 mm/10a and 4.7 mm/10a, respectively.The annual average contribution of surface runoff and groundwater flow to total runoff were 67.1% and 32.9%, respectively.However, the contribution rate of surface runoff and groundwater flow for the total runoff increase were 42.7% and 57.3%, respectively.This meant that the contribution of the increasing of groundwater flow to the increase of total runoff was larger than the increase of surface runoff.The change trend of soil water content kept relatively stable, with a rising trend of only 1.2 mm/10a.Over the past 50 years, the annual average precipitation and total runoff of the watershed have increased with climate change, with an increasing trend of 5.8 mm/10a and 8.2 mm/10a, respectively.According to the calculation results, the ET of the watershed accounted for 70.9% of the precipitation.The surface runoff and groundwater flow changes remained consistent with the total runoff, showing an increasing trend with the rates of 3.5 mm/10a and 4.7 mm/10a, respectively.The annual average contribution of surface runoff and groundwater flow to total runoff were 67.1% and 32.9%, respectively.However, the contribution rate of surface runoff and groundwater flow for the total runoff increase were 42.7% and 57.3%, respectively.This meant that the contribution of the increasing of groundwater flow to the increase of total runoff was larger than the increase of surface runoff.The change trend of soil water content kept relatively stable, with a rising trend of only 1.2 mm/10a.

Quantification of Water Balance Components at the Landscape Scale
The composition ratio and area percentage of water balance components at landscape scale are shown in Figure 6a.The alpine meadow occupied the largest area (51.1% of the watershed), and, as a consequence, accounted for the largest portion of each water balance component.The runoff generated from the alpine meadow contributed 55.5% to the total runoff of the watershed, with 31.9% from the surface runoff and 23.6% from the groundwater flow, respectively.The runoff coefficient of the alpine meadow was 0.33 (shown in Figure 6b) which equals to the runoff coefficient of the whole watershed.Considering the runoff contribution, runoff coefficient, and the area percentage, we could see that the alpine meadow was the most important vegetation type in the YLX watershed.Grassland accounted for 14.2% of the study area, which could lead to canopy interception and the capability of regulation and storage at the root zone; the water-yielding capacity was low with the lowest runoff coefficient which was 0.24.ET from the grassland landscape accounted for 16.0% of that in the whole watershed, while the groundwater flow from the grassland only accounted for 4.7% of the total groundwater flow of the watershed, which indicated that there was little extra water infiltrating into the underground.The runoff generated from the bare land (15.4% of the watershed) contributed 15.9% to the total runoff of the watershed, with 12.6% from the surface runoff and only 3.3% from the groundwater flow.The ET and soil water content formed in the bare land were relatively low, contributing 8.4% and 7.5% of the total ET and soil water content of the watershed, respectively.Because the bare land mainly is composed by alpine cold desert, which was in a cold environment, the soil of the alpine cold desert was thinner, and was unable to hold more soil moisture.So, the runoff coefficient of the bare land was relatively higher, equaling 0.46.There was an argument in the mountainous area of the inland river basin in semi-arid and arid regions that was whether the forest could generate runoff.In this study, the forest land included a large portion of shrub, so there was 16.9% of the total runoff of the watershed generated from the forest land (forest and shrub, 18.8% of the watershed) with 12.4% from the surface runoff and 4.5% from the groundwater flow.The runoff coefficient of the forest land was 0.32 as shown in the Figure 6b.The runoff coefficient of the glacier was 0.87; however, because the area percentage was very small, the contribution only occupied 3.2% of total runoff.

Quantification of Water Balance Components at the Landscape Scale
The composition ratio and area percentage of water balance components at landscape scale are shown in Figure 6a.The alpine meadow occupied the largest area (51.1% of the watershed), and, as a consequence, accounted for the largest portion of each water balance component.The runoff generated from the alpine meadow contributed 55.5% to the total runoff of the watershed, with 31.9% from the surface runoff and 23.6% from the groundwater flow, respectively.The runoff coefficient of the alpine meadow was 0.33 (shown in Figure 6b) which equals to the runoff coefficient of the whole watershed.Considering the runoff contribution, runoff coefficient, and the area percentage, we could see that the alpine meadow was the most important vegetation type in the YLX watershed.Grassland accounted for 14.2% of the study area, which could lead to canopy interception and the capability of regulation and storage at the root zone; the water-yielding capacity was low with the lowest runoff coefficient which was 0.24.ET from the grassland landscape accounted for 16.0% of that in the whole watershed, while the groundwater flow from the grassland only accounted for 4.7% of the total groundwater flow of the watershed, which indicated that there was little extra water infiltrating into the underground.The runoff generated from the bare land (15.4% of the watershed) contributed 15.9% to the total runoff of the watershed, with 12.6% from the surface runoff and only 3.3% from the groundwater flow.The ET and soil water content formed in the bare land were relatively low, contributing 8.4% and 7.5% of the total ET and soil water content of the watershed, respectively.Because the bare land mainly is composed by alpine cold desert, which was in a cold environment, the soil of the alpine cold desert was thinner, and was unable to hold more soil moisture.So, the runoff coefficient of the bare land was relatively higher, equaling 0.46.There was an argument in the mountainous area of the inland river basin in semi-arid and arid regions that was whether the forest could generate runoff.In this study, the forest land included a large portion of shrub, so there was 16.9% of the total runoff of the watershed generated from the forest land (forest and shrub, 18.8% of the watershed) with 12.4% from the surface runoff and 4.5% from the groundwater flow.The runoff coefficient of the forest land was 0.32 as shown in the Figure 6b.The runoff coefficient of the glacier was 0.87; however, because the area percentage was very small, the contribution only occupied 3.2% of total runoff.

Quantification of Water Balance Components and Changing Trends across Elevation Zones
Water balance analysis at different elevations zones is of critical importance for understanding the water balance of an alpine mountainous area in an inland river basin.Based on the vertical gradient of the landscape and the elevation, five elevation zones were divided in this paper, including zones of: 1637 to 2800 m, 2800 to 3500 m, 3500 to 4000 m, 4000 to 4500 m, and 4500 to 5062 m.
The composition ratio and area percentage of water balance components at the elevation zone scale was demonstrated in Figure 7a.The zone of 3500 to 4000 m had the largest area (40.2% of the watershed), and meanwhile, accounted for the largest portion of each water balance component.The

Quantification of Water Balance Components and Changing Trends across Elevation Zones
Water balance analysis at different elevations zones is of critical importance for understanding the water balance of an alpine mountainous area in an inland river basin.Based on the vertical gradient of the landscape and the elevation, five elevation zones were divided in this paper, including zones of: 1637 to 2800 m, 2800 to 3500 m, 3500 to 4000 m, 4000 to 4500 m, and 4500 to 5062 m.
The composition ratio and area percentage of water balance components at the elevation zone scale was demonstrated in Figure 7a.The zone of 3500 to 4000 m had the largest area (40.2% of the watershed), and meanwhile, accounted for the largest portion of each water balance component.The runoff generated from this elevation zone contributed 38.7% of the total runoff of the watershed with 23.2% from surface runoff and 15.5% from the groundwater flow.The runoff coefficient of this zone was 0.32 (shown in the Figure 7b) which was close to the runoff coefficient of the whole watershed (0.33).The zones of 1637 to 2800 m, 2800 to 3500 m, 4000 to 4500 m, and 4500 to 5062 m covered 5.5%, 26.6%, 25.4%, and 2.3% of the watershed, and contributed 2.6%, 21.7%, 32.0, and 5.1% of the total runoff, respectively.We could see that the zone of 3500 to 4500 m occupied a great part of the watershed area (65.6%), and contributed most of the total runoff (70.7%).The zone of 4500 to 5062 m was mainly covered by snow and glacier, and the zone of 1637 to 2800 m was the arid steppe zone and the portion of each water balance component in these two zones was lower than 5% as shown in Figure 7a.However, because of the melting water of the glacier and snow in zone of 4500 to 5062 m, the runoff coefficient was the highest with the value of 0.78.The zone of 4000 to 4500 m was located in a cold and humid environment, which means the proportions of ET and soil water content were relatively low.The dominant runoff component was surface runoff, with a relatively higher runoff coefficient (0.43).The zone of 2800 to 3500 m was warm and dry, covered by thick soil.The precipitation in this zone was primarily consumed by evapotranspiration and stored in the soil gaps as soil water, so the runoff coefficient was low (0.27).The zone of 1637 to 2800 m was in a drier and warmer environment, so there was little runoff generated with the runoff coefficient was 0.20.Wang et al. [42] used the stable isotope technique to trace the major source area of the mountainous runoff generation of the Heihe River and realized that the mountainous runoff was generated mostly at high altitudes.They revealed that 80.2% of the annual total mountainous runoff amount was generated at the area above 3600 m.The results indicated that we should pay much more attention to the impact of climate change in the high altitudes on the water balance components in mountainous watersheds of inland river basins.
Water 2016, 8, 472 11 of 16 runoff generated from this elevation zone contributed 38.7% of the total runoff of the watershed with 23.2% from surface runoff and 15.5% from the groundwater flow.The runoff coefficient of this zone was 0.32 (shown in the Figure 7b) which was close to the runoff coefficient of the whole watershed (0.33).The zones of 1637 to 2800 m, 2800 to 3500 m, 4000 to 4500 m, and 4500 to 5062 m covered 5.5%, 26.6%, 25.4%, and 2.3% of the watershed, and contributed 2.6%, 21.7%, 32.0, and 5.1% of the total runoff, respectively.We could see that the zone of 3500 to 4500m occupied a great part of the watershed area (65.6%), and contributed most of the total runoff (70.7%).The zone of 4500 to 5062 m was mainly covered by snow and glacier, and the zone of 1637 to 2800 m was the arid steppe zone and the portion of each water balance component in these two zones was lower than 5% as shown in Figure 7a.However, because of the melting water of the glacier and snow in zone of 4500 to 5062 m, the runoff coefficient was the highest with the value of 0.78.The zone of 4000 to 4500 m was located in a cold and humid environment, which means the proportions of ET and soil water content were relatively low.The dominant runoff component was surface runoff, with a relatively higher runoff coefficient (0.43).The zone of 2800 to 3500 m was warm and dry, covered by thick soil.The precipitation in this zone was primarily consumed by evapotranspiration and stored in the soil gaps as soil water, so the runoff coefficient was low (0.27).The zone of 1637 to 2800 m was in a drier and warmer environment, so there was little runoff generated with the runoff coefficient was 0.20.Wang et al. [42] used the stable isotope technique to trace the major source area of the mountainous runoff generation of the Heihe River and realized that the mountainous runoff was generated mostly at high altitudes.They revealed that 80.2% of the annual total mountainous runoff amount was generated at the area above 3600 m.The results indicated that we should pay much more attention to the impact of climate change in the high altitudes on the water balance components in mountainous watersheds of inland river basins.The changing trends in water balance components for different elevation zones were shown in Table 4.The runoff from the five elevation zones were all having an increasing trend with the largest (smallest) value was 8.8 mm/10a (6.9 mm/10a) in the elevation zone of 3500 to 4000 m (4500 to 5062 m), and that all of the runoff change trends were statistically significant.So, the total runoff change trend at the watershed scale was 8.2 mm/10a.The surface runoff and groundwater flow from the five elevation zones were all increasing during the last 50 years, and the groundwater flow was more obvious, indicated by the results at the elevation zone scale.Therefore, the increasing rate of groundwater flow was larger than that of the surface runoff at the watershed scale with the value of 4.7 mm/10a and 3.5 mm/10a, respectively.The trends in precipitation change at the elevation zone scale were all positive.There was a phenomenon that increasing rates were larger as elevation increased.However, a change in the trend of precipitation at the elevation zone scale was not significant.The changing trends in water balance components for different elevation zones were shown in Table 4.The runoff from the five elevation zones were all having an increasing trend with the largest (smallest) value was 8.8 mm/10a (6.9 mm/10a) in the elevation zone of 3500 to 4000 m (4500 to 5062 m), and that all of the runoff change trends were statistically significant.So, the total runoff change trend at the watershed scale was 8.2 mm/10a.The surface runoff and groundwater flow from the five elevation zones were all increasing during the last 50 years, and the groundwater flow was more obvious, indicated by the results at the elevation zone scale.Therefore, the increasing rate of groundwater flow was larger than that of the surface runoff at the watershed scale with the value of 4.7 mm/10a and 3.5 mm/10a, respectively.The trends in precipitation change at the elevation zone scale were all positive.There was a phenomenon that increasing rates were larger as elevation increased.However, a change in the trend of precipitation at the elevation zone scale was not significant.The annual average values of water balance components at the sub-watershed scale were different as shown in Figure 8.The annual average precipitation varied from 370 to 670 mm at the sub-watershed scale, decreasing from east to west on the space.The annual average precipitation was larger in the eastern sub-watersheds than that of the western sub-watersheds.The annual average sub-watershed ET varied from 275 to 540 mm with the ET of 384.2 mm at the watershed scale, as shown in Figure 8.The total runoff varied tremendously at the sub-watershed scale with the annual average value from 37 mm to 310 mm (the annual average value was 151.5 mm during 1964 to 2013 at the watershed scale) demonstrating that the sub-watersheds had different water yield capacity because of the different underlying surface.The spatial distribution pattern of the surface runoff was similar with that of the total runoff, which was larger in the central sub-watersheds and smaller in the east and west.The annual average sub-watershed surface runoff varied from 29 mm to 275 mm (101.6 mm at the watershed scale), as shown in Figure 8.The spatial distribution pattern of the groundwater flow was opposite to that of the surface runoff, showing that there was little ground water (the annual average value was smaller than 18 mm) flowing into the river in the central sub-watersheds and a considerable amount of ground water feeding the river in the western and eastern sub-watersheds with the largest value of 126 mm in the most western sub-watershed.The annual average groundwater flow of the whole watershed was 49.9 mm according to the summation of the sub-watershed scale results.The annual average soil water content varied from 10 mm to 165 mm at the sub-watershed scale (96.9 mm at the watershed scale), as shown in Figure 8.
Trends in the change of water balance components at the sub-watershed scale were shown in Figure 9.The increasing rate of the precipitation reduced from west to east, with the largest value of 1.23 mm/a in the most western sub-watershed.Only the three most eastern sub-watersheds' precipitation were decreasing during the past 50 years, and the precipitation of the other 40 sub-watersheds were increasing, which lead to an increasing rate of 5.8 mm/10a at the watershed scale.The change trend of ET at the sub-watershed scale varied between −1.52 and 0.76 mm/a.All except one runoff of sub-watersheds had increasing trend from 0.18 mm/a to 1.81 mm/a.Consequently, the total runoff at the watershed scale presented an increasing trend of 8.2 mm/10a.The increasing rate of the total runoff in the west part of the watershed was more significant compared to that of the central and eastern regions.Trends in the sub-watershed surface runoff varied from −0.06 to 0.80 mm/a.There was only one sub-watershed having a decreasing surface runoff with the trend of −0.06 mm/a, the others had an increasing trend.The change trend of the surface runoff at the whole watershed was 3.5 mm/10a during 1964 to 2013.The spatial distribution of the change trend of groundwater flow was very similar to that of the total runoff at the sub-watershed scale.The range of the change trend of the sub-watershed groundwater flow was −0.18 to 1.30 mm/a.Similarly, the change trend of groundwater flow of the whole watershed was increasing with the value of 4.7 mm/10a, based on the results of the sub-watershed scale.The change trend of soil water content at the sub-watershed scale varied between −0.32 and 0.86 mm/a (1.2 mm/10a at the watershed scale), respectively, as shown in Figure 9.

Discussion
Mountainous watersheds are the main water generating regions in inland river basins.Better understanding of the variation characteristics and the spatial distribution of the water balance components experiencing climate change is crucial for sustainable water resource management in semi-arid and arid inland river basins [43][44][45][46].Previous studies have facilitated our current understanding at specific scales.However, the multi-scales analysis is necessary to fully understand the interaction between the hydrological processes and environment.This study provided insights into quantities and change trends of the water balance components for the water resource region at different spatial scales (watershed scale, landscape scale, elevation zone scale, and sub-watershed scale), taking the upper reaches of the Heihe River Basin as a case.From the results, we could fill in the multi-scale knowledge gap of water balance of the mountainous watershed in inland river basins.This information is very useful for developing an overview of where the water resources are coming from and what changes have been taken place during the past half century.This study will provide a theoretical reference for the water resources management of the inland river basins of China [47,48].
There are several limitations in this study.First, the lack of field observation data-for example, soil moisture and actual evapotranspiration-hampers the validation of the simulation.Each landscape has different impacts on the hydrological process and in mountains, the water balance components vary in different landscapes [22].Moreover, spatial heterogeneity of soil has great impact on the dynamic processes of hydrological systems.The lack of observational data will influence the

Discussion
Mountainous watersheds are the main water generating regions in inland river basins.Better understanding of the variation characteristics and the spatial distribution of the water balance components experiencing climate change is crucial for sustainable water resource management in semi-arid and arid inland river basins [43][44][45][46].Previous studies have facilitated our current understanding at specific scales.However, the multi-scales analysis is necessary to fully understand the interaction between the hydrological processes and environment.This study provided insights into quantities and change trends of the water balance components for the water resource region at different spatial scales (watershed scale, landscape scale, elevation zone scale, and sub-watershed scale), taking the upper reaches of the Heihe River Basin as a case.From the results, we could fill in the multi-scale knowledge gap of water balance of the mountainous watershed in inland river basins.This information is very useful for developing an overview of where the water resources are coming from and what changes have been taken place during the past half century.This study will provide a theoretical reference for the water resources management of the inland river basins of China [47,48].
There are several limitations in this study.First, the lack of field observation data-for example, soil moisture and actual evapotranspiration-hampers the validation of the simulation.Each landscape has different impacts on the hydrological process and in mountains, the water balance components vary in different landscapes [22].Moreover, spatial heterogeneity of soil has great impact

Discussion
Mountainous watersheds are the main water generating regions in inland river basins.Better understanding of the variation characteristics and the spatial distribution of the water balance components experiencing climate change is crucial for sustainable water resource management in semi-arid and arid inland river basins [43][44][45][46].Previous studies have facilitated our current understanding at specific scales.However, the multi-scales analysis is necessary to fully understand the interaction between the hydrological processes and environment.This study provided insights into quantities and change trends of the water balance components for the water resource region at different spatial scales (watershed scale, landscape scale, elevation zone scale, and sub-watershed scale), taking the upper reaches of the Heihe River Basin as a case.From the results, we could fill in the multi-scale knowledge gap of water balance of the mountainous watershed in inland river basins.This information is very useful for developing an overview of where the water resources are coming from and what changes have been taken place during the past half century.This study will provide a theoretical reference for the water resources management of the inland river basins of China [47,48].
There are several limitations in this study.First, the lack of field observation data-for example, soil moisture and actual evapotranspiration-hampers the validation of the simulation.Each landscape has different impacts on the hydrological process and in mountains, the water balance components vary in different landscapes [22].Moreover, spatial heterogeneity of soil has great impact on the dynamic processes of hydrological systems.The lack of observational data will influence the simulation accuracy.Second, we neglect the impact of human activities, the change trends of water balance components are completely related to climatic factors such as precipitation and temperature.Although several studies have shown that the climate variability is mainly influenced the hydrological processes in mountainous watersheds of the inland river basins [29,41,48-50], somehow human activities will have some influence on the water cycle.Third, the SWAT model does not incorporate the soil freezing/thawing processes.The study area is located in a mountainous region within which there is distributed a wide range of permafrost and seasonally cold regions.The soil freezing/thawing processes will influence the simulation accuracy [27].

Figure 1 .
Figure 1.The study area, meteorological and hydrological stations in the YLX watershed.

Figure 1 .
Figure 1.The study area, meteorological and hydrological stations in the YLX watershed.

Figure 3 .
Figure 3. Monthly river discharge of the YLX hydrological station during the calibration period (a); and the validation period (b); the uncertainties of simulation (95PPU); and the monthly precipitation of the YLX watershed.

Figure 4 .
Figure 4. Correlation between the observed and simulated monthly river discharge of the YLX hydrological station during the calibration period (a) and the validation period (b).

Figure 3 .
Figure 3. Monthly river discharge of the YLX hydrological station during the calibration period (a); and the validation period (b); the uncertainties of simulation (95PPU); and the monthly precipitation of the YLX watershed.

Figure 3 .
Figure 3. Monthly river discharge of the YLX hydrological station during the calibration period (a); and the validation period (b); the uncertainties of simulation (95PPU); and the monthly precipitation of the YLX watershed.

Figure 4 .
Figure 4. Correlation between the observed and simulated monthly river discharge of the YLX hydrological station during the calibration period (a) and the validation period (b).

Figure 4 .
Figure 4. Correlation between the observed and simulated monthly river discharge of the YLX hydrological station during the calibration period (a) and the validation period (b).

Figure 5 .
Figure 5.The change trend of precipitation and air temperature (a); and the correlation of runoff coefficient and precipitation and air temperature (b); from 1964 to 2013.

Figure 5 .
Figure 5.The change trend of precipitation and air temperature (a); and the correlation of runoff coefficient and precipitation and air temperature (b); from 1964 to 2013.

Figure 6 .
Figure 6.The percentage of the water balance components and area (a); and the runoff coefficient (b); at the landscape scale.

Figure 6 .
Figure 6.The percentage of the water balance components and area (a); and the runoff coefficient (b); at the landscape scale.

Figure 7 .
Figure 7.The percentage of the water balance components and area (a); and the runoff coefficient (b); at the elevation zone scale.

Figure 7 .
Figure 7.The percentage of the water balance components and area (a); and the runoff coefficient (b); at the elevation zone scale.

Water 2016, 8 , 472 13 of 16 Figure 8 .
Figure 8.The annual average value of water balances components at the sub-watershed scale.

Figure 9 .
Figure 9.The change trend of water balances components at the sub-watershed scale during 1964 to 2013.

Figure 8 . 16 Figure 8 .
Figure 8.The annual average value of water balances components at the sub-watershed scale.

Figure 9 .
Figure 9.The change trend of water balances components at the sub-watershed scale during 1964 to 2013.

Figure 9 .
Figure 9.The change trend of water balances components at the sub-watershed scale during 1964 to 2013.

Table 1 .
Hydro-meteorological stations that were used in the SWAT simulation.

Table 2 .
Parameters related in the calibration procedure as well as the fitted values.
Notes: 1 r_: parameter value is multiplied by (1 + a given value) or relative change; v_: parameter value is replaced by given value or absolute change.

Table 3 .
Model performance based on monthly river discharge.

Table 3 .
Model performance based on monthly river discharge.

Table 4 .
The change trend of water balance components for different elevation zones during 1964 to 2013 (unit: mm/10a).Notes: * represents the significance at level of 0.1; ** represents the significance at level of 0.05.3.5.Quantification of Water Balance Components and the Change Trends at the Sub-Watershed Scale