Abstract

Glacier changes are driven by glacier melt, which in turn affects streamflow. This paper describes an accounting scheme for glacier area change distribution across elevation profiles for application in the glacier module of the Soil and Water Assessment Tool (SWAT) model. In addition to volume-area scaling relationship in the module, the paper introduced volume-length scaling relations to estimate changing glacier terminus and update glacier area changes between equilibrium line altitude (ELA) and the terminus. The improved scheme was used in the nested Urumqi Glacier No.1 catchment and Urumqi River Basin in Tienshan Mountains, China. Comparison of the simulated and observed data suggested that the new scheme accurately reproduced the length and area changes of Glacier No.1. The contributions of glacier melt and ice melt to runoff were estimated at 71% and 38% for Glacier No.1 Hydrological Station and 11.1% and 5.8% for Yingxiongqiao Hydrological Station, respectively. This helped to better interpret long-term monitored glacio-hydrological processes of Glacier No.1 and the variation of glacier melt contribution to streamflow at the catchment scale.

INTRODUCTION

Glacier melt is critical for streamflow in glacierized regions around the world. With climate change, the importance of the hydrological role of glaciers in glaciology, hydrology, and water resources has increased (Immerzeel et al. 2010). Glacio-hydrological models are powerful tools that are increasingly used to quantify the impact of glacier changes on watershed hydrology, including spatiotemporal projections of future water availability (Hock 2005; Bliss et al. 2014; La Frenierre & Mark 2014).

Glacier modeling is the basis for simulating hydrological processes in basins under substantial glacier cover. At basin scale, glacier effect is the sum of the effects of individual glaciers under different topographic and climatic conditions. Currently, individual glaciers are lumped into a single glacier unit in just one category of hydrological models and then the hydrological domain divided into glaciated and unglaciated areas, as in the Hydrologiska Byråns Vattenbalansavdelning (HBV) model (Akhtar et al. 2008), reducing the difficulties of dealing with a massive number of small glaciers. In another category, a single glacier is divided by grids in which individual glaciers and glacier patches are lumped into a bigger unit as in HBV-EC (Stahl et al. 2008; Jost et al. 2012) and Spatial Processes in Hydrology (SPHY) model (Lutz et al. 2014; Terink et al. 2015). Also using the grid approach, glacier mass balance and ice dynamics were added to the Glacier Evolution Runoff Model (GERM) (Huss et al. 2008; Farinotti et al. 2012) and the Variable Infiltration Capacity (VIC) model (Zhang et al. 2013; Zhao et al. 2015). Irrespective of the method used, uncertainties exist in glacier area or length estimation due to the decomposition individual glaciers into different patches belonging to different grid cells which thereby merges individual glaciers and/or glacier patches into a virtual glacier within a grid cell.

To consider individual glacier geometry, one category of distributed hydrological models uses the hydrologic response unit (HRU) concept. The HRU concept is used to account for spatial heterogeneities of land cover, soil and topography; as in Water Availability in Semi-Arid Environments (WASA) (Duethmann et al. 2015) and Soil and Water Assessment Tool (SWAT) (Neitsch et al. 2005). Luo et al. (2013) proposed the glacier hydrological response unit (GHRU) approach and incorporated it into SWAT to simulate hydrological processes in glacier-dominated watersheds. The GHRU approach treats glaciers individually in a sub-basin and uses elevation bands to account for altitudinal variations in precipitation and temperature.

Modeling glacier hydrology and transient runoff response as a result of glacier changes usually involves the modeling of glacier mass balance and dynamics in river basins. Due to its simplicity and lower data requirement, the combination of the degree-day method (Hock 2003; Liu et al. 2009) and volume-area (V-A) scaling relation (Chen & Ohmura 1990; Bahr et al. 2015) is commonly used in glacier mass balance and area-change calculations. This approach was also used by Luo et al. (2013) in the glacier-enhanced SWAT model. In this model, the V-A scaling converts the calculated volume-change (mass balance) into changes in glacier area, but cannot update glacier area change distribution across elevation profile. The glacier area may increase even above equilibrium line altitude (ELA) for a retreating glacier, which is usually not the case in practice. Some approaches assume glacier area to be constant and update glacier surface based on the accumulation area ratio (AAR) method (Horton et al. 2006). However, this method cannot reproduce transient response of glaciers to climate change (Li et al. 2015). Thus, Huss et al. (2008) proposed a Δh-parameterization, which is an empirical glacier-specific function based on observations, for calculating changes in glacier surface elevation and area. The approach is widely used in glacio-hydrological modeling (Huss et al. 2010; Duethmann et al. 2015; Li et al. 2015), but needs parameters that are difficult to determine in regions with scarce field data. Several methods are used to determine glacier area change with elevation bands. One such method assumes that glacier advance or retreat first occurs in the lowest elevation band (Radic et al. 2008; Moller & Schneider 2010). This assumption may ignore simultaneous retreat or advance of glaciers due to changes in length and width at lower ablation zone. Radic et al. (2008) assumed a constant-normalized glacier area-elevation distribution and glacier area loss or gain along glacier length under combined V-A and volume-length (V-L) scaling relations. Although glaciers shorten under such conditions, glacier area decreases in the highest elevation band, which hardly occurs especially for thick and extensive glaciers.

The objective of this paper was to overcome these issues by incorporating V-L relationship into SWAT glacier module (Luo et al. 2013) and redistributing glacier area between the terminus and ELA. First, the V-L scaling relation was introduced into the glacier module to determine glacier terminus elevation and the glacier area change distribution across altitudinal profile updated. Second, the improved approach was applied to the source region of Urumqi River Basin in Tienshan Mountains in China and the approach evaluated using detailed field data on Urumqi Glacier No.1 and Urumqi River Basin. Last, the streamflow components at Glacier No.1 and Urumqi River Basin were identified based on the improved approach to better interpret long-term monitored glacio-hydrological processes at Glacier No.1 and the catchment scale. The improved approach proved to be a more practical method of glacier change simulation and with an expected broad application in glacier hydrological studies.

MATERIALS AND METHODS

SWAT model

The SWAT model is a distributed, process-based hydrologic model for the simulation of different physical processes in a watershed with varying soil, land use, and management conditions over long periods of time (Gassman et al. 2007; Arnold et al. 2010). For modeling purposes, SWAT partitions a watershed into numerous sub-basins and uses ArcSWAT interface HRU to account for spatial heterogeneities of land cover, soil, and slope within the sub-basin. The hydrology simulation is split into two major divisions – the first part is the land phase of the hydrological cycle that controls the amount of water flow into the main channel in each sub-basin. The second part is the water or routing phase that defines the movement of water through the channel network of the watershed to the outlet. The sub-basins are divided into elevation bands to account for orographic effects on both precipitation and temperature. Then, snowmelt and rainfall runoff are calculated at each band. While snowmelt is controlled by air/snowpack temperature, melting rate, and areal coverage of snow, snowpack does not melt until snowpack temperature exceeds a threshold value. The snowmelt, which uses temperature-index (TI) approach, is calculated as a linear function of the difference between the average maximum air temperature of the snowpack and the base/threshold temperature for snowmelt. At present, SWAT can only simulate rainfall and snowmelt runoff processes.

The officially released SWAT versions do not simulate glacier hydrology. Luo et al. (2013) integrated a coupled glacier module into SWAT2005 (Neitsch et al. 2005) to simulate glacier dynamics and mass balance. The algorithm uses a dynamic GHRU approach for individual glacier simulations (Figure 1), where glacier area can increase or decrease in a GHRU.

Figure 1

A schematic of the SWAT-integrated GHRU approach as presented by Luo et al. (2013).

Figure 1

A schematic of the SWAT-integrated GHRU approach as presented by Luo et al. (2013).

The snow and glacier melt were calculated using the TI method (Neitsch et al. 2005), which is differently adopted for snow melt and glacier melt. Details of the glacier module are given in Luo et al. (2013) and calibration of TI parameters for snow and ice melt is described in the next sections.

