Quantifying the contribution of glacier runoff to streamflow in the upper Columbia River Basin, Canada

Glacier melt provides important contributions to streamflow in many mountainous regions. Hydrologic model calibration in glacier-fed catchments is difficult because errors in modelling snow accumulation can be offset by compensating errors in glacier melt. This problem is particularly severe in catchments with modest glacier cover, where goodness-of-fit statistics such as the Nash-Sutcliffe model efficiency may not be highly sensitive to the streamflow variance associated with glacier melt. While glacier mass balance measurements can be used to aid model calibration, they are absent for most catchments. We introduce the use of glacier volume change determined from repeated glacier mapping in a guided GLUE (generalized likelihood uncertainty estimation) procedure to calibrate a hydrologic model. This approach is applied to the Mica basin in the Canadian portion of the Columbia River Basin using the HBV-EC hydrologic model. Use of glacier volume change in the calibration procedure effectively reduced parameter uncertainty and helped to ensure that the model was accurately predicting glacier mass balance as well as streamflow. The seasonal and interannual variations in glacier melt contributions were assessed by running the calibrated model with historic glacier cover and also after converting all glacierized areas to alpine land cover in the model setup. Sensitivity of modelled streamflow to historic changes in glacier cover and to projected glacier changes for a climate warming scenario was assessed by comparing simulations using static glacier cover to simulations that accommodated dynamic changes in glacier area. Although glaciers in the Mica basin only cover 5 % of the watershed, glacier ice melt contributes up to 25 % and 35 % of streamflow in August and September, respectively. The mean annual contribution of ice melt to total streamflow varied between 3 and 9 % and averaged 6 %. Glacier ice melt is particularly important during warm, dry summers following winters with low snow accumulation and early snowpack depletion. Although the sensitivity of streamflow to historic glacier area changes is small and within parameter uncertainties, our results suggest that glacier area changes have to be accounted for in future projections of late summer streamflow. Our approach provides an effective and widely applicable method to calibrate hydrologic models in glacier fed catchments, as well as to quantify the magnitude and timing of glacier melt contributions to streamflow.


Introduction
In many mountainous regions, glacier melt makes significant contributions to streamflow, particularly in late summer during periods of warm, dry weather (Koboltschnig et al., 2008;Stahl and Moore, 2006;Verbunt et al., 2003;Zappa and Kan, 2007).Understanding the quantity and timing of these contributions is important for a range of purposes, including short-term and seasonal forecasting of reservoir inflows and long-term projections of the potential hydrologic effects of climate change.This knowledge is particularly critical given that these contributions are likely to decrease in the medium to longer term as glaciers retreat (Gurtz et al., 2003;Koboltschnig et al., 2008;Marshall et al., 2011;Stahl et al., 2008), with implications for both water resources management and aquatic ecology (Moore et al., 2009;Zappa and Kan, 2007).
Glacier contributions to streamflow have been reported in different ways.Huss (2011) reported the "melt from glacierized areas", which by definition also includes the snow melt component, while Stahl et al. (2008) reported only the glacier ice melt component as the relevant contribution of glaciers to streamflow because that is the component that diminishes as a direct result of glacier retreat.In catchments where glacier mass balance and snowline observations exist, a water balance approach can be used to estimate glacier ice melt contributions to streamflow (Schär et al., 2004;Young, 1982).Alternatively, empirical analysis of the contrasting responses of glacier-fed and unglacierized catchments can provide insight (Stahl and Moore, 2006;Schaefli and Huss, 2011).Deterministic hydrologic models can also be used to quantify glacier melt contributions to streamflow (Koboltschnig et al., 2008;Schaefli and Huss, 2011;Stahl et al., 2008).The use of models to quantify glacier melt contributions to streamflow requires adequate representation and parameterization of glacier ice melt processes.However, in cases where streamflow observations are the only available data for model calibration, an incorrect simulation of glacier ice melt can be offset by compensating errors in the simulation of snow accumulation and snow melt (Konz and Seibert, 2010;Schaefli and Huss, 2011;Stahl et al., 2008), resulting in "equifinality" -i.e. the existence of multiple parameter sets that provide adequate streamflow simulations despite differences in predictions of snow and ice processes.Equifinality introduces substantial uncertainty into model-based estimates of glacier melt contributions to streamflow.Problems associated with equifinality can be reduced by constraining a model with additional information.Previous studies that quantified the contribution of glacier melt to streamflow reduced equifinality by incorporating glacier mass balance data or equilibrium line altitudes, in addition to streamflow data (Konz and Seibert, 2010;Moore, 1993;Schaefli et al., 2005;Stahl et al., 2008).Unfortunately, glacier mass balance observations are sparse and typically unavailable for most catchments.
Most modelling studies that focused on glacier melt contributions to streamflow examined catchments with substantial glacier cover, typically in excess of 10 % of the catchment area (Koboltschnig et al., 2008;Schaefli and Huss, 2011;Stahl et al., 2008).However, Stahl and Moore (2006) found that the effects of glacier cover on late-summer streamflow can be detected in catchments with as little as 2 to 5 % glacier coverage.Equifinality may be especially problematic in large catchments with modest glacier cover (less than 10 %) given the relatively small variance in streamflow associated with glacier melt contributions.In a recent study focused on macro-scale catchments in Europe with low or modest glacier coverage, Huss (2011) stated that "[m]ass balance data for 50 glaciers in the Swiss Alps . . .[were] central to this study".Such extensive glacier mass balance data sets are uncommon outside Europe, limiting the geographic transferability of the approach to other larger catchments outside Europe.
Another challenge in modelling glacier melt contributions to streamflow is that they can be influenced by changes in glacier cover.In catchments with high glacier cover, several studies demonstrated that projected glacier changes over the next few decades need to be accounted for in hydrologic modelling to avoid biased predictions (Gurtz et al., 2003;Koboltschnig et al., 2008;Stahl et al., 2008).However, it is not clear from the current literature if accounting for future glacier changes is also necessary in large catchments with modest glacier cover.It is also uncertain how sensitive hydrologic simulations are to historic changes in glacier cover, particularly over recent decades.
The objective of this study was to develop an approach for estimating the magnitude and timing of glacier melt contributions to streamflow in large catchments with modest glacier cover and no mass balance observations, along with an assessment of uncertainty.This study used glacier volume and area changes to assist in calibration, which were derived from analyzing sequential digital elevation models (DEM) and maps of glacier cover.In addition, we address the sensitivity of modelled streamflow to historic changes in glacier cover, and to projected glacier changes for a typical climate warming scenario.

