The contribution of glacier melt to streamflow

One sixth of the world's population, and one quarter of its gross domestic product, resides in areas that rely on snow or glacier melt for a majority of its water supply (Barnett et al 2005), with melt from seasonal snow packs the dominant source. Glaciers contribute substantially to water resources (Hock 2005), especially in the High Asia region, which forms the headwaters of many of that continent's largest rivers. Despite concerns about glacier extent changes and implications for water supplies (Cruz et al 2007), the global contribution of glaciers to water supply is not well known (Armstrong 2010). While some previous studies have attempted to analyze the contribution of glacier discharge to water supply, most of these studies have been site-specific. For instance, Immerzeel et al (2010) found that the glacier melt portion of streamflow in five major Southeast Asian rivers contributes significantly to the food supply of at least 60 million people. Kaser et al (2010) made estimates of seasonal glacier melt globally, however they compared these estimates with precipitation, rather than discharge, hence their estimates are not directly applicable to water supply vulnerability. Armstrong (2010) summarizes the current state of understanding of glacier contributions to water supply as, '[p]revious assessments of the glacier melt impact on surface water supply have been primarily either highly qualitative or local in scale, and in some cases, appear to be simply incorrect'. Here, we produce global estimates of the monthly maximum contribution of glacier-derived discharge to streamflow, and track the fractional contribution of glacier-derived to total river discharge through the global stream network. We focus on the monthly maximum contribution and the associated population affected.

We define glacier discharge as discharge from glaciers whose source is perennial snow, firn or ice. We include all ice caps covering less than 50 000 km² and other permanent ice globally outside of Antarctica and Greenland. Our effort is directed towards bounding the contribution of glacier discharge to water supply; we do not attempt to project future changes. Nonetheless, our results are important for understanding climate change impacts, as they provide a basis for identifying regions that are at risk to future changes in glacier-derived discharge.

Determining the contribution of glacier discharge to streamflow globally is complicated by the absence of meteorological data for all but a very small number of the world's glaciers. We therefore attempted to bound the glacier contribution to streamflow using an energy balance approach applied to known glacierized areas, based on the Global Land Ice Measurements from Space (GLIMS) (Raup et al 2000, Armstrong et al 2005) and the Digital Chart of the World (DCW) (Defense Mapping Agency 1992) datasets. Because the GLIMS project is an ongoing effort, we supplemented the GLIMS data with DCW by creating the union of the two datasets for glacier area. The merger of GLIMS (data reference year 2010) with DCW (data reference year 1992) likely overestimates total glacier area because DCW is known to include some ephemeral snowfields. We accepted this bias given our interest in estimating an upper bound to glacier-derived river discharge and used the combined dataset to calculate glacierized fractions for each model grid cell.

Total discharge (S_T) and an upper bound for glacial discharge (S_Gu) were estimated in the following manner. To estimate S_T, we implemented the Variable Infiltration Capacity (VIC) land surface hydrological model (Liang et al 1994) at a spatial resolution of one-quarter degree. The VIC model is a semi-distributed, physically based land surface model, which solves the water and land surface energy balance for each individual model grid cell. Subgrid topographical features are represented through the use of multiple elevation zones within each model grid cell (up to five in this case). The VIC snow model represents the snowpack as two layers for purposes of thermal computations and uses an energy balance approach that accounts for refreezing of melt water in the pack, the role of vegetation where present and a snow age dependent albedo (Andreadis et al 2009). The model has been widely applied to water and energy budget studies, streamflow and drought forecasting, and climate change assessments (Maurer et al 2001, Nijssen et al 2001a, Christensen and Lettenmaier 2007, Koster et al 2010).

Daily precipitation, minimum and maximum temperatures, and wind speed are the primary forcings for the VIC model. Other forcings, such as downward solar and longwave radiation and humidity, are estimated by the model from these variables using the Thornton and Running (1999) and Tennessee Valley Authority (1970) algorithms for downward solar and longwave radiation, respectively, as implemented by Maurer et al (2002). Precipitation was taken from the Tropical Rainfall Measurement Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) for the 50°S–50°N latitude band (Huffman et al 2007). Outside of this band, precipitation was taken from Sheffield et al (2006). Temperature and wind data were taken from the National Centers for Environmental Prediction (NCEP) Reanalysis 2 dataset (Kanamitsu et al 2002). All data were interpolated to a one-quarter degree grid. As part of the interpolation, the daily temperature minima and maxima were adjusted using a pseudo-adiabatic lapse rate of 6.5 °C/1000 m based on elevations from the Global Land One-kilometer Base Elevation Project (GLOBE) digital elevation model (DEM) (GLOBE Task Team et al 1999). Global soil parameters were derived as described in Nijssen et al (2001b). Soil depths and baseflow parameters were taken from Nijssen et al (2001a). Global vegetation parameters were derived as in Su et al (2005).

