Introduction

The terrestrial cryosphere, including seasonal snow cover and glaciers, has reduced in extent and volume in recent decades (WGMS 2017). This is considered to be primarily a response to increasing global temperatures, with secondary effects associated with local precipitation patterns (Marzeion et al. 2014). Across High Mountain Asia (HMA), glacier trends show a mixed signal, with slight positive mass balance trend in the Kunlun, negative mass balances in the eastern Himalaya, and somewhat stable glacier systems in the Pamir Alay and Karakorum (Brun et al. 2017). This is due in part to higher-elevation regions of HMA remaining near or below freezing during much of the year, even in the presence of a warmer climate (Bolch et al. 2012). In light of these changing environmental conditions, societally relevant questions arise regarding the sustainability of melt-sourced river flow and groundwater utilized by people, agriculture, and industry downstream. Assessing the impact of a changing cryosphere on water resources requires quantifying the individual contributions of snow melt and ice melt to runoff, because each has a discretely different role in the timing and amount of water arriving downstream.

Seasonal snowpacks renew in winter and melt in the spring or early summer, providing the bulk of runoff in the early high flow season. In contrast, glacier ice melt water is sourced from a long-term reservoir that has been built over decades to centuries. Ice melt supplies streamflow later in the summer after the snowpack has been exhausted, and it can act as a streamflow buffer during low snow years or drought. The long-term influence of climate on each of these water sources also varies, necessitating that they be treated separately. This segregated approach will likely help to better anticipate areas of future water stress and water conflict, enabling early adaptation and mitigation strategies.

Most attempts at assessing changes to the cryosphere typically have combined seasonal snow and glacier melt contributions together in downstream hydrology. Notably, some indirect attempts have been made to separate these two sources. Methods include using a depletion curve to estimate a percentage of runoff at the end of the melt season to be the glacier ice melt contribution (Immerzeel et al. 2009), or estimating runoff from glacier ice melt as all melt generated within the glacierized grid cell, including seasonal snow that also exists in the Randolph Glacier Inventory (RGI) (Arendt et al. 2012) outlines (Lutz et al. 2014). Hydrograph separation has also been pursued, which identifies glacier melt as the primary input towards river discharge after all seasonal snow below the estimated regional snow line elevation (SLE) has been depleted (Mukhopadhyay and Khan 2014). While these methods are logical, they fall short of actually and independently quantifying the specific contributions of glacier ice melt and seasonal snow melt.

Here we use a spatially consistent suite of remote sensing data and reanalysis input to assess the geographically and seasonally varying roles of ice and snow melt to river flow across five major river basins of HMA: Syr Darya, Amu Darya, Indus, Ganges, and Brahmaputra (northwest to southeast, Fig. 1). These rivers cover a range of geographic and climatic settings, from arid Central Asia to the monsoonal eastern Himalaya. This assessment is systematic in its consistent methodology and use of data sources across this large region. This work was conducted in the context of formal capacity-building collaboration with 11 partner institutions in 8 countries across HMA (Bhutan, Nepal, India, Pakistan, Afghanistan, Kyrgyzstan, Tajikistan, Kazakhstan) as part of the CHARIS (Contributions to High Asia Runoff from Ice and Snow) project. Partners provided hydrological, meteorological, glacier mass balance, and hydrochemistry data, and also participated in training workshops covering a range of glaciological and geospatial skills. This collaboration is essential to the development of a thorough and systematic regional-scale assessment of the individual contributions of seasonal snow melt and glacier ice melt to water resources of the region, and for the success of long-term water resource research in the region.

Fig. 1
figure 1

Study area outlines of the five major river basins and the calibration basins used in this study. The major basins were obtained from the Global Runoff Data Centre (GRDC), with the Indus basin outline modified to eliminate a lobe in the northeast of the basin (Khan et al. 2014) which matches that used by Pakistan’s Water and Power Development Authority (WAPDA). Also shown are 20% snow probability and persistent snow/ice, for the period of record (2001–2014) derived from Moderate Resolution Imaging Spectroradiometer (MODIS) products: MOD10A1 (snow cover, light blue) and MODICE (MODIS Persistent Ice, minimum snow and ice, dark blue) at 500 m resolution. Snow probability refers to the likelihood of snow cover on any given day over the period of record. Pie charts represent estimated mean annual contributions of snow, ice, and rainfall to runoff in the subject basins above 2000 m. Background shaded relief image from Natural Earth