Luo et al. (2013) did not describe the scheme for updating glacier area change distribution across elevation bands; and the scheme used in the code to update glacier area distribution across elevation bands is as follows:

  1. An elevation distribution function is adopted to distribute glacier area across elevation. This function normalizes glacier height, which is expressed as the (H−Hter)/(Htop−Hter); where H is the elevation and the subscripts top and ter indicate the top and terminus of the glacier.

  2. The glacier area is updated using the V-A relationship that is, in turn, used to update ELA based on simulated glacier mass balance across the elevation.

  3. The relationship between ELA and glacier area and the maximum and minimum elevations adopted by Xie & Liu (2010) is used to update the minimum elevation of the glacier. The relationship is acquired from observation data of glaciers in Tienshan Mountains. When glacier area is updated, ELA is accordingly updated and glacier tongue elevation is eventually derived from the empirical relationships discussed above. Finally, altitudinal distribution of the updated glacier area is determined through the normalized glacier area-elevation curve.

  4. Finally, the area distribution across elevation is updated using the distribution ratio profile described in (1).

The weakness of this scheme is that glacier area can increase even above ELA for retreating glaciers, which is usually not the case in practice. This work was undertaken to overcome this weakness by integrating the V-L relationship into the SWAT glacier module (Luo et al. 2013). The new scheme performs as the following functions:

  1. Glacier length-terminus elevation relationship is determined.

  2. Glacier volume and ELA are updated using glacier mass balance simulation.

  3. Glacier area change is calculated using V-A relationship.

  4. Glacier length change is calculated using the updated glacier volume via V-L scaling relationship.

  5. Glacier terminus elevation is next updated using the length-elevation relation.

  6. Glacier area change is redistributed between terminus and ELA using negative mass change weight in ablation zone.

Glacier volume-area scaling relation

Glacier V-A scaling relation is defined as (Chen & Ohmura 1990; Stahl et al. 2008; Bahr et al. 2015):  
formula
(1)
where V is glacier volume (m·km2); A is glacier area (km2); and ca and γ are empirical constants.
Volume change (ΔV) in year y (which is generally glacier mass balance for the period from the 1st September to the 31st August the following year) of an individual glacier is calculated as:  
formula
(2)
where bi is net mass balance in the ith elevation band (m w.e.); Ai is glacier area in the ith elevation band (m2); ρ is the bulk density of ice, generally taken as 0.9 kg/m3; i is the sequential number of elevation bands starting from the top of the glacier; and N is the total number of elevation bands. Glacier volume (Vy) in the successive year is updated as:  
formula
(3)
where y is the successive year. The glacier area is then updated as:  
formula
(4)
Glacier area change in the preceding year of mass balance is given as:  
formula
(5)
Then, glacier length is derived from the V-L scaling relation (Bahr 1997; Radic et al. 2008) as:  
formula
(6)
where L is glacier length (km); and vl and q are constants.
Glacier length is next updated as:  
formula
(7)
 
formula
(8)
If glacier surface is a uniform plane, change in glacier terminus elevation is estimated as:  
formula
(9)
where ΔHter is change in glacier terminus elevation (km) and α is the apparent inclination of glacier surface plane, given as:  
formula
(10)
where Htop is top elevation of glacier; Hter0 is glacier terminus elevation at the start of the simulation; and L0 is the glacier length at the start of the simulation. The terminus elevation change is derived from the change in glacier length.

Glacier area change with elevation band

In the GHRU approach, a new scheme for updating glacier area change across altitudinal profile is proposed as follows. The glacier area is assumed to change across the entire ablation zone, and no change is assumed to occur in the accumulation zone of a glacier. The accumulation and ablation zones are defined by ELA, which is determined by the altitude at which the calculated annual cumulative mass balance across the elevation band is equal to zero. Thus, glacier area change in the ith elevation band of the ablation zone is calculated for the preceding year based on the mass balance of the ablation zone as:  
formula
(11)
where M is the total number of elevation bands in the ablation zone. The approach described above is finally incorporated into the SWAT model.

Runoff component

At GHRU scale, runoff consists of water yield from glaciated and unglaciated areas. Here, runoff from the glaciated area is regarded as glacier melt, which comprises ice melt, snowmelt, and direct rainfall-runoff on the ice. Ice melt and snowmelt both are calculated using the classical degree-day method, with two degree-day factors (June 21 and December 21) as representative seasonal change patterns. During summer, rainfall can directly generate runoff on the surface of the ice. The different components derived from different origins are expected to occur on distinct temporal scales.

Following Lutz et al. (2014), the streamflow consists of rainfall runoff, snow melt and glacier melt runoff, and baseflow. Baseflow is released from groundwater storage, which is estimated using the two-reservoir method (Luo et al. 2012).

In this study, GHRU is a unique case of SWAT-based HRU that is used to simulate hydrological processes on glaciers. The GHRU-enhanced SWAT model is applied to the source region of Urumqi River, which has 4% glacier cover and long-term observatory glacier, and to Urumqi Glacier No. 1. The simulations of Glacier No. 1 Basin (Glacier No. 1 GHRU) were assessed against observational data for capacity and performance of the GHRU approach at glacier scale. Also, the simulation of Urumqi River Basin was assessed for capacity and performance of GHRU at basin scale.

Study area

Urumqi River Basin is in the Eastern Tienshan Mountains in North China (Figure 2). The river originates from the north slope of Tiager Peak, a characteristic continental temperate arid region with westerly winds. Yingxiongqiao Hydrological Station (YHS) is located at 43°22′N, 87°12′E and an elevation of 1,920 m. The drainage area of the basin above YHS is 923.4 km2, with snow/glacier melt, rainfall, and groundwater as the primary sources of water. Glacier area in this basin was 38 km2 (4%) in 1964. Field-observed data showed that the rate of shrinkage of the glacier area in 1960–2007 was 37%.

Figure 2

A map showing the location of Urumqi Glacier No.1 Basin and Urumqi River Basin in China. The numbers depict the sequence of the sub-basins.

Figure 2

A map showing the location of Urumqi Glacier No.1 Basin and Urumqi River Basin in China. The numbers depict the sequence of the sub-basins.

Urumqi Glacier No. 1 (with glacier area ratio of 58% in 1964 in the headwater region of Urumqi River) is China's only glacier cataloged in the World Glacier Monitoring Service (WGMS) network. Glacier No. 1 is a typical continental valley glacier that is representative of glaciers in inland Central Asia. Shrinkage due to increasing temperature in Tienshan Mountains reduced Glacier No. 1 into two independent glaciers in 1993, now called the East and West Branches of Glacier No. 1. Since then, the West Branch has been retreating a little faster than the East Branch due to different response to climate change and retreat rates (Xu et al. 2011). Mass balance observations have been done every year since 1959, with interruptions during 1967–1979. Mass balance is calculated on contour maps of accumulation and ablation, using data from permanent stake network on the glacier and snow pits (Cui et al. 2013). The area of Glacier No. 1 was mapped in 1962, 1964, 1973, 1981, 1986, 2001, 2006, and 2010. The volume of Glacier No. 1 was measured in 1964, 1981, 1986, 2001, and 2006. In 1964, the area was 1.94 km2 and the length was 2.2 km. Drainage area of Glacier No. 1 is 3.34 km2, and discharge from this drainage area is monitored at Glacier No. 1 Hydrological Station (GHS) (Figure 2).

The terminus elevation of Glacier No. 1 has been measured each year since 1980, except for the period 1982–1985. The distance of glacier terminus retreat was 139.72 m along the flow line in 1962–1993, with an average retreat rate of 4.5 m/a. From 1993 to 2010, the East Branch retreated by 59.6 m, at an average rate of 3.5 m/a. For the same period, the West Branch retreated by 101.6 m, at an average rate of 6.0 m/a (Zhang et al. 2014). The terminus elevation of Glacier No. 1 was 3,730 m in 1964, while it was 3,766 m for the East Branch and 3,832 m for the West Branch in 2007.

Model setup

The delineation of sub-basins and GHRUs/HRUs in the GHRU-enhanced SWAT model was based on digital elevation model (DEM), the glacier map, land cover and soil map.

DEM

Using the 90-m-resolution DEM obtained from CGIAR Consortium for Spatial Information (CGIAR-CSI), the headwater region of Urumqi River above YHS was delineated into 21 sub-basins. Also, the reaches and main channel of the river were generated from the DEM (see Figure 2). The ten elevation bands of each sub-basin were next delineated in ArcSWAT.

Glacier map