The VIC model was initially run at a daily time step for a total of 18 yr to initialize the model's soil water storages, and to create a snow 'reservoir' in the areas defined as glaciers. The modified snowpack was located in the highest elevation band (mountain glaciers) unless the elevation of the lowest band was less than 200 m (valley or tidewater glaciers). For each glacierized grid cell the highest or lowest elevation band area was adjusted to match the areal coverage of the glacier area determined with GLIMS and DCW. Following the 18 yr spin-up, model simulations were analyzed for the period from January 1998 to December 2006 at a three-hourly time step.

To distinguish glacier melt from the calculation of seasonal snowmelt and other sources of runoff, which are included in our estimate of S_T, we separately estimated the maximum possible amount of melt (S_Gu). This calculation was performed over the fraction of each grid cell that was glacier covered according to the GLIMS and DCW sources and used the energy balance terms calculated by the VIC model over the glacierized fraction of each grid cell. For these areas, we assumed that all available energy was converted to melt. The energy available for melt in our calculation of S_Gu is the sum of net radiation and turbulent heat exchange. Compared to the other components, advected heat is generally small (Lachapelle 1959) and was ignored. Downward shortwave and longwave radiation, as well as latent and sensible heat fluxes were taken from the VIC simulations. We estimated emitted longwave radiation using an emissivity of 1.0 and an assumed surface temperature of 0 °C. Reflected shortwave radiation and hence net radiation is particularly sensitive to albedo, which can vary widely for glaciers and is highly variable over small areas and from year to year (Barnett et al 2005, Hock 2005, Cruz et al 2007). Albedo is highest for newly fallen snow and decreases as snow melts or is metamorphosed into glacier ice. For the energy balance model (S_Gu), we used a single effective albedo, an assumption that is common in snow models (Andreadis et al 2009). We compiled 100 glacier ice albedo estimates from 14 publications, which ranged from 0.03 to 0.85 with a mean of 0.39, a standard deviation of 0.20, and a median of 0.35 (see supplementary materials available at stacks.iop.org/ERL/7/034029/mmedia). The values fall into three categories: those stated in the text or in a table in published articles or reports, those derived from figures, and those averaged from a stated minimum and maximum. In the case of continuous time series of albedo values, the average value was taken. The highest values are from glaciers covered partially or fully with recently fallen snow; the lowest values are from glaciers with substantial debris, black carbon and/or biomass growth. Given our interest in estimating a maximum contribution, we used the lower 25th percentile of the 100 values (0.24). We chose not to use the lowest value, which likely is associated with debris covered glaciers, for which other thermal factors such as insulation by the debris cover offset the effects of a lower albedo on melt.

The upstream area for each location in the river network was derived from the flow direction network of Wu et al (2011). Spatially aggregated, monthly runoff (S_T) and monthly glacier melt (S_Gu) were calculated for each location by summing the monthly grid cell estimates over the contributing area. All glacier runoff was assumed to leave the glacier immediately as river discharge. The contribution of glacier discharge to total river discharge was then calculated as S_Gu/S_T. From these estimates, we identified the maximum monthly contribution, which generally occurs in summer when seasonal snowmelt from the non-glacierized part of a river basin is low, glacier melt is high, and other non-glacier sources of runoff are low (aside perhaps from tropical glaciers). We excluded the winter months from the analysis for North America, because low flows during the winter for glacierized areas resulted in spuriously high ratios in some cases.

The model domain (figure 1) consisted of all river basins with glaciers in their headwaters and for which river discharge could therefore be affected by glacier melt. This area was determined based on the glacier coverage and the flow network. To assess populations potentially at risk to changes in glacier discharge, we overlaid the Gridded Population of the World 2010 dataset (CIESIN 2005) on areas for which the maximum monthly contribution of glacier melt exceeded a certain threshold.