Full size image

Study region and basins

The study area encompasses the five glacierized mountain ranges, the Himalaya, Karakoram, Hindu Kush, Pamir, and Tien Shan, which contain the headwaters of five major rivers of High Asia: the Ganges, Brahmaputra, Indus, Amu Darya, and Syr Darya, totaling roughly three million square kilometers (Fig. 1). The region covers a wide range of geographic, area-elevation distributions and climatic conditions (Figs. 1 and online resource Fig. S1), which significantly affect the role of melt water in relation to other inputs in each basin. For example, the Brahmaputra has about half the total area of the Ganges basin but nearly three times the area above 2000 m (Table 1).

Table 1 CHARIS major basins, above indicated gauge location and above 2000 m. Basin delineation based on shapefiles from the Global Runoff Data Centre (GRDC), with the outline of the Indus modified to remove a lobe in the northeast portion of the Upper Indus, as identified by Khan et al. (2014). Basin areas include internal endorheic basins. Persistent snow and ice area estimated from MODICE (2 Strikes, 2001–2014). Relative percent contributions calculated with respect to sum of Rain (Aphrodite V1101R2), and modeled melt from: snow on land (SOL), snow on ice (SOI), and exposed glacier ice (EGI), for elevations above 2000 m
Full size table

The dry western mountain ranges such as the Karakoram, Pamir, and to a somewhat lesser extent the Tien Shan are characterized by little precipitation during the dry summer season, causing the hydrology to be dependent on spring and summer snow melt followed by glacier ice melt bolstering inputs later in the summer (Hagg and Mayer 2016). In contrast, in the humid central and eastern Himalayan areas, discharge is dominated by monsoon precipitation with glacier ice melt contributions augmenting monsoon inputs. Therefore, in many areas of the Himalaya glacier ice melt serves as a flow buffer during drought rather than a major source of river flow (Thayyen and Gergan 2011).

To capture this west to east variability in climate and glacier regimes, we identified a headwater catchment for each major basin to use for calibration (Fig. 1). These basins serve to calibrate the melt model across the different climatic settings. Calibration basins were chosen because they (a) have discharge records spanning several years during the MODIS satellite period, (b) include snow and ice melt inputs, and (c) are of computationally efficient size.

Glacier ice melt and seasonal snow melt modeling methodology

Melt modeling framework

We estimate snow and ice melt volume as the product of melt depth, and snow or ice covered area, following the approach of Rango and Martinec (1995). Snow and ice melt are estimated using a temperature index (TI) melt model, also referred to in the literature as a degree-day melt model (Hock 2003). TI models are based on the strong empirical relationship between snow or ice melt and positive air temperatures. They are well suited to applications with limited availability of meteorological forcing data, especially when melt estimates are required at daily or longer time-scales and for semi-distributed levels of spatial aggregation (e.g., elevation bands). The approach has been used widely in hydrological modeling, ice dynamic modeling, and climate sensitivity studies (Hock 2003). However, because TI models use only near-surface air temperature to estimate melt, they do not account for the influence of topographic shading, slope or aspect on melt rate. Nevertheless, comparisons of TI models with full energy balance models indicate that the simple empirically based models perform as well as or, in some cases, better than more complex models when accurate meteorological information is not available (Réveillet et al. 2018; Magnusson et al. 2015).

