Temperature-index models are popular tools for glacier melt-modeling due to their minimal data requirements and generally favorable performance. We examine the effects of temperature forcing provenance and extrapolation on the performance of one such model applied to a small glacier in the Saint Elias Mountains of northwestern Canada. The model is forced with air temperatures recorded (a) on two glaciers, (b) at two nearby ice-free locations, and (c) by two low-elevation valley stations. We extrapolate these temperatures using constant lapse rates and assess model performance by comparing measured and modeled cumulative summer ablation at a network of stakes over five melt seasons. When the model is calibrated individually for each temperature forcing and lapse rate, the variation in model performance is modest relative to inter-annual variations associated with melt-season conditions and calibration data quality. Despite <30% variation in estimated summer ablation arising from the combined influences of temperature forcing and lapse rate, the resulting variations in estimated annual mass balance can be significant (>100% in some cases). While model parameters calibrated in this way suffer from error compensation and exhibit equifinality, the lapse rates associated with minimum model error exhibit inter-annual variation that can be related to prevailing meteorological conditions. When the model is instead calibrated at the point scale without employing a lapse rate, and the resulting parameters are paired with an arbitrary temperature forcing, lapse rates associated with minimum model error vary widely between forcing types and years. Low-elevation stations distal from the study site sometimes outperform the calibration station, but the prescribed lapse rate becomes critical in this case. With either calibration method, lapse rates that minimize model error for the valley stations are generally steeper than the measured environmental lapse rates.
Introduction
In light of projected climate changes and the uncertainties associated with estimates of attendant sea level rise (e.g. Lemke et al., 2007) and runoff from glacierized catchments (e.g. Moore et al., 2008; Immerzeel et al., 2010; Kaser et al., 2010; Pellicciotti et al., 2010), regional models of climate-driven glacier change are needed. Point-scale energy-balance studies of glacier melt provide useful insights into the detailed physical mechanisms linking glaciers and climate (e.g. Oerlemans and Klok, 2002; Pellicciotti et al., 2005), while spatially distributed melt models extend these relationships to larger scales (e.g. Brock et al., 2000b; Hock and Holmgren, 2005, MacDougall and Flowers, 2011; Petersen and Pellicciotti, 2011). To evaluate the impact of glacier melt on sea level and regional water resources, consideration of multiple glaciers and regional-scale glacier changes are necessary (e.g. Raper and Braithwaite, 2006; Arendt et al., 2008; Radić et al., 2013).
The relative scarcity of detailed meteorological data for input to regional glacier melt models presents a challenge. In the absence of the necessary data required for physically based energy balance models, temperature-index models of snow- and ice-melt have often been used. Classical temperature-index (or degree-day) models are based on an empirical relationship between melt and air temperature. This approach has been in use since the 19th century (Finsterwalder and Schunk, 1887). The popularity of degree-day models owes to their minimal data requirements, their conceptual and computational simplicity, and their generally good performance. Ohmura (2001) attributed this performance to the strong correlation between air temperature and incoming longwave radiation, one of the leading heat sources in the surface energy balance. The major drawbacks of the classical degree-day model are its decreasing accuracy with increasing temporal resolution, its inability to accurately model spatial melt variability (e.g. Hock, 2003), its tendency to drift away from physically based models in long-duration simulations (Hock et al., 2007), and its questionable spatial and temporal transferability (e.g. MacDougall et al., 2011).
Various modifications of the classical degree-day approach have been introduced in an attempt to address these limitations, most involving the introduction of terms related explicitly to radiation. Two notable examples are the models of Hock (1999) and Pellicciotti et al. (2005). Hock's (1999) model modifies the degree-day factor with a term dependent on the potential direct radiation, whereas the model of Pellicciotti et al. (2005) separates the radiation and temperature terms. Both of these “enhanced” temperature-index models have demonstrated improved performance over the classical degree-day model, especially at sub-daily resolution, and have been widely implemented to simulate glacier melt and runoff (e.g. Huss et al., 2007; Zhang et al., 2007; Ragettli and Pellicciotti, 2012). They also exhibit notably improved spatial and temporal transferability over the classical model (e.g. Carenzo et al., 2009; MacDougall et al., 2011).
Temperature-index models have been forced both with temperature records collected above the glacier surface (e.g. Braithwaite, 1995; Hock, 1999; Carenzo et al., 2009) and with those collected from surrounding ice-free locations (e.g. Jóhannesson et al., 1995; Huss et al., 2008; Shea and Moore, 2010). While physically based energy-balance models benefit from measurements of air temperature on the glaciers themselves (e.g. Gu
mundsson et al., 2006), it is common to find that temperatures measured outside the glacier boundary layer lead to better performance of temperature-index models under many circumstances (e.g. Lang, 1968; Lang and Braun, 1990; Gumundsson et al., 2009; Wheler, 2009). Due to the ability of a non-melting surface to be heated, temperature records collected in ice-free regions carry information about shortwave radiation, an important component of the energy balance (e.g. Ohmura, 2001). This explains at least some of the improvement when using temperatures measured distal from the glacierized region. Other studies have documented unrelated performance improvement with increasing distance between study glacier and meteorological station. Letréguilly (1988) found improved correlation between glacier mass balance and meteorological quantities with distance, which may have been the result of differences in data quality among stations. Data quality for a particular application (see Whitfield, 2012) is affected by station siting, among other factors. Stations sited in mountainous areas, as compared to valleys, may preferentially sample conditions more representative of the free atmosphere and be affected by local topography (e.g. Lang and Braun, 1990).
In addition to selecting an appropriate temperature forcing, extrapolation of temperature from station to study site must be considered. Although a standard environmental lapse rate of -6.5 K km-1 is often assumed and has been measured in mountainous areas (e.g. Blöschl, 1991), the use of a constant and uniform lapse rate in melt modeling has been increasingly called into question (e.g. Marshall et al., 2007; Gardner and Sharp, 2009; Petersen and Pellicciotti, 2011). Lapse rates are known to vary on many timescales, from daily to seasonal (see Minder et al., 2010, and references therein), and are well documented to be less steep over melting surfaces (within the boundary layer) (e.g. Greuell and Böhm, 1998; Marshall et al., 2007) and affected by the development of katabatic winds (e.g. Shea and Moore, 2010; Petersen and Pellicciotti, 2011; Jiskoot and Mueller, 2012). Though the variance of lapse rates has been shown to increase with wind speed (e.g. Blöschl, 1991), wind speeds below a certain threshold are generally not strongly correlated with glacier melt (e.g. Lang, 1968), in part because incoming longwave and net shortwave radiation are independent of atmospheric movement (e.g. Ohmura, 2001). Katabatic winds are an important exception to the general statement above, in that they tend to flatten the lapse rate across the glacier surface by transporting cold air to lower elevations. Katabatic winds develop with sufficient fetch and in the presence of diurnal heating, and result in reduction of outgoing longwave radiation due to surface cooling as well as increased turbulent energy fluxes (sensible and latent) due to the elevated windspeeds (e.g. Shea and Moore, 2010; Jiskoot and Mueller, 2012).
In this study we aim to identify what combinations of temperature forcing and lapse rate yield the best estimates of glacier melt, rather than to maximize the fidelity of air temperature extrapolation to the glacier surface directly. As demonstrated in previous studies, the temperature records with the most information content may not be those that produce the most accurate melt estimates (e.g. Blöschl, 1991). We use an enhanced temperature-index model applied over five melt seasons to a glacier on the continental side of the Saint Elias Mountains, Yukon, Canada (Fig. 1). We force the model with six temperature records: two from locations on glacier surfaces, two from adjacent ice-free mountain areas, and two from low elevation valleys. Cumulative summer ablation, measured at a network of stakes across the glacier, is used to calibrate the model with two different strategies. Though we have documented lapse rate variability on many timescales within our study region, here we restrict our focus to constant lapse rates and experiment with values ranging from -3.5 to -10.0 K km-1. This broad range includes the values most commonly adopted for glacier melt modeling, as well as the mean values we have measured during the melt season both on and off the study glacier.
Study Area
The Saint Elias Mountains, located in southeastern Alaska, southwestern Yukon, and northwestern British Columbia (Fig. 1), are characterized by high relief and host an unusual concentration of surge-type glaciers (e.g. Clarke and Holdsworth, 2002). This region provides a unique opportunity to study the relationship between glaciers and climate in the context of these two characteristics. Our primary field site is situated on the continental side of the Saint Elias Mountains, in the Donjek Range (60°50′N, 139°10′W; see Fig. 1, part b). The study glacier is polythermal (Wilson et al., 2013) and has a history of surging (Johnson and Kasper, 1992). It is ∼5 km long and 5.3 km2 in area, spans an elevation range of 1970–2960 m above sea level (a.s.l.) with an equilibrium line altitude of ∼2550 m a.s.l., and is situated to the southeast of the Donjek Range crest in a valley with predominantly southerly exposure. Energy-balance modeling suggests a radiation-dominated ablation regime, with net radiation supplying ∼80% of the melt energy in 2008 and 2009 (MacDougall and Flowers, 2011). Previous studies of this glacier have addressed the proglacial sedimentary record (Johnson and Kasper, 1992; Johnson, 1997), glacier dynamics (De Paoli and Flowers, 2009; Flowers et al., 2011), glacier thermal structure (Wilson et al., 2013), and energy-balance modeling (Wheler and Flowers, 2011; MacDougall and Flowers, 2011). For consistency with previous work, we hereafter refer to the study site as “South Glacier” (Fig. 1, part c).
Field Measurements and Derived Model Inputs
AIR TEMPERATURE AND LAPSE RATES
The Simon Fraser University Glaciology Group (SFUGG) maintains automatic weather stations (AWSs) year-round at four locations in the Donjek Range (Fig. 1, part b): at ∼2300 m a.s.l. in the ablation areas of South Glacier (SG; since July 2006) and North Glacier (NG; since May 2007), and at ∼2200 m a.s.l. in two ice-free locations in the Canada Creek (CC) and Duke River (DR) headwaters (since August 2006). Temperature is measured using a Campbell Scientific HMP45C212 TRH probe, for which the manufacturer reports ±0.1 °C thermistor absolute accuracy and interchangeability, and linearity of ±0.09 °C. Instruments are identical among stations and are installed at a nominal height of 2 m. Temperature sensors are shielded but not actively ventilated and are installed on the north side of the tripod mast. Data are recorded at 5-minute intervals using a Campbell CR1000 data logger and are averaged to hourly means for this study.
FIGURE 1.
Study site. (a) Saint Elias Mountains showing study area (box) and locations of Environment Canada weather stations at Burwash Landing (BL) and Haines Junction (HJ). (b) Detail of Donjek Range study area with AWS locations labeled: South Glacier (SG), North Glacier (NG), and ice-free locations near the headwaters of Canada Creek (CC) and Duke River (DR). (c) Contour map of South Glacier with locations of ablation stakes and AWS.
Two additional AWSs are operated by Environment Canada near the study area (Fig. 1, part a): one ∼60 km to the north in Burwash Landing (BL; 805 m a.s.l.) and one ∼90 km to the east in Haines Junction (HJ; 595 m a.s.l.). These stations meet World Meteorological Organization (WMO) standards for siting and exposure. Haines Junction (#2100630) is a Meteorological Service of Canada (MSC) Autostation and is inspected and calibrated every 6 months. Burwash Landing is a regional climate station and is inspected and calibrated annually, though records from three different Burwash Landing stations had to be concatenated for this study: “Burwash A” (#2100182; 806.2 m a.s.l.) for dates up to 3 July 2008, “Burwash AWOS” (#2100185; 805.3 m a.s.l.) for 3 July 2008–30 November 2011, and “Burwash” (#2100181; 805.3 m a.s.l.) for dates after 30 November 2011. We group the AWSs in this study by their geographic settings and describe them as: (1) “glacier” (SG, NG), (2) “mountain” (CC, DR), and (3) “valley” stations (BL, HJ). All stations used in this study are situated below the mean elevation of the study glacier.
Temperature forcings for the enhanced temperature-index model are derived from the six AWSs (Table 1) using a constant lapse rate. We test values from -3.5 to -10.0 K km-1 in increments of 0.5 K km-1. A value of -6.5 K km-1 has been adopted in many previous studies and is often cited as a “standard” value (see Minder et al., 2010, for a critical discussion). Other studies have documented improved model performance using lapse rates less steep than -6.5 K km-1 over glacier surfaces (e.g. Gardner and Sharp, 2009), consistent with the influence of the glacier boundary layer (e.g. Shea and Moore, 2010). We have measured a value of -6.4 K km-1 across the surface of South Glacier from May to September 2011 using four temperature microloggers (Onset HOBO H08-032-08) mounted on stakes drilled into the ice from elevations of 2050 to 2670 m a.s.l. We calculate mean summer lapse rates between South Glacier (AWS on ice) and the low elevation (ice-free) valley stations, Burwash Landing and Haines Junction, for melting conditions during 2007–2012 to be -6.5 and -6.0 K km-1, respectively. Although lapse rates vary on diurnal to inter-annual timescales, determination of time-variable lapse rates generally requires information from more than one station. In adopting a constant lapse rate, we mimic the typical situation of having only one meteorological station in the region.
WINTER BALANCE
The winter balance, taken as the water-equivalent snowpack depth at the end of the accumulation season, forms the initial condition for the enhanced temperature-index model. To determine the winter balance, snow depths were probed at 14–18 stake locations plus up to 16 other sites on South Glacier, and several snow pits were excavated across a range of elevations in May of 2007, 2008, 2009, 2011, and 2012 (see Table 2). Spring field work is intended to coincide with the end of the accumulation season (and hence peak snow depth); however, due to inter-annual variability this timing is only approximate. The winter balance at each stake is computed from these snow depth measurements and the integrated snow-pit densities. Poles above the highest pit and below the lowest pit use the integrated snow densities from the highest and lowest pits, respectively; poles between the highest and lowest pits use an average integrated density. For melt-model calibration runs, we use the winter balances calculated at the poles themselves to initialize the model. For the purpose of calculating glacier-wide mass balance, the winter balance bw is estimated across the glacier surface using the following linear regression:
with glacier surface elevation z (m), slope m (°), curvature c (m-1), and constants β0-β3. Slope and curvature are calculated numerically on a 30 m digital elevation model. The addition of curvature to the regression improves the fit in some years, but not all. See Table 3 for regression parameters and Table 4 for resulting winter balances.
TABLE 1
Automatic weather stations (AWS) operated by the Simon Fraser University Glaciology Group (SFUGG) and Environment Canada (EC) used in this study. Elevations correspond to AWS locations.
SUMMER ACCUMULATION
Summer snow events can cause large reductions in short-term melt rates by increasing the glacier surface albedo and reducing the surface roughness length (e.g. Brock et al., 2000a). Accumulation events on South Glacier are identified using an ultrasonic depth gauge record (SR50 Sonic Ranger, ±0.4%) taken from the South Glacier AWS. The record is processed to daily averages and accumulation identified as a decrease in distance-to-surface. Finer-than-daily resolution for the accumulation record is precluded due to noise in the data. Sub-daily accumulation events that are offset by a greater amount of daily melt therefore go undetected.
We estimate summer accumulation across the glacier surface by extrapolating the accumulation measured by the ultrasonic depth gauge using a constant accumulation—elevation lapse rate of 2 cm w.e. km-1. This value is the average of three values, each determined by regression of snow accumulation on elevation as measured at South Glacier on three different days in July between 2006 and 2011. In each case, snow depth and density were measured at stake locations during or immediately after an accumulation event. Adopting accumulation—elevation lapse rates of 1 and 3 cm w.e. km-1, rather than the reference value of 2 cm w.e. km-1, leads to differences in total estimated summer accumulation of up to 25% of the reference value. This is a potentially significant source of uncertainty, not just for estimates of summer accumulation, but for estimates of ablation and therefore net balance. The overall importance of summer accumulation is evident when compared with winter balances (Table 4), with both quantities being on the order of decimeters: 0.17–0.33 m w.e. for summer accumulation compared to 0.34–0.57 m w.e. for winter balance. The comparable magnitudes of summer accumulation and winter balance arise, in some cases, due to significant accumulation events in May and June after the onset of melt.
TABLE 2
Dates of stake measurements (observations) and melt-season simulations, along with number of stakes measured and used in model initialization (# spring stakes), number of stakes measured in mid-summer (# summer stakes), and number of stakes with records spanning full-season simulation period (# stakes with full record). Site visits were planned to roughly bracket the melt season, except in the case of 2009. Simulation start dates are dictated by the dates of the early-season observations. Model tuning is performed for the period defined by the early- and late-season observation dates, though the simulation period may be longer.
TABLE 3
Winter-balance regression parameters determined from Equation (1).
SUMMER ABLATION
Cumulative summer ablation at up to 18 stake locations on South Glacier is derived from measured stake heights and measured or assumed surface density. An estimate of uncertainty is made for each stake height measurement using a 30 cm × 30 cm PVC plate placed at the base of the stake. Stake measurements are made from the plate to the top of the stake, while surface height variability under the plate is estimated by eye as a measure of local uncertainty. This is combined with uncertainty in surface density to estimate the uncertainty in cumulative surface ablation at each stake. Uncertainty associated with the representativeness of point ablation stake measurements within the digital elevation model (DEM) grid cells is not quantified. We assume ice and firn densities of 900 and 550 kg m-3, respectively, though all calculations have been repeated with an assumed firn density of 800 kg m-3 for comparison. We measure snow density in the field, but where data are missing we assign a value of 200 kg m-3 based on the mean value of measured new-snow densities at our field site.
Stake height and snow density were measured near the beginning of the melt season in all years included in this study, and near the end of the melt season in all years but 2009 (Table 2). Stakes are also measured several times in the middle of the melt season, most often in July. The cumulative summer surface ablation across the stake network is used for calibration of the enhanced temperature-index model.
DEMS AND POTENTIAL DIRECT SOLAR RADIATION
A DEM of the glacier surface was constructed from real-time kinematic (RTK) Global Positioning System (GPS) measurements and digitized map contours for the surrounding topography (De Paoli, 2009; Wheler, 2009). Potential direct solar radiation is computed at half-hourly intervals using the ArcGIS Solar Analyst tool and knowledge of the surrounding terrain (Fu and Rich, 2000). The potential radiation values are summed to obtain hourly values.
Enhanced Temperature-Index Model
We adopt the temperature-index model of Hock (1999) because of its widespread use within the community and its minimal input requirements. Simulated hourly melt rates Ms in each DEM grid cell are calculated as follows:
with n = 4 (h day-1), MF a temperature melt factor similar to the degree-day factor (mm w.e. day-1 °C-1), asnow/ice the radiation melt factor for snow or ice (mm w.e. h-1 °C -1 m2 W-1), I the potential direct clear-sky solar radiation (W m-2), and T the air temperature (°C). The value of I varies in time and space due to the combined effects of the position of the sun, altitude, surface slope, surface aspect, and surrounding topography. Hourly, as opposed to daily, model time-steps have the advantage of resolving diurnal cycles in temperature and radiation, which can lead to positive values of melt energy even on days with sub-freezing mean temperatures. A constant firn-line elevation of 2450 m a.s.l. is assigned based on field observations. Above this firn line, snow is assumed to be arbitrarily deep and no snow-to-ice transition can occur. When a summer snowfall event occurs, the formerly ice-exposed grid cells take on the snow radiation melt factor until the snow is ablated.
TABLE 4
Glacier-wide mass balance components for South Glacier in m w.e.: winter balance (Bw) estimated by extrapolation of stake measurements, summer accumulation (Cs) estimated by extrapolation of processed ultrasonic depth gauge record assuming an accumulation—elevation lapse rate of 2 cm w.e. km-1 (values in parentheses show range with lapse rates of 1–3 cm w.e. km-1), summer ablation (As) estimated by model tuning with MAE (see also Fig. 3, bottom row), and annual mass balance Bn = Bw + Cs + As. Terminology and notation taken from Cogley et al. (2011). The range of As represents simulations forced with lapse rates from -4.5 K km-1 to -7.0 K km-1 for glacier (SG, NG), mountain (CC, DR), and valley (HJ, BL) AWS temperature records.
MODEL TUNING AND EVALUATION
We adopt two different approaches to model calibration in this study. In the first (henceforth “CAL1”), model parameters are derived for each year, temperature forcing, and lapse rate by minimizing the discrepancy between simulated and measured cumulative melt at the ablation stakes. This minimization is carried out using two different error metrics: the mean absolute error (MAE),
and the relative root mean squared error (RMSE), where Msi and Mmi are the simulated and measured cumulative melt, respectively, at the ith ablation stake,
is the mean measured melt at all ablation stakes, and N is the number of ablation stakes. MAE and RMSE are normalized here by
following Hock (1999) and are also known as the relative error and the relative standard deviation, respectively.A brute force method is used to calibrate the model wherein a range for each parameter is defined (e.g., asnow/ice = 0–30 × 10-3 mm w.e. h °C -1 m2 W-1, MF = 0–10 mm w.e. day-1 °C) and sampled using increments of 0.5 × 10-3 mm w.e. h °C -1 m2 W for the radiation factors and 0.25 mm w.e. day-1 °C -1 for the melt factor. Rather than attempting to define a restricted and physically plausible range of parameter values (based, for example, on previous literature), we simply sample a large range and require only that the parameter values be positive. Each combination of parameter values is tested, and the “tuned” model defined by the parameters that minimize either MAE or RMSE. The precise range of parameters tested varies such that the minimum error does not lie along the upper edge of parameter space (i.e. highest value of MF or asnow/ice). The parameter set associated with minimum MAE may be different than that associated with minimum RMSE, thus there can be more than one tuned model. This calibration exercise is repeated using stake data from all years simultaneously to determine a “master” set of model parameters. In making use of all stake data, the tuning method described above uses information collected over a broad region of the glacier. In calibrating the model for each forcing, no assumption is made about the quality of individual temperature records. The disadvantage of this method is that it invites error compensation by allowing simultaneous variation of the temperature forcing, lapse rate, and model parameters.
In the second approach to model calibration (henceforth “CAL2”), parameters are derived for individual years using only the local temperature forcing (SG) and a single co-located melt record. Our ablation data restrict us to the years 2008, 2009, and 2012 for CAL2. The parameter-space sweep and error minimization are done exactly as described above for CAL1. The advantage of CAL2 is that no lapse rate is required in the calibration. The disadvantage is that it relies on single records of temperature and ablation for each year, thus assuming that the calibrated parameters are intrinsic to the site (rather than a function of the temperature forcing) and representative of the glacier as a whole.
Results and Discussion
CAL1: MODELED AND MEASURED CUMULATIVE ABLATION AT STAKE LOCATIONS
Comparisons between cumulative ablation modeled and measured at all stake locations for South Glacier are shown in Figure 2 for each year, each temperature forcing, and lapse rates of -4.5 to -7.0 K km-1, for CAL1 models tuned using MAE. Visual inspection of the different panels in Figure 2 reveals similar patterns of model—data mismatch across all forcings in a given year. This mismatch is quantified as MAE in Figure 3, where modeled glacier-wide summer ablation is also shown for each year, forcing, and lapse rate. Results using RMSE for model tuning are visually similar. For the simulations shown in Figure 2, the inter-annual variation of MAE (Fig. 3) is generally greater than the variation of MAE as a function of temperature forcing or of the lapse rates presented. This result corroborates what has been implicitly shown in other studies: that the location of the temperature station has little bearing on simulation of cumulative ablation (e.g. Jóhannesson et al., 1995), and goes some way toward explaining the success of temperature-index models even when forced with temperature records collected tens of kilometers from the study glaciers (e.g. Huss et al., 2008).
In Figure 3 the lowest values of MAE (≤0.12) are achieved for summers 2009 and 2012, while MAEs for 2011 are all higher than 0.14; MAEs for 2007 and 2008 are intermediate, between 0.11 and 0.17. Using RMSE rather than MAE as a tuning metric produces the same qualitative pattern of errors between years. Although it might be reasonable to expect, this inter-annual variation is not strictly related to the magnitude of cumulative summer ablation (Fig. 3, bottom row). For example, the lowest errors are found for the years of highest (2009) and lowest (2012) cumulative ablation.
Though inter-annual variation in MAE dwarfs that due to forcing type (glacier, mountain, valley) (Fig. 3, top row), there is still some variation in MAE with forcing. Glacier stations do not necessarily produce the lowest errors, as illustrated in 2011 where the local forcing, SG, produces higher values of MAE than the valley and mountain stations for some lapse rates (Fig. 3, top row). An inspection of this temperature record does not reveal any obvious explanation, but the correlation between SG and NG temperature records is slightly lower for 2011 than for other years (e.g. R = 0.914 for 2011 versus R = 0.925 for 2008).
FIGURE 2.
Comparison of modeled and measured cumulative summer ablation As at South Glacier stake locations (each symbol represents an individual stake). Each column represents a different melt season (left to right: 2007, 2008, 2009, 2011, 2012) and each row a different temperature forcing (top to bottom: SG, NG, CC, DR, BL, HJ). The top two, middle two, and bottom two rows represent forcings from “glacier” (circles), “mountain” (triangles), and “valley” (squares) AWSs, respectively. Lapse rates are indicated by symbol shading according to the legends at right. Modeled values plotted here were tuned by minimizing mean absolute error (MAE). Propagated uncertainties in measured ablation are not plotted, in the interest of clarity.
The variation in MAE produced by the range of lapse rates shown in Figure 3 is generally greater than that produced by the different forcings. The largest variations in MAE as a function of lapse rate usually occur with the valley stations, an intuitive result considering the significant difference in elevation between these stations and the study site (see Table 1). To further explore the dependence of MAE on lapse rate, we tabulate the lapse rates that produce minimum MAE for each forcing and year (Table 5). When averaged over years and forcings, it becomes clear that inter-annual variations in lapse rate are greater than variations with forcing type using the CAL1 approach: 2009 is characterized by the shallowest (least negative) values, 2008 and 2011 by the steepest (most negative) values, and 2007 and 2012 by intermediate values. Averaged across all years, MAE is minimized with steeper lapse rates for valley stations (-6.7, -7.0 K km-1) compared to glacier stations (-5.8, -6.2 K km-1). The mean optimal lapse rate for mountain station CC lies between those of the two valley stations, while that for mountain station DR lies between those of the two glacier stations.
FIGURE 3.
Mean absolute error (MAE) between modeled and measured cumulative summer ablation for South Glacier stake network (top row) (see Fig. 2) and glacier-wide cumulative summer ablation estimated with tuned model (bottom row). Each symbol represents a simulation for the entire stake network (top row) or the entire glacier (bottom row). Each column represents a different melt season (left to right: 2007, 2008, 2009, 2011, 2012). Each cluster of symbols within each panel represents a different temperature forcing plotted in order of increasing distance from the study site: South Glacier (SG), North Glacier (NG), Canada Creek (CC), Duke River (DR), Burwash Landing (BL), and Haines Junction (HJ). Forcing type is coded with symbols: “glacier” (circles), “mountain” (triangles), and “valley” (squares), while lapse rates are indicated by symbol shading according to the legend (in panels at right).
Using RMSE, rather than MAE, to tune the model yields similar but not identical results to those described above (not shown). The implication of the choice of lapse rate for the calculated glacier-wide summer ablation, in many cases, is modest and independent of whether MAE or RMSE is used for tuning in the CAL1 approach.
In Figure 3, model error is dominated by inter-annual variation over variation with forcing or lapse rate. We speculate that this is rooted at least partly in data quality, which is itself indirectly a function of mass-balance conditions. For example, the high MAEs in 2011 are related in part to stakes being buried or lost, and snow depth measurements being made in the vicinity of stake locations rather than precisely at the stakes themselves. In contrast, stakes not found in spring 2012 were immediately replaced and measured. Data quantity may also be a factor, with 2011 having the fewest stakes (12) found during the spring measurement period and the fewest stakes (6) with complete melt-season records of any year (see Table 2). Interestingly, 2009 is the year with the greatest number of stakes (15) having a complete melt-season record, and also the year with the lowest model MAE. Corroboration of this explanation for MAE variability would require further (e.g. cross-validation) analysis.
CAL1: MODELED GLACIER-WIDE ABLATION
In estimating cumulative glacier-wide summer ablation (Fig. 3, Table 4), the model must be initialized with extrapolated (rather than measured) winter balance as described above. The modeled cumulative summer ablation typically varies by less than ∼0.2 m w.e. as a function of lapse rate between -4.5 and -7.0 K km-1, for a given temperature forcing type (Table 4); valley and mountain stations produce variations of 0.08–0.21 m w.e. and 0.05–0.17 m w.e., respectively, while glacier stations produce variations of 0.04–0.15 m w.e. in all years but 2009 (0.24 m w.e. variation).
Estimated cumulative summer ablation varies over a similar range when lapse rate is held fixed and temperature forcing is varied. For a given lapse rate between -4.5 and -7.0 K km-1 in a given year, this variation with temperature forcing is 0.06–0.26 m w.e. (not separated out in Table 4). The collective variations in cumulative summer ablation illustrated in Figure 3 (bottom row) can lead to significant variations in estimated net balance (Table 4). The combined influences of lapse rate (between -4.5 and -7.0 K km-1) and forcing type yield estimated net balances that differ by up to 0.34 m w.e. (-0.61 to -0.95 m w.e. for 2007); 2011 is particularly poor in a relative sense, with net balance estimates from -0.20 to -0.46 m w.e.
Note that the differences above arise in the presence of uniformly calculated values of winter balance and summer accumulation (Bw and Cs in Table 4). Using 2011 for illustration, had we assumed precipitation lapse rates of 1–3 cm w.e. km-1 (rather than the reference value of 2 cm w.e. km-1), estimated cumulative ablation with valley station forcings would range from -0.94 to -1.34 m w.e. (rather than -0.94 to -1.13 m w.e.) and the corresponding net balances from -0.21 to -0.73 m w.e. (rather than -0.27 to -0.46 m w.e.).
FIGURE 4.
Comparison of model parameters MF (melt factor, top row), asnow (radiation factor for snow, second row), aice (radiation factor for ice, third row), and modeled glacier-wide cumulative summer ablation As (last two rows) for models tuned to South Glacier data using relative root mean square error (RMSE) versus mean absolute error (MAE). Each symbol represents a simulation for the entire stake network (top three rows) or the entire glacier (bottom two rows), and columns 1–5 each represent a different melt season (left to right: 2007, 2008, 2009, 2011, 2012). The right-most column (boxed) uses data from all years simultaneously to tune model parameters; the last row (boxed) shows glacier-wide cumulative summer ablation for each year modeled with the resulting parameters. Forcing type is coded with symbols: “glacier” (circles), “mountain” (triangles), and “valley” (squares), while lapse rates are indicated by symbol shading (legends in column 1 apply to all panels). Units given on y-axes also apply to x-axes.
In the absence of an impractically dense stake network, some methodology is required to estimate glacier-wide mass balance from a series of point measurements. The modeling procedure used here could be considered an informed interpolation/extrapolation technique, highlighting the uncertainty in values that may be reported as “observed” glacier-wide mass balance.
CAL1: MODEL PARAMETER VARIABILITY
The transferability of a model in space or time is partially a function of the applicability and stability of its parameters (e.g. Saltelli et al., 2004; Carenzo et al., 2009). The common problem of equifinality (e.g. Beven, 1993) afflicts our modeling procedure and results in different combinations of model parameters yielding equally good results (Fig. 4) as has been noted elsewhere (e.g. Hock, 1999; Schuler et al., 2007; Carturan et al., 2012). Parameters in Figure 4 vary with year, temperature forcing (SG, NG, CC, DR, BL, HJ), lapse rate (-4.5 to -7.0 K km-1), and tuning metric (MAE, RMSE). When considering all five melt seasons, inter-annual parameter variation stands out and could be a result of inter-annual differences in the surface energy balance (e.g. MacDougall and Flowers, 2011), combined with parameter equifinality. The approach used in CAL1 invites error compensation, thus making it difficult to ascribe meaning to the calibrated parameters.
One way to minimize parameter variation is to tune the model with data from multiple years. This approach not only increases the volume of calibration data, but provides a wider sampling of surface energy-balance conditions. Model parameters tuned with data from all five melt seasons (last column, Fig. 4) exhibit less variability, particularly for asnow and aice. Values of MF still span a wide range as a result of permitting lapse rate and model parameters to vary together in CAL1. Despite the general scatter of MF values, however, values of all three model parameters are tightly clustered for both glacier stations (SG, NG) (circles in last column of Fig. 4). Glacier-wide cumulative summer ablation simulated with parameters calibrated using data from all years varies more than that simulated with parameters calibrated individually for each year (compare bottom two rows of Fig. 4), though the variation is less pronounced if considering only one forcing type (glacier, mountain, valley).
For individual glacier and mountain stations, lapse rates produce little of the parameter variation shown in Figure 4. Lapse rate does play a significant role, however, in the variation of MF using valley stations, steeper lapse rates requiring higher values of MF to achieve the same cumulative melt (see squares especially in 2009–2012 in top row of Fig. 4). This is consistent with the greater influence of lapse rate on forcings derived from low elevations. The forcing itself is responsible for some of the parameter variation seen in Figure 4, but there is no clear and persistent pattern in the parameter values as a function of forcing type (glacier, mountain, valley).
Differences in parameter values for models tuned with MAE versus RMSE are pronounced in 2008 for all stations, and in 2011 for mountain stations. Higher melt factors (MF) and correspondingly lower radiation factors (asnow/ice) are obtained for models tuned with MAE compared to RMSE (Fig. 4); RMSE tuning even produces MF = 0 in many cases in 2008. When MF = 0, shortwave radiation plays a stronger role in determining spatial variations in ablation as compared to temperature. The low winter balance in 2008, combined with the dependence of the winter-balance distribution on surface curvature (β3 in Table 3), would lead to a significant role for albedo in determining the spatial variations in cumulative ablation for 2008. An ablation field strongly influenced by albedo variations is likely to exhibit significant spatial heterogeneity, giving rise to high variance in the data. Because RMSE weights large errors (including outliers) more heavily than MAE, model tuning using RMSE versus MAE is most likely to differ for data sets with high variance. Year 2008 provides a stark example of parameter equifinality whose cause remains difficult to pinpoint. However, the variation of modeled cumulative summer ablation with forcing and lapse rate in 2008 is unusually low (see tightly clustered symbols in Fig. 4, fourth row, column 2008).
Note that the results above apply strictly to the enhanced temperature-index model used here, which has a radiation component. We have found valley stations to outperform both glacier and mountain stations when using the classical degree-day model, owing to the higher amplitude variations in air temperature recorded by the valley stations (Wheler, 2009). This result has been obtained in other environmental settings (e.g. Gumundsson et al., 2009) and has a theoretical basis: air temperatures measured over non-melting surfaces contain information related to incoming longwave and shortwave radiation, as well as sensible heat, thus capturing the most important terms in the surface energy balance (e.g. Ohmura, 2001). It also suggests that our mountain stations may be situated within the regional glacier boundary layer (Wheler, 2009; Shea and Moore, 2010), which seems plausible given the southwesterly synoptic flow in the region (Moore et al., 2002).
The inter-annual variation in tuned model parameters cautions against confident application of the model to melt seasons for which calibration data are not available. This result contrasts with the general conclusions of some previous work demonstrating good spatial and temporal melt-model transferability (e.g. Carenzo et al., 2009; Shea et al., 2009). However, differences in model formulation, the provenance and extrapolation of air temperatures, the quality, quantity, nature, and processing of mass balance data, parameter tuning procedures, the physical environment, and the spatial and temporal scales over which studies are conducted may explain some of the differences in study outcomes. Our results also suggest that a more robust model tuning procedure than we have used here is warranted, if the intention is to apply the model without recalibration to other sites or years. Such a procedure might involve using additional data metrics for calibration (e.g. snowline retreat; Schuler et al., 2007) and/or employing a more sophisticated optimization scheme (e.g. Heynen et al., 2013). Here we have chosen not to tune the model to output from another model (e.g. hourly ablation rates from an energy-balance model) (c.f. MacDougall et al., 2011), though there is precedent for this strategy (e.g. Pellicciotti et al., 2005) as it can be used to improve the temporal model performance. At the very least, making use of data from multiple melt seasons in the model calibration leads to some improvement in parameter stability.
CAL2: MODEL ERROR AS A FUNCTION OF FORCING AND LAPSE RATE
Error compensation should be reduced in CAL2 compared to CAL1 by employing a single set of model parameters for each year, leaving temperature forcing and lapse rate as the only variables. Model errors are generally higher for CAL2 than CAL1, as is expected without parameter recalibration (Fig. 5). Although the calibration station SG produces similar errors in all three years (MAE < 0.2), it only performs best of all stations in one year (2008). Stations NG, HJ, and BL outperform SG in 2009, while NG, CC, and BL outperform SG in 2012. The highest errors are produced by NG in 2008 (MAE ∼ 0.3) and CC in 2009 (MAE > 0.35).
With the exception of Canada Creek (CC) in 2008 and 2009, glacier and mountain forcings produce a relatively flat response to MAE with lapse rate, as expected for stations close in elevation to the study site. By contrast, the choice of lapse rate is critical with the low-elevation valley stations (HJ, BL). The lapse rates that minimize MAE in CAL2 are somewhat consistent between forcing types in a given year, but vary markedly between years. For example, lapse rates between -5.5 and -6.0 K km-1 minimize MAE with glacier forcings in 2009, and those between -7.5 and -8.0 K km-1 minimize MAE with valley forcings; in 2012 the situation is reversed: -8.0 and -9.5 K km-1 minimize MAE with glacier forcings, and -6.0 and -6.5 K km-1 minimize MAE with valley forcings. The use of a single temperature station and melt record in the calibration makes CAL2 vulnerable to errors in these records. Any such errors will be expressed in the lapse rates in CAL2, rather than being partitioned between lapse rates and model parameters as in CAL1.
FIGURE 5.
Mean absolute error (MAE) between modeled and measured cumulative summer ablation for South Glacier stake network for 2008 (left), 2009 (center), and 2012 (right) as a function of prescribed lapse rate. Forcing type is coded with line type: bold solid (glacier: SG, NG), bold dashed (mountain: CC, DR), and fine solid (valley: BL, HJ). Legend in center panel applies to all panels. Model is calibrated with local (SG) forcing and a single melt record collected at SG AWS location in all simulations shown.
INTERPRETATION OF TUNED LAPSE RATES
Lapse rates that minimize model error using valley-station forcings (Table 5 for CAL1, Fig. 5 for CAL2) can be compared to those calculated directly between valley stations (HJ, BL) and the South Glacier AWS (Fig. 6). In all but one case, the lapse rates that minimize model error are steeper (more negative) than the measured lapse rates averaged over the model calibration period (compare values in Table 5 for CAL1 and Fig. 5 for CAL2 to measured values in Fig. 6). The model calibration (for both CAL1 and CAL2) thus generally demands lower extrapolated temperatures than would be predicted with measured lapse rates between the valley stations (HJ, BL) and the study site (SG). This result would be consistent with systematically overestimated air temperatures at the SG AWS, perhaps as a result of inadequate sensor ventilation.
Despite the systematic differences between measured lapse rates and those that minimize model error, the inter-annual variation in tuned lapse rates for valley stations with CAL2 follows that of the measured variations, with 2009 having the steepest values, followed by 2008, and then by 2012. This pattern of variation also holds for mountain stations with CAL2 but not for glacier stations, for which the shallowest lapse rates are found in 2009 as in CAL1. Year 2009 with CAL1 is the exception to the pattern described above, with measured lapse rates being steeper than those that minimize model error for valley stations. The model calibration therefore imposes higher extrapolated air temperatures than would otherwise be predicted at the study site with the measured lapse rates. This hints at a requirement for additional melt energy in 2009 for CAL1.
We speculate that the inter-annual variation in lapse rates that minimizes model error in both the CAL1 and CAL2 approaches, though different from one another, reflects real differences in prevailing meteorological conditions. Figure 7 presents mean values of meteorological variables measured at the South Glacier AWS for the calibration and simulation periods in each year (see Table 2 for dates). Year 2009 stands out as being characterized by high mean air temperature, high wind speed, high radiation, and low relative humidity. This warm, dry, and windy year saw the greatest cumulative ablation of all years and is associated with the least negative lapse rates for CAL1 across all forcings (Table 5) and for CAL2 with glacier forcings. Years 2008 and 2011, by contrast, were the two coolest, with moderate mean wind speeds and moderate values of radiation; these two years are associated with the steepest lapse rates for CAL1 (Table 5). The structure of the inter-annual lapse rate variation for CAL1 across all forcings resembles the overall structure of the mean SG air temperature in Figure 7, part a.
TABLE 5
Lapse rates associated with minimum MAE for all CAL1 simulations (all values in K km-1) along with mean values by year (last row) and forcing (last column). Note that Figures 2, 3, and 4 present results only for lapse rates from -4.5 K km-1 to -7.0 K km-1, while the range of lapse rates tabulated here is broadened to capture the minimum model error in each case.
FIGURE 6.
Mean monthly lapse rates between South Glacier AWS (SG) and valley stations Burwash Landing (BL) and Haines Junction (HJ) calculated for melting conditions (T > 0 °C) at SG. Mean values for the model calibration period are listed for each year.
Examination of hourly lapse rate variations measured across South Glacier over one melt season (unpublished data, Simon Fraser University Glaciology Group) shows a strong diurnal variation in lapse rate, with the most negative values occurring in the late afternoon while (downslope) wind speeds are at a maximum (not shown). This observation suggests katabatic winds may play a role in the energy-balance regime of this glacier, as has been found elsewhere (e.g. Shea and Moore, 2010), and may explain some of the lapse rate variations with forcing type in Table 5 (CAL1), as well as the anomalous lapse rates of 2009. Although usually associated with surface cooling, katabatic winds enhance turbulent energy exchange at the surface that can lead to enhanced ablation (e.g. Jiskoot and Mueller, 2012). The high mean air temperatures and elevated wind speeds in 2009 (Fig. 7) may indicate a role for katabatic flows in contributing to the high rates of summer ablation, and therefore the need for higher extrapolated air temperatures in the model for this year across all forcings for CAL1 and with glacier forcings for CAL2.
Summary and Conclusions
We have examined the effect of temperature forcing provenance and lapse rate on the performance of an enhanced temperature-index model using six different temperature stations, a large range of lapse rates, and two different approaches to parameter calibration. The temperature stations span a range of elevations and represent both ice-free and ice-covered areas. The lapse rates extend well beyond the range of mean measured values during melt-season conditions in the study area, both over ice and over ice-free terrain. The purpose of this study is to inform the choice of temperature forcing and lapse rate in the context of empirical melt modeling, rather than to examine the most faithful means of reproducing air temperatures measured over glaciers. The latter question is a separate one that is particularly relevant to energy-balance modeling.
We find that:
(1) The source of the temperature forcing has only a modest impact on the simulated cumulative ablation, provided that model parameters are calibrated individually with each forcing. If model parameters are calibrated with one forcing and then used with other forcings without recalibration, model error generally increases and distal stations sometimes outperform the calibration station.
(2)The derived model parameter values exhibit inter-annual variation, and in two years of five, vary significantly with the error metric used for model tuning (MAE or RMSE). Temperature-forcing provenance and lapse rate contribute to parameter variation, but not in a clear and systematic way. Equifinality is evident in the model tuning, and points to the need for a more robust parameter optimization scheme. Tuning with data from all years simultaneously produces more consistent parameter values, particularly within a given forcing type, with some compromise in the fidelity of simulated cumulative ablation for individual years.
(3) The influence of lapse rate is enhanced when using temperatures from low-elevation (“valley”) stations to drive the melt-model, but modest for stations located near, and at similar elevations to, the study site (“glacier” and “mountain” stations). When model parameters are recalibrated for each forcing, lapse rates that minimize model error average to -5.8 and -6.3 K km-1 for glacier stations, and -6.7 and -7.0 K km-1 for valley stations. When model parameters are used without recalibration, the choice of lapse rate becomes critical for low-elevation stations. Inter-annual variation of model lapse rates is greater than variation with forcing type, and appears related to prevailing surface and meteorological conditions.
(4) Despite the modest variations in simulated cumulative ablation introduced individually by temperature-forcing provenance and lapse rate, together these factors can produce differences in calculated glacier-wide ablation that may be significant. In any given year, the estimated glacier-wide cumulative summer ablation varied by 0.11–0.35 m w.e. or by 11–26%, depending on the temperature forcing and lapse rate employed. This uncertainty can lead to large variations in the estimated annual net balance, for example, from -0.20 to -0.46 cm w.e. in 2011.
FIGURE 7.
Meteorological variables measured at South Glacier AWS (SG) averaged over calibration (cal) and simulation (sim) periods. (a) Air temperature. (b) Relative humidity. (c) Wind speed. (d) Incoming shortwave radiation. (e) Net all-wave radiation.
The results of this study suggest that temperatures collected within a few kilometers of the site of interest do not necessarily lead to clear improvements in the simulation of cumulative summer ablation, compared to temperatures collected within 100 km at much lower elevations. This, combined with the modest variation of modeled ablation as a function of the lapse rates tested here, points to some flexibility in choosing a forcing for empirical melt models provided adequate calibration data exist. Our results caution against adopting off-the-shelf model parameters and using them with arbitrary temperature records and prescribed lapse rates, particularly if the temperature records come from elevations much different than the study site.
The conclusions above hold for our subarctic study site and its elevation range, but may not necessarily generalize to other environmental settings. Furthermore, they apply specifically to the simulation of cumulative summer ablation. Simulation of hourly melt rates (not examined here) may lead to different conclusions related to the value of various temperature forcings (e.g. Wheler, 2009) and the treatment of lapse rate (e.g. Petersen and Pellicciotti, 2011). While we have focused our attention on constant lapse rates here, other studies have highlighted the potentially important effects of spatially (e.g. Minder et al., 2010) and temporally variable lapse rates (e.g. Gardner and Sharp, 2009) on melt-model performance. Devising means of including these variations, without generating onerous increases in data requirements, is a productive direction for future research.
Acknowledgments
We are grateful for funding from the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Foundation of Innovation, the Canada Research Chairs program, the Northern Scientific Training Program, the Polar Continental Shelf Program, Simon Fraser University (SFU), and the SFU Community Trust Endowment Fund through the Climate Change Impacts Research Consortium. Permission to conduct this research was granted by the Kluane First Nation, Parks Canada, and the Yukon Territorial Government. Support from the Kluane Lake Research Station (KLRS) and Kluane National Park and Reserve is greatly appreciated. We are indebted to Andy Williams, Sian Williams, Lance Goodwin (KLRS), Doug Makkonen and Steve Soubliere (Trans North Helicopters) for logistical support, and to many others for field assistance including Laurent Mingo, Alex Jarosch, Patrick Belliveau, Flavien Beaud, Nat Wilson, Mauro Werder, and Christian Schoof. Discussion with Valentina Radić, and critical comments from Hester Jiskoot, Martin Funk, and one anonymous reviewer led to an improved manuscript.
References Cited
mundsson
,
H. Björnsson
,
F. Pálsson
, and
H. Haraldsson
,
2006: Energy balance of Brúarjökull and circumstances leading to the August 2004 floods in the river Jökla, N-Vatnajökull.
Jökull
, 55: 1–18. Google Scholar
mundsson
,
H. Björnsson
,
F. Pálsson
, and
H. Haraldsson
,
2009: Comparison of energy balance and degree-day models of summer ablation on the Langjökull ice cap, SW-Iceland.
Jökull
, 59: 1–18. Google Scholar