Zoom In Zoom Out Reset image size

Figure 2 shows the areas for which the maximum monthly fraction of glacier-derived discharge was greater than 5%, 10%, 25% and 50%. The largest inferred fractions are for mid to late summer. The largest areas identified are streams with headwaters in the Himalayan Mountains, coastal Alaska, the southern Andes, Iceland and the Alps. A number of smaller areas were also identified. Table 1 shows the affected populations by region. The results show that for the 5% threshold, the population potentially affected is about 600 million (8.9% of global population). For 10%, 25% and 50% thresholds the estimated populations are about 370 million (5.4%), 140 million (2.1%) and 120 million (1.8%), respectively.

Zoom In Zoom Out Reset image size

Two points stand out from the results. First, while the total domain for which there is any inferred signature of glacier discharge is quite large (46 200 000 km² or about 35% of total global land area exclusive of Greenland and Antarctica), the area for which the glacier contribution exceeds even the lowest threshold (5%) for the maximum melt month is much smaller (8 190 000 km² or about 6.2% of global land area). Second, the largest populations at risk are in Asia. The Indus stands out as a major river basin in which a large glacier contribution during part of the year combines with a high population density. This basin is directly responsible for most of the population above the 25% and 50% thresholds. Outside of Asia, both contributing area and affected population are small for maximum monthly thresholds greater than 25%.

To evaluate our procedure for bounding glacier-derived river discharge, we compared results of our glacier discharge estimation procedure to published values for the Gulkana and Wolverine glaciers, two United States Geological Survey (USGS) Benchmark glaciers in Alaska (Fountain et al 1997). Both glaciers have long-term records of mass balance, precipitation and streamflow near the glacier terminus. Wolverine Glacier is located in a maritime environment, whereas Gulkana Glacier is in a continental environment (Bitz and Battisti 1999). Gulkana Glacier covers 19.6 km² in a drainage basin of 31.6 km² (stream gauge location about one km downstream of the glacier terminus). Wolverine Glacier is similar in size (16.8 km²) and lies in a 24.6 km² drainage basin, the gauge for which is about 150 m downstream of its terminus.

From 1998 to 2006, the month of maximum discharge for the Gulkana glacier was July, with an average recorded discharge of 10.4 m³ s⁻¹ (880 mm/month). Precipitation during this month was 113 mm. Since part of this precipitation contributes to evapotranspiration, almost all of the July discharge originates from melt. The observed discharge compares well with our predicted July average glacier discharge of 9.7 m³ s⁻¹. Maximum potential discharge occurs in June for Wolverine Glacier with values of 6.3 m³ s⁻¹ (660 mm/month) (recorded) and 5.2 m³ s⁻¹ (predicted). Precipitation during this month was 36 mm. Given our simplifying assumptions, this agreement is encouraging, because at the small scale our estimates compare well with the observations for these two glaciers.

Table 2 contains further comparisons between our estimates of the fractional contribution of glacier discharge and other published values. Where published values are given for river basins that constitute more than one grid cell, the grid cell corresponding to the basin outlet was used. As a result, comparisons are most relevant for point references. The references that use water balance approaches derive glacier contribution indirectly by subtracting estimated water balance components from observed streamflow. Table 2 includes direct comparisons of like quantities depending on what was reported in the literature (for example, annual estimates or estimates for a particular month) as well as our estimated maximum monthly glacier contribution fraction and the month in which it occurred.

Taken as a group our estimates of glacier-derived discharge from table 2 are lower than published values, and the differences are larger than for the well-instrumented Gulkana and Wolverine Glaciers. Among the possible reasons are that some of the published values are known to include seasonal snowmelt (these values are noted). For some reported values the location of the measurements is poorly defined in the publication and our model simulations are likely to include a larger area given the resolution of our model. In addition, the reported values show considerable interannual variability, which is not reflected in our multi-year average. For example, the annual glacier contribution fraction for the Rio Santa ranges from 12%–20% to 40% depending on the averaging period (see table 2). When comparing the monthly contribution for specific months (bottom half of table 2), our overall maximum often compares better with the previously reported value than the value for the same month, reflecting possible shortcomings in our method for aggregating runoff over an upstream area without regard for travel times within the glacier and in the channel, although this effect is likely small since we focus on monthly averages. Because we do not represent storage of liquid water within the glacier, modeled glacier discharge may enter the stream network too soon (Stenborg 1970, Tangborn et al 1975). This could put glacier runoff out of phase with non-glacier runoff simulated by the hydrology model. Also, as glacier discharge moves downstream it is compared to total streamflow that may be subject to diversions by water engineering works and agriculture (Kaser et al 2010); these effects are not represented in the hydrology model.