The melt model was run over the 2001–2014 period coincident with the MODIS record. This 14-year period includes both high and low flow years. Near surface air temperatures are from ERA-Interim reanalysis. Areas of seasonal snow cover and permanent snow and ice are derived from MODIS using the MODIS Snow Covered Area and Grain Size (MODSCAG) algorithm (Painter et al. 2009). To account for spatial variability in climate, the five major basins are separated into sub-basins. For each sub-basin, daily snow and ice melt volumes are calculated for 100 m elevation bands and aggregated to give basin total runoff. For each of the five major basins, total runoff is the sum of runoff from the sub-basins. The model is calibrated on monthly and annual runoff for the five calibration sub-basins (hereafter calibration basins) (Fig. 1). The TI model, snow and ice mapping, generation of daily temperature forcing, other model inputs, and calibration procedures are described in more detail in the following sections.

Model input

Elevation data

A fundamental requirement for implementing any snow and ice melt model is the elevation distribution of each type of surface, which can be extracted from a digital elevation model (DEM), represented here as gridded data. In our model, the DEM is used to (1) derive the altitudinal distribution, or hypsometry, of the various surface types; (2) downscale air temperatures; and (3) define the boundaries of the drainage basin using watershed analysis tools. The horizontal and vertical accuracy requirements differ somewhat for these objectives, the last being the most demanding. Flow direction depends on the first derivative of the DEM and flow divides can be sensitive to the location and magnitude of slopes.

We use a void-filled version of the Shuttle Radar Topography Mission (SRTM) DEM produced under the NASA MEaSUREs Program, identified as “SRTMGL3” (NASA Jet Propulsion Laboratory 2013). Data were obtained from NASA’s Land Processes Distributed Active Archive Center (LPDAAC) in 1° × 1° tiles and then mosaicked for the study region. In SRTMGL3, the void filling was achieved using elevation data primarily from the ASTER Global Digital Elevation Model (GDEM) v2 data (NASA Jet Propulsion Laboratory 2009). The ASTER GDEM2 data brought in spurious data attributable to unscreened cloud contamination. We identified these errors by thresholding on deviations from a median filter. These were filled by interpolation, exclusive of locations where rapid elevation changes, such as at peaks and ridges, triggered the median filter. Grids were re-projected from the native geographic coordinate system to the MODIS sinusoidal grid at 500 m resolution for delineation of major basin boundaries and melt modeling. The 90 m version of the DEM, re-projected to Albers Equal Area Conic, was used to delineate boundaries for calibration basins.

Near-surface air temperature

Near-surface air temperatures are necessary to estimate snow and ice melt. However, observations of air temperatures in the mountains of High Asia (and in other regions) are sparse, especially at elevations with seasonal snow cover and glaciers. Temperatures must either be extrapolated from valley stations or derived from air temperatures from global or regional numerical weather prediction models, such as atmospheric reanalyses. Reanalyses are well suited to deriving temperatures for the large extent of the study region but are not without their problems. Near-surface temperature fields can have biases on the order of 10 °C when compared with temperatures measured at surface meteorological stations. These biases vary both in space and with time (online resource Fig. S2).

We use a quasi-physical methodology following the approach of Jarosch et al. (2012) to downscale air temperatures from the European Reanalysis (ERA)-Interim (Dee et al. 2011) to the 500 m by 500 m MODIS grid resolution. In brief, the downscaling method uses analysis fields of upper-atmosphere air temperatures to derive a reference level temperature (taken here as 0 m a.s.l.) and lapse rate for each ERA-Interim grid cell, using linear regression. Derived lapse rates and reference temperatures are then interpolated to 500 m grid cells using bilinear interpolation. Finally, near-surface air temperatures are calculated from interpolated reference temperature, lapse rate, and grid-cell elevation. Anomaly correlations and root mean squared error between daily downscaled temperatures and measured temperatures were calculated for winter (DJF), spring (MAM), summer (JJA), and autumn (SON) for 13 meteorological stations in the Upper Indus Basin, spanning elevations from 1479 to 4440 m, for the period 2001–2010 (online resource Fig. S3a–d). Correlation coefficients are between 0.5 and 0.9 for spring, summer, and fall when the bulk of snow and ice melt occur, but between 0.3 and 0.8 in winter. Root mean squared errors (RMSE) are between 3 and 7 °C.

Daily snow and ice maps