The glacier map was from the First Glacier Inventory of China released by the Cold and Arid Regions Science Data Center (Liu et al. 2000). The inventory records area, length, and volume of individual glaciers. There was a total of 124 glaciers above YHS in Urumqi River Basin, with a total area of 38 km2 in 1964. The maximum, minimum, and mean glacier areas were 3.28 km2, 0.02 km2, and 0.31 km2, respectively. The recorded glacier area in Tienshan Mountains was mapped from aerial photos and topographic maps and glacier volume in Urumqi River Basin calculated from empirical relations between measured glacier thickness and area (Liu & Ding 1986). The ca and γ constants of the V-A scaling relation were derived from regression analysis based on measured data for the 124 glaciers in Tienshan Mountains. Similarly, vl and q constants of the V-L scaling relation were derived from regression analysis based on the records. The estimated constants γ, ca, q, and vl for Urumqi River Basin were 1.383, 41.6 km−2.15, 2.534, and 8.96 km−6.6, respectively. The same values for the basin were adopted for the γ and q constants for Glacier No. 1, as the two parameters are always constant for a given region. However, ca and vl change from one glacier to another (Bahr et al. 2015). To get an accurate initial volume and length of Glacier No. 1 based on the area in 1964, the constants ca and vl were set at 41.6 and 14.4, respectively. The relative errors of calculated to observed volumes of Glacier No. 1 were within the relatively tight bounds of ±5% (Figure 3(a)). The relative errors in the calculated lengths, compared to the observed lengths, were in the range of −19.0 to +2.1% (Figure 3(b)).

Figure 3

Comparisons of estimated (a) volume by V-A scaling relation and (b) length by V-L scaling relation with the measured values of Glacier No. 1 (the labels are relative errors).

Figure 3

Comparisons of estimated (a) volume by V-A scaling relation and (b) length by V-L scaling relation with the measured values of Glacier No. 1 (the labels are relative errors).

Land cover

