Introduction
The magnitude of glacier ablation depends on the energy balance at the glacier surface. For remote glacier locations, the large amounts of data required to determine this balance are often unavailable. For this reason, degree-day models (DDMs) based on the strong observed relationship between glacier melt and near-surface air temperatures are often used to estimate surface melt. Most studies that employ DDMs must downscale or interpolate climate model output and/or station observations to the area of interest. For spatially distributed models, temperature fields must be generated to reflect local conditions, which often include complex topography. This is most commonly done through the use of digital elevation models and a constant free air lapse rate (change in temperature with elevation) that is often assumed to be the moist adiabatic lapse rate (6–7°Ckm–1). The problem with this approach is that temperature lapse rates measured over melting glacier surfaces have been found to be systematically less steep than the free air values (Reference Greuell and BöhmGreuell and Böhm, 1998; Reference Braun and HockBraun and Hock, 2004; Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and StephensHanna and others, 2005; Reference Marshall, Sharp, Burgess and AnslowMarshall and others, 2007). These deviations are mainly the result of temperatures over melting surfaces being fixed at 0°C, while free atmospheric temperatures generally increase with decreasing elevation. This can generate strong temperature gradients between the surface and overlying atmosphere that increase with decreasing elevation when air temperature is above 0°C. The strong temperature gradients drive a sensible heat flux that cools the near-surface atmosphere at lower elevations and effectively reduces the magnitude of the near-surface temperature lapse rates (Reference Greuell and BöhmGreuell and Böhm, 1998).
Recent investigations into near-surface temperature lapse rates in the Canadian High Arctic show that lapse rates vary on both seasonal and daily timescales (Reference Marshall, Sharp, Burgess and AnslowMarshall and others, 2007). Gardner and others (unpublished results) found that during the summer months (JJA) lapse rates vary systematically with daily mean air temperatures in the lower (i.e. 750 mbar) troposphere. They concluded that much of the observed covariability was due to an increase in the climate sensitivity of screen-level temperatures with increasing elevation (cf. Reference Denby, Greuell and OerlemansDenby and others, 2002), where ‘climate sensitivity’ is defined as the change in near-surface temperature relative to a change in the free atmospheric temperature.
In this study, we examine the implications of including variable near-surface lapse rates (VLR) in DDM simulations of the surface mass balance of Devon Ice Cap, Nunavut, Canada. Surface mass-balance estimates are determined for the years 1980–2006 and are evaluated against 23 years (1980–2003) of surface mass-balance measurements from the northwest sector of the ice cap (Reference KoernerKoerner, 1970). These results are also compared with estimates derived using the assumption that model temperatures can be downscaled using a constant near-surface lapse rate equal to either the measured summer (JJA) mean lapse rate (GLR) or the moist adiabatic lapse rate (MALR).
Site and Data Description
Devon Ice Cap is a large (14 000km2; Reference Dowdeswell, Benham, Gorman, Burgess and M. SharpDowdeswell and others, 2004) dome-shaped ice mass covering much of the eastern part of Devon Island, Nunavut (Fig. 1). From 1961 to the present, the Geological Survey of Canada (GSC) has made physical measurements of the net surface mass balance b n along a 40 km long stake network following a path from the glacier summit (1880ma.s.l.) to the terminus of the Sverdrup outlet glacier in the northwest sector of the ice cap (Fig 1; Reference KoernerKoerner, 2005). Previous estimates of surface mass balance (mmw.e. a–1) for the main ice cap (excluding the southwest arm) vary greatly between study periods: –76±50 (1961–66 (Reference KoernerKoerner, 1970)), –30±60 (1963–2000 (Reference Mair, Burgess and SharpMair and others, 2005)) and +28±70 (1996 (Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others, 2007)).
Fig. 1. Devon Ice Cap shown with region boundaries and transects along which field measurements have been taken.
Continuous meteorological measurements have been made on the ice cap since 1997, with as many as 7 automatic weather stations and 13 temperature sensors operating along 2 transects in 2004 (Fig. 1). Near-surface temperature lapse rates derived from these transects show strong linear relationships (r = 0.72) with summit elevation temperatures (Fig. 2).
Fig. 2. 1998–2004 average daily summer (JJA) lapse rates β for the northwest transects plotted against respective mean daily summit air temperatures.
Mass-Balance Model
Glacier melt
The Degree-Day glacier Mass Balance model (DDM) determines the daily quantity of melt mt (mmw.e.) as a function of the mean daily air temperature Tt (˚C) using a factor of proportionality referred to as the degree-day factor DDF (mm d–1 ˚C):
Owing to differences in the surface energy-balance characteristics of snow and ice (including albedo, shortwave penetration, thermal conductivity and surface roughness), separate values are used for snow (DDFs) and ice (DDFi). Reported factors vary greatly from study to study, depending on local variations in the surface energy balance (albedo, incoming solar radiation, surface roughness, wind speed and air temperature), how melt was determined (measured or modeled), the time period over which the study was conducted, and how air temperatures were determined (daily or hourly mean). Reference HockHock (2003) provides a comprehensive review of the temperature-index melt-modeling method. For this study we adopt the values DDFs = 3.3 mmd–1 °C and DDF〉 = 8.2 mmd–1 °C, which are the mean values determined by Reference BraithwaiteBraithwaite (1995) at two locations (880 and 790ma.s.l.) on the Greenland ice sheet. Reference BraithwaiteBraithwaite (1995) obtained these values through the use of an energy-balance model over a 6 year period and they agree well with factors determined from stake measurements made at the same locations. These are the same factors used by Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others (2007) to model melt over Devon Ice Cap for 1996, and compare with DDFs = 4 and DDF〉 = 6–14mmd–1 °C used by Reference Mair, Burgess and SharpMair and others (2005) to model melt over the ice cap for the years 1963–2000.
Surface air temperature and downscaling
Daily mean surface air temperatures were downscaled to the entire ice cap from summit elevation temperatures using both modeled mean daily and assumed constant lapse rates and a 3000 ×3000 digital elevation model (DEM; Canada3D: available from Natural Resources Canada). Daily mean 2 m air temperatures from the North American Regional Reanalysis (NARR; Reference MesingerMesinger and others, 2006) were extracted for the model gridcell centered on the summit of the ice cap and were corrected for a mean summer (JJA) temperature offset of +4°C relative to 8 years (1997–2004) of air-temperature measurements from an on-glacier weather station located near the summit of the ice cap. The 4°C offset in the NARR air temperature is partly due to the coarse resolution of the NARR model topography, which models the maximum elevation of Devon Ice Cap at an elevation 400 m lower than the true summit elevation. After adjusting for the mean difference, measured and modeled summer daily mean summit air temperatures show good agreement (r = 0.86) except for systematic underestimations of maximum and overestimation of minimum summer summit temperatures. To correct for this bias, we add or subtract corrections drawn randomly from a normal distribution with a range of 0–2°C to all NARR summer temperatures above –4.5°C and below –7.5°C, respectively (Fig. 3). Surface air temperatures were generated for both ablation season days and non-ablation season days in order to determine mean annual and winter temperatures needed for the calculation of maximum meltwater refreeze amounts, which are discussed below.
Fig. 3. NARR bias-adjusted summer (JJA) daily summit air temperatures (T) plotted against measured summit temperatures (1880 m). The two solid lines show the linear regression relationship between NARR and measured summit temperatures before (grey: r = 0.86) and after (black: r = 0.85) the temperatures have been adjusted for overly warm temperatures on relatively cool days and overly cool temperatures on relatively warm days.
To investigate the impact on mass balance of using different near-surface temperature lapse rates β to downscale air temperatures, three separate surface air-temperature fields were generated. The first was generated using a near-surface lapse rate equal to the moist adiabatic lapse rate 
the second using mean measured lapse rates 
for all days during the ablation season (15 May–15 September) and all winter days (16 September–14 May: 
and the third using a modeled variable daily lapse rate β
VLR(°Ckm–1) during the ablation season and 
in winter. β
VLR was estimated using the lapse-rate equation developed from over 10 years of screen-level air-temperature measurements made over Devon Ice Cap (Gardner and others, unpublished results):
where T
750 is the 3 day running mean of the daily anomaly (summer mean removed) in NARR 750 mbar temperature. T
750 was derived from a six-gridcell domain (64×86 km) centered over the ice cap. To illustrate the range in variability of β
VLR, mean daily summit temperatures are plotted with βVLR for the highest (1986: Fig. 4a) and lowest (2001; Fig. 4b) measured mass-balance years in the study period. On 26 July 2001, when temperatures peaked at 5°C, βVLR = –3°C km–1. On this day, downscaling summit temperatures with β
VLR results in sea-level temperatures that are 4°C and 7°C cooler than temperatures downscaled using 
and 
respectively. This represents ∽30–55mmw.e. less ice melt for this day.
Fig. 4. NARR summit elevation temperatures plotted with variable modeled daily lapse rates β VLR for representatively (a) cold and(b) warm years. Modeled lapse rates are less steep when temperatures are anomalously warm, and steeper when temperatures are anomalously cold.
The three model runs forced with mean daily air-temperature fields downscaled using β
VLR, and are referred to throughout this paper as the VLR, MMLR and MALR runs, respectively.
Precipitation
Daily precipitation was derived from the NARR and corrected for a mean over-precipitation bias of 100mmw.e. a–1, which was determined from the comparison of NARR mean annual precipitation at the summit with the net mass balance measured at the same location for the years 1979–2003 (personal communication from R.M. Koerner). Taking point values at the center of the NARR gridcells, the bias-adjusted coarse-resolution (3×32 km), NARR daily precipitation fields were linearly interpolated based on horizontal distance to match the resolution of the DEM (∽200×1000 m). The phase of the daily precipitation Pt (mmw.e.) was determined at each gridcell based on the estimated mean daily gridcell air temperature Tt :
Here we have chosen the critical surface temperature for phase change as 0°C (Reference Bassford, Siegert and DowdeswellBassford and others, 2006). Other studies use values as high as 2°C (Reference OerlemansOerlemans, 1991). Using a critical value of 2°C would increase net mass-balance estimates of Devon Ice Cap on average by 40 ±14 mmw.e. relative to the value used here. At the summit of the ice cap, the standard error of the interpolated precipitation was determined as 40 mmw.e. a–1.
Refreezing of meltwater
All liquid water (melt or rain) generated within any gridcell in a given year is assumed to refreeze in that gridcell until the volume of refrozen water F = F max, where F max is the maximum amount of refreeze that can occur in that cell. Once the liquid water exceeds Fmax, all additional rain or melt generated in that year is assumed to leave the glacier as runoff. We use the term ‘refreeze’ to refer to the total amount of meltwater and rain that refreezes on the glacier during a given mass-balance year and can re-melt later in that year and run off; the term ‘internal accumulation’ refers to the amount of water that refreezes, does not run off, and remains on the glacier into the following mass-balance year. There is no communication between gridcells, so meltwater only affects refreeze and internal accumulation quantities of the cell in which that meltwater is generated.
F max was determined at each gridcell using the method proposed by Reference Wright, Wadham, Siegert, Luckman, Kohler and NuttallWright and others (2007). This method determines F max based on a simplified estimate of the total amount of energy required to change the temperature profile of the upper 10–15m of the glacier from an end-of-winter temperature profile to a post-refreezing profile:
where c 〉 is the specific heat capacity of ice at 0°C (2097 J kg–1 °C–1), d 〉 is the maximum depth to which the annual temperature cycle penetrates, L f is the latent heat of fusion for water (333.5 kJ kg–1), T a is the mean annual air temperature and T w is the mean winter temperature (here taken as the mean temperature outside of the ablation season). On Devon Ice Cap the annual temperature cycle penetrates to ∽10–15m depth. Using this value for d 〉 results in F max values which are much too large (>1000mmw.e.). We instead set d 〉 = 1 m, which gives mean F max values at the summit of the ice cap that are ∽60% (130 mm) of the annual snowfall, the value often used for the upper threshold for meltwater refreezing on the Greenland ice sheet (Reference ReehReeh, 1991), and is supported by field measurements made near the equilibrium-line altitude of the Greenland ice sheet (Reference Braithwaite, Laternser and PfefferBraithwaite and others, 1994).
Results
Mass-balance components
Average (1980–2006) annual snowfall, rainfall, refreeze and melt generated from the VLR model run are shown in Figure 5. Similar spatial patterns in the four mass-balance components are seen in all three model runs, but magnitudes vary as all components are dependent on the surface air temperature. Averaged over the entire ice cap (including the southwest arm) for the period of study, mean annual snowfall (156±11mmw.e. a–1), rainfall (43 ±10mmw.e. a–1) and meltwater refreeze (90±7mmw.e. a–1) vary little between model runs, with slightly lower (20mmw.e. a–1) snowfall and slightly higher (20mmw.e. a–1) rainfall amounts for the MALR run due to a higher frequency of above-zero temperatures in that run. There are, however, large differences in melt quantities, with mean annual estimates of 730±250, 960±350 and 1840±440mmw.e. a–1 for the VLR, MMLR and MALR model runs, respectively. Because the model runs differ greatly in both the melt amounts and the altitudes at which melt occurs, the three model runs produce very different internal accumulation distributions (Fig. 6) and thicknesses (8, 12 and 3mmw.e. a–1 for the VLR, MMLR and MALR model runs, respectively).
Fig. 5. Mean annual snow (a), rain (b), meltwater refreeze (c) and melt (d) for Devon Ice Cap averaged over the years 1980–2006. Values determined from the VLR model run (see text for details).
Fig. 6. Annual mean internal accumulation formation estimated from the (a) VLR, (b) MMLR and (c) MALR model runs.
Comparison with measurements (northwest sector)
Modeled point-location net mass-balance (b n) estimates were extracted from the model output along the GSC’s northwest transect and compared with their respective stake measurements at 100m elevation intervals from 150 to 1750 ma.s.l. for the period of available data (1980–2003; personal communication from R.M. Koerner) (Fig. 7). All model estimates are similar and agree well with the stake measurements at elevations above 1450 m, but they diverge with decreasing elevation as estimated melt quantities increase. At 150 m, all model runs significantly overestimate melt quantities relative to measurements by as much as 2800mm for the MALR run. Overall, the VLR and MMLR model runs significantly outperform the MALR model run, with 23 year mean ± standard errors averaged across all elevations of +25±260 and –140±340mmw.e. a–1, respectively, compared with –880±920mmw.e. a–1 for the MALR. Using the sector boundaries outlined by Reference KoernerKoerner (1970), cumulative annual net mass balances B n for the three model runs were estimated for the northwest sector of the ice cap and are shown together with an estimate derived from stake measurements in Figure 8. B n estimates vary greatly between model runs with mean ± standard errors of 9±120, –113±150, –577±223mmw.e. a–1 for the VLR, MMLR and MALR runs relative to measured values. The VLR model run produces B n that are close to measured values, with a 1980–2003 cumulative difference of only 0.2 mw.e.
Fig. 7. Average (1980–2001) annual point-location net mass-balance (b n) model estimates extracted at 100m intervals along the northwest transect shown with respective mass-balance stake network measurements.
Fig. 8. Modeled cumulative surface mass balance B n for the northwest sector of Devon Ice Cap plotted with estimates derived from measurements.
Net mass balance of Devon Ice Cap
Annual B n values were calculated for the main portion of the ice cap (southwest arm excluded: 12 100 km2 (Reference Dowdeswell, Benham, Gorman, Burgess and M. SharpDowdeswell and others, 2004)) and are shown in Figure 9. The MALR model run produces mass-balance values that are, on average, 4 times more negative than values determined from the VLR model run and ∽2.5 times more negative than values determined from the MMLR run. During more positive mass-balance years (e.g. 1986 and 2004) there is good agreement between the VLR and MMLR model runs, but the differences are larger during more negative years (e.g. 2001 and 2006). This is because the total annual positive degree-days (PDDs) for the two runs are most similar when lower-tropospheric temperatures are seasonally low, and most different when temperatures are seasonally warm (i.e. during high melt years). For the whole ice cap (southwest arm included), B n values are 24–41% more negative than for the main ice cap because the entire southwest arm lies below 900m and often lies completely within the ablation zone (Fig. 10).
Fig. 9. Net mass balance B n for the main Devon Ice Cap (southwest arm excluded) for all three model runs.
Fig. 10. Model run estimates of mean annual net mass balance B n for the main (excluding southwest arm) and whole Devon Ice Cap.
Using perturbed historical accumulation measurements and a constant lapse rate of –4.6°Ckm–1 (mean 2001 lapse-rate measure over Devon Ice Cap) to extrapolate sea-level temperature data from the International Arctic Buoy Program over different sectors of the ice cap, Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others (2007) estimated the 1996 net mass balance for the main ice cap to be 28 mmw.e. This value is 42, 140 and 730 mmw.e. larger than 1996 estimates computed using the VLR, MMLR and MALR, respectively. We attribute most of the difference between the VLR model run and this estimate to the slightly higher (12mmw.e.) snowfall estimate used by Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others (2007), which results in both greater net accumulation at higher elevations and less melt at lower elevations, since snow has a lower DDF than ice. The differences between the estimate of Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others (2007) and the estimates generated using MMLR and MALR are attributed to an overestimation of surface air temperatures in the MMLR and MALR simulations (much higher melt and slightly lower snowfall) due to the use of overly steep lapse rates.
Based on our comparisons with independent model estimates of the 1996 net glacier mass balance for the main ice cap and 24 years of measurements made along the northwest stake transect, we take the VLR model results as the best estimate of surface mass balance. We estimate that the main ice cap lost an average of 333 ±120mmw.e. a–1 over the period 1980–2006, where the error associated with this estimate is taken to be the standard error of the modeled, relative to measured, net mass balance for the northwest sector.
Discussion and Conclusion
The use of different air-temperature lapse rates in the downscaling of surface air temperatures over Devon Ice Cap has a large impact on estimates of glacier mass balance produced with the degree-day method. Using the moist adiabatic lapse rate (–6.5°Ckm–1) to downscale summit elevation temperatures over the ice cap generates net mass-balance estimates that are 4 times more negative than estimates made using a modeled variable lapse rate and ∽2.5 times more negative than estimates made using a mean measured summer lapse rate. Over the main ice cap, mass-balance estimates made from temperature fields downscaled with the mean measured summer lapse rate overestimate melt by ∽180±100mmw.e. a–1 on average relative to estimates made from temperature fields downscaled using a modeled variable lapse rate. During more positive mass-balance years, model results generated using the two methods agree closely, but they diverge by as much as 550 mmw.e. a–1 during the most negative mass-balance years.
When the DDM is forced with surface air-temperature fields downscaled using modeled variable daily lapse rates, it produces estimates of net surface mass balance that agree well (9±120 mm) with the 23 years of measurements made along the northwest transect. This same method produces a 1996 estimate of net mass balance for the main ice cap that agrees closely (42mmw.e. higher) with the previous 1996 estimate made by Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others (2007). However, forcing the model with these surface air temperatures only results in a 25% (80mmw.e. a–1) improvement in the standard error of point-location net mass-balance estimates relative to estimates made from temperature fields downscaled using a mean summer lapse rate. Like all DDMs, model results are highly sensitive to the choice of degree-day
factors. While we have chosen not to tune the model using the degree-day factors, model root-mean-square errors similar to those achieved when the DDM is forced with a variable lapse rate could be achieved by adjusting the degree-day factor for ice to ∽4.5–6.5 and ∽2–3.5mmw.e. d–1 when employing the mean measured summer lapse rate and the moist adiabatic lapse rate, respectively. When comparing against previously published values (Reference HockHock, 2003), a degree-day factor < 5.0mmw.e. d–1 for ice is unrealistic. Therefore, the model can only be tuned with a realistic degree-day factor for ice to produce comparable results to those achieved using the untuned variable lapse rate when the mean measured summer lapse rate is employed. The root-mean-square error across all three models is relatively insensitive to degree-day factor for snow between 2 and 6 mmw.e. d–1.
Our best estimate of the net surface mass balance of the main ice cap over the years 1980–2006 is –330 ±120mmw.e. a–1, which is much lower than previous estimates of –76±50 (1961–66 (Reference KoernerKoerner 1970)), –130±60 (1963–2000 (Reference Mair, Burgess and SharpMair and others, 2005)) and 28±70mmw.e. a–1 (1996 (Reference Shepherd, Du, Benham, Dowdeswell and MorrisShepherd and others, 2007)). The more negative model estimate can be attributed to differences in mean summer conditions between the different time periods for which estimates have been derived, and to large uncertainties in estimated accumulation rates due to large spatial variations in snowfall (Reference KoernerKoerner, 1966), which are poorly reproduced by the 32×32 km resolution NARR output. In particular, total accumulation amounts measured by Reference KoernerKoerner (1966) were found to be 2–6 times larger for the southeast (largest sector) than for the northwest sector of the ice cap and are nearly equivalent to the regridded NARR model output.
In conclusion, DDMs are highly sensitive to the choice of lapse rate when models are forced with downscaled temperature fields. For the DDM used in this study, use of a variable daily lapse rate estimated from lower-tropospheric (750 mbar) temperatures to downscale surface air temperatures gives significantly better mass-balance estimates than a constant lapse rate equal to either the summer mean or the moist adiabatic lapse rate.
Acknowledgements
This work was supported by the Natural Sciences and Engineering Research Council (NSERC Canada) (through Discovery, Northern Supplement and Equipment Grants to M. Sharp, and an Alexander Graham Bell Canada Graduate Scholarship to A.S. Gardner), the Canadian Foundation for Climate and Atmospheric Sciences through the Polar Climate Stability Network, Environment Canada through the Cryosphere System in Canada Project (CRYSYS), the Northern Scientific Training Program (Indian and Northern Affairs Canada), the Institute for Geophysical Research and Circumpolar Boreal Alberta Research Grant fund (University of Alberta) and by the Alberta Ingenuity Fund (through scholarships to A.S. Gardner). We thank the late R.M Koerner for making available surface mass-balance and temperature datasets.