Time series of maps of daily snow and ice cover are derived from a combination of MODSCAG fractional snow covered area (fSCA) and grain size (Painter et al. 2009), and a subsidiary product MODICE (Painter et al. 2012a). MODSCAG fSCA have been cloud masked and temporal data gaps filled (Dozier et al. 2008; Rittger et al. 2012, 2016). Four surface type classes are mapped: (1) bare ground or snow-free debris covered glacier ice, (2) snow on land (SOL), (3) snow on glacier ice (SOI), and (4) exposed glacier ice (EGI). The three snow and ice surface classes are the basis for separating runoff into snow and ice melt components. As an initial step, maps of persistent snow and ice are generated from maps of minimum exposed snow and ice for each year between 2001 and 2014 produced by a modified version of the MODICE algorithm. MODICE uses daily maps of MODSCAG fSCA to identify pixels where fSCA drops below 15% at any point in the year as “ice-free.” The remaining snow and ice covered pixels represent the minimum extent of snow and ice for that year. Pixels with persistent snow and ice are defined as those MODICE pixels classified as snow or ice in five or more years of the 14-year record. This “one-strike” classification of “ice free” pixels can be problematic in areas of high topographic relief and geolocation uncertainty because shading can cause MODSCAG to return spuriously low fSCA. To solve this problem, we require that fSCA must drop below 15% for at least two dates before a pixel is classified as “ice-free”: the “two-strike” method. MODICE maps provide a constant basemap that is used to distinguish SOL and SOI surface classes. Persistent snow and ice pixels are classified as snow or ice in a following step.

Figure 2 demonstrates the workflow of the classification process once persistent snow and ice have been identified. For each day, MODSCAG maps are used to identify bare ground or debris covered ice pixels where fSCA is less than 15%. All other pixels are classified as snow or ice covered. Snow covered pixels that are not coincident with persistent snow and ice pixels are classified as SOL. SOI and EGI pixels are then classified based on thresholding MODSCAG grain sizes (Painter et al. 2009) with the expectation that grain size retrievals are smaller for snow and larger for ice. Grain size thresholds were determined by iteratively comparing maps of snow and ice produced from MODSCAG grain size using different thresholds with maps of snow and glacier ice derived from suitable monthly scenes of 30 m Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) collected of the Hunza basin in 2013 using a semi-automated band-ratio technique (online resource Fig. S4). Three alternative products in addition to MODSCAG grain sizes were also tested against the 30 m Landsat 8 maps: grain size from the MODIS Dust Radiative Forcing in Snow (MODDRFS) (Painter et al. 2012b), albedo from the MODIS BRDF/Albedo product (MCD43) (Schaaf et al. 2011), and albedo from MOD10A1 (Klein and Stroeve 2002). Binary metrics (precision, recall and F) for the four products are presented in Table 2. MODSCAG grain sizes were best able to differentiate snow from ice. Other products scored lower in binary metrics, for example MOD10A1 never successfully maps snow or ice on some of the exposed glacier tongues (online resource Fig. S5).

Fig. 2
figure 2

Schematic showing classification of surface type

Full size image
Table 2 Binary statistics calculated by comparing each Moderate Resolution Imaging Spectroradiometer (MODIS) product to Landsat 8 for detecting snow on ice versus exposed glacier ice averaged over all pixels shown in Fig. 5. MODIS Snow Covered Area and Grain size (MODSCAG) performs the best in both recall and precision producing the highest F statistic
Full size table

Melt modeling implementation