Study area
The study focused on Mica basin, a major tributary to the Canadian portion of the Columbia River Basin in British Columbia (Fig. 1).Mica basin has a drainage area of 20 742 km 2 , with elevation ranging from 579 m above sea

Data
Data from five climate stations within or just outside Mica basin were available for modelling (Fig. 1).FLK (1985FLK ( -1993)).Errors associated with air temperature measurements, instrumental error and the possibility of bias associated with site characteristics are believed to be small relative to spatial and temporal variability in Mica basin.
Streamflow data used in this study are daily inflows to Kinbasket Reservoir, computed by BC Hydro from a water balance based on the rates of release through the dam and changes in water level.Although evaporation from the reservoir is not included in the computed inflows, estimates based on reservoir area and potential evaporation indicate it should not exceed about 1 % of inflow.BC Hydro calculates an estimate of the average daily inflow error as the root mean squared error (RMSE) between quality-controlled and raw inflow data for all their reservoirs (F.Weber, BC Hydro, personal communication, 2011).The RMSE for daily Mica inflow observations is approximately 10 %.For monthly and annual data, where daily random errors are likely to cancel out, RMSE was not calculated but should be smaller than the daily error and primarily reflect the presence of systematic errors.
Snow water equivalent (SWE) data for three snow pillows located at the FLK, MOL, and RGR climate stations (Fig. 1) were available from 1995 onwards.Glacier coverages were derived from Landsat Thematic Mapper scenes for 2005 and 2000 and from high altitude aerial photography for 1985 (Bolch et al., 2010).Glacier volume loss was calculated from digital elevation models (DEM) derived from the 1999 Shuttle Radar Topography Mission (SRTM) and from aerial photographs taken between 1982 and 1988, which have a median weighted date for Mica basin of 1985 (Schiefer et al., 2007).The estimated ice volume loss from 1985-1999 was 7.75 km 3 .Taking mapping uncertainty into account, ice volume loss from 1985-1999 lies between 6 and 9 km 3 .
The geodetic estimate of the rate of thickness change for the Columbia region of British Columbia was −0.53 m yr −1 (Schiefer et al., 2007), while the rate specifically for Mica basin was −0.43 m yr −1 .This rate of thinning is slightly less than that indicated by in situ measurements of mass balance at Peyto Glacier, Alberta, which averaged approximately −0.6 m yr −1 between 1966 and 2005 (Statistics Canada, last access: 21 October 2011).Peyto Glacier is located to the east of Mica basin in the Rocky Mountains, which receives less snowfall than Mica basin.Therefore, the difference between the rates of mass loss for Mica basin and Peyto Glacier is consistent with the differences in the climatic settings.