Our energy balance approach does not account for all of the complex processes occurring within and upon glaciers and their drainage basins. In our model for estimating glacier melt, all glaciers are uniform, flat slabs with full exposure (no terrain blockage) to solar radiation and a fixed, relatively low albedo; hence for south facing glaciers (northern hemisphere) we likely underestimate net radiation, and hence discharge. However, these effects are likely greatest for small headwater glaciers, and should average out for larger and/or multiple glaciers. The assumption of a glacier surface (and implicitly, subsurface) temperature of 0 °C during the melt period generally will result in an overestimation of discharge (Pellicciotti et al 2009) because refreezing within the glacier is not represented (Pfeffer et al 1998). Complex local glacierized environments may introduce errors in turbulent heat flux estimates. However, apart from maritime glaciers, net radiation is the main contributor to melt in both high-energy conditions (clear and warm) and low-energy conditions (cloudy and cool) for the vast majority of glaciers (Brock et al 2000, Datt et al 2008), and the contribution of turbulent fluxes to the overall glacier energy balance is generally small, especially since latent and sensible heat are usually of opposite sign. Advected energy is not included in our energy balance model but in most cases is small (Lachapelle 1959).

Biases in model-simulated streamflow are another source of error. Model-simulated streamflow generally agrees well with observations when the model forcings (especially precipitation) are well known (Maurer et al 2002). The direction of simulation errors varies from watershed to watershed, and can be considered essentially random on a global basis.

The contribution of non-seasonal (i.e., permanent, or at least long-term with respect to our observation period) loss of glacier ice is another potential source of error. Our method computes total discharge, whether or not the glacier is in balance. However, we do assume that the glacier area is fixed through the period of record. For most glaciers however, changes in area over the 9 yr analysis period are small relative to the glacier area at the beginning of the period. For receding glaciers, our approach will bias our estimates upward somewhat.

As the climate warms, our estimates of glacier contribution fractions may increase as mass loss accelerates, but ultimately will be more than canceled by loss of glacier area. It is important to note, however, that reduction in the fraction of glacier discharge does not imply reduction of streamflow, which depends on characteristics of precipitation and snow accumulation and ablation on non-glacier fed portions of the river basins in question. Although not explicitly examined in our analysis, glacier-derived runoff generally is less variable on an interannual basis than non-glacier-derived runoff, so regardless of the direction and magnitude of changes in total runoff that might accompany glacier recession, variability is likely to increase.

Our analysis is intended to provide, at least approximately, a worldwide upper-bound estimate of the glacier discharge contribution to river discharge and populations potentially at risk to glacier retreat. For all glacier melt contribution thresholds that were evaluated, more than 90% of the at risk population lives in Asia. We find that no more than 8.9% of the global population lives in river basins that depend on seasonal glacier discharge for at least 5% of river discharge in the peak-melt month; no more than 5.4% rely on 10% of glacier discharge in the peak-melt month, no more than 2.1% rely on 25% and less than 1.8% rely on 50%. In general, our estimates of the glacier discharge contribution to streamflow are lower than previously published values. While our estimates of the glacier melt contribution to streamflow are no substitution for detailed studies of glacier dynamics in any specific geographic setting, in aggregate, they provide a methodologically consistent global identification of regions in which glacier melt contributes significantly to streamflow during at least part of the year. While the number of at risk people at each of the threshold levels identified in table 1 is subject to considerable uncertainty, there is much less doubt regarding the relative importance of glacier melt to total streamflow in the areas identified in figure 2. As such, our results identify regions and populations that are at risk to future changes in glacier-derived discharge.

The research on which this paper is based was supported in part by Grant No. NNX10AP90G from the National Aeronautics and Space Administration to the University of Washington. The data reported in this paper are archived at ftp://ftp.hydro.washington.edu/pub/HYDRO/data/published/Schaneretal_2012_Global/.