The melt model calculates melt from SOI, SOL and EGI surfaces (“Daily snow and ice maps” section) as the product of the melt depth and area of each surface. Melt depth is the product of near-surface air temperature and a degree day factor (DDF, mm °C−1 day−1), when temperature is greater than 0 °C. We assume no melt when temperatures are below freezing. Separate DDFs are used for snow (SOI and SOL) and ice (EGI) surfaces. Snow and ice DDFs vary seasonally, approximated by a sine function with a maximum at the northern hemisphere summer solstice (June 21), following the approach used in SRM SNOW-17 (Anderson 2006). Each snow and ice DDF functions are defined by two parameters, a minimum and maximum. The TI-melt model is implemented for a major basin by dividing the major basin into sub-basins. The melt model is only run for sub-basins that have snow cover during the 2001 to 2014 period. For a selected sub-basin, the SRTMGL3 DEM at 500 m is then used to divide the basin into 100 m elevation bands. Time series of daily SOL, SOI, and EGI area are calculated for each elevation band (“Daily snow and ice maps” section). Time series of downscaled near-surface air temperature (“Near-surface air temperature” section) are also compiled for each elevation band. Basin runoff is the sum of melt from SOI, SOL, and EGI surfaces, aggregated to monthly totals, and basin monthly rainfall. Evapotranspiration is removed from this sum for the calibration process described next, but not removed in summaries for the five major basins (e.g., Figs. 1, 3, 4, and 5). Monthly basin rainfall is estimated from Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation (APHRODITE) V1101R1 APHRO_MA data set (Yatagai et al. 2012). Monthly basin ET is estimated from MOD16 (Running and Qiaozhen Mu - University of Montana and MODAPS SIPS - NASA 2015).

Fig. 3
figure 3

Major basin contributions to melt by elevation bands (bottom x-axis), overlaid with the elevational distribution of glacier area (black line, top x-axis, based on Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM)). Note the variation in x-axis scale between basins

Full size image
Fig. 4
figure 4

Mean monthly hydrologic inputs summed across elevations greater than 2000 m in the five basins (note that y-axis range varies). Large extractions for irrigation at lower elevations and non-natural flow patterns due to reservoir operations in some basins prevent meaningful comparisons with streamflow at lower elevation pour points. The volume contribution by basin at each elevation band is provided in Online Resource Table S2

Full size image
Fig. 5
figure 5

Annual percent contribution of modeled melt and rainfall from APHRODITE (Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation) over the period 2001–2014 for each of the five major basins

Full size image

The model is calibrated against measured streamflow for five calibration basins (Fig. 1.) using a simulated annealing (SA) algorithm (Henderson et al. 2003). For each calibration basin, the SA algorithm is used to explore the large potential solution space and identify a set of best snow and ice, minimum and maximum DDFs by finding a global minima for the cost function defined by equally weighted components of annual volume difference and monthly RMSE. We explored 100 cycles and 50 potential trials per cycle, for a possible 5000 combinations, plus 1 initial set, of the four parameter DDFs. At each step in the process, randomly generated sets of potential new parameters were chosen, with the following constraints: (1) 0.0 < DDF ≤ 60 mm °C day−1, (2) minimum snow (ice) DDF ≤ maximum snow (ice) DDF, and( 3) minimum snow DDF ≤ minimum ice DDF. The best 20 sets of DDFs (defined as generating the lowest calibration cost) are used to produce an ensemble melt estimate for all years with available snow cover and temperature data (2001–2014). Ensemble results are summarized to estimate potential variability in modeled melt volumes. The best model DDFs for each calibration basin are then used for all sub-basins in the respective large basins.

Results

Contributions to runoff above 2000 m in the five major basins from melt from SOL, SOI, and EGI, and rain are summarized in Table 1. We limit our analysis to above 2000 m because contributions to runoff below 2000 m are dominated by rain. The Ganges is distinct from the other four basins because rain exceeds melt contributions to runoff, with rain contributing 52% of runoff and combined melt contribution about 48%. In the four other basins, melt from SOL dominates runoff, contributing more than 65%. In contrast to the Ganges, rain contributes between 23 and 26% in these other basins. The Indus and Brahmaputra basins have the largest SOI contribution (6% and 7%, respectively). Melt from EGI is small for all basins, especially for the Ganges where EGI contributions are less than 1%. Of note is the Amu Darya, in which EGI contributions are 8%, exceeding the contribution from SOI. The higher percentage contributions to runoff from EGI in the arid, western basins are due to the combination of ice melt volumes relative to the lower overall water inputs to the basin (Table S2).