The HBV-EC hydrologic model
The HBV-EC model is a Canadian variant of the HBV-96 model (Lindstrom et al., 1997).It has been incorporated into the EnSim Hydrologic modelling environment (now known as Green Kenue) (Canadian Hydraulics Centre, 2010).The ability of HBV-EC to provide accurate predictions of streamflow in British Columbia's mountain catchments was demonstrated in an intercomparison study of watershed models for operational river forecasting (Cunderlik et al., 2010;Fleming et al., 2010).The model algorithms have been described in detail by Hamilton et al. (2000), Canadian Hydraulics Centre (2010) and Stahl et al. (2008).Key features are presented below.
To minimize computational effort, HBV-EC is based on the concept of grouped response units (GRUs), which contain grid cells having similar elevation, aspect, slope, and land cover.HBV-EC has the capability to model four landcover types: open, forest, glacier and water.To represent lateral climate gradients, HBV-EC allows for subdividing a basin into different climate zones, each of which is associated with a single climate station and a unique parameter set.Water draining from non-glacier GRUs is routed through two lumped reservoirs representing "fast" and "slow" responses.To predict the discharge for a given time step, HBV-EC sums output from the two non-glacier reservoirs and the reservoirs associated with glacier GRUs (see below).
The temperature-index-based snow melt algorithm from HBV-96 was adapted by Hamilton et al. (2000) to account for the effects of slope, s, aspect, a, and forest cover.In HBV-EC, daily snowmelt (M) (mm day −1 ) is calculated from daily mean air temperature (T air ) ( • C) as follows: where C 0 is a base melt factor (mm day −1 • C −1 ) that varies sinusoidally between a minimum value (C min ) at the winter solstice to a maximum value at the summer solstice (C min + C) to account for seasonal variations in solar radiation, and C is the increase in melt factor between winter and summer solstices.The melt ratio for forests (MRF) ranges between 0 and 1 and reduces melt rates under forests compared to melt at open sites.The coefficient AM controls the sensitivity of melt rates to slope s and aspect a and thus mimics the effects of spatial variations in solar radiation.For glacier GRUs, melt is computed as for an open site (MRF = 1) until the previous winter's snow accumulation has ablated.At that point, glacier melt is computed by multiplying open site melt by the coefficient MRG, which typically ranges between 1 and 2, to represent the reduction in surface albedo.
Storage and drainage of meltwater and rain for each glacier GRU are modelled using linear reservoirs.The outflow coefficient (KG) for each GRU depends on snow depth, ranging from a low value (KG min ) when the GRU has deep snow cover to a maximum value (KG min + KG) when the GRU is snow-free (Stahl et al., 2008).This representation accounts for seasonal changes in the efficiency of the glacier drainage system.Glaciers in HBV-EC cannot vary in area or volume during a model run without stopping and restarting the simulation.The net mass balance for each GRU is calculated from time series of SWE and glacier ice melt for each glacier GRU.The total mass balance for the Mica basin is calculated from area-weighted net mass balances from each elevation band.For more details see Stahl et al. (2008).

Calibration and testing
The model was calibrated for the period 1985-1999, the same period for which the glacier volume loss was calculated.Calibration runs were split into two time periods, each with a five-year spin-up period to ensure that storages in the model, in particular the slow reservoir storage, equilibrated with the forcing data.Simulations for the first period, 1985-1992, used the 1985 glacier coverage, while the second period, 1992-1999, was based on glacier coverage from 2000.The updating of glacier area during the model calibration was done to be consistent with the updating in a long term simulation to assess the importance of glacier area updating described in Sect.2.7, below.Glacier net mass balance (b n ) for the entire basin was derived from net mass balances for each GRU and compared to geodetically calculated glacier volume loss.
The period 2000-2007, with glacier cover based on data from 2005, was used as an independent test period.Model predictions were compared to observed streamflow data and SWE data from the three snow pillow sites (FLK, MOL, RGR, in Fig. 1).Although HBV-EC does not explicitly represent temporal variability in precipitation and temperature lapse rates, our model setup does account for some degree of seasonal variability in vertical climatic gradients by delineating Mica basin into five climate zones -partly based on elevation -and forcing each climate zone with a separate climate station.Since HBV-EC does not predict state variables for a specific location but only for each GRU, observed SWE data were compared with simulated SWE from the GRUs in which the snow pillows are located.The SWE data were not used in model calibration, and thus represent an independent test of the model.