The land cover map was obtained from Global Land Cover Map (GlobCover 2009, 300 m resolution) (http://due.esrin.esa.int/globcover) using the UN Land Cover Classification System (LCCS). The main land cover types included mixed forest (FRST), grassland (RNGE), and barren land (SWRN); accounting respectively for 21.7%, 39.3%, and 30.8% of the basin.

Soil map

The Harmonized World Soil Database V1.2 (HWSD, available at http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/index.html?sb=1), with a scale of 1:1,000,000 for China, was used to generate the soil layer. In this study, the classification of HWSD soil map units was based on FAO-90 soil classification system. Each soil unit was defined in terms of type, texture, depth, bulk density, organic carbon content, and electrical conductivity; which are the soil variables in the SWAT soil database. The soil in the study area was classified into calcic chernozems (6.2%), glacier/permanent snow (10.6%), haplic greyzems (4.6%), gelic leptosols (44.6%), and mollic leptosols (34.1%).

GHRUs and HRUs generation

The drainage area of each GHRU consisted of glaciated and unglaciated portions. The drainage area boundaries were delineated from DEM and overlaid with the glacier map, generating 124 GHRUs in the Urumqi River Basin. Each GHRU was also delineated into ten elevation bands based on the DEM, and glaciated and unglaciated areas calculated on ArcGIS platform using the surface analysis function. Based on 1964 data, Glacier No. 1 ranked second in terms of area, with a glaciated area of 1.94 km2. The drainage area of Glacier No. 1 GHRU was 3.34 km2, monitored at Glacier No. 1 Hydrological Station (GHS). The altitudinal distribution of Glacier No. 1 GHRU is plotted in Figure 4.

Figure 4

Altitudinal distribution of glaciated and unglaciated areas of Glacier No. 1 Basin in 1964.

Figure 4

Altitudinal distribution of glaciated and unglaciated areas of Glacier No. 1 Basin in 1964.

The HRUs were generated from land cover and soil maps, and a total of 150 HRUs was generated for 21 sub-basins in Urumqi River Basin.

Meteorological data

The GHRU-enhanced SWAT model was driven by daily meteorological data from Daxigou Meteorological Station (DMS), located at an elevation of 3,543.8 m in Urumqi River Basin. Average annual air temperature in the basin for 1959–2006 is −5.1 °C, with a maximum of 0 °C and a minimum of −9.0 °C. The average temperature in summer (from June to August) is 4.2 °C. Then the mean annual precipitation is 453 mm, 67% of which is in summer. Yang et al. (1988) noted that measured precipitation over Glacier No. 1 is underestimated by ∼30% and proposed a correction factor of 1.3 (i.e., multiplying DMS data by 1.3).

Model calibration and validation

The 1964 glacier map was used as the initial distribution of glaciers in the basin, which was the starting time of the model simulation. Annual discharge at GHS for 1980–2006 was adapted after Li et al. (2010b) and used to evaluate simulated runoff for Glacier No. 1 Basin. Observations of daily discharge at YHS for 1964–1989 were obtained from China's Hydrological Yearbook compiled for the Xinjiang autonomous prefecture. The data series was split into two segments, the one for 1964–1974 used for calibration and the other for 1980–1989 used for validation.

In this study, Glacier No. 1 data (including observed glacier area, length, and mass balance) were used to calibrate glacier parameters. The data, including glacier area and terminus elevation in the Second Glacier Inventory of China (Guo et al. 2014), were compared with the data of corresponding glaciers in the First Glacier Inventory of China. The observed changes in glacier area and terminus elevation were used to calibrate glacier hydrological processes in the model.

Table 1 lists the most sensitive parameters for the GHRU-enhanced SWAT model based on other studies in Tienshan Mountains (Luo et al. 2013; Gan et al. 2015; Ma et al. 2015; Wang et al. 2015).

Table 1

List of sensitive parameters and calibration values of GHRU-enhanced SWAT model for the Urumqi River Basin

Variable Description Unit Literature values Literature Calibration values 
PLAPS Precipitation lapse rates mm km−1 220 (annual) Yang et al. (1998)  10–160 (monthly) 
TLAPS Temperature lapse rates °C km−1 −2.7–6.0 Cui et al. (2013)  −2.7–6.0 
GMFMX Ice melt factor on June 21 mm °C−1 day−1 10.1–14.3 Liu et al. (1998); Cui et al. (2013)  2.5–17.4 
GMFMN Ice melt factor on December 21 mm °C−1 day−1 0.98–5.1  1.6–10.8 
SMFMX Snow melt factor on June 21 mm °C−1 day−1 3.1 Liu et al. (1998)  3.5 
SMFMN Snowmelt factor on December 21 mm °C−1 day−1   2.5 
ALPHA_BF Baseflow alpha factor  0–1  0.02–0.8 
RCHRG_DP Deep aquifer percolation factor  0–1 Luo et al. (2012)  0.4 
CH_K2 Channel effective hydraulic conductivity mm/hr 0–300  10–50 
Variable Description Unit Literature values Literature Calibration values 
PLAPS Precipitation lapse rates mm km−1 220 (annual) Yang et al. (1998)  10–160 (monthly) 
TLAPS Temperature lapse rates °C km−1 −2.7–6.0 Cui et al. (2013)  −2.7–6.0 
GMFMX Ice melt factor on June 21 mm °C−1 day−1 10.1–14.3 Liu et al. (1998); Cui et al. (2013)  2.5–17.4 
GMFMN Ice melt factor on December 21 mm °C−1 day−1 0.98–5.1  1.6–10.8 
SMFMX Snow melt factor on June 21 mm °C−1 day−1 3.1 Liu et al. (1998)  3.5 
SMFMN Snowmelt factor on December 21 mm °C−1 day−1   2.5 
ALPHA_BF Baseflow alpha factor  0–1  0.02–0.8 
RCHRG_DP Deep aquifer percolation factor  0–1 Luo et al. (2012)  0.4 
CH_K2 Channel effective hydraulic conductivity mm/hr 0–300  10–50 

Studies have noted that temperature lapse rates (TLAPS) and precipitation lapse rates (PLAPS) in the headwater region of Urumqi River vary during the year (Cui et al. 2013). Cui et al. (2013) put TLAPS for the months from January to December respectively as −3.0, −3.2, −4.1, −5.3, −6.0, −5.3, −6.0, −5.7, −5.8, −4.2, −3.3, and −2.7 °C km−1, which values were used in this study.

Monthly PLAPS were obtained using the trial and error method. Luo et al. (2013) used annual average value of 45 mm km−1 for Manas River Basin, which underestimated streamflow in the basin. As Manas River Basin is ∼250 km from Urumqi River Basin, the monthly PLAPS was manually calibrated relative to the mass balance and area change of Glacier No. 1 and relative to discharge data at GHS and YHS at distance interval of 10 mm km−1. The calibrated PLAPS for the months from January to December were 10, 10, 20, 30, 70, 120, 160, 110, 40, 20, 20, and 10 mm km−1, respectively. Based on the values, average annual precipitation over Glacier No. 1 and Urumqi River Basin was estimated at ∼590 and ∼520 mm, respectively. Precipitation falls as liquid, solid, or a mixture of both. The model classifies precipitation as rain or freezing rain/snow using temperature boundary (SFTMP). The default value of 1.0 °C for SFTMP (Fontaine et al. 2002) was used in this study.

Snow and ice melt were both calculated by different degree-day factors in the GHRU approach (Luo et al. 2013; Gan et al. 2015; Ma et al. 2015). The degree-day factors for both ice and snow were also calibrated. As snowmelt factor for Glacier No. 1 is ∼3.1 mm °C−1 day−1 (Liu et al. 1998), the maximum (SMFMX) and minimum (SMFMN) values occurring on summer and winter solstices were set at 3.5 mm °C−1 day−1 and 2.5 mm °C−1 day−1, respectively. The GHRU defined ice melt factor allows seasonal variation with maximum (GMFMX) and minimum (GMFMN) values at summer and winter solstices (Luo et al. 2013). Calibration was first done using observed glacier area and mass balance of Glacier No.1 and streamflow records at the outlet of Glacier No.1. Initial values of ice melt were based on literature values in Table 1. The final calibrated GMFMX and GMFMN of Glacier No.1 were 5.0 mm °C−1 day−1 and 7.9 mm °C−1 day−1, respectively. For individual glaciers, it was assumed that ice melt factor did not change with location on the ice surface. As the mass balance of other glaciers was not observed, the GMFMX and GMFMN were simply multiplied by a ratio of the values obtained at Glacier No.1 from the area changes between 2007 and 1964.

The Nash–Sutcliffe efficiency (NSE), the percent bias (PBIAS), and coefficient of determination (R2) were used to evaluate the model performance. The evaluation of streamflow simulation adopted the ratings suggested by Moriasi et al. (2007). The NSE and PBIAS are calculated as:  
formula
(12)
 
formula
(13)
where Qiobs is the ith streamflow observation; Qisim is the ith simulated streamflow value; Qmean is the mean observed streamflow; and k is the total number of observations.

RESULTS

Model calibration and validation

Glacier No. 1

Streamflow: The simulated average annual streamflow of GHS for 1980–2006 was ∼2.3 × 106 m3 and the PBIAS value at GHS for 1980–2006 was 0.1%. Moreover, the observed and simulated annual streamflows all increased at 0.4 and 0.6 × 108 m3/10 yr, with R2 = 0.53, representing a satisfactory simulation result (Figure 5). However, in some years, the simulation errors were much larger, indicating some accuracy concerns about the year-to-year deviations of the simulated discharge from the observed discharge. However, the calculated evaluation indices from long-term data series and GHRU approach suggested that the simulated streamflow was reliable at glacier scale.

Figure 5

Comparison of simulated to measured annual runoff at the Glacier No. 1 Hydrological Station.

Figure 5

Comparison of simulated to measured annual runoff at the Glacier No. 1 Hydrological Station.

Glacier area, length and mass balance: Unique to Glacier No. 1, the simulation well depicted glacier retreat in the last decades (Figure 6(a)). The simulated area shrinkage of Glacier No. 1 was 0.31 km2 (16.0%) in 1964–2006, compared with the observed 0.27 km2 (13.9%) in the same period. The range of annual relative error was between −3.5% and +1.1%, with a mean of −1.4%.

Figure 6

Comparisons of simulated to measured values of (a) area, (b) length, (c) annual mass balance, and (d) cumulative mass balance of Glacier No. 1.

Figure 6

Comparisons of simulated to measured values of (a) area, (b) length, (c) annual mass balance, and (d) cumulative mass balance of Glacier No. 1.

The V-L scaling relation was used to calculate glacier length change and thus change in glacier terminus elevation. The simulated glacier length retreat before 1993 agreed well with the observed one (Figure 6(b)). After 1993, the simulated values were closer for the West Branch than for the East Branch. The range of deviation of the simulated annual length retreat from the observed length was −6.8 to +7.6 m for the West Branch and −5.1 to +10.7 m for the East Branch. Correspondingly, the range of annual change in terminus elevation was −2.4 to +2.7 m for the West Branch and −1.8 to +3.2 m for the East Branch. The model-simulated glacier length was only the ‘virtual length’ of the glacier branches.

The simulation captured the patterns of inter-annual variations in glacier mass balance (Figure 6(c)), with simulated and observed values in fair agreement in terms of cumulative mass balance (Figure 6(d)). Although glacier mass loss was overestimated for 1970–1995, the mean difference between the observed and simulated mass balance for the period of study was 0.028 (−0.733, 0.550) m w.e., with a mean relative error of 5% (−507%, +671%).

Urumqi River Basin

Streamflow at hydrological stations: At basin scale, NSE for the simulated streamflow at YHS for the calibration and validation periods was greater than 0.65, which was ‘good’ (Table 2). The high NSE values suggested that the hydrological simulation reliably captured the seasonal patterns of discharge. The hydrographs also showed very good agreement between daily simulated and observed values (Figures 7(a) and 7(b)). Also, correlation analysis suggested that the model accurately simulated the 93–94% monthly discharge (Figures 7(c) and 7(d)). The PBIAS values for the two periods were both within 10%, rated as ‘very good’ (Table 2). However, the model overestimated or underestimated some daily flash floods (Figures 7(a) and 7(b)). This over/underestimation could be due to the lapse rates used, which were based on meteorological data from a single station to derive precipitation and temperature values for the entire basin. The approach neglected orographic precipitation and assumed that precipitation at the station meant precipitation everywhere in the basin and vice versa. Furthermore, since the study area was the headwater region of Urumqi River and was relatively less influenced by anthropogenic activities, the simulation did not consider the impact of such activities on the hydrological process.

Table 2

Evaluation statistics for simulated daily discharge at Yingxiongqiao Hydrological Station

Period Data series NSE Rating PBIAS (%) Rating R2 
Calibration 1965–1974 0.73 Good 1.6 Very good 0.79 
Validation 1980–1989 0.72 Good 2.3 Very good 0.77 
Period Data series NSE Rating PBIAS (%) Rating R2 
Calibration 1965–1974 0.73 Good 1.6 Very good 0.79 
Validation 1980–1989 0.72 Good 2.3 Very good 0.77 
Figure 7

Comparisons of simulated and measured discharges at Yingxiongqiao Hydrological Station in Urumqi River Basin: (a) daily discharge for calibration period, (b) daily discharge for validation period, (c) monthly discharge for calibration period, and (d) monthly discharge for validation period.

Figure 7

Comparisons of simulated and measured discharges at Yingxiongqiao Hydrological Station in Urumqi River Basin: (a) daily discharge for calibration period, (b) daily discharge for validation period, (c) monthly discharge for calibration period, and (d) monthly discharge for validation period.

Glacier area and terminus elevation change in Urumqi River basin: The GHRU approach in combination with V-A and V-L scaling relations was used to determine the rates of change in size of individual glaciers at basin scale. The area retreat ratios of the individual glaciers were calculated for 2007 relative to 1964 based on glacier simulations results driven by the glacier inventories. Correlation analysis suggested that the ratios of the simulated glacier retreat agreed well with glacier inventory retreat (Figure 8). Some small glaciers disappeared in the inventory, but persisted until the end of the simulation, with some 10% or so of the glacier areas remaining. The simulation showed 33% decline in total glacier area, which was comparable to the 37% reduction derived from the two phases of glacier inventory and other studies (Li et al. 2010a; Wang et al. 2011).

Figure 8

Comparison of simulated area reduction of glaciers to that derived from two phases of glacier inventory in Urumqi River Basin.

Figure 8

Comparison of simulated area reduction of glaciers to that derived from two phases of glacier inventory in Urumqi River Basin.

The basin-wide glaciers were classified into four groups based on glacier size (Table 3). Based on glacier retreat processes since 1964, the area of one of the five glaciers in group 4 (G4) was 3.28 km2 and those of the others were less than 2.0 km2. Smaller glaciers, especially with area less than 0.25 km2, shrank faster than larger ones (Figure 9(a)). Recent studies showed that small glaciers have larger variabilities both in magnitude and area changes direction due to their rapid response to climate change (Kriegel et al. 2013; Unger-Shayesteh et al. 2013).

Table 3

Glacier change statistics for the Urumqi River Basin study area

Glacier group Size (km2Area in 1964 Observed area change (2007–1964)
 
Simulated area change (2006–1964)
 
km2 km2 km2 
G1 Size <0.25 7.82 −3.91 −50 −3.68 −47 
G2 0.5 > Size ≥ 0.25 9.93 −3.97 −40 −3.67 −37 
G3 1 > Size ≥ 0.5 11.40 −3.76 −33 −3.19 −28 
G4 Size ≥ 1 8.80 −2.29 −26 −2.02 −23 
Total  37.95 −13.93 −37 −12.56 −33 
Glacier group Size (km2Area in 1964 Observed area change (2007–1964)
 
Simulated area change (2006–1964)
 
km2 km2 km2 
G1 Size <0.25 7.82 −3.91 −50 −3.68 −47 
G2 0.5 > Size ≥ 0.25 9.93 −3.97 −40 −3.67 −37 
G3 1 > Size ≥ 0.5 11.40 −3.76 −33 −3.19 −28 
G4 Size ≥ 1 8.80 −2.29 −26 −2.02 −23 
Total  37.95 −13.93 −37 −12.56 −33 
Figure 9

Different sizes of simulated glacier area shrinkage (a) and observed and simulated glacier terminus elevation change (b) in the Urumqi River Basin.

Figure 9

Different sizes of simulated glacier area shrinkage (a) and observed and simulated glacier terminus elevation change (b) in the Urumqi River Basin.

Also, terminus elevation change calculated from change in glacier length was simulated. The terminus elevation changes in individual glacier were different and the simulated glacier amounts within different elevation bands were close to observed values (Figure 9(b)). During 1964–2006, ∼70% of the glaciers retreated by 0–200 m, but some 10% slightly advanced, which was mainly for glaciers in high-altitude regions. Moreover, only about 5% of the glaciers had extensive length retreat in excess of 300 m. Therefore, the terminus elevation of glaciers in Urumqi River Basin mostly increased by 0–200 m.

In general, the simulated results based on the combination of the GHRU approach and the V-A and V-L scaling relations provide changes in glacier area, length, and terminus elevation for each glacier. The parameters related to the glacier can be calibrated more accurately using this approach. The approach could therefore control the simulated changes in glacial morphology based on observed data and, in turn, determine glacier response to climate change at basin scale.

Simulated GRHU-based runoff components

Glacier No.1 GHRU runoff components

GHRU-simulated runoff consists of glacier melt from the glaciated area and then snowmelt and water yield from the unglaciated area. The GHRU-simulated glacier melt contributed 71% to runoff at GHS during 1964–2006. Ice melt, snowmelt on ice, and rainfall-runoff on ice accounted, respectively, for 54%, 30%, and 16% of the glacier melt. Seasonally, runoff hydrographs for unglaciated area, snowmelt on ice, and rainfall runoff on ice were similar (Figure 10). The runoff components started rising in April and dropped in July through October. There was ice melt concentration in July through September, accounting for 93% of the annual volume. While the ratio of ice melt to glacier melt was 54%, the hydrograph in Figure 10 showed that it was critical for smooth discharge at GHS.

Figure 10

Simulated intra-annual runoff components at Glacier No. 1 Hydrological Station for the period 1964–2006.

Figure 10

Simulated intra-annual runoff components at Glacier No. 1 Hydrological Station for the period 1964–2006.

Runoff component hydrographs in Urumqi River Basin

The GHRU approach was to estimate meltwater for each individual glacier in Urumqi River Basin. The contribution ratios of glacier melt to river discharge varied significantly from glacier to glacier. The contribution to river runoff of the largest glacier was ∼1%, with Glacier No. 1 ranking second (0.54%). The glacier contribution ratio was linearly correlated with glacier area (R2 ≥ 0.97), a near-overlay of the 1:1 slope and a near-zero Y-intercept. In total, glacier melt contributed 11.1% to discharge into the Urumqi River in the period 1964–2006, close to the 10.0% proportion reported by Yang (1991).

Discharge at YHS included two components: surface runoff and baseflow. Part of the rainfall, snowmelt and glacier melt left HRU or GHRU as surface runoff, while the other part infiltrated into the soil profile and later discharged as saturated lateral flow. Additionally, part of the infiltration water recharged the aquifer systems in the basin and was later released from storage as baseflow.

Urumqi River is fed mainly by rainfall runoff (Table 4). The basin precipitation is concentrated more in summer, accounting for 59.6% of the annual amount. The estimated shares of glacier melt, snowmelt, rainfall runoff, and baseflow were 11.1%, 10.6%, 37.8%, and 40.5%, respectively. Rainfall-runoff accounted for 63.5% of the annual surface runoff. Surface runoff was mainly concentrated in summer, accounting for 73.8% of the annual volume. Then, rainfall-runoff accounted for 70% of the summer surface runoff.

Table 4

Statistics of simulated seasonal precipitation, temperature and runoff components at Yingxiongqiao Hydrological Station for the period 1964–2006

Season Temperature (°C) Precipitation (mm) Discharge (mm) Runoff component (%)
 
Surface runoff
 
Baseflow 
Glacier melt Snow melt Rainfall 
Spring −1.9 111.5 24.1 6.2 32.6 26.0 35.2 
Summer 8.0 308.9 206.7 12.6 10.0 45.8 31.6 
Autumn −1.8 61.2 38.3 9.4 2.9 12.6 75.1 
Winter −12.4 36.4 10.8 0.0 0.0 0.0 100.0 
Annual −2.0 518.0 279.9 11.1 10.6 37.8 40.5 
Season Temperature (°C) Precipitation (mm) Discharge (mm) Runoff component (%)
 
Surface runoff
 
Baseflow 
Glacier melt Snow melt Rainfall 
Spring −1.9 111.5 24.1 6.2 32.6 26.0 35.2 
Summer 8.0 308.9 206.7 12.6 10.0 45.8 31.6 
Autumn −1.8 61.2 38.3 9.4 2.9 12.6 75.1 
Winter −12.4 36.4 10.8 0.0 0.0 0.0 100.0 
Annual −2.0 518.0 279.9 11.1 10.6 37.8 40.5 

Due to the small glacier coverage, glacier melt had little effect on intra-annual discharge hydrograph (Figure 11). Although summer glacier melt accounted for 83.6% of the annual volume, summer glacier melt generated only 12.6% of summer runoff. Moreover, glacier melt included ice melt, supra-glacial snowmelt and rainfall-runoff on ice. The contributions of these three components to glacier melt runoff were 49%, 34% and 17%, respectively. Then ice melt contributed only 5.8% to total annual runoff, and glacier had minimal effect on streamflow of the basin.

Figure 11

Simulated intra-annual runoff components at Yingxiongqiao Hydrological Station for the period 1964–2006.

Figure 11

Simulated intra-annual runoff components at Yingxiongqiao Hydrological Station for the period 1964–2006.

DISCUSSION

V-A and V-L scaling relations

In this study, the combination of V-A and V-L scaling relations in the GHRU approach was designed to determine changes in glacier geometry as well as glacier influence on hydrological processes. Glacier terminus is essential for the simulation of ice melt, which is usually stronger in the glacier tongue than in other parts of the glacier ablation zone (Vincent & Six 2013). Glacier area simulation is also critical for calculating mass balance in accumulation and ablation zones as well as for glacier melt (Zhang et al. 2007; Yang et al. 2011; Engelhardt et al. 2014). The two power-law relationships were initially established purely on an empirical basis (Bahr 1997; Bahr et al. 1997). From detailed review of the physical processes of V-A and V-L scaling relations, Bahr et al. (2014) noted that exponent γ of V-A scaling relation was constant (1.375 for glaciers), and that although ca was constant for any given glacier, it changed from region to region. In practice, these constants are derived from field observations. Bahr et al. (1997) verified the physical basis of V-A scaling relationship from 144 glacier measurements and obtained a constant of γ = 1.36. Liu et al. (2003) obtained a constant of ca = 40.0 and γ = 1.35 for glaciers in Tienshan Mountains, China. Based on the HBV-EC model, Stahl et al. (2008) used ca = 28.5 and γ = 1.35 to estimate glacier area change in Bridge River Basin in Southern Chilcotin Mountains of Canada. From the geometry of over 300 glaciers, Bahr (1997) noted that glacier volume varied with length and that the exponent q = 2.0–2.1. Radic et al. (2008) used q = 2.2 to test the sensitivity of glacier volume evolution.

The V-A scaling relation was proposed for application in dense glacier regions and not for individual glaciers. This is because very large errors were possible when applied to single glaciers for which the ca scaling parameter (Bahr et al. 2015) and the cl parameter in V-L scaling relation were not previously established. This study aimed to generate glacio-hydrological processes at basin scale using the GHRU approach. Because of the availability of long-term observation data, Glacier No. 1 was used as a case in point to test GHRU performance. The application of the relations to Glacier No. 1 was carefully assessed in terms of simulated length bias which increased after the glacier separation in 1993. However, the annual change in glacier terminus elevation was within ±5 m, which was far less than the elevation band height (usually within 50–100 m or more). Although the separation of Glacier No.1 had a significant effect on the mass balance dynamics and response to climate change for the two branches (Xu et al. 2011), the model treated the two separated glaciers as one glacier. Thus, the model could not simulate glacier disintegration, a major weakness of the GHRU approach. There is the need for further improvement on this element of glacier simulation because glacier disintegration is a common phenomenon during continuous shrinkage. The simulated glacier mass balance was also underestimated after the separation (Figure 6(c)). The simulated annual glacier length showed a considerable variance with the observed value (Figure 6(b)). This was partly due to the use of V-A and V-L scaling relations in modeling glacier changes, which assumes no lag in γ and q constants in response to climate change (Moller & Schneider 2010; Bahr et al. 2015). Nevertheless, the agreement between model results and observations for Glacier No. 1 was good over the long run. The overall performance of the model simulation indicated that the combination of V-A and V-L scaling relations robustly simulated glacier geometrical change and glacier melt at basin scales (Radic et al. 2008; Lutz et al. 2014).

Glacier influence at different scales

As the GHRU approach was built on individual glaciers, the modified SWAT model determines glacier influence on hydrological processes at different scales (e.g., sub-basin scale) and on grouped glaciers of different sizes. Most importantly, the model was designed to analyze the response of glaciers and runoff to climate change. Previous studies have highlighted the importance of glaciers in compensating for precipitation shortfalls and buffering river runoff during dry summers (Jansson et al. 2003; Sorg et al. 2014).

Based on the intra-annual hydrograph at GHS and YHS, it was noted in this study that ice melt more significantly influenced runoff delay at glacier scale than at basin scale. At the inter-annual scale for Glacier No. 1, total runoff at GHS varied significantly from year to year within the range of 348–929 mm (Figure 12(a)). The range of annual ice melt was 48–491 mm and with an average of 248 mm, contributing 13.9–58.7% to GHS discharge with a mean of 36.8%. For the basin, the range of annual runoff was 167–452 mm, with an average of 280 mm (Figure 12(b)). However, ice melt was only a minor component of the total runoff, which was 2.0–12.7% and with a mean of 5.8%. Although precipitation and temperature both continued to increase since 1964, runoff at the two scales responded differently to the two meteorological elements. Runoff at GHS increased by 6.53 mm/a, far higher than the increase in precipitation (Figure 12(a)). However, runoff at YHS increased by 2.57 mm/a, less than that of the basin-wide precipitation (Figure 12(b)). Evapotranspiration increased as well at the two scales, but not as much as precipitation. Ice melt at GHS increased by 3.04 mm/a, while the linear rate was negative at YHS, indicating that the discrepancy in inter-annual runoff at the two scales was due mainly to ice melt. Total runoff from Glacier No. 1 increased significantly due to the increase in both temperature and precipitation. Li et al. (2010b) also noted that increase in glacier runoff coincides with increase in temperature and precipitation. Moreover, runoff in river basins due mainly to ice melt is largely affected by temperature. Wang et al. (2015, 2016) reported that total runoff in river basins with large glacier area ratios and ice melt contributions to runoff are mostly affected by temperature and that it is mainly influenced by precipitation in basins with low ice melt. For the more strongly glacierized Sari-Djaz catchment, the dominant role of increasing temperature was the leading cause of discharge changes (Duethmann et al. 2015). Then, for the Kakshaal catchment (which is covered by less glacier), the major part of the increase in water input is due to increased precipitation (Duethmann et al. 2015). The large difference between these areas was mainly due to glacier coverage ratios at different scales.

Figure 12

Plots of simulated inter-annual runoff components at (a) Glacier No. 1 Hydrological Station and (b) Yingxiongqiao Hydrological Station.

Figure 12

Plots of simulated inter-annual runoff components at (a) Glacier No. 1 Hydrological Station and (b) Yingxiongqiao Hydrological Station.

Uncertainties

Although the simulation of the basin runoff using GHRU-enhanced SWAT model had good results, there were some uncertainties, due mainly to the input data and parameters. One parameterization limitation of the glacio-hydrological model was the lack of sufficient data in Urumqi River Basin. This had a compensation effect on precipitation and temperature (Luo et al. 2013; Baraer et al. 2015) and some bias contributions by glacier melt to river runoff. Fortunately, weather data at DMS, Glacier No.1 mass balance and runoff at GHS were available for the calibration of Glacier No.1 parameters. Two glacier inventories (1960–2007) were used to calibrate each glacier area change. However, the use of data from only one weather station to drive the model was a possible source of uncertainty in the hydrological processes at the basin scale. Figure 13 shows the plots of sensitivity of PLAPS and TLAPS to the simulated basin-wide streamflow and glacier area change by changing the calibrated values at regular steps. It was noted that while change in PLAPS significantly influenced NSE and PBIAS of simulated streamflow, change in TLAPS mainly affected PBIAS. The model performance was at least ‘satisfactory’ when the PLAPS varied within ±20% or TLAPS varied within ±4 °C. Because PBIAS range was large (±25%), ranging from ‘satisfactory’ to ‘very good,’ the simulated total runoff varied in a broad range too. Also the inner ratios of the runoff component changed relatively more substantively, especially glacier melt proportion. This indicated that a tight PBIAS boundary was essential for glacio-hydrological simulation in glacierized basins. In this study, calculated PBIAS values were both within ±3% for the calibration and validation periods. This suggested that PBIAS had a minor effect on determining glacier melt contribution to total runoff.

Figure 13

Sensitivities of PLAPS and TLAPS to simulated streamflow at Yingxiongqiao Hydrological Station (a and b) and relative change in glacier area in Urumqi River Basin (c and d).

Figure 13

Sensitivities of PLAPS and TLAPS to simulated streamflow at Yingxiongqiao Hydrological Station (a and b) and relative change in glacier area in Urumqi River Basin (c and d).

As the glaciers were sensitive to climate change and change in TLAPS and PLAPS significantly affected glacier change simulations, the best values of PLAPS and TLAPS for the basin was determined by the glacier area change. It is evident in Figures 13(c) and 13(d) that uncertainty in TLAPS had more effect on glacier area change than that in PLAPS. When PLAPS increased or decreased by 10% or 20%, glacier area retreat increased or decreased by 7% or 15% during the 2000s. When again TLAPS increased or decreased by 1 °C, glacier area retreat decreased or increased by 12% or 14% during the 2000s. For the same period, when TLAPS increased or decreased by 2 °C, glacier area retreat decreased or increased by 22% or 30%. Therefore, glacier area change was an important index for PLAPS and TLAPS determination in glacierized river basins.

Also, the V-A and V-L approach used in this study had some uncertainties. However, the advantage of this approach was that it needed less data to parameterize at the cost of neglecting surface elevation change. Neglecting glacier surface elevation change caused some uncertainties in long-term simulations. This is important in the investigation of long-term glacier change under warming climatic conditions. It, however, must be pointed out, that this approach required extensive field data to parameterize the model, which were usually not available in our case. Incorporation of surface elevation change into the current approach was an interesting idea that needs more attention in further research studies.

CONCLUSIONS

In this study, the GHRU approach was updated with V-A and V-L scaling relations and used to determine glacier terminus elevation and realize simulating glacier retreat and advance in the whole ablation zone. Based on observed values for Glacier No. 1, changes in glacier area, length, and runoff were simulated and the model parameters validated at glacier scale. Both NSE and PBIAS indices of streamflow at the hydrological stations suggested good GHRU-enhanced SWAT model performance in reproducing runoff at the basin scale. Based on the improved GHRU approach, changes in glacier area and terminus elevation were also well simulated at basin scale.

Based on the simulated results of GHRU-enhanced SWAT model, streamflow components at GHS and YHS were determined. The contributions of glacier melt and ice melt to runoff were, respectively, 71% and 38% at GHS (glacier area ratio = 58%) and, respectively, 11.1% and 5.8% at YHS in Urumqi River Basin (glacier area ratio = 4%). At both scales, the contributions of glacier to streamflow were significantly different. The GHRU-enhanced SWAT model better determined glacier influence on hydrology in glacier-covered basins under climate change.

Uncertainties in the simulation of the glacio-hydrological processes for specific glaciers arose from the parameterization of V-A and V-L scaling relations, which were based on glacier groups rather than individual glaciers. Nevertheless, the regional parameterization of the scaling relations met the object of the study, which was the development of a glacio-hydrological approach for application at basin or regional scale as against an approach for application at glacier scale. At basin scale, uncertainties in the simulated glacier changes via the GHRU approach could only be significant for simulations driven by scarce data on glacier change. Therefore, more observation data were also critical for accurate simulation of glacial-hydrologic modeling.

ACKNOWLEDGEMENTS

This study was funded by Chinese Academy of Sciences (No. KZZD-EW-12), Natural Science Foundation of China (Grant No. 41130641 and 41741025), and China Postdoctoral Science Foundation (No. 2016M60116).

REFERENCES

REFERENCES
Arnold
,
J. G.
,
Allen
,
P. M.
,
Volk
,
M.
,
Williams
,
J. R.
&
Bosch
,
D. D.
2010
Assessment of different representations of spatial variability on SWAT model performance
.
Trans. Asabe
53
(
5
),
1433
1443
.
Bahr
,
D. B.
1997
Width and length scaling of glaciers
.
J. Glaciol.
43
(
145
),
557
562
.
Bahr
,
D. B.
,
Meier
,
M. F.
&
Peckham
,
S. D.
1997
The physical basis of glacier volume-area scaling
.
J. Geophys. Res.-Sol. Earth
102
(
B9
),
20355
20362
.
Bahr
,
D. B.
,
Pfeffer
,
W. T.
&
Kaser
,
G.
2014
Glacier volume estimation as an ill-posed inversion
.
J. Glaciol.
60
(
223
),
922
934
.
Bahr
,
D. B.
,
Pfeffer
,
W. T.
&
Kaser
,
G.
2015
A review of volume-area scaling of glaciers
.
Rev. Geophys.
53
(
1
),
95
140
.
Baraer
,
M.
,
McKenzie
,
J.
,
Mark
,
B. G.
,
Gordon
,
R.
,
Bury
,
J.
,
Condom
,
T.
,
Gomez
,
J.
,
Knox
,
S.
&
Fortner
,
S. K.
2015
Contribution of groundwater to the outflow from ungauged glacierized catchments: a multi-site study in the tropical Cordillera Blanca, Peru
.
Hydrol. Process.
29
(
11
),
2561
2581
.
Bliss
,
A.
,
Hock
,
R.
&
Radic
,
V.
2014
Global response of glacier runoff to twenty-first century climate change
.
J. Geophys. Res-Earth
119
(
4
),
717
730
.
Chen
,
J.
&
Ohmura
,
A.
1990
Estimation of Alpine glacier water resources and their change since the 1870s
.
IAHS Publ.
193
,
127
135
.
Duethmann
,
D.
,
Bolch
,
T.
,
Farinotti
,
D.
,
Kriegel
,
D.
,
Vorogushyn
,
S.
,
Merz
,
B.
,
Pieczonka
,
T.
,
Jiang
,
T.
,
Su
,
B. D.
&
Guntner
,
A.
2015
Attribution of streamflow trends in snow and glacier melt-dominated catchments of the Tarim River, Central Asia
.
Water Resour. Res.
51
(
6
),
4727
4750
.
Engelhardt
,
M.
,
Schuler
,
T. V.
&
Andreassen
,
L. M.
2014
Contribution of snow and glacier melt to discharge for highly glacierised catchments in Norway
.
Hydrol. Earth Syst. Sci.
18
(
2
),
511
523
.
Farinotti
,
D.
,
Usselmann
,
S.
,
Huss
,
M.
,
Bauder
,
A.
&
Funk
,
M.
2012
Runoff evolution in the Swiss Alps: projections for selected high-alpine catchments based on ENSEMBLES scenarios
.
Hydrol. Process.
26
(
13
),
1909
1924
.
Fontaine
,
T. A.
,
Cruickshank
,
T. S.
,
Arnold
,
J. G.
&
Hotchkiss
,
R. H.
2002
Development of a snowfall-snowmelt routine for mountainous terrain for the soil water assessment tool (SWAT)
.
J. Hydrol.
262
(
1–4
),
209
223
.
Gassman
,
P. W.
,
Reyes
,
M. R.
,
Green
,
C. H.
&
Arnold
,
J. G.
2007
The soil and water assessment tool: historical development, applications, and future research directions
.
Trans. Asabe
50
(
4
),
1211
1250
.
Guo
,
W. Q.
,
Xu
,
J. L.
,
Liu
,
S. Y.
,
Shangguan
,
D. H.
,
Wu
,
L. Z.
,
Yao
,
X. J.
,
Zhao
,
J. D.
,
Liu
,
Q.
,
Jiang
,
Z. L.
,
Wei
,
J. F.
,
Bao
,
W. J.
,
Yu
,
P. C.
,
Ding
,
L. F.
,
Li
,
G.
,
Li
,
P.
,
Ge
,
C. M.
&
Wang
,
Y.
2014
The Second Glacier Inventory Dataset of China (Version 1.0)
.
Cold and Arid Regions Science Data Center at Lanzhou
.
Doi:10.3972/glacier.001.2013.db
.
Horton
,
P.
,
Schaefli
,
B.
,
Mezghani
,
A.
,
Hingray
,
B.
&
Musy
,
A.
2006
Assessment of climate-change impacts on alpine discharge regimes with climate model uncertainty
.
Hydrol. Process
20
(
10
),
2091
2109
.
Huss
,
M.
,
Farinotti
,
D.
,
Bauder
,
A.
&
Funk
,
M.
2008
Modelling runoff from highly glacierized alpine drainage basins in a changing climate
.
Hydrol. Process.
22
(
19
),
3888
3902
.
Huss
,
M.
,
Jouvet
,
G.
,
Farinotti
,
D.
&
Bauder
,
A.
2010
Future high-mountain hydrology: a new parameterization of glacier retreat
.
Hydrol. Earth Syst, Sci.
14
(
5
),
815
829
.
Immerzeel
,
W. W.
,
van Beek
,
L. P. H.
&
Bierkens
,
M. F. P.
2010
Climate change will affect the Asian water towers
.
Science
328
(
5984
),
1382
1385
.
Jansson
,
P.
,
Hock
,
R.
&
Schneider
,
T.
2003
The concept of glacier storage: a review
.
J. Hydrol.
282
(
1–4
),
116
129
.
Jost
,
G.
,
Moore
,
R. D.
,
Menounos
,
B.
&
Wheate
,
R.
2012
Quantifying the contribution of glacier runoff to streamflow in the upper Columbia River Basin, Canada
.
Hydrol. Earth Syst. Sci.
16
(
3
),
849
860
.
Kriegel
,
D.
,
Mayer
,
C.
,
Hagg
,
W.
,
Vorogushyn
,
S.
,
Duethmann
,
D.
,
Gafurov
,
A.
&
Farinotti
,
D.
2013
Changes in glacierisation, climate and runoff in the second half of the 20th century in the Naryn basin, Central Asia
.
Global Planet Change
110
,
51
61
.
Li
,
Z. Q.
,
Li
,
K. M.
&
Wang
,
L.
2010a
Study on recent glacier changes and their impact on water resources in Xinjiang, northwestern China
.
Quatern. Sci.
30
(
1
),
96
106
.
Li
,
Z. Q.
,
Wang
,
W. B.
,
Zhang
,
M. J.
,
Wang
,
F. T.
&
Li
,
H. L.
2010b
Observed changes in streamflow at the headwaters of the Urumqi River, eastern Tianshan, central Asia
.
Hydrol. Process.
24
(
2
),
217
224
.
Liu
,
C. H.
&
Ding
,
L. F.
1986
The newly progress of glacier inventory in Tianshan Mountains
.
J. Glaciol. Geocryol.
8
(
2
),
167
170
.
Liu
,
S. Y.
,
Ding
,
Y. J.
,
Wang
,
N. L.
&
Xie
,
Z. C.
1998
Mass balance sensitivity to climate change of the Glacier No.1 at the Urumqi River Head, Tianshan Mts
.
J. Glaciol. Geocryol.
20
(
1
),
9
13
.
Liu
,
C. H.
,
Shi
,
Y. F.
,
Wang
,
Z. T.
&
Xie
,
Z. C.
2000
Glacier resources and their distributive characteristics in China – A review on Chinese glacier inventory
.
J. Glaciol. Geocryol.
22
(
2
),
106
112
.
Liu
,
S. Y.
,
Zhang
,
Y.
,
Zhang
,
Y. S.
&
Ding
,
Y. J.
2009
Estimation of glacier runoff and future trends in the Yangtze River source region, China
.
J. Glaciol.
55
(
190
),
353
362
.
Luo
,
Y.
,
Arnold
,
J.
,
Allen
,
P.
&
Chen
,
X.
2012
Baseflow simulation using SWAT model in an inland river basin in Tianshan Mountains, Northwest China
.
Hydrol Earth Syst. Sci.
16
(
4
),
1259
1267
.
Lutz
,
A. F.
,
Immerzeel
,
W. W.
,
Shrestha
,
A. B.
&
Bierkens
,
M. F. P.
2014
Consistent increase in High Asia's runoff due to increasing glacier melt and precipitation
.
Nat. Clim. Change
4
(
7
),
587
592
.
Ma
,
C. K.
,
Sun
,
L.
,
Liu
,
S. Y.
,
Shao
,
M. A.
&
Luo
,
Y.
2015
Impact of climate change on the streamflow in the glacierized Chu River Basin, Central Asia
.
J. Arid Land
7
(
4
),
501
513
.
Moriasi
,
D. N.
,
Arnold
,
J. G.
,
Van Liew
,
M. W.
,
Bingner
,
R. L.
,
Harmel
,
R. D.
&
Veith
,
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. Asabe
50
(
3
),
885
900
.
Neitsch
,
S. L.
,
Arnold
,
J. G.
,
Kiniry
,
J. R.
&
Williams
,
J. R.
2005
Soil and Water Assessment Tool Theoretical Documentation, Version 2005
.
Radic
,
V.
,
Hock
,
R.
&
Oerlemans
,
J.
2008
Analysis of scaling methods in deriving future volume evolutions of valley glaciers
.
J. Glaciol.
54
(
187
),
601
612
.
Sorg
,
A.
,
Mosello
,
B.
,
Shalpykova
,
G.
,
Allan
,
A.
,
Clarvis
,
M. H.
&
Stoffel
,
M.
2014
Coping with changing water resources: the case of the Syr Darya river basin in Central Asia
.
Environ. Sci. Policy
43
,
68
77
.
Stahl
,
K.
,
Moore
,
R. D.
,
Shea
,
J. M.
,
Hutchinson
,
D.
&
Cannon
,
A. J.
2008
Coupled modelling of glacier and streamflow response to future climate scenarios
.
Water Resour. Res.
44
(
2
),
W02422
.
Terink
,
W.
,
Lutz
,
A. F.
,
Simons
,
G. W. H.
,
Immerzeel
,
W. W.
&
Droogers
,
P.
2015
SPHY v2.0: spatial processes in HY drology
.
Geoscientific Model Development
8
,
2009
2034
.
doi:10.5194/gmd-8-2009-2015
.
Unger-Shayesteh
,
K.
,
Vorogushyn
,
S.
,
Farinotti
,
D.
,
Gafurov
,
A.
,
Duethmann
,
D.
,
Mandychev
,
A.
&
Merz
,
B.
2013
What do we know about past changes in the water cycle of Central Asian headwaters?
A Review. Global Planet Change
110
,
4
25
.
Wang
,
S. J.
,
Zhang
,
M. J.
,
Li
,
Z. Q.
,
Wang
,
F. T.
,
Li
,
H. L.
,
Li
,
Y. J.
&
Huang
,
X. Y.
2011
Glacier area variation and climate change in the Chinese Tianshan Mountains since 1960
.
J. Geogr. Sci.
21
(
2
),
263
273
.
Wang
,
X. L.
,
Luo
,
Y.
,
Sun
,
L.
,
He
,
C. S.
,
Zhang
,
Y. Q.
&
Liu
,
S. Y.
2016
Attribution of runoff decline in the Amudarya river in Central Asia during 1951–2007
.
J. Hydrometeorol.
17
(
5
),
160309141121003
.
Xie
,
Z. C.
&
Liu
,
C. H.
2010
Introduction to Glaciology
.
Shanghai Science Popularization Press
,
Shanghai
(in Chinese).
Xu
,
X. K.
,
Pan
,
B. L.
,
Hu
,
E.
,
Li
,
Y. J.
&
Liang
,
Y. H.
2011
Responses of two branches of Glacier No. 1 to climate change from 1993 to 2005, Tianshan, China
.
Quatern. Int.
236
,
143
150
.
Yang
,
Z. N.
1991
Glacier Water Resources in China
.
Gansu Science and Technology Press
,
Lanzhou
.
Yang
,
D. Q.
,
Jiang
,
T.
,
Zhang
,
Y. S.
&
Kang
,
E. S.
1988
Analysis and correction of errors in precipitation measurement at the head of Urumqi River,Tianshan
.
J. Glaciol. Geocryol.
10
(
4
),
384
399
.
Yang
,
W.
,
Guo
,
X. F.
,
Yao
,
T. D.
,
Yang
,
K.
,
Zhao
,
L.
,
Li
,
S. H.
&
Zhu
,
M. L.
2011
Summertime surface energy budget and ablation modeling in the ablation zone of a maritime Tibetan glacier
.
J. Geophys. Res-Atmos.
116
,
D14116
.
Zhang
,
L.
,
Su
,
F.
,
Yang
,
D.
,
Hao
,
Z.
&
Tong
,
K.
2013
Discharge regime and simulation for the upstream of major rivers over Tibetan Plateau
.
J. Geophys. Res. Atmos.
118
,
8500
8518
.
doi:10.1002/jgrd.50665
.
Zhang
,
G. F.
,
Li
,
Z. Q.
,
Wang
,
W. B.
&
Wang
,
W. D.
2014
Rapid decrease of observed mass balance in the Urumqi Glacier No. 1, Tianshan Mountains, central Asia
.
Quatern. Int.
349
,
135
141
.
Zhao
,
Q. D.
,
Zhang
,
S. Q.
,
Ding
,
Y. J.
,
Wang
,
J.
,
Han
,
H. D.
,
Xu
,
J. L.
,
Zhao
,
C. C.
,
Guo
,
W. Q.
&
Shangguan
,
D. H.
2015
Modeling hydrologic response to climate change and shrinking glaciers in the highly glacierized Kunma Like River Catchment, Central Tian Shan
.
J. Hydrometeorol.
16
,
2383
2402
.
doi:10.1175/jhm-d-14-0231.1
.