Figure 3 shows the contributions of melt (km3) SOL, SOI, and EGI by elevation, along with the distribution of glacier area by elevation (km2). Melt from SOL dominates at all elevations, except for the 5000 m to 6000 m band in the Indus, where combined melt from SOI and EGI is the dominant contribution. The Indus basin has the largest glacierized area at high elevations. While melt from EGI noticeably contributes to inputs at the highest elevations, combined melt from SOL and SOI is still larger than melt from EGI in all elevation bands and all basins.

Seasonal patterns of melt contributions are shown in Fig. 4. July is the month with the highest runoff in the Indus, Ganges, and Brahmaputra. In the Amu Darya and Syr Darya, May and April, respectively, are the months with highest runoff. Again, the Ganges stands out because it is the only basin where rain contributions exceed melt water contributions in July, August and September. In the other four basins, the combined melt water contributions (SOL + SOI + EGI) exceed the contributions from rain in all months above 2000 m. Melt contributions from SOI and EGI peak in July and August for all basins except the Brahmaputra where the peak is in June and July. In the Syr Darya, Amu Darya and Indus the SOI + EGI peak occurs after much of the SOL melt has occurred. In the Ganges and Brahmaputra, months with the highest runoff from SOI and EGI coincide with peak flow, reflecting the compounding effect of the monsoon and melt on flow in these basins. Only in the Amu Darya, does melt from EGI and SOI make a significant contribution, especially in August and September when these components contribute close to half of runoff in those months. EGI and SOI also make significant contributions to runoff in the Indus but these contributions remain less than contributions from SOL or rain, even in July and August when contributions from EGI and SOI are at their highest. These contrasting seasonal patterns of runoff reflect arid climate conditions of the Syr Darya and Amu Darya basins and the monsoonal climates of the Ganges and Brahmaputra. The Indus, while drier than the monsoonal Ganges and Brahmaputra, still receives some monsoon rainfall, especially in the eastern tributaries, for example the Chenab, Jhelum, and Sutlej. These contrasting seasonal patterns as well as the varied relative importance of snow and ice melt to basin runoff across HMA’s east-west gradient have been noted in other studies (e.g., Immerzeel and Bierkens 2012; Mankin et al. 2015).

Year to year variations in runoff components expressed as percent of total estimated contributions are shown in Fig. 5. SOL and rain exhibit similar variability between basins. Except for the Ganges, variations of SOL are larger than variations of rain. In the Ganges, both the magnitude and variability of SOL and rain are similar. SOL and rain in this basin appear to be inversely correlated, suggesting that melt from SOL may compensate for years with lower rain contributions and vice versa. Both variability and magnitude of SOI and EGI are smaller than the magnitudes of other components. It is likely that variations in these components only influence water supply at higher elevations or in months where these components make up a larger proportion of runoff.

Discussion

Comparison of the results presented here with estimates of runoff contributions presented in other studies is difficult because each study has used different definitions for glacier ice and snow melt, and has calculated runoff for different basin definitions. Immerzeel et al. (2010) estimated glacier ice melt contributes about 40% of meltwater runoff in the Indus basin above 2000 m. This contribution is much larger than our estimated 3% contribution from EGI relative to the much larger snowmelt (SOI + SOL) contributions, 73% (Table 1). Other published estimates of glacier melt contribution for the Indus basin are for either small head water basins (21%) (Mukhopadhyay and Khan 2014) or for the larger Upper Indus basin defined by the stream gauge at Besham Qila (32%) (e.g., Immerzeel et al. 2009; Lutz et al. 2014). Using a simple physical distributed process model based on area-altitude hypsometry, Yu et al. (2013) estimate the glacier ice contribution to the Indus to be 18%. The same study suggests snowmelt likely generates most of the remaining Indus River inputs since the majority of river flow is realized during the summer melt season. While direct comparison is not possible, it is worth noting that all of these published estimates of glacier melt contribution are considerably larger than the combined contribution from EGI and SOI presented here. Additional melt in other studies may be from snow on land nearby the glaciers that is not separated out. As our results indicate, the additional contribution from SOL generally overwhelms ice melt contributions as a much larger relative proportion of inputs to the basin. The nature of other methods prevents full separation of snow melt from the ice melt, with some quantity of snow melt implicitly lumped into the ice melt calculation. The other major basins in our study (Ganges, Brahmaputra, Amu, and Syr) have received less attention from the community than the Indus.