A "guided" GLUE approach to address parameter uncertainty
A common approach to address uncertainty in model predictions is to generate random samples from the usually highdimensional parameter space and subsequently to pick the best performing parameter sets according to one or multiple criteria (e.g.Konz and Seibert, 2010;Stahl et al., 2008, for glacier related applications).However, in a high-dimensional parameter space, random sampling with even thousands of model runs does not guarantee that the "best" parameter combinations are found.Without prior knowledge of how well the "best" possible solution performs, the modeller will usually relax criteria in order to obtain enough acceptable parameter sets, with the possible result that criteria for acceptable parameter sets are more relaxed than necessary.Particularly in large catchments with moderate glacier cover, this approach could result in high uncertainties for glacier ice melt estimates.To ensure that the final ensemble parameter set contains solutions that perform similarly to the "best" possible solution(s) within a parameter space, we modified the Generalized Likelihood Uncertainty Estimation (GLUE) approach outlined in Beven and Freer (2001) and Freer et al. (1996) to an approach that can be described as a "guided" GLUE approach (Fig. 2).
Our calibration procedure starts with finding a benchmark parameter set by maximizing the Nash-Sutcliffe efficiency (E) or, in terms of GLUE, the generalized likelihood measure.This was done with genoud (Mebane and Sekhon, 2011), an optimization algorithm in R (R Development Core Team, 2011) that combines evolutionary algorithm methods with a steepest gradient descent algorithm.A large negative number was returned for parameter sets that did not fulfill the multiple criteria listed in Fig. 2 to ensure that the optimization algorithms not only maximized E but also searched for solutions that met the additional criteria.In a second step, a Latin Hypercube Search (LHS) with 10 000 model runs was performed.Latin hypercube designs are most often used in high-dimensional problems, where it is important to sample efficiently from distributions of input parameters.Parameter sets from the 10 000 model runs were constrained by criteria given in Fig. 2. If no parameter sets with Nash-Sutcliffe efficiencies greater than the benchmark efficiencies minus a threshold were found, all parameter sets were rejected.There are two ways to proceed when no parameter sets are found by the LHS: either increase the sample size or adjust the prior parameter distributions (decrease the ranges).Increasing the sample size is the favored solution because it should lead to a more diverse set of parameters.However, the number of model runs is limited by computational power (even with multiple CPUs it would take weeks for Mica basin).Given time constraints, we chose to adjust the prior parameter distributions.With adjusted (narrowed) parameter ranges, the LHS was repeated until enough parameter sets (∼20-30) were found that fulfilled all criteria (i.e."behavioral" parameter sets in GLUE terminology).The parameter ranges for model calibration and uncertainty analysis were based on default values provided in the HBV-EC manual (Canadian Hydraulics Centre, 2010), values reported in previous studies (Hamilton et al., 2000;Stahl et al., 2008), the authors' experience with applying HBV-EC on other catchments, and by visually testing the influence of parameters on the simulated hydrograph.With the modest glacier cover in Mica, a visual inspection of simulated hydrographs provided more information on the sensitivity of modelled streamflow to the various glacier parameters than a single goodness of fit measure such as E. Prior parameter distributions for LHS were assumed uniform at all stages.

