Modeling melt from glaciers is crucial to assessing regional hydrology and eustatic sea level rise. The transferability of such models in space and time has been widely assumed but rarely tested. To investigate melt model transferability, a distributed energy-balance melt model (DEBM) is applied to two small glaciers of opposing aspects that are 10 km apart in the Donjek Range of the St. Elias Mountains, Yukon Territory, Canada. An analysis is conducted in four stages to assess the transferability of the DEBM in space and time: 1) locally derived model parameter values and meteorological forcing variables are used to assess model skill; 2) model parameter values are transferred between glacier sites and between years of study; 3) measured meteorological forcing variables are transferred between glaciers using locally derived parameter values; 4) both model parameter values and measured meteorological forcing variables are transferred from one glacier site to the other, treating the second glacier site as an extension of the first. The model parameters are transferable in time to within a <10% uncertainty in the calculated surface ablation over most or all of a melt season. Transferring model parameters or meteorological forcing variables in space creates large errors in modeled ablation. If select quantities (ice albedo, initial snow depth, and summer snowfall) are retained at their locally measured values, model transferability can be improved to achieve ≤15% uncertainty in the calculated surface ablation.
Melt models are valuable tools for assessing the impact of climate change on glacier mass balance, changes in regional hydrology, and eustatic sea level rise (e.g., Hock 2005; Lemke et al. 2007). Glacier melt models have been applied over vast regions—even globally—by calibrating the models to well-studied glaciers and applying them unchanged over much larger areas (e.g., De Woul and Hock 2005; Oerlemans et al. 2005; Raper and Braithwaite 2006; Schneeberger et al. 2003). Both physically based melt models, which are derived from a surface energy balance, and empirically based melt models, which correlate melt with temperature and in some instances radiation, have been developed for glacierized regions (e.g., Hock 2005, 2003; Pellicciotti et al. 2005). In practice, physically based melt models employ both measured and empirical parameters, while empirical melt models use empirical parameters alone. Empirical melt models have been applied at regional scales in which transferability in space and time must be assumed, despite evidence that such transferability is not universally possible (Hock et al. 2007; Hock 2003). Two studies have specifically explored the transferability of melt model parameters between glaciers within the same mountain range. Carenzo et al. (2009) explored the transferability of an enhanced temperature-index melt model (Pellicciotti et al. 2005) within the Swiss–Italian Alps. This study found the model to be highly transferable between three glaciers in three summers, except when the model was calibrated separately for overcast and clear conditions. Wheler (2009) explored the transferability of the Hock (1999) enhanced temperature-index melt model within the Donjek Range of the St. Elias Mountains. He found a small reduction in model skill when transferring the model in space and high model transferability in time.
Energy-balance melt models are the most physically justified form of melt model (Hock 2005) and have been developed for several of the best-studied glaciers in the world (e.g., Anderson et al. 2010; Anslow et al. 2008; Arnold et al. 1996; Hock and Holmgren 2005; Klok and Oerlemans 2002). Distributed energy-balance melt models (DEBMs) either extrapolate or parameterize each of the components of the surface energy balance across a grid that covers the glacier surface and compute the energy available for melt in each grid cell. Despite the advantage of being able to model the effects of climate change on individual components of the energy balance (Hock et al. 2007), these models have not been used for regional studies because of their large data requirements and perceived complexity (Hock and Holmgren 2005; Pellicciotti et al. 2005). Here we take the first steps to applying a DEBM regionally by assessing the transferability of model parameters and meteorological forcing variables between two small glaciers in the Donjek Range of the St. Elias Mountains.
2. Study site
The St. Elias Mountains, located in northwestern North America, are characterized by extreme topographic gradients (Clarke and Holdsworth 2002) and host one of the largest glacierized regions outside of the Arctic or Antarctic (Arendt et al. 2008). During the latter decades of the twentieth century, glaciers in southeastern Alaska and the Coast Mountains contributed more to sea level rise than any other glacierized region (Kaser et al. 2006; Lemke et al. 2007). Within this region, the greatest single contribution came from the glaciers of the St. Elias Mountains (Berthier et al. 2010). The Donjek Range is located in the southwestern Yukon Territory of Canada—at the eastern edge of the St. Elias Mountains, just beyond the contiguous ice fields (Fig. 1a). The range is separated from the Gulf of Alaska by the highest peaks of the St. Elias Mountains and therefore, despite being less than 100 km from the ocean, experiences a continental climate (L’Heureux et al. 2004). This study is conducted on two unnamed mountain glaciers 10 km apart in the Donjek Range (Fig. 1b). One glacier has a predominantly southerly aspect and a surge-type dynamic regime (De Paoli and Flowers 2009) and is referred to as “South Glacier” (Fig. 1c). The other has a northwesterly aspect and is referred to as “North Glacier” (Fig. 1d). South Glacier is thought to have a polythermal structure (De Paoli and Flowers 2009) similar to that of Storglaciären in northern Sweden (Pettersson et al. 2004). The thermal regime of North Glacier is also presumed to be polythermal but has not been studied. The two glaciers have been studied in detail from 2007 to 2009 with a full complement of instruments in the 2008 and 2009 field seasons.
a. Field methods
Parallel sets of meteorological and mass-balance instrumentation were deployed on each glacier. An automatic weather station (AWS) is located in the central ablation zones of each glacier at ∼2300 m above mean sea level (AMSL). The AWSs are instrumented to measure typical meteorological variables as well as net radiation and incoming and reflected shortwave radiation (Table 1). An ultrasonic depth gauge (USDG), which measures surface lowering, is deployed several meters from each AWS. In addition to the continuous point record from the USDG, mass balance is measured at an array of 18 (17) ablation stakes deployed on South (North) Glacier (Figs. 1c,d). Stake heights are measured at least twice per year, at the beginning and end of the melt season. Snow pits are excavated near the AWSs and in the accumulation zones of the glaciers each May in order to evaluate the density and temperature structure of the snowpack. The albedometers were deployed on the AWSs only during the “summer season,” 5 May–11 September 2008 and 9 May–27 July 2009 on South Glacier and 8 May–14 September 2008 and 18 May–3 August 2009 on North Glacier. These dates define the maximum period for which the DEBM can be run and the length of time between them defines the summer season for our purposes.
2) Data processing
Meteorological data are gap filled using linear interpolation and averaged to hourly values. No gap in the three-year record is longer than 30 min for deployed and undamaged instruments. The net radiation record contains spurious five-minute readings lower than −300 W m−2. These readings are removed and replaced with linearly interpolated values. The wind monitor on the North Glacier AWS was destroyed in winter of 2008 and was not replaced until June 2008. The wind regime of North Glacier is dominated by katabatic flow during periods of positive air temperatures. To fill the gap in the wind record, wind speeds were simulated using the katabatic coefficient parameterization, developed by Klok and Oerlemans (2002), which relates wind speed to glacier temperature for positive glacier temperatures. The best-fit parameterization reproduces the existing wind record during positive temperatures with a correlation coefficient R of 0.56. Using this wind parameterization may be a significant source of error in the first 41 days of the North Glacier 2008 simulation. The AWS on North Glacier was tilted 6° to the south when visited in August of 2009. This tilt produces an error in the measurement of incoming direct shortwave radiation that was corrected using the relationship between the slope and aspect of a plane and the incident radiation on that plane [Oke 1987; Eq. (2)] and by assuming a linear drift from horizontal. This tilt would, under the most detrimental condition (summer solstice at noon), produce a 7% overestimation of direct shortwave radiation. At most times the error is much smaller.
The USDG record of surface lowering is converted into estimates of accumulation and ablation at the AWS locations. Daily averages are taken to remove noise in the records (Brazenec 2005). Positive differences in the distance to surface are taken as accumulation events and negative differences as ablation events. Fresh snow density is taken as 200 kg m−3 (Paterson 1994) (field measurements of fresh snow density range from 128 to 252 kg m−3). Water-equivalent ablation from an ablation event is calculated using the density of either the May snow pit density at the AWS, a linear combination of snow pit density and cumulative fresh snow density, or ice density (900 kg m−3) if the spring snowpack has been removed. The ablation stake measurements are used to estimate mass balance following the method of Østrem and Brugman (1991). Ablation stakes were surveyed at the beginning and end of the summer season on North Glacier and at weekly to monthly intervals during the summer season on South Glacier. During each survey the height of the stake, depth of snow, and snow density are recorded. An estimate of uncertainty is made for each height measurement using a 30 × 30 cm2 disk placed at the base of the stake. Changes in stake height, snow density, and snow depth (when a snow-to-ice transition occurs) are used to estimate the change in mass at each stake. Snowfall recorded at the AWS during the summer season is used to estimate summer surface accumulation and to infer ablation missed between surveys. Only the estimates of cumulative surface ablation for the full summer season at the stake locations are used as validation data for this study. No attempt was made to interpolate ablation from the stake locations to the rest of the glacier [see MacDougall (2010) for additional details].
3) Initial snow depth and snowfall events
The glacier-wide initial snow depth in snow water equivalent (SWE) is calculated by interpolating measured May snow depths from the ablation stake locations using linear regressions on slope and elevation for South Glacier and linear regression on elevation for North Glacier (Wheler 2009; Table 2). The relationship between SWE and slope and SWE and elevation is significant on both glaciers; however, on North Glacier the relationship with respect to slope and elevation does not yield improvements over the relationship with respect to elevation alone. Snow pit density is used to estimate water equivalent (w.e.) snow depth at the ablation stake locations. The area-averaged initial snow depth is the estimated winter accumulation. Snowfall events are recorded at the AWS by the USDG. These events are extrapolated to the rest of the glacier using an empirical precipitation lapse rate (Γp; Table 3), a temperature melt–freeze threshold (1°C), and the rainfall records from the AWS. The temperature melt–freeze threshold discriminates between solid and liquid precipitation and is treated as a constant for simplicity. The empirical precipitation lapse rate was derived from field measurements of summer snowfall events at different elevations at each study site (Wheler 2009).
4) Temperature lapse rates
Three (four) Onset HOBO microloggers are deployed on South (North) Glacier at various elevations (Figs. 1c,d). These microloggers measure temperature and relative humidity. Temperature records for each of the microloggers are averaged and a linear regression on elevation is used to estimate the temperature lapse rates (ΓT). The lapse rate for South Glacier 2008 and South Glacier 2009 is −6.0 K km−1 and for North Glacier 2008 is −5.3 K km−1. The temperature lapse rate for North Glacier 2009 could not be derived because of sensor malfunction and is therefore assumed to be the same as North Glacier in 2008.
5) Aerodynamic roughness length
In an attempt to limit the need for model tuning, we measured the aerodynamic roughness length of the glacier surface—the height above the ground at which wind speed reaches zero (Oke 1987). This parameter is needed for the estimation of turbulent heat fluxes and is frequently used to tune DEBMs (e.g., Hock and Holmgren 2005; Anslow et al. 2008). Roughness length was measured using the microtopographic approach pioneered by Munro (1989), wherein the aerodynamic roughness is related directly to the physical roughness of a surface profile. Instead of measuring the distance from a pole or taught string to the surface to create a profile (e.g., Brock et al. 2006; Munro 1989), the method of Rees (1998) was implemented. In this method, a topographic profile is created by inserting a dark backdrop into the surface and taking a digital photograph of the backdrop–snow boundary. This image is then processed into a surface profile. The backdrop insertion method works only for soft surfaces. For ice and hard-firn surfaces, a similar but new method to capture a surface profile is used: a dark string is laid across the surface such that it is in contact with the surface at all points. A meter stick is placed horizontally behind this string such that the string is equidistant from the meter stick as viewed from above to minimize parallax. A surface profile is extracted from a digital photograph of the string and meter stick (MacDougall 2010). The string method produces surface roughnesses ∼30% lower than those produced by the original method of Munro (1989) at the same locations. As aerodynamic roughness length relates to turbulent fluxes logarithmically, the difference in magnitude between the two methods creates only a small bias. We measure a mean snow roughness length of 0.4 mm and mean ice roughness length of 0.8 mm. These values are comparable in magnitude to those found for high-latitude glaciers (Brock et al. 2006).
b. Modeling methods
The DEBM created for this study is similar to those created by Hock and Holmgren (2005) and Anslow et al. (2008). Innovations in this DEBM are 1) a new simple subsurface model for calculating surface temperature and subsurface heat flux, and 2) the incorporation of a time-evolving parameterization for snow aerodynamic roughness length. The surface energy balance of ice or snow is written as
where Sin is incoming shortwave radiation, α is the albedo of the ice or snow surface, Lin is the incoming longwave radiation, and Lout is the outgoing longwave radiation. Here QH is the sensible heat flux, the energy exchanged between the glacier and the atmosphere; QL is the latent heat flux, the heat transferred to or from the glacier through sublimation, deposition, evaporation, or condensation; Qg is the heat transferred to and from the glacier subsurface when the ice or snow of the subsurface changes temperature; QR is the heat flux from rain, the sensible heat released when rain is cooled to the freezing point; and QM is the energy available to melt ice or snow. The sensible heat transferred to the glacier from the cooling of rainfall to the freezing point (QR) is disregarded, as rainfall totals are small in our field area.
1) Shortwave radiation
To extrapolate measured incoming shortwave radiation to each digital elevation map grid cell, Sin is broken into direct and diffuse components using the empirical relationship of Collares-Pereira and Rabl (1979), following Hock and Holmgren (2005). Direct shortwave radiation is only incident on grid cells that are unshaded by surrounded terrain. At each time step, topographic shading is calculated using the traverse of the sun through the celestial hemisphere and the local horizon at each grid cell. Direct shortwave radiation is corrected from the measured horizontal to the slope and aspect of each grid cell following Oke (1987):
where I is the direct shortwave radiation, Ie is the direct shortwave radiation effective on each grid, θ is the incident angle of the sun on the surface, and θh is the incident angle of the sun on a horizontal surface.
Diffuse radiation is assumed to originate from all parts of the sky equally. The two sources of diffuse shortwave radiation are that incident from the unobscured fraction of the sky and that reflected by surrounding terrain. Diffuse shortwave radiation is therefore calculated as
where Do is the diffuse radiation from an unobscured sky, Vf is the sky-view factor (the fraction of the sky that is unobscured by topography), and αter is the albedo of the surrounding terrain. The albedo of the surrounding terrain was estimated in the field by measuring the albedo of rock slabs that had fallen from nearby slopes. The effect of snow cover on surrounding slopes is ignored, as snow is rarely present on surrounding slopes during the summer season. The measured shortwave radiation at the AWS already includes the effects of surrounding terrain, therefore Do must be backed out of Eq. (3).
If the AWS is shaded by surrounding topography, all measured incoming shortwave radiation is assumed to be diffuse. The ratio of diffuse to incoming shortwave radiation is assumed to stay constant until the next time step at which the AWS is unshaded. This ratio is used to estimate direct shortwave radiation for the period of AWS shading (Hock and Holmgren 2005).
Many parameterizations have been proposed for modeling the evolution of snow albedo throughout a melt season (e.g., Arnold et al. 1996; Brock et al. 2000; Hock and Holmgren 2005; Oerlemans and Knap 1998). We use a slightly altered version of the albedo parameterization developed by Hock and Holmgren (2005) in the DEBM (Fig. 2). The parameterization is written as
where αt−1 is the albedo at the previous time step, αt is the albedo at the current time step, nd is the number of days since the last snowfall, Ps is the measured rate of snowfall, T is temperature, and a1 to a4 are constants that must be found through calibration (Hock and Holmgren 2005). The original parameterization proposed by Hock and Holmgren (2005) included an alteration of the albedo based on a cloudiness factor that has not been successful outside of Storglaciären (R. Hock 2009, personal communication). When this parameterization was applied to our study sites, snow albedo in certain grid cells falls to absurdly low values, prompting us to impose a minimum value of 0.66 for snow albedo (αs-lim) (Paterson 1994). Firn is treated in a similar fashion as snow. When a snow–firn transition occurs, an instantaneous change in albedo of 0.03 is applied to account for the accumulation of debris in firn (αstof), after which the firn is treated the same as snow, following Hock and Holmgren (2005). Firn is given a different lower albedo limit than snow of 0.56 (αf-lim) (Paterson 1994). A constant elevation firn line was assigned based on field observations (Ef). No attempt was made to account for firn–ice transitions, as the depth of the firn at our study sites is unknown. Ice is assumed to have a constant albedo (αi) (e.g., Hock and Holmgren 2005; Oerlemans and Knap 1998), taken as the mean value measured at the AWS location, with a constant elevation-dependent albedo decrease applied to account for greater debris cover at lower elevations (dαi/dZ). A lower limit to ice albedo is imposed and set equal to the lowest measured values of ice albedo at each glacier (αi-lim). Rates of ice darkening were estimated using photographs taken during ablation stake surveys (Corripio 2004).
3) Outgoing longwave radiation and subsurface heat flux
Outgoing longwave radiation is calculated from the temperature of the ice or snow surface (Ts) according to the Stefan–Boltzmann relationship
where ε is the emissivity of the surface and σ = 5.67 × 10−8 J s−1 m−2 K−4 is the Stefan–Boltzmann constant. Following convention in energy-balance modeling, the emissivity ε of ice and snow is taken as 1.0 (e.g., Arnold et al. 1996; Hock and Holmgren 2005). The temperature of the surface and the subsurface heat flux are calculated using a simple subsurface scheme, where the subsurface flux Qg is taken as the residual of the energy balance when the energy balance is negative. The subsurface flux is forced into a thin subsurface layer such that
where Δt is the model time step in seconds, ρs is the surface material density, cs is the specific heat capacity of ice, and h is the thickness of the subsurface layer. The depth of the subsurface was chosen as h = 0.10 m to maximize model stability while maintaining a shallow enough depth to produce physically plausible changes in surface temperature. For this scheme to work during periods of extended subfreezing temperatures, a second passive layer is added to the model. This second layer takes up cold content if the surface temperature drops below a defined threshold (Tsub = −30°C). When the residual of the energy balance is positive, the energy warms the subsurface layer but the surface is prevented from melting until all of the cold content in the passive layer is balanced by the subsurface heat flux, thus enforcing the conservation of energy. A priori knowledge of Ts is necessary to calculate Lout and vice versa; therefore, an iterative loop initialized with a surface temperature of 0°C is required to solve the energy balance.
4) Incoming longwave radiation
Incoming longwave radiation is computed as the residual of the radiative energy balance at the AWS location and is assumed to be constant over the entire glacier, following Hock and Holmgren (2005). The radiative balance is
where QN is the measured net radiation and Sref is reflected shortwave radiation. Incoming longwave radiation emitted from surrounding terrain Lterrain was calculated following Anslow et al. (2008):
where εter is the emissivity of surrounding terrain taken as 0.95 (Anslow et al. 2008), Ta is air temperature, and Lsky is the incoming longwave from an unobscured sky. Although this parameterization is an obvious oversimplification, it has been shown to work well (Greuell et al. 1997). Quantity Lsky must be backed out of Eq. (9) as the AWS is affected by surrounding topography:
where all values are for the AWS location.
5) Turbulent heat fluxes
Sensible and latent heat fluxes are calculated using the bulk aerodynamic approach as in other recent DEBM studies (e.g., Anderson et al. 2010; Anslow et al. 2008; Brock et al. 2000; Hock and Holmgren 2005). Despite glacier boundary layers violating some of the assumptions that are used to derive the bulk aerodynamic equations, the method has been shown to work well in glacierized environments (Denby and Greuell 2000). Sensible heat flux QH is calculated using this approach as
where ρa is the density of air, cp is the heat capacity of air, k = 0.4 is the von Kármán constant, uz is wind speed, Tc is the air temperature relative to the freezing point and adjusted for air pressure, z is the height at which the measurements are made above the ground, zo is the aerodynamic roughness length of the surface, zoT is the roughness length for temperature, L is the Obukhov length, and ΨM and ΨH are stability constants. Latent heat flux QL is similarly calculated as
where ez is the vapor pressure at height z, es is vapor pressure at the surface, Lυ is the latent heat of vaporization (as in Anslow et al. 2008), pc is the atmospheric pressure, and zoe is the roughness length for humidity. The equation for L depends on QH such that an iterative loop is necessary to calculate these energy fluxes (Hock and Holmgren 2005). The stability constants ΨM and ΨH are calculated for stable conditions using the nonlinear stability functions of Beljaars and Holtslag (1991), following Hock and Holmgren (2005) and Forrer and Rotach (1997). For the infrequent unstable conditions, the commonly used Businger–Dyer expressions are used (Paulson 1970).
6) Aerodynamic roughness length z0
The aerodynamic roughness length (z0) used in the bulk aerodynamic approach is often used as a tuning constant in DEBMs (e.g., Anslow et al. 2008; Hock and Holmgren 2005). For this study we instead attempt to parameterize this quantity. Brock et al. (2006) made the first attempt to create a time-evolving parameterization of the roughness length of snow, relating the natural logarithm of snow roughness length to the logarithm of the sum of positive degree days since the last snowfall event:
where Pdd is the base 10 logarithm of the sum of daily maximum temperatures since the last snowfall event and b1 to b4 are empirical constants. A snowfall threshold is defined above which the sum is reset to zero (Zthr = 0.01 m). We calculate the empirical constants using a least squares fit between the parameterization and snow roughness lengths measured during the 2009 field season. Insufficient data existed to produce independent fits for North Glacier and South Glacier, so the datasets were combined to provide one parameter set for the entire study area. The 2009 parameter values were used for the 2008 simulations. Brock et al. (2006) failed to find a time-evolving roughness length parameterization for ice. We therefore use the mean value of the measured roughness of ice for each glacier (zoi). Firn is treated simply as very old snow. The roughness lengths for temperature zoT and water vapor zoe were taken as two orders of magnitude smaller than the roughness length zo, following Hock and Holmgren (2005).
7) Model calibration
The model is calibrated using the record of albedo at the AWS to solve for the four empirical albedo parameters (a1to4) separately for each glacier and for both summers. The snow roughness length evolution parameters are calibrated to the roughness length measurements. The mass-balance data are not used to calibrate the model. For consistency with previous studies, an experiment was carried out wherein the model was tuned by calibrating the roughness length of snow and ice using the mass-balance data. The experiment resulted in a 1% reduction in root-mean-square error (RMSE) between the measured and simulated ablation, indicating that tuning the model with mass-balance data would not significantly improve model performance.
4. Model validation and control runs
Four control runs are carried out to assess model skill and create references from which to judge model performance in the transfer experiments. The model is run for both glaciers for the summers of 2008 and 2009 using the locally derived parameter values (Tables 3 and 4) and forced with locally measured meteorological variables. The output from these runs is compared to the measured ablation at the stake locations (Figs. 3a–d) and the daily USDG-estimated ablation (Fig. 4). Error in the simulations relative to the ablation stake measurements is quantified using two statistics: mean percentage error (MPE), a measure of model bias,
where Ms,i is the simulated ablation in grid cell i, Mm,i is the measured melt at the stake located in the grid cell i, Mm is the mean of the measured melt, and n is the number of ablation stakes surveyed. For assessing model performance, we found it useful to compute RMSE both as a relative error (RMSEp expressed in %) and as an absolute error (RMSEa expressed in m w.e.).
The comparison of modeled and measured ablation at the stake locations (Figs. 3a–d) shows that in three of the control experiments, the model underestimates ablation by 10%–11%; for North Glacier in 2008, the model overestimates ablation by 18%. The relative (absolute) RMSE ranges from 14% to 30% (0.11–0.29 m) with North Glacier 2008 having the highest relative RMSE and South Glacier 2008 having the highest absolute RMSE. The simulations for summer 2009 produce a smaller RMSE than those for summer 2008 but a similar magnitude bias. The comparison of modeled and measured ablation at the AWS location shows that the simulated ablation tracks the USDG record within 20 cm w.e. (Figs. 4e–h). Ablation for South Glacier 2008 is underestimated by the model at the beginning of the season and between days 210 and 240 but is close to the USDG record in midsummer. Modeled ablation for South Glacier 2009 closely tracks the record except near the snow–ice transition date and near the end of the season. Ablation for North Glacier 2008 is overestimated by the model for most of the summer after the USDG began functioning with periods of convergence just after the sensor begins working and in the days before the late summer snowfall. The North Glacier 2009 simulation closely tracks the USDG record except during a period in midsummer centered around the snow–ice transition date. Overall model performance is similar to that achieved using enhanced temperature-index models (Wheler 2009), a result that is consistent with other efforts at physically based melt modeling (Hock 2005).
Areally averaged energy-balance components (Table 5) show that both glaciers are dominated by net radiation, as expected in a continental climate, with the radiative heat fluxes contributing between 60% and 83% of the melt energy. Areally averaged surface mass-balance components (Fig. 5, first column) show that both glaciers had negative net mass balances for the 2008 balance year (−0.18 m w.e. for South Glacier and −0.07 m w.e. for North Glacier) and negative net mass balances at the time the monitoring program was discontinued in late July of 2009 (−0.23 m w.e. for South Glacier and −0.25 m w.e. for North Glacier). Figures 4a–d display the spatially distributed ablation calculated for both glaciers for both years. These maps of modeled ablation are used as references against which to compare the results of the transfer experiments.
Model sensitivity to parameter values and uncertainties in meteorological forcing were assessed by individually perturbing each of the 21 parameters (Tables 3 and 4) and each of the 10 meteorological forcing variables while holding all other quantities at their respective control values [following Anslow et al. (2008) with the exception that we have omitted the Monte Carlo analysis]. Parameter values were perturbed by ±10% from their control values. The results of these perturbations (Fig. 6) show that the model is most sensitive to albedo parameters, particularly the lower limit of snow albedo and the a2 rate parameter, and firn-line elevation. An important caveat when interpreting Fig. 6 is that the ±10% perturbation does not account for the plausible range of values that each parameter can adopt. For example, the roughness length of ice can vary over many orders of magnitude (Brock et al. 2006), while the plausible range of the lower limit of snow albedo is likely much smaller than the 0.59–0.73 range examined. The meteorological forcing variables were varied within their measurement error (Table 1) following Anslow et al. (2008). The model was most sensitive to uncertainties in measured wind speed (42% variation in relative RMSE) and estimated initial snow depth (38% variation in relative RMSE).
5. Transferability experiment design
Model transferability is the ability of a model calibrated for one time and location to produce realistic results at another time and/or location. We describe transferability here both in terms of model parameter values and meteorological forcing variables. The former is the ability of parameters calibrated or measured for one time or location to describe another, while the latter is the ability of meteorological variables measured at one site to force a model at another site. Ideally, a model will have both high parameter transferability and high forcing transferability, such that a model calibrated at one well-studied site can be used successfully at surrounding sites.
Parameter transferability is assessed in three experiments composed of four tests each (Fig. 7a) examining parameter transferability in time, in space, and in space and time together. In each of these experiments, the parameter values derived for one glacier and year are used in place of those locally derived for the other glacier and/or year. That is, the 11 parameter values displayed in Table 4 for a specific glacier and year were substituted with values from another column of Table 4. The meteorological forcing variables are locally measured for these experiments.
Transferability of the meteorological forcing is assessed in four tests (Fig. 7b), wherein the variables measured at one glacier are used to drive the model for the other glacier for the same time period. No meteorological variable transfer is conducted in time because such a transfer is illogical. The meteorological variables are transferred to the other glacier in the same way that the meteorological variables are extended from the AWS to the each grid cell on a single glacier in the control runs. This approach is no different than extrapolating AWS measurements across a large glacier. The AWS is not assumed to be in any location other than its own physical location. The parameters are held at their locally derived values. For the purpose of this experiment initial snow depth is grouped with meteorological variables, although it is an initial condition related to the integrated local wintertime precipitation. Initial snow depth is transferred between sites by imposing the elevation and slope regression parameters derived for one glacier on the other. In most practical applications, parameter values and meteorological variables would be transferred jointly in space. The primary value of the individual parameter and meteorological transfers in space is therefore to partition the sources of error in model transferability.
Four joint-transferability tests are conducted wherein both the parameters and meteorological variables associated with one glacier are transferred to the other for the same year (Fig. 7c). This experiment effectively treats one glacier as an extension of the other. A final experiment is conducted wherein all but several key parameter values and meteorological variables are transferred from one glacier to the other. This experiment tests whether model transferability can be improved with only a small number of locally derived quantities.
To distinguish the many simulations, each is given a unique code in which “P” indicates parameters, “M” indicates meteorological variables, “S” indicates South Glacier, “N” indicates North Glacier, and years are represented by their last two digits. Quantities that have been transferred are shown in bold. For example, the code S08-PN08-MN08 represents the simulation for South Glacier 2008 using parameter values and meteorological variables associated with North Glacier in 2008.
a. Parameter transfer in time
The comparison between simulated and measured ablation at the stake locations (Figs. 3e–h) shows that temporal parameter transfer tests produce about the same magnitude of error as the control runs for South Glacier and higher errors than the control runs for North Glacier. The deviations in spatially distributed cumulative surface ablation between the temporal parameter transfer tests and control runs (Figs. 8a–d) are small for both glaciers in both years. Relative to the control, North Glacier 2008 shows a region of overestimated ablation near the firn–ice transition line, while North Glacier 2009 shows a region of underestimated ablation in the same location. The absolute differences in surface ablation and mass balance calculated in these tests compared to the control run are relatively small (TP in Table 6), ranging from 1 to 4 cm w.e. The relative similarity in mass-balance components between these tests and the control runs is seen in Fig. 5 (cf. column C to TP).
b. Parameter transfer in space and time
Transferring parameters in space or in space and time together produces similar results with significantly higher errors than in the control runs (Figs. 3i–p). MPE, a measure of model bias, is as high as 55% and relative (absolute) RMSE as high as 126% (0.53 m) in these tests. The deviation in spatially distributed ablation between these tests and the control runs (Figs. 8e–l) is large in the ablation zones of the glaciers. On South Glacier, ablation is underestimated in the ablation zone with respect to the control runs; the reverse is seen for the ablation zone on North Glacier. Deviations tend also to occur near the firn–ice transition region of each glacier. Differences in mass-balance components relative to the control runs are larger than those for the temporal parameter transfer experiment; the magnitude of the differences varies by glacier and year (Fig. 5, compare columns SP, TP, and C). Absolute differences in areally averaged modeled ablation between the control runs and the spatial parameter transfer tests range from 1 to 13 cm w.e. and differences in mass balance from 0 to 10 cm w.e. (SP and SPT in Table 6).
Transferring parameters in space or in space and time produces large differences in the modeled distribution of ablation and, consequently, in the calculated mass balance. The contrasting average ice albedo—lower for South Glacier than North Glacier (Table 4)—may explain part of this low transferability.
c. Meteorological variable transfer
The comparison between ablation simulated in the meteorological variable transfer experiment and measured at the stake locations (Figs. 9e–h) shows that the South Glacier tests, interestingly, produce smaller errors than the control runs; the North Glacier tests produce much larger errors than the control runs.
The differences in spatially distributed ablation between these tests and the control runs (Figs. 10a–d) show that ablation is overestimated in the ablation zone of South Glacier and underestimated in the ablation zone of North Glacier with respect to the control run. The comparison of mass-balance components between these tests and the control runs (Fig. 5, compare columns C and M) shows large differences in all of the mass-balance components, even resulting in the opposite sign for the net mass balance of North Glacier. The absolute value of the differences in areally averaged ablation ranges from 8 to 14 cm w.e. in this experiment, not significantly larger than the upper bound for the spatial parameter transferability experiments. The absolute value of the differences in areally averaged mass balance ranges from 22 to 50 cm w.e., much greater than that for the parameter transfer experiments (M in Table 6). Transferring meteorological variables (including the water equivalent initial snow depth) creates large errors in our ablation calculations, suggesting that forcing DEBMs with meteorological variables and accumulation derived at a remote site in the same region can produce unreliable output.
d. Parameter and meteorological transfer
Figures 9i–l compare measured ablation at the stake locations to ablation simulated in the joint-transferability experiment, in which both parameter values and meteorological variables were swapped. In three of four cases, the errors are much larger than in their respective control runs; for South Glacier 2009 the errors are roughly similar. The error magnitudes for the North Glacier tests are similar to, but slightly less than, those for the meteorological variable transfer experiment.
The differences in spatially distributed ablation between these tests and the control runs (Figs. 10e–h) show patterns of overestimation and underestimation for both glaciers in 2008 and for South Glacier in 2009. The pattern for North Glacier in 2009 exhibits uniform underestimation in the ablation zone with respect to the control runs. The mass-balance components (Fig. 5, compare columns M and MP) are similar to those for the meteorological transfer experiment. The absolute value of differences in ablation between these tests and the control runs ranges from 8 to 18 cm w.e. and in mass balance ranges from 12 to 55 cm w.e. (MP in Table 6).
Transferring both parameter values and meteorological variables fails to produce a better result than transferring meteorological variables alone, cautioning against the widespread application of site-specific meteorological variables (including accumulation) for accurate melt modeling in our study area.
e. Selective parameter and meteorological transfer
Given that joint transfer of parameter values and meteorological variables arguably does not produce accurate estimates of mass balance, we investigate whether retaining several locally derived quantities improves model performance. Based on the variability of parameters between study sites and years (Table 4) and the results of our sensitivity tests, we identify three potentially important quantities and fix them at their locally derived values: ice albedo, initial snow depth, and summer snowfall rate.
The comparison between the ablation simulated in this experiment and that measured at the stake locations (Figs. 9m–p) shows errors larger than in the control runs, but smaller than in all of the other transferability experiments for North Glacier and all but the temporal parameter transferability experiment for South Glacier.
The differences in spatially distributed ablation between these tests and the control runs (Figs. 10i–l) show regions of overestimation and underestimation for all of the tests with respect to the control runs, but of lower magnitude than in the joint-transfer tests. The comparison of the mass-balance components (Fig. 5, compare column C to MP − αi, swe, Ps) shows that these tests produce results reasonably similar to the control runs. The absolute differences in ablation and mass balance between these tests and the control runs are identical and range from 0 to 9 cm w.e. (MP − αi, swe, Ps in Table 6).
This experiment demonstrates that ablation can be replicated to within 9 cm w.e. of the values estimated in the control runs at our sites using a DEBM with transferred parameter values and meteorological variables, provided local values of summer snowfall, initial snow depth, and ice albedo are known.
The areally averaged energy-balance components calculated by our model can be compared between the two glacier sites (Table 5). South Glacier receives more net shortwave radiation than North Glacier and turbulent fluxes are larger on North Glacier than South Glacier. This pattern is expected considering the glaciers’ respective aspects. Compared to other glaciers where DEBMs have been implemented, the Donjek Range study glaciers derive a relatively high proportion of melt energy from net radiation: net radiation constitutes 83% of the melt energy for South Glacier in 2008, 81% for South Glacier in 2009, 60% of the melt energy for North Glacier in 2008, and 76% for North Glacier 2009. These values are comparable to Haut Glacier d’Arolla in summer 1990 (82% net radiation) (Arnold et al. 1996), and South Cascade Glacier in summer 2004 (62% net radiation) and summer 2005 (68% net radiation). The relative contribution of net radiation to melt energy is higher for our study glaciers than that found for Storglaciären in the summers of 1993 (35% net radiation) and 1994 (57% net radiation) (Hock and Holmgren 2005) and than that found for Brewster Glacier between 2004 and 2008 (45% net radiation) (Anderson et al. 2010). We make these comparisons in areally averaged fluxes from different studies, bearing in mind that each study covers a different length of time and employs a different model design.
The high transferability of parameters in time suggests that this model could be run for many years at these sites without recalibration, opening the possibility of using regional climate model output to project future ablation. High transferability of both model parameters and meteorological variables is desirable for regional melt modeling, where it is convenient to treat all glaciers in a range as extensions of the study glacier. Our experiments have shown that extending the DEBM from one to many glaciers in the Donjek Range will probably lead to large errors in calculated ablation, unless ice albedo, initial snow depth, and summer snowfall are known for each glacier basin. Ice albedo can be derived from satellite remote sensing (Klok et al. 2003), but remotely deriving initial snow depth and summer snowfall presents a greater challenge. Snow depth can be derived from satellite remote sensing, even in complex terrain (Shi and Dozier 2000). However, spatial resolution, data availability at the right time of year, and uncertainties in the derived snow depth itself remain challenges in implementing remotely sensed snow depth in high-resolution modeling. Recent advances in orographic snow redistribution models (Dadić et al. 2010) suggest that it might be possible to use reanalysis precipitation data, coupled to a model of wind redistribution of snow, to produce realistic snow distributions for each basin of interest. The computational intensity of such calculations remains a significant obstacle. Existing physical and statistical orographic precipitation models have been developed and used in complex terrain (e.g., Kunz and Kottmeier 2006; Schuler et al. 2008); however, these models are unable to account for the transport of snow from valley walls to glaciers.
One of the ultimate aims of our research program in the Donjek Range is to create reliable regional-scale glacier melt models. In considering how best to extend the DEBM to other glaciers, one should consider both the variability of climate and the variability of glacier characteristics in the region. The primary differences between the glaciers are their locations on opposing sides of a topographic divide and that South Glacier exhibits a surge-type dynamic regime. North and South Glaciers have similar elevation-area distributions. Comparing the mean values of the measured meteorological variables shows that there is less difference between the sites in terms of barometric pressure, temperature, and incoming radiation than the measurement uncertainty for each of those variables. The only major meteorological difference between the sites is that North Glacier has a lesser initial snow depth than South Glacier. Summer snowfall and relative humidity were lower on North Glacier in 2008, but not in 2009 when snowfall totals were small (3 cm w.e. for South Glacier versus 5 cm w.e. for North Glacier in 2009).
Parameter values, particularly those related to albedo, vary much more substantially between the two sites than the meteorological variables. The difference in ice albedo between the two sites is curious, as ice albedo is expected to be similar for glaciers with commensurate debris inputs and dust deposition. This difference in ice albedo may be a product of South Glacier’s surge-type dynamic regime; anecdotally, field personnel have witnessed debris-rich water emerging at the surface from moulins and fractures on South Glacier. Much of this debris appears to remain on the ice, presumably lowering the ice albedo.
The most important difference between the two sites appears to be related to precipitation and orography. We hypothesize that the study glaciers are typical of their respective sides of the range, except perhaps in ice albedo. This hypothesis could be explored by deploying additional AWSs on the respective sides of the range and by extracting ice albedo values from satellite data for the other glaciers in the range.
Temperature-index models and simplified energy-balance models, by their nature, implicitly take account of all of the components of the energy balance with a small number of parameters (Hock 2005). Intuitively, this implies that DEBMs, having a strong physical basis, should exhibit higher transferability in space and time. Two studies that have examined the transferability of enhanced temperature-index models (Wheler 2009; Carenzo et al. 2009) found these models to be highly transferable in space and time. On the surface, the success of these simple models and the limitations on DEBM transferability found in this study might appear to undermine the case for using physically based models for regional melt modeling. It is important to recall, however, that neither Wheler (2009) nor Carenzo et al. (2009) attempted to transfer initial snow depth from one site to another. This is a serious caveat, considering how important knowledge of this quantity appears to be for DEBM model transferability. Also important to recall is that the results of simplified melt models have been shown to drift away from those of physically based models when ablation is projected decades into the future (Hock et al. 2007). Simplified models may therefore be a better approach to regional melt modeling over short periods of time or in locations where insufficient data exist to calibrate a DEBM.
DEBMs depend on parameterizations of complex processes to estimate the evolution of snow albedo and aerodynamic roughness length. They also must frequently make oversimplified assumptions about meteorological quantities—for example, that wind fields are uniform in space. Replacing key parameterizations with physically based submodels should be a priority. Advances have recently been made in the physical modeling of snow albedo (Gardner and Sharp 2010) and efforts should be made to integrate physical albedo models into DEBMs.
A distributed energy-balance melt model has been constructed for two glaciers in the Donjek Range of the St. Elias Mountains in the Canadian Yukon. The model was calibrated separately for the two glaciers for the summers of 2008 and 2009 using albedo and aerodynamic roughness length data. A series of experiments was performed to assess the spatial and temporal transferability of model parameters and meteorological variables. We find the model parameters to be transferable between years to within a bound of 4 cm w.e. (9%) in the total calculated surface ablation; model parameters are transferable in space to within a bound of 13 cm w.e. (28%). Transferring meteorological variables, including accumulation, produces large differences in the calculated ablation, as does transferring model parameters and meteorological variables jointly. The joint parameter and meteorological variable transferability can be greatly improved (to within a bound of 9 cm w.e., 15%), if locally measured ice albedo, initial snow depth, and summer snowfall are retained. Our results suggest that caution must be exercised when transferring melt models in space. If one does impose a melt model unaltered on surrounding glaciers, errors of ∼30% in total calculated surface ablation should be expected in study areas such as ours. Techniques must be developed to remotely determine certain glaciological quantities before regional glacier melt models can achieve high reliability in complex topography. Further study is required to determine whether similar limitations to regional model transferability exist for other glacierized regions.
We are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Foundation of Innovation (CFI), the Canada Research Chairs (CRC) program, the Northern Scientific Training Program (NSTP), and Simon Fraser University for funding. AHM is grateful for his support from NSERC PGS-M. 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 Stephen Soubliere (Trans North Helicopters) for logistical support, and to P. Belliveau, B. Wheler, G. Rosenkjær, N. Roux, A. Jarosch, L. Mingo, J. Logher, F. Anslow, C. Schoof, S. Heðinsdóttir, and A. Rushmere for field assistance. We thank an anonymous reviewer and F. Pellicciotti for their careful reviews.
Corresponding author address: Gwenn Flowers, Department of Earth Sciences, Simon Fraser University, 8888 University Drive, Burnaby BC V5A 1S6, Canada. Email: firstname.lastname@example.org