Our results (Figs. 1, 3, 4, and 5) demonstrate the varying importance of ice melt and snow melt across elevations, seasons, and regions. Glacier ice melt plays a noticeable role in river flow only above 3000 m in the Amu Darya and Indus basins. While these are highly agricultural regions with economies heavily reliant on consistent irrigation water, most agricultural activity occurs below 3000 m. Across HMA in general, seasonal snow melt is the more dominant melt water input to river flow above 3000 m. In terms of cryospheric water resource vulnerabilities, changes to snowpack accumulation, extent, and melt timing are poised to have a profound impact on the current and future hydrologic regimes of most of HMA. The high dependence of both runoff and groundwater on melt inputs (Hill et al. 2017) make the Aral Sea basin especially vulnerable to future water stress in a warmer world.

Other cryospheric reservoirs, such as rock glaciers and permafrost, may be sources of runoff but their contribution has not been estimated. Rock glaciers and permafrost have been shown to provide relatively substantial contributions at the headwater scale in mountain basins (Cowie et al. 2017), but their relevance is anticipated to be trivial over the regional domain of this study given the magnitude of other inputs. Contributions of melt from debris covered glaciers have also not been included in our estimates. Thick debris cover on glaciers in the Karakoram and Himalaya generally reduces melt but melt rates are not zero (e.g., Mihalcea et al. 2006). Our estimates, therefore, likely underestimate melt contribution.

Uncertainties imposed by air temperature

The downscaling approach used in this study, and described in the “Near-surface air temperature” section, is applied to relatively coarse resolution output from the ERA-Interim reanalysis. Jarosch et al. (2012) downscaled temperature from a mesoscale model, with a better representation of topography. However, whichever resolution or type of model is used it is clear that the approach is limited by the quality of the model. Most output fields from reanalyses have biases. Biases may result for a number of reasons: coarse model resolution, necessarily simplified representation of topography, or limited observations to constrain models. A fundamental issue, relevant to the downscaling approach used here, is that the temperature lapse rate derived from the upper-air data is a free air lapse rate, which is different from gradients in near-surface air temperature. Near-surface air temperatures are not only influenced by upper-level air masses but also by advection of air masses in the boundary layer and exchange of moisture and energy between the surface and overlying air (Stull 1989). These boundary layer processes are not included in our downscaling approach.

Where meteorological stations are available, few are located above the snowline or on glacier surfaces. Snow and ice will influence overlying air masses in a way different from off-glacier and off-snow locations, especially during melting when surface fluxes are towards snow/ice surfaces (Oke 1995). The near-surface boundary layer above melting glacier surfaces often has a strong, stable inversion layer and is accompanied by a persistent katabatic wind, especially under clear conditions (Van Den Broeke 1997), which influences how the glacier micro-climate is coupled to meteorological conditions in surrounding areas. A limited survey of air temperature measured on the Morteratsch glacier, Switzerland, and air temperatures measured at nearby meteorological stations at a valley location, Samedan, and a high elevation location, Corvatsch, showed that air temperatures at the glacier station were better correlated with the high elevation station (Oerlemans 2001). This was especially true during the winter, when the valley station experienced much colder temperatures than either the glacier or the high elevation station. These conditions are likely to occur in HMA. Unfortunately, high elevation data are not generally available in this region. It is clear that care needs to be taken in developing bias correction methods.

Uncertainties in melt model

Snow cover uncertainties in MODSCAG were assessed in both Painter et al. (2009) and Rittger et al. (2012). They found root mean squared errors of 5% and 10% and recall and precision statistics above 0.9 for all regions tested. Painter et al. (2009) also found that grain size uncertainties in fully snow-covered pixels were on the order of 4% but fractionally covered pixels were not assessed.