Assessing the sensitivity of streamflow to glacier area changes
To assess the sensitivity of streamflow simulations to historic glacier area changes, we compared simulations for the latter part of the test period using the earliest glacier coverage (1985) to contrast with the simulations based on the 2005 coverage (which was used for the latter part of the test period).We also assessed the sensitivity of simulated streamflow to projected changes in glacier area under a typical climate warming scenario by comparing ensemble simulations using a static (observed) glacier cover from 2005 throughout the 100 yr simulation period to simulations with a dynamic glacier cover.Forcings for HBV-EC were based on simulations from one global climate model (GCM), the Canadian CGCM3.1-T47, with the A1B emissions scenario, which represents a mid-range warming scenario (Nakicenovic et al., 2000).Daily output from the GCM was downscaled for input to HBV-EC using the TreeGen algorithm (Stahl et al., 2008).Projected changes in glacier cover under the A1B emission scenario, also derived from CGCM3.1-T47, were simulated using the UBC Regional Glaciation Model, a physically based, spatially distributed model of glacier dynamics (G.Clarke, The University of British Columbia, personal communication, 2011; Clarke et al., 2011).At the time of writing, publications describing the Regional Glaciation Model and its application to the Columbia Mountains are in preparation, so we cannot provide details of model output here.The objective here was not to present comprehensive projections of future changes to inflow, but rather to help identify if and when glacier area updating is of relevance in a large basin with modest glacier cover; hypothetical glacier decreases could have been used as an alternative.Model output indicates that glacier area will decrease by 83 % of the glacier extent in 2000 by the end of this century.

Modelling the contributions of glacier ice melt to streamflow
An estimate of the contribution of glacier ice melt to discharge and the associated uncertainty was calculated as the difference between streamflow simulations with and without glaciers for each ensemble member.In the no-glacier runs, all glacier cover was converted to open land cover to account for the fact that snowmelt and rainfall runoff from the areas currently covered by glaciers would occur even if the glaciers completely disappeared.

Model calibration and uncertainty analysis
The benchmark parameter set obtained by the combined evolutionary-steepest gradient optimization matched observed streamflow data with E of 0.93 for the calibration period (1985-1999) and 0.95 for the test period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007).A first 10 000 run LHS within the initial parameter ranges (parameter range step 1 in Fig. 2) found no acceptable parameter sets that met all criteria.Although 28 parameter sets had E > 0.91, all of these parameter sets were rejected because none fulfilled all of the additional criteria.In the absence of prior knowledge of the benchmark E, the common procedure would now have been to find acceptable solutions with relaxed criteria.However, the benchmark parameter set indicates that there are better performing solutions within the initial parameter space; 10 000 runs are too few to sample the parameter space for acceptable solutions.A second LHS with adjusted parameter ranges found 17 acceptable parameter sets, but histograms indicated that two parameters in the acceptable parameter sets were predominantly sampled near a range boundary and therefore a third LHS with slightly refined parameter ranges was performed.From the 10 000 model runs in the third LHS, 705 parameter sets met the final criterion of E > 0.92, but only 23 of these also met the additional criteria (Fig. 2).
The calibrated parameters with the highest correlations to glacier net mass balance (b n ) and E are temperature lapse rate (T lapse ), melt factor at winter solstice (C min ), and precipitation lapse rate (P lapse ) (Table 1).Other important parameters are the ratio between melt rates for glacier ice and seasonal snow (MRG) and the increase of melt factor between winter and summer solstice ( C). MRG and C are both correlated with b n at all steps during the uncertainty analysis, but are correlated with E only in the first LHS with wide parameter ranges.The routing parameters with the highest correlation with E are the fast reservoir release coefficient (K F ) and the fraction of runoff directed to the slow reservoir (FRAC).The exponent to adjust the linearity of the release rate of the fast reservoir, α, has little influence on E. Glacier reservoir coefficients and the melt ratio for forest (MRF) show weak correlation to both glacier volume change and E.  Step 1 Step 2 Step 2 Step 3 Step A wide range of modelled glacier volume changes can lead to values of E close to the benchmark (Fig. 3).Results from the first LHS suggest that equifinal parameter solutions are possible with glacier volume losses ranging from 5 to 40 km 3 .This point underlines the advantage of using observed glacier volume changes to constrain model parameters, particularly in large basins with modest glacier coverage like Mica.Note that the second LHS gives higher maximum values of E because of the greater sampling density within the restricted parameter space and not necessarily because the glacier volume loss is close to the observed.More intense sampling within the parameter space that leads to higher glacier volume losses could possibly have led to higher E at higher glacier volume losses as well.Combinations that lead to increases in glacier volume yielded a substantial decrease in E.

