Introduction
Dynamic glacier and ice-sheet models are used to investigate climate–cryosphere relationships for both past and future climatic conditions (Reference Hostetler and ClarkHostetler and Clark, 2000; Reference Otto-Bliesner, Marshall, Overpeck, Miller and HuOtto-Bliesner and others, 2006), with surface forcing provided by mass-balance models. While the energy budget approach provides the most physically rigorous basis for modeling glacier ablation, temperature-index or degree-day models have proven to be effective and robust for many glaciological and hydrological applications, particularly in light of their modest data requirements (Reference MooreMoore, 1993; Reference Jóhannesson, Laumann and KennettJóhannesson and others, 1995; Reference Ritz, Fabre and LetréguillyRitz and others, 1996; Reference HockHock, 1999; Reference Marshall and ClarkeMarshall and Clarke, 1999; Reference OhmuraOhmura, 2001). Given the advances in spatial modeling of daily meteorological data in mountain regions, particularly air temperature (Reference Thornton, Running and WhiteThornton and others, 1997; Reference Stahl, Moore, Floyer, Asplin and McKendryStahl and others, 2006), it is feasible to apply degree-day models to estimate historic variations of glacier melt over broad regions using standard climate data (Reference Braithwaite and RaperBraithwaite and Raper, 2007). Even where more complex models are preferred, or when new data sources such as re-analysis fields are used, degree-day models can provide a baseline for judging the value of more complex approaches for different variables and timescales (e.g. diurnal streamflow vs seasonal mass balance; Reference HockHock, 1999). Furthermore, degree-day models provide a simple basis for estimating mass-balance sensitivity to climate (e.g. Reference Braithwaite, Zhang and RaperBraithwaite and others, 2003; Reference De Woul and HockDe Woul and Hock, 2005).
A challenge in applying temperature-index models is the specification of snow- and ice-melt factors, which vary among observational sites (e.g. Reference HockHock, 2003) and between observation periods (Reference Rango and MartinecRango and Martinec, 1995) when calculated using in situ measurements of temperature and snow- and/or ice-melt. Several dynamic ice-sheet models have used values of 3 and 8 mm d−1 °C−1 for snow- and ice-melt factors, respectively (Reference Ritz, Fabre and LetréguillyRitz and others, 1996; Reference Marshall and ClarkeMarshall and Clarke, 1999; Reference Tarasov and PeltierTarasov and Peltier, 1999; Reference Casal, Kutzbach and ThompsonCasal and others, 2004), based on observational data (Reference HockHock, 2003). However, the validity of these values for use with (1) regional air-temperature fields and (2) seasonally integrated positive degree-days has not been established.
If melt factors determined from in situ measurements are to be used, gridded and/or extrapolated regional temperature data must be reconciled with air temperatures measured within the glacier boundary layer. For example, Reference Braithwaite, Zhang and RaperBraithwaite and others (2003) incorporated a ‘cooling effect’ to account for the differences between air temperatures measured above the glacier and those extrapolated from regional stations. The other issue is the temporal stability of melt factors, which vary throughout the melt season due to changes in albedo, insolation and the relative contributions of surface energy-balance fluxes (Reference Rango and MartinecRango and Martinec, 1995; Reference Singh and SinghSingh and Singh, 2001). Thus, melt factors measured at discrete periods during the course of the melt season may not represent the season as a whole.
The objective of this study is to introduce a method for estimating bulk seasonal melt factors suitable for application in regional mass-balance models driven by gridded/interpolated air-temperature fields, with a geographical focus on western Canada. The analysis addresses the following questions: (1) how variable are the melt factors among glaciers; (2) at a given glacier, how much do the melt factors vary from year to year; (3) how sensitive are the fitted melt factors to the lapse rate assumed in the temperature interpolation; (4) what impact does the assumption of the presence/absence of firn have on fitted melt factors; and (5) can melt factors determined at one long-term mass-balance site generate accurate predictions of summer balance at other sites? The derived values are also used to calculate static mass-balance sensitivities to summer temperature increases.
Study area and Methods
Study area
The southwestern Cordillera of Canada are comprised of numerous mountain chains including the southern Coast Mountains, the Columbia Mountains and the Canadian Rocky Mountains. Annual precipitation is delivered mainly in winter by low-pressure systems approaching from the southwest, producing strong west–east precipitation gradients and local distributions determined by windward and leeward slopes. Total annual precipitation averages 1230 mm at Whistler, British Columbia and 480 mm at Golden, British Columbia (Environment Canada, http://www.climate. weatheroffice.ec.gc.ca; Fig. 1) while mean net winter balances from 1965 to 1995 range between 1.75 m at Place Glacier and 1.19 m at Peyto Glacier (Reference DyurgerovDyurgerov, 2002). January (July) mean temperatures are −3.0°C (13.9°C) at Whistler (658 m a.s.l.) and −9.7°C (17.1°C) at Golden (780m a.s.l.), indicating that eastern sites will have a higher degree of continentality. Net glacier mass balances throughout the region have been predominantly negative since massbalance measurements began in 1965/66.
Mass-balance data
Mass-balance records at nine sites in southwest Canada (Fig. 1; Table 1) were obtained from published sources (Reference DyurgerovDyurgerov, 2002) and government reports (Reference Mokievsky-Zubok, Ommanney and PowerMokievsky-Zubok and others, 1985). Reported specific mass-balance data (bj ,z ) for year j represent the winter (b w) or summer (b s) balance (in mm w.e.) at elevation z. For calculating positive degree-days, we assumed that reported specific mass-balance values represent the value for the midpoint of a given elevation band z. For convenience, summer balance is treated as a positive value.
The period spanned by the mass-balance data includes substantial climatic variability. The period samples both cold and warm phase regimes of the Pacific Decadal Oscillation, as well as warm and cold phases of the El Niño–Southern Oscillation (Reference Moore and DemuthMoore and Demuth, 2001).
Initial analysis of mass-balance data indicated that some data filtering was necessary. Reported values that remained constant with elevation or followed fixed gradients (e.g. 500 mm w.e. (100 m)−1) were removed from the analysis. These values were assumed to be estimated due to (1) loss of data associated with snow burial or meltout of ablation poles or (2) weather or logistical constraints preventing observation. We also removed observations where summer balance equaled zero (n = 69), as well as specific balance observations for Sykora Glacier that were duplicated from observations for Bridge Glacier (n = 55). The two glaciers were reported as one unit from 1977 to 1982 (Reference Mokievsky-Zubok, Ommanney and PowerMokievsky-Zubok and others, 1985). As a consequence of this data filtering, only individual specific mass-balance values were analyzed and no attempt was made to estimate mean specific summer balances.
Air-temperature data
Near-surface air temperatures were reconstructed for each glacier site using an inverse-distance-squared interpolation of federal and provincial climate data sources (Reference Stahl, Moore, Floyer, Asplin and McKendryStahl and others, 2006). To estimate average daily temperatures at the elevation bands of each site, we assumed a constant lapse rate of 6.0°C km−1, which is consistent with previous melt-modeling studies in temperate locations (Reference MooreMoore, 1993; Reference Jóhannesson, Laumann and KennettJóhannesson and others, 1995; Reference VincentVincent, 2002).
Mass-balance sites used in this study are substantially higher in elevation than the climate stations (cf. Reference Stahl, Moore, Floyer, Asplin and McKendryStahl and others, 2006, fig. 1). However, low-elevation climate data are more likely to reflect variations in insolation and thermal characteristics of the regional air mass than climate data collected from within the glacier boundary layer, and thus may provide a more suitable positive degree-day (PDD) estimate (Reference Lang, Kraijenhoff and MollLang, 1986; Reference VincentVincent, 2002). In order that calculated melt factors were compatible with regional or gridded temperature fields, we did not consider glacier cooling effects.
The dates of mass-balance collection varied between sites and years, but were not reported for many observations. For the observations for which the measurement dates are known, the mean dates are 20 May and 28 September. To calculate PDDs, we assumed a nominal ablation season extending from 15 May to 30 September.
Derivation of melt factors
Separate melt factors for snow and ice can be fitted using an approach described below. The model was applied using PDD sums computed between 15 May and 30 September with a lapse rate of 6.0°C km−1.
Below the glacier equilibrium-line altitude (ELA), specific summer balance is the sum of snow (M s) and ice (M i) melt (in mm), if melt from summer snowfall events is assumed to be negligible and ice melt occurs once the winter snow cover has been removed:
Seasonal snow and ice melt can be calculated as the product of melt factors (k s and k i) and the accumulated PDD for each surface (PDDs and PDDi):
where δ and ∊ are error terms.
Assuming that mid-summer snowfalls are small compared to the total winter balance, seasonal snowmelt (M s) will be approximately equivalent to the amount of snow initially available for melting (or b w):
Combining Equation (5) with Equation (2) yields
Inserting Equation (6) into Equation (4) and combining with Equations (1) and (2) yields an expression for the summer balance below the ELA:
which has the form of a multiple linear regression through the origin.
Above the ELA,
Equation (8) has the form of a simple linear regression through the origin. At the ELA,
Equations (7–9) suggest a piecewise linear regression model of the form of Equation (7) for b w < k s PDD, and of the form of Equation (8) for b w ≥ k s PDD. Melt factors for snow and ice were optimized for each site through non-linear least-squares regression using the ‘nlinfit’ function in Matlab. Summer balance b s was the dependent variable and b w and PDD were the independent variables. For each glacier, models were fitted using values for all years and elevation bands to generate one value for each of k s and k i. In addition, melt factors were estimated for each year for Peyto and Place Glaciers, which have the longest periods of record.
Sensitivity analyses of fitted melt factors
There are two sources of uncertainty in the calculation of PDD sums: (1) the length of the ablation season over which to calculate the sums (which should be based on the dates of mass-balance observations) and (2) the lapse rate used to adjust interpolated air temperature for elevation. To account for the former, four different ablation season lengths (ABL; Table 2) were used in the calculation of accumulated PDDs. Sensitivity of fitted melt factors to the method of temperature extrapolation was examined by calculating PDD totals using (1) a constant lapse rate of 6.5°Ckm−1 and (2) mean monthly lapse rates which varied between 3.0 and 8.0°C km−1, reflecting surface observations (Reference Stahl, Moore, Floyer, Asplin and McKendryStahl and others, 2006).
Another source of uncertainty in the fitted melt factors stems from the fact that the model does not distinguish between ice melt and firn melt when mass balances are negative. Firn is a winter accumulation of snow which persists through more than one melt season, giving it a lower albedo than fresh snow and a higher albedo than ice. While most degree-day models do not explicitly account for the presence or absence of firn, it is possible that the presence of firn might result in a lower fitted melt factor for ice. The piecewise model used in this study assumes that once the snow accumulations of the current winter are removed, ice melt begins. If firn were present, this would result in a reduced melt rate and thus a lower melt factor.
To evaluate the sensitivity of the fitted melt factors to the presence of firn, the surface type was estimated for each elevation band from cumulative and annual net balances (b n). Observations where cumulative b n was positive and annual b n was negative were assumed to indicate that firn melting had occurred, and these data were flagged and removed from the statistical analysis. In addition, data from the first measurement year were removed if subsequent years demonstrated continuous firn cover. Piecewise regressions were then repeated at each site using the dataset filtered for firn melt.
Testing regional predictions of seasonal ablation
Specific summer balances at all sites were predicted using PDD j ,z calculated from regional temperature interpolations and two sets of melt factors. The first set was derived from Place Glacier and was applied to assess the accuracy of summer balance predictions based on melt factors from a single long-term mass-balance site. The melt factors for Place Glacier were partly chosen because of the length of the mass-balance record, but also for its proximity to the majority of the sites studied. A second set of melt factors (k s = 3.0 mm°C−1 d−1, k i = 8.0 mm°C−1 d−1) was chosen to represent parameters assumed in several dynamic ice-sheet models (Reference Ritz, Fabre and LetréguillyRitz and others, 1996; Reference Marshall and ClarkeMarshall and Clarke, 1999; Reference Casal, Kutzbach and ThompsonCasal and others, 2004). These values are within the range of melt factors calculated for glacier sites using melt observations and in situ temperature measurements (Reference Braithwaite and ZhangBraithwaite and Zhang, 2000; Reference HockHock, 2003).
Seasonal snow- and ice-melt totals were estimated following Equations (2) and (4) by first using k s and PDD j ,z to deplete the observed winter balance b w. Remaining positive degree-days were directed toward ice melt using k i, and snow- and ice-melt totals (in mm w.e.) were summed to obtain an estimate of b s (Equation (1)).
Static summer mass-balance sensitivity
The static summer mass-balance sensitivity (ST ) to changes in temperature (Reference Oerlemans and FortuinOerlemans and Fortuin, 1992; Reference Braithwaite and ZhangBraithwaite and Zhang, 1999; Reference Braithwaite, Zhang and RaperBraithwaite and others, 2003; Reference De Woul and HockDe Woul and Hock, 2005) was determined by prescribing a daily 1 K increase, recalculating PDD totals and then re-estimating seasonal melt totals. The approach used here assumes no change in winter precipitation, and does not include the effect of summer snowfalls or a lengthened ablation season. Temperature sensitivity for each site was calculated as the mean difference between original and perturbed melt estimates, averaged over all data points.
Results
Fitted melt factors
Fitted melt factors from the reference model run (15 May–30 September and a lapse rate of 6.0°C km−1) range from 2.31 to 3.61 for k s and 3.61 to 5.57 for k i (Table 3). Except for Zavisha Glacier (where R 2 is 0.37), model fits are reasonably strong with R 2 up to 0.90. Calculated values of k s and k i are statistically different, with the exception of one outlying ice-melt factor derived from the weakest regression (Zavisha Glacier). Peyto Glacier has the strongest fit as well as the lowest k s and the highest k i.
Snow- and ice-melt factors are relatively consistent among sites, with mean (standard deviation) k s =3.04 (0.38) and k i = 4.59 (0.59). The mean values are similar to the optimized melt factors derived from a mass-balance and streamflow-calibrated hydrological model at Bridge Glacier (Reference Stahl, Moore, Shea, Hutchinson and CannonStahl and others, 2008), where k s = 3.09 and k i = 4.70 on 21 June. They also fall within the seasonal range of streamflow calibrated melt factors for the Lillooet River basin, located west of Place Glacier and north of Helm Glacier (Reference MooreMoore, 1993).
Melt factors calculated for Peyto Glacier (Fig. 2; Table 4) vary modestly from year to year, with low standard deviations and high correlation coefficients. Ice-melt factors exhibit higher variability between years than snowmelt factors. At Place Glacier, melt factors generally demonstrate a similar magnitude of interannual variability. However, the 1978 season shows k s to be higher than k i, suggesting a limit to this approach. In particular, fitting melt factors to individual years greatly reduces the degrees of freedom, and correspondingly increases the error in the fitted melt factors.
Degree-day and melt factor sensitivity
Changing the ablation season length produced a non-linear response in accumulated positive degree-days with elevation for all lapse rate scenarios. Figure 3 illustrates two cases of this sensitivity for Bridge and Place Glaciers. At the higher elevations, PDD are nearly equal for all ablation season lengths. This is because mean daily temperatures are substantially greater than 0°C only on the warmest days (typically in June, July and August). At lower elevations, there is a larger spread between PDD totals calculated for different ablation seasons, but averaged over all elevations for Place Glacier there is only a 2% difference between the model run (15 May–30 September) and ABL3 (1 May–30 September). Extending the ablation season by 2months (1 April–30 October) increased accumulated PDD by an average of 4%, but only PDD accumulated between the dates of winter and summer mass-balance collection are of interest. Therefore, computed PDD are relatively robust to uncertainty in the dates of mass-balance observations. Examples of the sensitivity of PDD sums to lapse rate model are shown in Figure 4. For all sites and elevations, PDD calculated using varying monthly lapse rates are similar to those calculated using a fixed lapse rate of 6.0°Ckm−1. For the examples shown, a fixed lapse rate of 6.5°C km−1 decreased PDD accumulated between 15 May and 30 September by 13% for Bridge Glacier (1978) and by 15% for Place Glacier (1987) relative to the totals calculated using 6.0°C km−1.
Sensitivity analyses for the 12 PDD models (four ablation seasons and three lapse rate scenarios) for Peyto and Place Glacier melt factors are presented in Table 5. Melt factors are nearly equivalent for the variable lapse rate and the 6.0°C km−1 fixed lapse rate, while a greater lapse rate results in higher fitted melt factors. Varying the time period over which PDD are calculated produces relatively small changes in calculated k s, whereas k i exhibits a stronger sensitivity to assumed ablation season length. However, even for k i, the difference between the model run and the most likely range of ablation season lengths (ABL3, ABL4) is no more than 10% of the base case value. Accounting for the presence or absence of firn resulted in small and non-systematic changes to the fitted melt factors. Of the original 1029 mass-balance observations, 387 were removed due to the possible presence of firn and the melt factors were recalculated. At Peyto Glacier k i increased by 4%, but the mean k i for all sites increased by only 0.1% following the removal of mass-balance observations containing firn. Aerial photographs and the strongly negative cumulative mass balance of these sites suggest that firn extents are generally quite low in this region. However, melt factors for ice will be underestimated at sites where firn is more extensive, and in such cases mass-balance observations from the firn area should be excluded from the regression analysis.
Regional melt estimation
Melt factors assumed in several modeling studies (k s = 3.0, ki = 8.0) generally over-predicted ablation. Melt factors fitted for Place Glacier (k s = 2.59, k i = 4.51) provided more accurate predictions of summer balance at most sites, although there was a tendency to under-predict at several sites including Bridge and Woolsey Glaciers (Fig. 5). Both sets of melt factors overestimated b s for low summer ablation totals, and summer melt was poorly estimated at Helm and Zavisha Glaciers. The former was more accurately estimated using the higher ice-melt factor.
Discussion
Fitted ice-melt factors calculated in this study are lower than those observed in situ at other locations using on-glacier air-temperature measurements (Reference BraithwaiteBraithwaite, 1995; Reference HockHock, 1999). One possible explanation for this result is the fact that the methods presented here do not account for summer snowfalls, which would bias (reduce) the fitted melt factors. However, this source of bias is believed to be small at the lower elevations of the study glaciers. Summer precipitation in southern British Columbia is much lower than winter snow accumulation, and average summer (JJA) temperatures at the lower elevations are 9°C, indicating that most summer precipitation will fall as rain. Summer snowfalls may have a stronger influence at the highest elevations.
The most likely cause for the difference in ice-melt factors found in this study and those calculated from in situ observations is that screen-level temperatures measured within the glacier boundary layer are affected by sensible heat exchange and katabatic flow (Reference Van den BroekeVan den Broeke, 1997; Reference Greuell and BöhmGreuell and Böhm, 1998; Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004). Temperatures observed within the boundary layer are lower relative to air temperatures measured over non-glacier surfaces at the same elevation. Furthermore, this cooling over the glacier is greatest at the lowest elevations because of stronger katabatic flow and enhanced sensible heat exchange. The use of on-glacier screen-level temperatures for calculating melt factors therefore produces a higher ice-melt factor relative to a melt factor calculated from extrapolated valley-bottom temperatures, which are not influenced by the glacier boundary layer. The fitted snowmelt factors are probably similar to in situ observations because the estimated values of k s are strongly influenced by observations of snowmelt above the ELA, where boundary-layer development is weakest.
One major implication of this research is that mass-balance models using the degree-day approach should use degree-day factors appropriate to the air-temperature data used to drive the model. In particular, if regionally interpolated air-temperature fields are used, the degree-day factors should be consistent with the lapse rates used to adjust air temperature for elevation.
Melt factors calculated annually for two long-term datasets show modest interannual variability, with coefficients of variation ranging from 14 to 25%. Interannual variability in melt factors arises from sampling variability and also from variations in weather conditions that influence the relative contributions of different heat fluxes. Sensitivity to variations in heat fluxes appears to be greater for melting of ice than snow, judging by the greater variability of k i. In principle, a process-based model should be able to account explicitly for this variability. An important task for future study will be to assess whether this potentially superior performance can be realized for regional prediction of glacier melt, especially in situations where the driving meteorological variables must be estimated from a regional station network.
Summer melt was underestimated at both Helm and Zavisha Glaciers, the smallest glaciers in the sample. The fitted melt factors at Helm Glacier were the highest (for k s) and second highest (for k i) among the glaciers studied. The reasons for these high melt factors are not obvious, and detailed field investigations would be required to estimate the energy fluxes responsible. The snowmelt factor was higher at Zavisha Glacier than that estimated for Place Glacier, while the ice-melt factor was lower. However, these estimates are highly uncertain due to the poor fit (R 2 = 0.37) and should be interpreted with caution.
At several sites, low values of specific summer balance are overestimated using melt factors from various dynamic ice model studies (k s = 3.0, k i = 8.0) and values derived from Place Glacier (k s = 2.59, k i = 4.51) (Fig. 5). These errors might result from summer snowfalls which had not been accounted for. This would reduce melt totals by temporarily increasing the surface albedo. Alternatively, nocturnal refreezing would increase the cold content of the snowpack, and this effect would also be greatest at the highest elevations. A lower snowmelt factor might therefore be required to accurately model high-elevation melt at both daily and seasonal timescales.
Conclusions
A piecewise linear regression model provided a straightforward approach for estimating melt factors using regionally interpolated air temperature and specific mass-balance data, yielding values that were broadly consistent over southwest Canada. Fitted ice-melt factors from this study are lower than values calculated at other temperate sites using in situ temperature observations, which may reflect the effects of katabatic flow and surface heat exchange on air temperatures in the glacier boundary layer. Calculated melt factors are relatively insensitive to the assumed ablation season length, but somewhat more sensitive to the choice of lapse rate. Removing observations where firn was present resulted in negligible changes to the fitted melt factors at these sites. However, firn will likely have more influence on melt factors at locations where it is more extensive. Snowmelt factors calculated annually for two long-term datasets were relatively stable, despite the broad range of climatic conditions that occurred during the study period.
Melt factors fitted to mass-balance data from a single site (Place Glacier) provided reasonable predictions at most other sites representing both maritime and continental climates, although there was a tendency for under-prediction at several sites. The combination of regionally interpolated air temperatures and a degree-day model appears capable of generating first-order estimates of regional summer balance, which can provide a benchmark against which to judge the predictive ability of more complex (e.g. energy-balance) models applied at a regional scale. In particular, a more process-based model should, in principle, be able to account for the factors that produce the inter-glacier and interannual variations in melt factors. Mass-balance sensitivity analyses indicate that a temperature increase of 1 K will increase ablation in the region by 0.51 m w.e. a−1 on average.
Specific mass-balance data reported by elevation band represent an under-utilized source of information for calculating seasonally integrated melt factors for snow and ice. The method introduced here can be easily duplicated at other mass-balance sites where appropriate estimates of accumulated PDD are available, and they can also be used to evaluate melt factors by year or by elevation, given a sufficient number of observations. A compilation of melt factors estimated from a range of sites might reveal systematic geographic or climatic variations in melt factors, with important consequences for global cryospheric models and estimates of future rates of sea-level rise.
Summer mass-balance sensitivity to a 1 K temperature increase
Calculated ST for glaciers in western Canada range between −0.43 and −0.56 m w.e. a−1 K−1 (Table 3), which is within the range of values reported by other studies (Reference Oerlemans and FortuinOerlemans and Fortuin, 1992; Reference De Woul and HockDe Woul and Hock, 2005; Reference Braithwaite and RaperBraithwaite and Raper, 2007). These values are relatively low given the temperate location and moderate accumulation rates (average accumulation at the ELA of Place Glacier is 1.90m w.e.), but they do not include changes in the precipitation type due to increasing temperature. Maritime sites (Helm and Place Glaciers; Fig. 1) tend to have higher sensitivities than continental sites (Peyto Glacier), which is consistent with previous studies.
Acknowledgements
This study has been supported by funding from the Natural Sciences and Engineering Research Council through a scholarship to J.M.S. and a Discovery Grant to R.D.M., the Government of Canada through the Climate Change Impacts and Adaptation Research Program (grants A676 and A875) and the Canadian Foundation for Climate and Atmospheric Science through its support of the Western Canada Cryospheric Network. G.K.C. Clarke, R. Braithwaite, M. Citterio, R. Hock and an anonymous reviewer provided useful comments on earlier versions of this work.