As with all temperature-index melt modeling, uncertainties in derived degree-day factors are influenced by biases in all calibration inputs, including temperature biases; choice of daily or hourly time steps; uncertainties in snow and ice area mapping; estimates of basin rainfall and evapotranspiration; the size of elevation bands used in the aggregation from raster inputs to drainage-scale areas of SOL, SOI, and EGI; and uncertainties in measured streamflow. Lastly, the calibration optimization function is mathematically underdetermined, so multiple potential values of the snow and ice DDFs can potentially satisfy the optimization. However, we explored a large range of DDFs for snow and ice. Based on modeled melt volumes by surface type and calibration basin from the 20 best simulated annealing models, we found that uncertainty in modeled melt volume is derived largely from interannual variability in area of snow on land (Fig. 6).

Fig. 6
figure 6

Uncertainty in modeled melt volumes by surface type and calibration basin, derived from 20 best simulated annealing models, 2001–2014. Spread in modeled melt volume is derived largely from interannual variability in area of snow on land

Full size image

Snow and ice melt factors provide an indication of the sensitivity of our estimates of melt to temperature. For example, snow melt factors in the Hunza (Indus calibration basin) are between 10.41 and 10.48 mm K−1 day−1 and ice melt factors between 55.86 and 56.94 mm K−1 day−1 (Table S3). This relates to an additional or decreased melt volume of ± 10 mm day−1 and ± 55 mm day−1 with a change in temperature of ± 1 °C. The combined effect of temperature biases and melt factors could have a considerable impact on estimates of runoff. However, care should be taken with such an interpretation. The calibration procedure using observed discharge, in effect, also corrects for biases in temperature, along with biases and uncertainty in other inputs and the model itself, by tuning melt factors to give the best estimate of runoff. For example, a low melt factor for snow could compensate for both high air temperatures or an overestimate of snow covered area. This compensating effect is demonstrated by the large melt factor obtained in calibration, especially for ice, which is larger than values reported in the literature (e.g., Hock 2003).

Conclusion

We quantified the separate contributions to downstream water resources from glacier ice melt and seasonal snow melt across five major High Asia basins (Fig. 1) using a new method to separate these sources at a daily time step. Overall, glacier ice melt is an important water source for western High Asia (Syr Darya, Amu Darya, and Indus basins) primarily during the summer (June–August), but it is less significant in eastern regions (Ganges and Brahmaputra basins) that are dominated by monsoon rainfall during the summer months. Snow melt contributes approximately 10 to 20 times more water than glacier ice melt across the entire study region. Total snow melt contribution across the full basins above 2000 m is between 65 and 72% for the Syr Darya, Amu Darya, Indus, and Brahmaputra, and 43% for the Ganges.

This study’s estimates of water volumes coming from glacier ice and seasonal snow melt provide valuable insight into the specific contributions to streamflow, which allows for development of conceptual models of how water contributions may change in the future. More certain assessments of future and long-term availability of water resources in the HMA region, especially in the higher elevation regions of Central Asia, require quantifying the non-renewable portion of EGI (e.g., ice melt that is not replenished by seasonal snow accumulating on the glacier surface). This necessitates a closer analysis and modeling of glacier mass balance trends that can then be used to inform hydrologic modeling. Changes in precipitation trends (e.g., amount, rain vs. snow) also need to be considered to make accurate water resource projections.

Examining the elevation-specific contributions to runoff in association with the actual location of water users within a basin (population, agriculture, and industry) will also further the climate adaptation relevance of water scarcity studies that currently only consider the physical hydrologic processes and source waters. Socio-hydrologic research that incorporates cross-disciplinary factors will likely be well served through developing valuable in-country and cross-boundary partnerships like those established in the CHARIS project.

The research relationships built across the region during this work were critical for acquiring difficult-to-access data sets and supporting remote field work. The work has also advanced related research within the countries that are most vulnerable to the changing hydrologic regimes of High Asia. The long-term improvement of future modeling efforts depends on the vested interest of these parties, potentially leading to better data, more boots on the ground, and advocacy for a denser network of climate observations. Fundamentally, this facilitates a better understanding of the shifting relationship between the cryosphere and water resources, with local research communities contributing to data-driven management and policy decisions.