Model testing
Model testing on streamflow for the period 2000-2007, using the observed glacier extents from 2005, yielded an efficiency of 0.95 for the best model (Fig. 4a), a slightly better performance than during the 1985-1999 calibration period (E = 0.93).The goodness of fit of the best parameter set derived by the Latin hypercube search is essentially the same as the fit obtained by the combined evolutionary-steepest gradient optimization.
All 23 behavioral parameter sets reproduce the seasonal peak flows as well as low flows, but have difficulty with modelling intense rainfall events, especially during autumn (Fig. 4b).This is not surprising, since one of the two reservoirs (the slow reservoir) is primarily used to model the low flows during winter, and the single fast reservoir cannot simultaneously represent runoff generation due to melt and rainfall given the differences in their spatial patterns and nonlinearity.Since this model weakness only appears to affect rainfall-generated daily peak flows, it should not detract from the estimation of glacier melt contributions to streamflow, especially over monthly or longer time scales.
Despite the difference in spatial scales associated with modelled and observed SWE, SWE predicted by HBV-EC shows reasonable agreement with observations, with linear regressions between predicted and observed having R 2 of 0.82, 0.76, and 0.86 for the Molson Creek, Floe Lake, and Mount Revelstoke snowpillows, respectively (Fig. 5).The model accurately predicts the timing of the onset of snowmelt as well as the rate of decrease of SWE during the ablation period at all three snow pillow sites.However, the model tends to underestimate peak SWE.For some years this underprediction is within the expected error of SWE measurements (5 % according to Gray and Male, 1981).Snow pillows tend to overestimate SWE due to snow creep, which puts additional load on the pillows (Gray and Male, 1981).However, there are some station-years in which the underestimation is too large to be simply attributed to measurement errors in the snow pillows (e.g.1996-1997 at Floe Lake).The timing of these types of errors is not consistent among stations.For example, in the water year 1996-1997, peak SWE was reproduced reasonably accurately at Molson Creek and Mount Revelstoke, but strongly underpredicted at Floe  Lake.Inconsistent variations in gauge catch efficiency could explain at least some of this underprediction.This inconsistent pattern of errors could also partly reflect the inherent variability in precipitation patterns from year to year, which are probably not fully represented through the use of fixed vertical gradients in each climate zone.

Historic contributions of glacier ice melt to streamflow
The mean annual contribution of ice melt to total streamflow varied between 3 and 9 % and averaged 6 % (Fig. 6).Trend analysis revealed no significant increase or decrease of the annual contributions of glacier ice melt with time.For annual, August and September flows, the uncertainty bounds between runs with and without glaciers are close but do not overlap for most years (Figs.6 to 8).Contributions of glacier ice melt to discharge at Mica dam dominantly occur in August and September (Figs.7  and 8).Mean August streamflow, calculated from the ensemble mean, would be up to 25 % lower if there were no glaciers, although the interannual variation of contributions is relatively high (standard deviation = 7 % of the simulated mean flow with glaciers) (Fig. 7).The relative contribution of glaciers is highest in September, when ice melt can provide up to 35 % of the discharge.September contributions of ice melt are also less variable over time, with a standard deviation of 5 % of the simulated mean flow with glaciers (Fig. 8).
Figure 9 presents the mean and range of ensemble predictions for simulations with glaciers to simulations where glaciers have been removed and replaced by open land cover for two years with contrasting hydroclimatic conditions.Glacier runoff is particularly important in years with early snowmelt such as 1998, the year with the highest modelled ice melt (Fig. 9), when glaciers can contribute to more than 20 % of the flow for periods of more than two months.In years with late snow melt, such as the year 2000, glaciers have a minor effect on discharge.In July, some years have a higher discharge in the no-glacier scenario because the glacier routing routine stores water in the early part of the melt season and releases it later.The glacier reservoir can lag flows from a few days up to several weeks, depending on the parameter values.This type of seasonal storage effect has been documented in previous studies of glacier hydrology (e.g.Stenborg, 1970).

Sensitivity of streamflow to glacier area changes
Figure 10 presents the sensitivity of streamflow to historical glacier area changes by comparing streamflow simulations using the earliest glacier coverage (1985) with the simulations based on the 2000 glacier coverage for 1998, the year with the highest historical glacier icemelt.Although late summer flows are lower for simulations based on the smaller glacier coverage in 2000, the sensitivity of streamflow to historical glacier area changes is small and within the variation associated with parameter uncertainty.
Averaged over the Mica basin, the climate scenario based on the A1B emission scenario generates a warming of 3.0 • C and an 11 % increase in precipitation.Figure 11 shows that the effect of glacier retreat on predicted streamflow can be significant, depending on the stage and rate of deglaciation.Differences between the simulations for static and dynamic glacier cover start to emerge from parameter uncertainty around 2060, when glacier area decreased by 40 % relative to the glacier area in 2000.

