Combined use of volume-area and volume-length scaling relationships in glacio-hydrological simulation

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 the 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.

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.

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−H_ter)/(H_top−H_ter); 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 the ablation zone.

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.

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 km², with snow/glacier melt, rainfall, and groundwater as the primary sources of water. Glacier area in this basin was 38 km² (4%) in 1964. Field-observed data showed that the rate of shrinkage of the glacier area in 1960–2007 was 37%.

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 km² and the length was 2.2 km. Drainage area of Glacier No. 1 is 3.34 km², 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.

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

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.

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 km² in 1964. The maximum, minimum, and mean glacier areas were 3.28 km², 0.02 km², and 0.31 km², 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 c_a 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, v_l and q constants of the V-L scaling relation were derived from regression analysis based on the records. The estimated constants γ, c_a, q, and v_l 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, c_a and v_l 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 c_a and v_l 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)).

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.

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%).

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 km². The drainage area of Glacier No. 1 GHRU was 3.34 km², monitored at Glacier No. 1 Hydrological Station (GHS). The altitudinal distribution of Glacier No. 1 GHRU is plotted in Figure 4.

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.

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).

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).

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⁻¹, 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⁻¹ 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⁻¹. 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⁻¹, 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⁻¹ day⁻¹ (Liu et al. 1998), the maximum (SMFMX) and minimum (SMFMN) values occurring on summer and winter solstices were set at 3.5 mm °C⁻¹ day⁻¹ and 2.5 mm °C⁻¹ day⁻¹, 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⁻¹ day⁻¹ and 7.9 mm °C⁻¹ day⁻¹, 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.

Streamflow: The simulated average annual streamflow of GHS for 1980–2006 was ∼2.3 × 10⁶ m³ 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 × 10⁸ m³/10 yr, with R² = 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.

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 km² (16.0%) in 1964–2006, compared with the observed 0.27 km² (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%.

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%).

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 (Figure 7(a) and 7(b)). Also, correlation analysis suggested that the model accurately simulated the 93–94% monthly discharge (Figure 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 (Figure 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.

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).

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 km² and those of the others were less than 2.0 km². Smaller glaciers, especially with area less than 0.25 km², 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).

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.

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 shows that it was critical for smooth discharge at GHS.

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 (R² ≥ 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.

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.

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 c_a 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 c_a = 40.0 and γ = 1.35 for glaciers in Tienshan Mountains, China. Based on the HBV-EC model, Stahl et al. (2008) used c_a = 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 c_a scaling parameter (Bahr et al. 2015) and the c_l 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).

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.

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.

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 Figure 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.

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.

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).