Discussion
The Nash-Sutcliffe efficiency E exhibits an optimum value associated with negative glacier volume change; i.e. parameter sets that incorrectly predict positive mass balance do a worse job of predicting streamflow.In contrast to other studies (e.g.Stahl et al., 2008), the optimal value of E can be achieved by a wider range of glacier volume changes (−5 to −40 km 3 ), which reflects the modest glacier cover in the large Mica basin and the associated lower sensitivity of streamflow to glacier melt, relative to studies of more heavily glacierized catchments.However, even if there had been a narrow peak of E over glacier volume loss (like in Stahl et al., 2008), Schaefli and Huss (2011) showed that one has to be careful not to infer glacier volume loss from such a relation since, due to potential model structural errors, the "real" glacier volume loss does not need to coincide with the glacier volume loss that gives the best model performance.Because Stahl et al. (2008) had winter mass balance measurements, they were able to fix the climate parameters, T lapse and P lapse , at an initial step during model calibration, separately from the calibration using streamflow data.In contrast, our approach does not allow climate related parameters to be calibrated independently from the streamflow simulations.Hence, a greater amount of uncertainty in these parameters (wider parameter ranges) is propagated through to our streamflow predictions.
The guided GLUE approach for model calibration clearly demonstrated the value of glacier volume change for constraining model parameters that control snow accumulation and glacier mass balance.Schaefli and Huss (2011) also combined the classical GLUE method with a global optimization algorithm, although with a different goal: they used the global optimization algorithm to find model structural errors while we used global optimization to gauge how well the model can potentially perform to aid the decision to accept or reject parameter sets.Schaefli and Huss (2011) concluded that seasonal mass balance data are a prerequisite for constraining model parameters.However, our results indicate that even glacier-wide volume change spanning several years can constrain model parameters such that the glacier ice contribution to streamflow can be quantified with reduced uncertainty.Given the general lack of mass balance data worldwide, our approach should prove useful for assisting in model calibration, particularly in large basins with modest glacier cover, where goodness-of-fit indices like the model efficiency are less sensitive to the streamflow variability related to glacier contributions.If a hydrologic model will be used to make future projections of the effects of climate change, it is imperative that a model be able to simulate glacier mass balance with reasonable accuracy, not just streamflow.
Given that repeated DEM mapping may only be available over periods of a decade or more, the approach applied here will require relatively long calibration periods -15 yr in the case of Mica basin.The future projections (Fig. 11) demonstrate that changes in glacier area over decadal time scales can potentially have a significant influence on summer streamflow.Therefore, under conditions of rapid glacier retreat, it may be necessary to represent the effects of glacier shrinkage during the calibration period, even in large catchments like Mica basin with modest glacier cover.Current long term planning data sets for the Columbia Region, used for example in hydroelectric operations planning, either do not account for glaciers (Hamlet et al., 2010) or assume static glaciers (Schnorbus et al., 2011;Bürger et al., 2011).Our results suggest that, for climate change impact assessments where glaciers are projected to recede substantially, the effects of glacier recession on streamflow have to be considered even in basins with modest glacier cover (less than 10 %).
It is awkward to incorporate automated glacier area updates during model calibration and during long-term future projections using a model like HBV-EC, which represents a basin using static land cover.In our study, these updates were accomplished using rather complicated scripting.Despite the efficient structure of HBV-EC, which employs grouped response units and lumped reservoirs, the calibration process took substantial processing time (one week for 10 000 model runs on five CPUs, two weeks for the evolutionary optimization on one CPU).This challenge is not unique to the HBV-EC model, since most existing model codes that we are aware of do not allow for changes in land cover during a simulation run.One solution would be to develop a new model code that can accept updated land cover information without having to stop and restart execution.
Glacier contributions to Mica basin streamflow are greatest in August and September.Although minor in terms of long-term average flows, glacier ice melt is especially important during relatively warm, dry weather in summers following a winter with low snow accumulation and early snowpack depletion.These conditions can be critical from both water supply and ecological perspectives.Therefore, water managers and aquatic ecologists need to appreciate the hydrologic significance of glacier melt, even in large basins with moderate glacier cover.In a large basin such as Mica, glacier ice melt can contribute up to 25 % to 35 % of late summer streamflow.

Conclusions
Use of glacier volume change in the calibration procedure effectively reduced parameter uncertainty and helped to ensure that the model was accurately predicting glacier mass balance as well as streamflow.This approach should be widely useful for quantifying glacier contributions to streamflow in glacier-fed catchments where mass balance observations are lacking.One drawback to the approach is that the calibration period must span the interval between glacier maps, in this case 15 yr.Because glacier cover can change significantly over decadal and longer time periods -historically and in future projections -approaches should be developed to accommodate glacier cover changes during model simulations.
Although glaciers only cover 5 % of the Mica basin, they contributed up to 25 % of mean August flow and 35 % of mean September flow in the historic period.These contributions are particularly important during periods of warm, dry weather following winters with low accumulation and early snowpack depletion.Glacier retreat over the twentyfirst century could therefore have significant implications for streamflow during critical late-summer periods.
Acknowledgements.This work was supported financially by the Canadian Foundation for Climate and Atmospheric Science through its support of the Western Canadian Crysopheric Network, a contract from BC Hydro Generation Resource Management, and a Natural Sciences and Engineering Council of Canada Discovery Grant to RDM.Tobias Bolch and Erik Schiefer contributed to the development of historic glacier masks and digital elevation models.Sean Fleming, Frank Weber and Scott Weston of BC Hydro Generation Resource Management assisted with access to data and project administration.Sean Fleming and Frank Weber also provided detailed reviews of earlier versions of this work.Garry Clarke, Faron Anslow, Alex Jarosch and Valentina Radic conducted the simulations of future glacier response under climatic warming scenarios.Alex Cannon assisted with downscaling GCM output.The anonymous reviewers provided constructive comments that helped us improve the ms.

Fig. 1 .
Fig. 1.Location of the Mica basin and locations of climate stations used to force the hydrological model.

FindFinal
Fig. 2.Flow chart illustrating the "guided" GLUE approach for calibration and uncertainty analysis.E is the Nash-Sutcliffe efficiency.

Fig. 3 .
Fig. 3. Nash-Sutcliffe efficiency (E) plotted against simulated glacier volume change for 10 000 model runs in the initial Latin Hypercube Search (black) and for 10 000 model runs in a Latin Hypercube Search with adjusted prior parameter distributions (blue).Red dots indicate acceptable parameter combinations.

Fig. 4 .
Fig. 4. Observed and simulated discharge for the test period (2000-2007).(a) observed and simulated discharge predicted with the best performing (E) parameter set; (b) observed and the ensemble of simulated discharge.

Fig. 5 .
Fig. 5. Ensemble simulated and observed snow water equivalents for three snow pillow sites.Simulations are for the GRU that corresponded to the snow pillow sites, and show the 5-95 % quantile range.

Fig. 6 .
Fig. 6.The effect of glaciers on mean annual discharge shown by comparing simulations with and without glaciers including uncertainty limits (5-95 % quantile range).

Fig. 7 .
Fig. 7.The effect of glaciers on mean August discharge shown by comparing simulations with and without glaciers including uncertainty limits (5-95 % quantile range).

Fig. 8 .
Fig. 8.The effect of glaciers on mean September discharge shown by comparing simulations with and without glaciers including uncertainty limits (5-95 % quantile range).

Fig. 9 .
Fig. 9.The effect of glaciers on discharge shown by comparing simulations with and without glaciers for the year with the highest modelled ice melt (1998) and the year with the lowest ice melt (2000).

Fig. 10 .
Fig. 10.Illustration of the sensitivity of simulated streamflow to glacier changes during the calibration period, based on a comparison of ensemble simulations for the year 1998 using the glacier coverages from 1985 and 2000.

Fig. 11 .
Fig. 11.Ensemble simulations of mean August streamflow for a climate scenario based on the A1B emissions scenario using a static glacier cover (based on the observed glacier cover in 2005) and a dynamic glacier cover.

Table 1 .
Description of calibrated model parameters, benchmark parameter values, parameter ranges for LHS, and correlation between each parameter and glacier net mass balance (b n ) and Nash-Sutcliffe efficiency (E).