Translator Disclaimer
1 May 2009 Modeling Ablation on Place Glacier, British Columbia, from Glacier and Off-glacier Data Sets
Author Affiliations +

The results of two simulations of hourly ablation, from late July to September 2002, at a site on the Place Glacier are described. First, ablation is modeled from a data set collected at the glacier site; second, from a data set collected off-glacier at a site below the glacier tongue. The glacier data set simulations, based as they are on global and reflected short-wave radiation measurements, a net long-wave radiation model and the bulk turbulent heat transfer approach, provides a reasonably good simulation of cumulative ablation, amounting to almost 1.2 m during the experimental period, mostly due to melting out of the preceding winter snowpack. Average melt component flux densities due to net short-wave radiation, net long-wave radiation, sensible heat and latent heat due to moisture exchange are 93.1, −22.6, 14.4 and 3.91 W m−2, respectively, during this time. The glacier site data are also used to fit a snow albedo model to albedo measurements for old and new snow cover. The albedo model is then used in the second simulation, which is based on global radiation measurements, a similar net long-wave radiation model, and a heat transfer approach in which turbulent mixing due to katabatic and geostrophic flow is parameterized from the off-glacier temperature data. The second simulation scheme performs best if the katabatic component of the parameterization scheme is suppressed because the wind regime at the glacier site appears to be intermittent, but the scheme itself is based upon the idea of a continuously flowing glacier wind.


A distributed model of glacier surface energy exchange provides the key to realistic simulation of snowline retreat during the summer melt period. This in turn is crucial to the effective routing of surface meltwater input to the glacier hydrologic system because the rules of storage and delay are different for snow cover than they are for ice (Hannah and Gurnell, 2001). At the heart of the matter lies the ability to effectively model ablation anywhere on the glacier surface, such that the effects upon energy transfer of transition from snow to ice cover are properly taken into account. This may be done in situ by using automatic weather station (AWS) data obtained on the glacier surface. In distributed modeling this is generally done with data collected at an AWS located off the glacier, where it is easier to find a stable site for maintaining a year-round measurement program.

At Place Glacier, the location of this study, instruments located at a glacier AWS site are used to measure ablation in situ and to collect data from which to simulate the ablation there, using micrometeorological theory to model turbulent heat transfer contributions to surface melt (e.g., Munro, 1990). Also, measurements made at an off-glacier AWS site are used to estimate summer snow accumulation and to provide data for ablation simulation at the same glacier AWS site, using a snow albedo model and an air temperature parameterization of turbulent transfer (e.g., Klok and Oerlemans, 2002). The data used for the simulations were collected from 20 July to 13 September 2002, an ablation period during which the melting out of the old winter snowpack was followed by fresh snow on ice events.

The goal of the study is to see how effective the off-glacier simulations are in comparison to the glacier based simulations. This is done by comparing cumulative ablation from each simulation to total ablation measured at the glacier site during the study period. To account for possible deviation between the two simulations, the off-glacier simulation comparison is done in steps, such that an element of the glacier-based simulation is replaced by an element of the off-glacier simulation, until all elements of the simulation depend on off-glacier data.

Study Site and Data


The Place Glacier (50°26′N, 122°36′W), located in the Coast Range of British Columbia, covers a small area of approximately 4 km2. The glacier originates in an accumulation area that descends from ∼2610 m a.s.l. to an ablation zone, which terminates at ∼1850 m a.s.l. (Fig. 1). In contrast to many other glaciers that have been studied, the area of the accumulation zone is relatively small, taking up little more than a third of the glacier area. In fact, cold air drainage from the accumulation zone of the glacier would be expected to diverge in two directions: mostly toward the northwest, along the main tongue of the glacier, and southeastward, over the smaller tongue of the Joffre Glacier. As would be expected of a Coast Range glacier, winter accumulation can generate a deep snowpack (Østrem and Brugman, 1991), thus delaying the exposure of glacier ice well into the summer ablation period.

Figure 1

Map of the Place Glacier and adjacent ice areas (white line boundary), showing 1 km grid of UTM Zone 10 (lower left corner 526,000 m E., 5,584,000 m N., NAD83), likely directions for cold air drainage (broken line arrows), AWS sites (+ elevations) and 1991 micrometeorological site.



The off-glacier AWS was installed in July 2001 at 1840 m a.s.l., approximately 0.5 km northwest of the glacier terminus (Fig. 1). Data collected at this station throughout the year include hourly averages of air temperature, Ta, relative humidity fraction, rh, global radiation, K↓, wind speed, u, wind direction, udir, and total precipitation (P) for each hour (Table 1). The following year, the glacier AWS was installed in an open area (sky view factor 0.98) near the upper end of the ablation zone, at approximately 2044 m a.s.l. (Fig. 1). Hourly data obtained here include surface change, δh, due to ice and snow ablation (or accumulation) and averages of Ta, rh, u, K↓, and K↑, the reflected short-wave radiation (Table 1).

Table 1

Instruments deployed at AWS sites.


The glacier AWS functioned until the spring of 2003, when it was crushed by snow loading. In 2004, its air temperature and humidity sensor were relocated to a ridge above the base AWS, at 1901 m a.s.l. (Fig. 1). Also noted in Figure 1 is the site of a micrometeorological experiment conducted in 1991 from which wind speed and direction data were reviewed for the discussion section of this paper.

It is a straightforward matter to obtain vapor pressure, by applying the relative humidity fraction to the saturation vapor pressure at air temperature, using any one of a number of empirical relationships for that purpose (e.g., Oerlemans, 2000). However, inspection of the relative humidity records from the two stations showed that while data from the glacier station would maximize at ∼100%, those from the base station would achieve values no greater than ∼92%, thus casting suspicion on the results. Therefore, the sensors were brought together at the base AWS, for calibration against ventilated wet-bulb psychrometer vapor pressure values. In the case of the off-glacier AWS, it was found that a 1.23 multiplier and a −125 Pa offset were required to correct vapor pressures obtained from temperature and relative humidity data, while the equivalent corrections for the glacier AWS were 1.1 and −50 Pa. The corrections are likely to be specific to the sensors used in this study, so application elsewhere is not advised unless confirmed by other calibrations.

The water equivalent of snowfall at the glacier AWS site was estimated from the off-glacier AWS precipitation gauge record, setting the temperature threshold for snowfall at 1.5 °C (Klok and Oerlemans, 2002). Following Cheng (2001), noise in the gauge record was filtered by first taking differences between two moving 24-hour sums of the data, each one hour apart, to minimize diurnal variation. Then, expressing the results as hourly values, three passes of a moving 15-hour median filter were used to remove the remaining negative values. Gauge catch error was addressed by fitting second order polynomials to U.S. Army Core of Engineers corrections for the Alter shield (Anonymous, 1956), such that for wind speeds scaled to gauge height, of 4 and 8 m s−1, respective multipliers of ∼1.6 and 2 applied to rain, ∼2.4 and 3.5 to snow. Using air temperature data from each AWS to separate snowfall from a total precipitation of ∼131 mm during the period of this experiment, total snowfall was estimated to be 4.83 mm at the off-glacier site, 30.7 mm at the glacier site.

A Hamming filter (Hamming, 1989) was applied to the acoustic sounder record of the glacier AWS, choosing after some trial and error a 5-hour moving window for the task. Given the uncertainties associated with precipitation gauge corrections, it is interesting to take Σ(+δh) from the filtered acoustic sounder record during the summer snowfall period and to treat it as a snowfall thickness record. This period spanned one week, yielding Σ(+δh)  =  222 mm, which requires a multiplier of 0.138 to convert it to the 30.7 mm w.e. estimated from the gauge. This is a value that could reasonably apply to the specific density of fresh snow (Goodison et al., 1981), thus bolstering confidence in the precipitation estimates.


Allowing for summer accumulation and assuming the effect of settling within the snowpack to be negligible, the readily observed effect of ablation anywhere on the glacier surface is δh over time period, δt:

in which Ps is the snowfall rate, QM the melt energy input rate, Lf the latent heat of fusion, and ρfs, ρs, and ρi are, respectively, the densities of fresh snow, old snow, and glacier ice. QM is computed as the sum of the heat fluxes due to net short-wave radiation, K*, net long-wave radiation, L*, sensible heat transfer, QH, and latent heat transfer due to water vapor, QE:
which involves a modeling approach. The details of the approach differ according to whether glacier or off-glacier AWS data are being used to estimate melt energy inputs at the glacier site.


Glacier AWS data allow a combination of measurements and modeling to be used in estimating the net radiation terms:

in which the short-wave radiation measurements are as defined for Table 1, L↓ is incoming long-wave radiation modeled according to Equation 10 below, and σ is the Stefan-Boltzmann constant. Ts, the surface temperature, is assumed to be at 273.15 K.

The QH and QE transfers incorporate measurements of air and surface vapor pressures, ea and es, as well as air and surface temperatures, Ta and Ts:

wherein ρ is the air density, cp the specific heat of air at constant pressure, ϵ the gram molecular weight ratio of water vapor to air, p the atmospheric pressure, and Lv the latent heat of vaporization. Equation 4 is the bulk transfer approach, applicable to a glacier AWS, where wind speed, uz, scaled to the temperature-humidity measurement height z (e.g., Oerlemans and Klok, 2002), and the bulk transfer coefficients, CH,E, control heat conductance through the glacier boundary-layer.

The bulk transfer coefficients depend upon the roughness scaling lengths for wind speed, zo, temperature, zH, and humidity, zE:

in which k is von Kármán's constant and ψ is a stability correction, the value of which is iteratively obtained by way of Monin-Obukhov stability scaling (Munro, 2004). Taking a fixed value of zo  =  1.21 mm from roughness element description (e.g., Munro, 1989), dynamic zH,E values follow from the scheme proposed and recently reviewed by Andreas (2002).


Using the off-glacier AWS data, the net radiation components are calculated as

in which K′↓ is modeled global radiation and αS is the surface albedo. In addition to a model for K′↓, a suitable model is required for αS in order to allow response to changing snow cover conditions. The L↓ model remains the same, except that off-glacier rather than glacier AWS data are used in Equation 10.

Turbulent transfers using off-glacier AWS data are modeled from the parameterization scheme developed by Oerlemans and Grisogono (2002):

in which two heat conductance values, Kb and Kkat act in parallel, thus replacing the bulk transfer coefficients and uz in Equation 4. The katabatic boundary-layer conductance, Kkat, is parameterized from off-glacier AWS temperature data, using the normal lapse rate to adjust it to the altitude of the glacier AWS. A background geostrophic conductance, Kb, is obtained by adjustment, either to fit cumulative δh from Equation 1 to δh measurements (e.g., Klok and Oerlemans, 2002), or to cause Equation 7 to converge upon Equation 4 (e.g., Munro, 2004).

The katabatic boundary-layer conductance follows from Oerlemans and Grisogono (2002), where the strength of turbulent transfer is parameterized according to the temperature difference between the air mass and the glacier surface for Ta > 0 °C, using the altitude adjusted off-glacier AWS data:

where κ  =  0.0004 is an empirically derived constant, g is acceleration due to gravity, and Pr  =  5 is the turbulent Prandtl number (Oerlemans and Grisogono, 2002). Following Munro (2004) the potential temperature lapse rate, γ, varies sinusoidally between 0.0015 and 0.0085 °C m−1over one day, the maximum occurring at mid-afternoon.


Following Munro and Young (1982), K′↓ is comprised of a direct component, S, which depends upon hourly cloud cover fraction, n, and a diffuse component, D, which in addition to n depends upon a surface average reflectivity, αS, cloud absorptivity, an  =  0.18, cloud base albedo, αb  =  0.5, and an hourly cloud top albedo, αn, the daily mean of which is close to αb:

in which clear sky direct, So, and diffuse, Do, radiation may be estimated from a suitable model, in this case that which is described in Munro and Young (1982). Lacking observer estimates of hourly n, and noting that n is non-linear in Equation 9c, we iteratively extract trial values for the cloud cover fraction, n*, from the glacier and off-glacier K↓ records between 0800 and 1600 each day, such that a suitable choice of n* at each site causes K′↓ to converge on the measurements.

This is in contrast to the approach taken by Greuell et al. (1997) who, using a combination of their own cloud observations and data from Sauberer (1955), related cloud cover to atmospheric transmissivity (Fig. 2). A similar type of relationship, using n* as the basis for cloud cover, indicated that although n* estimates of n appeared to be reasonable in comparison to Greuell et al. near the extremes of the cloud cover range, they were comparatively low toward the middle. Divergence between their work and ours is consistent with Suckling and Hay (1977), who noted the tendency for observers to overestimate cloud cover fraction and for sunshine records to underestimate n.

Figure 2

Pyrometrically derived transmissivity and cloud cover scatter plot using values obtained from base (open symbols) and glacier (closed symbols) AWS data, with line plot showing the relationship of Greuell et al. (1997, their equation 4).


Therefore, to the extent that the pyranometer behaves like a sunshine recorder, as suggested by Equation 9b, an iteratively derived cloud cover fraction would be too low between the limits of its range. Munro and Young (1982) used a weighted combination of observer and sunshine record cloud estimates to develop Equation 9, giving the greater weight to observed cloud cover. The effect of their approach is represented by setting n  =  n*0.75 for long-wave radiation modeling. The result is overlapping sets of glacier and off-glacier AWS cloud fraction values that are well described by Greuell et al. (1997) throughout the range of values plotted (Fig. 2), so n  =  n*0.75 is used here.

Long-wave radiation depends upon hourly n from 0800 to 1600, but average n for this period is used to model L↓ from 1600 to 0800 of the following day. The L↓ model itself follows from Klok and Oerlemans (2002):

where the terms within square parentheses constitute effective atmospheric emissivity (Greuell et al., 1997; Konzelmann et al., 1994). The clear sky emissivity, ϵa, may be estimated from a number of empirical relationships (e.g., Oke, 1987), though we follow Klok and Oerlemans (2002). Also, as is usually done for mountainous terrain (e.g., Dozier and Frew, 1990; Klok and Oerlemans, 2002; Munro and Young, 1982), terrain influence corrections are made to both K↓ and L↓.


Surface albedo was modeled with reference to albedo measurements made at the glacier AWS site: αS  =  K↓/K↑. They yield the hourly pattern plotted in Figure 3, wherein values toward sunrise and sunset rise above the mid-day minimum. Daily mean αS were obtained by defining albedo according to the daily sums of the radiation terms. They occupy positions (not shown in Fig. 3) within the daytime range, closer to the minimum. It is toward the daily mean that modeling is directed.

Figure 3

Hourly αS measurements (broken line plots) and modeled αS using data from glacier (horizontal bars) and off-glacier (open circles) AWS sites. Also shown are water equivalent of old and new snow cover (heavy line plot), and model αS for days 235 to 244 using d*  =  32 mm (light line plot). Inset is the RMSE map, with ‘+’ at the minimum RMSE location.


Several albedo models can be considered, each of which has its strengths and weaknesses (Brock et al., 2000). We have chosen the αS model of Oerlemans and Knap (1998):

in which αs(i) is snow albedo on the ith day after a snowfall event, αfirn  =  0.53 the firn albedo limit, αfs  =  0.9 the fresh snow albedo, and t*  =  21.9 days, taken from Oerlemans and Knap, is a time scale for snowpack aging. Then
in which the specific mass of the snow cover, MS(i), follows from the cumulative effect of the processes stated in Equation 1 on snowpack thickness:
where do is the initial thickness. Then, following Klok and Oerlemans (2002), MS(i)  =  ds(i) ρs/fs and M*  =  d*ρs/fs is an optical depth scale for which d*  =  32 mm and density is expressed as a specific gravity. Eventually, snow cover ablation results in MS  =  0, resulting in αS  =  αice for which the value of 0.34 used by these authors is adopted here.

As noted in Brock et al. (2000) there are disadvantages to the use of this model, such as higher than measured albedo values when old snow cover thins out to expose the ice (e.g., Fig. 3 of Klok and Oerlemans, 2002; Fig. 3 of this paper). Also, errors in αS will be propagated in MS(i), which is a problem for iterative modeling because this affects the calculation of MS(i) through QM in Equation 1. Therefore Brock et al. (2000) recommended an air temperature–based parameterization of albedo that has since been applied to distributed modeling of Storglaciären (Hock and Holmgren, 2005). Nevertheless, we use Equation 11 here because the modeling proceeds sequentially, thus minimizing error propagation through MS(i). Furthermore, it appeals to curiosity about whether higher than measured albedo values during snowpack thinning can be corrected by a suitable choice of optical depth, as discussed below.

Values for ρs and ρfs were not measured in the field, so optimal values were estimated according to a two-way minimization of the root-mean-square error (RMSE) between modeled and observed daily mean αS (Fig. 3 inset), yielding ρs  =  0.568, ρfs  =  0.161. In so doing, Equation 11 was initialized with do  =  1115 mm, the snowpack thickness first measured in the field, and by setting i  =  31 to allow for aging of snow cover prior to when measurements began. The one modification to the model was to allow at least a fivefold increase in d* for old snow, the better to capture the pattern of albedo decline before day 245, thus avoiding the higher than measured αS that would otherwise be the result during the final 10 days of old snow cover (Fig. 3, days 235–244). Despite noticeable disagreement on some days the pattern of change throughout the experimental period is well captured by modeling with glacier AWS data, yielding a RMSE  =  0.034 to compare with the RMSE  =  0.050 that would be the outcome for an old snow d*  =  32 mm. Also good is the αS pattern modeled from non-glacier AWS data.

The snow densities seem to be reasonable relative to other findings. For example, mid-April old snow densities of ∼0.3–0.4 at Place Glacier are known to increase to 0.5 or more by mid-June due to settling and meltwater percolation (Østrem and Brugman, 1991), so a value of ∼0.57 by midsummer is plausible. The fresh snowfall density is close to the value of 0.138 suggested from the precipitation gauge and sonic sounder comparison stated above. In terms of work on other glaciers, ρfs  =  0.161 is not much different from the value of 0.18 used by Klok and Oerlemans (2002) to convert snow depths on the Morteratschgletscher to water equivalents.

More problematic is the better than fivefold increase in d* for old snow, in this case to 170 mm. Tentatively, we suggest that this compensates for the effect of water building up in the base of the snow cover as melting proceeds, the effect being to increase the transmissivity of snow layers adjacent to the ice, thus enhancing the effect of αice as the snow cover thins. That said, there is no reason to suppose that d*  =  170 mm is generally applicable, thus inviting further investigation.

Results and Discussion


Incoming radiation comparisons for the two sites are shown in Figure 4. Despite the fact that the stations are less than 3 km apart there is noticeable scatter in the comparison of global radiation measurements (Fig. 4a), most likely as a result of spatial variation in partial cloud cover. At low solar elevations there could be scatter due to differing times of topographic shading as well as occasional scatter from snowfall accumulation on the glacier AWS pyranometer dome. In bulk statistical terms, however, both the comparison of mean values and the slope of the regression are consistent with the expectation that global radiation at the lower, less exposed off-glacier AWS site should be somewhat less than that measured at the glacier site. The effect of constraining the regression to pass through the origin has only a minor effect on the slope, raising it to 0.94 while maintaining the value of r2.

Figure 4

Hourly off-glacier and glacier AWS incoming radiation values for (a) global radiation measurements and (b) long-wave radiation estimates. The intersection of the broken lines in (b) denotes long-wave radiation at 315 W m−2.


Cloud cover cannot be used to interpret the differences between means in Figure 4b; in fact, they differ by only ∼3 W m−2. The implications of the regression parameters are only important well outside the data range plotted here; in fact, a constrained regression results in a slope that is very close to one. The scatter is substantial over this range, as one would expect from using two different sets of n estimates obtained from the pyranometer record of each AWS, as well as the temperature records specific to each station. The interesting point about Figure 4b, however, is that many L↓ values exceed 315 W m−2, which is the approximate value of blackbody radiation from melting ice. This indicates that L* can contribute to melting at times when cloud cover and air temperature are sufficiently high. At this well exposed site, where a view factor ∼0.98 applies, the effect of topography is to increase L↓ by ∼1 W m−2.

Comparisons of hourly temperature and humidity at the sites (Fig. 5) show that there is close agreement between the two sets of humidity data. In fact the contrast between means over the study period, 709 Pa off-glacier compared to 695 Pa on the glacier, is consistent with the expected atmospheric pressure change as elevation rises from 1840 to 2044 m. Contrasts between the temperature data sets are manifest in the diurnal temperature range, where the smaller range for the glacier AWS is consistent with its location over a melting surface. Both air temperature and vapor pressure tend to be above their respective ice point values of 0 °C and 6.11 hPa, except during summer snowfall events, though drops in vapor pressure below ice point value do occur briefly elsewhere in the record.

Figure 5

Hourly temperature and humidity series for the AWS sites, with glacier AWS values shown as dots, off-glacier AWS values as light line plots and precipitation as bars. Inset shows estimated and measured glacier temperature comparisons.


Peyto Glacier estimates of glacier air temperature from off-glacier AWS temperatures, adjusted according to the normal lapse rate, yielded values that were approximately 20% higher above ice point than those measured at a glacier AWS (Munro, 2004). The same procedure for the Place Glacier yields glacier Ta estimates that are only slightly more, even potentially less than measured according to the slope of the regression between them (Fig. 5, inset). This could mean that the glacier cooling effect is being felt at the off-glacier site, a thought supported by subsequent comparison of air temperatures measured at 1901 and 1840 m. The comparison indicates that if Ta at 1901 m, adjusted according to the normal lapse rate were used to estimate Ta at 1840 m, the result would be ∼0.7 °C warmer than Ta recorded at the off-glacier AWS. Therefore, it may be that the off-glacier AWS record is influenced by the glacier wind, thus making it less representative of the regional air mass than initially thought.

Also worthy to note are the contrasting summer wind speed regimes at the two AWS sites and the predominance of glacier wind directions in the off-glacier AWS record (Fig. 6). At first, the summer wind speed pattern of the glacier AWS was thought to reflect sensor malfunction because much of the record indicated wind speeds less than the anemometer stall speed. But AWS comparisons during the winter, when air temperatures are below 0 °C and regional winds dominate, show vigorous wind speeds for both data sets. This suggests that although air temperatures above 0 °C needed for glacier wind development prevail during the summer at Place Glacier, they do not necessarily lead to persistent katabatic flow. Instead, the summer pattern appears to reflect an intermittent glacier wind regime, one that is not easily amenable to current modeling approaches. Furthermore, it offers another interpretation of the small contrast between off-glacier and glacier mean air temperatures: that the cooling effect of this glacier could be rather weak.

Figure 6

Summer and winter hourly wind speeds for off-glacier (lines) and glacier (dots lines) AWS sites. Also plotted (broken lines) are the off-glacier wind directions, the 120° to 180° range encompassing glacier wind direction.


Bearing these considerations in mind, off-glacier AWS data were prepared for modeling ablation at the glacier site, first by adjusting Ta to 2044 m, using a normal lapse rate of −0.0065 °C m−1, then by adjusting vapor pressure according to the expected ratio (0.98) of air pressure at the two sites. No elevation adjustment was made to off-glacier global radiation data, its sole purpose being to generate cloud estimates for radiation modeling. Also, for modeling purposes, the 1.5 °C temperature threshold for glacier snowfall was determined by the adjusted off-glacier Ta, the effect of which was to raise the total snowfall estimate at the glacier site to 31.7 mm.


The surface energy balance components of QM (Fig. 7) show that the turbulent transfer contributions are small relative to the net radiation terms on most days. Values plotted in Figure 7, averaged over the whole study period, are 93.1 W m−2 for K*, −22.6 W m−2 for L*, 14.4 W m−2 for QH, and 3.91 W m−2 for QE. The value for QH is no more than half what has been reported elsewhere from bulk transfer estimates over melting ice and snow (e.g., Munro, 1990; Oerlemans, 2000; Pohl et al., 2006). Together, the turbulent transfers amount to less than the estimated long-wave radiation loss from the site. Also, they depend on wind speed data, so hourly QH and QE show the same discontinuous behavior noted in the wind speed record (Fig. 5). The expected consequence of this is extreme stability, such that the effect of applying the stability correction scheme described in Munro (2004) is minimal. The impression gained from these results is one of summer melt that is dominated by radiative transfer. Nevertheless, the turbulent transfers are included in modeling cumulative δh from Equation 1 because they still account for approximately 15% of the ablation.

Figure 7

Surface energy exchange components computed from hourly glacier AWS data.


The ability to capture the cumulative sonic sensor δh with Equation 1, using glacier AWS data to model QM, appears to be good near the beginning and end of the old snow decay period, less so in the middle of this period, where failure to account for density variation within the snowpack may explain the result (Fig. 8). Substantial disagreement is seen after the old snow is gone, when the step change expected from the first fresh snow events of days 245 and 246 does not stand out in the sensor record, resulting in a modeled offset to cumulative δh after day 246 of approximately 0.1 m. Nevertheless, modeled change after that time is broadly similar to measured change, even to the extent that the change expected from snowfall on days 251 and 252 appears to be the same for both. A comparison between simulation and measurements yields a RMSE  =  87 mm, 57% of which is systematic, according to the statistical measures of Willmott et al. (1985). Better agreement can be forced by injecting a 0.1 m step into the measurement record, beginning with day 246, in which case the RMSE  =  38 mm, 31% of which is systematic (Table 1).

Figure 8

Measured (dots) and modeled cumulative δh, using data from glacier AWS (heavy line) and off-glacier AWS (light line), with snowfall (bar graph). Upward displacement of measurements, δh ∼0.1 m, compensates for initial lack of sonic sounder response to snowfall.


The question remains as to why the sonic sensor record should initially be so unresponsive to fresh snow accumulation. There was no observer on site to witness events. But one may speculate from the low albedo recording on day 244 that the events entailed melting out of the old snowpack, followed by fresh snowfall that accounts for albedo rise through days 245 and 246 (Fig. 3). In such circumstances there could have been water on the surface which initially absorbed or melted snow on contact, thus compromising snow layer thickening. After day 246, when absorption and freezing had firmed the surface, snow accumulation was fully effective in thickening the snow cover.

Some support for this speculation can be gained by looking at daily Ps/Σ(+δh) values which are 0.395 and 0.302 for the first two days, when Ps was estimated to total 12.1 mm w.e., but 0.147, 0.020, and 0.117 on the last three days, when the total Ps w.e. estimate was 17.6 mm. The smallest Ps/Σ(+δh) values, 0.026 and 0.018, apply to days 247 and 248, when less than 1 mm w.e. snowfall was estimated for each day. The fact that Ps/Σ(+δh) for the first two days stands well above the rest is consistent with the idea of initial loss from accumulation due to melt and absorption on contact, a process which is not included in the model.


Off-glacier data were introduced into the ablation model in a series of steps, each step constituting a replacement of a glacier-based model element with one from the off-glacier approach. Because the sonic sensor appears to have missed the early new snow events, δh  =  0.1 m at the start of day 246 was applied to the measurement record, as stated above. Thus the experiment is to see if off-glacier data and procedures can do as well as glacier data and procedures. This is done with distinction between the systematic (RMSES) and unsystematic (RMSEU) parts of the RMSE (Table 2), as outlined in Willmott et al. (1985), stating the statistics in mm. Included are the slope and intercept of the best fit line to modeled and observed Σ(δh).

Table 2

Glacier and off-glacier AWS ablation model outcomes.


The modeling steps were as follows, the first three involving straightforward change of data or procedure, the last two minimization of RMSE error between modeled and observed Σ(δh):

  1. replacement of glacier Ta with lapse rate adjusted off-glacier Ta to distinguish the snowfall in the precipitation record;

  2. replacement of glacier incoming radiation terms with off-glacier modeled values, a key change being the use of off-glacier n and Ta to model K↓ and L↓;

  3. replacement of K↑ measurements with the albedo model, using off-glacier modeled K↓ and L↓, but retaining glacier AWS QH and QE to estimate MS in Equation 11;

  4. use of Equation 7 rather than Equation 4 to model QH and QE, thus replacing the bulk transfer procedure with the full katabatic parameterization scheme; and

  5. use of Kb alone rather than Kkat + Kb to model turbulent transfer, a step that follows from the outcome of the preceding step.

The only effect of changing to off-glacier Ta is to increase the glacier PS estimate by 1 mm, the effect upon Σ(δh) being too small to generate more than a 0.002 m change from the glacier AWS result. Fortuitously, the RMSE measures are noticeably smaller (Table 2, step 1). On average, the replacement of the incoming radiation terms lowers K↓ by 1 W m−2, raises L↓ by 3 W m−2, and the RMSE, now at its highest value of 41 mm, has equal parts of RMSES and RMSEU (Table 2, step 2). Despite the 1.2 W m−2 decrease in K* and 2.7 W m−2 increase in L* that follows from this, the change in Σ(δh) from the glacier AWS result is still well within 0.01 m.

Incorporation of the albedo model raises K* to 95.8 W m−2 and Σ(δh) to −1.228 m, the most extreme of the ablation simulations, but with a reduced RMSE (Table 2, step 3). A comparison of the off-glacier modeled albedo pattern to that modeled from glacier AWS data shows that the rise in average K* does not always signify reduction in off-glacier modeled albedo (Fig. 3). Off-glacier albedo outcomes are lower than those from the glacier AWS near the end of the old snow cover period and on some days of new snow cover, but they are higher on other days of new snow, as would accord with the increased snowfall outcome of step 1. Nevertheless, most of the change from Figure 3 is toward lower albedo, so it appears that the effect on MS of the rise in L↓ outweighed the relatively small reduction in K↓.

Following Klok and Oerlemans (2002), the QH and QE parameterization entailed adjusting Kb in Equation 7 to obtain agreement between measured and modeled ablation, using the criterion of minimal RMSE (Table 2, step 4). The parameterized fluxes are clearly smaller than the glacier AWS values that they have replaced and, because this affects MS in the albedo model, K* has changed slightly. Error measures have changed slightly as well, but the outcome was achieved by setting Kb  =  −0.0051 m s−1 and constraining QH + QE to 0 for Kkat + Kb < 0.

To use a negative Kb to achieve a fit is to adopt the physically untenable position that the geostrophic momentum flux is upward from the ground. Therefore, optimization was done again, this time setting Kkat  =  0 and Kb  =  0.0026 m s−1 to minimize the RMSE (Table 2, step 5). Now QH and QE are closer to the glacier AWS values than they were before and the RMSE measures are reduced to the levels that apply to step 1. Furthermore, the hourly QH pattern is closer to that of the bulk transfer estimates but well short of agreement (Fig. 9), as further discussed below. The key point, however, is that Kb is now positive and its order of magnitude is similar to that of values used elsewhere (Klok and Oerlemans, 2002; Munro, 2004).

Figure 9

A 10-day sample of sensible heat flux model results from Equation 4a, using glacier AWS data (dots) compared to Equation 7a results for steps 4 (light line) and 5 (heavy line), using off-glacier AWS data. Inset shows a sample of wind directions at the 1991 site.



One should not place undue emphasis upon the relative sizes of the RMSE measures for glacier and non-glacier AWS model outcomes, or the size of RMSES relative to RMSEU. Although the RMSE for the final off-glacier model result is slightly smaller than that of the glacier model result, the former benefits from optimization of the turbulent transfer terms according to the RMSE but the latter does not. This has the effect of allowing the choice of QH and QE to absorb the effects of errors in other aspects of the model. Had the results of the first three steps in changing to the off-glacier model not turned out as well as they did, departure of the turbulent transfer values from those of the glacier AWS could have been greater than indicated in Table 2.

The RMSES can be reduced by constrained regression, the effect upon the outcome in Table 2, step 5 being a slope of ∼1.01 and a RMSES  =  4 mm. In fact, constraining the intercept to zero for each outcome listed in Table 2 resulted in a slope close to one. The fact that each of the outcomes has a substantial intercept stems from the obvious divergence between measured and modeled ablation that is seen in the middle of the old snow decay period (Fig. 8). Perhaps this reflects the density variation within the snowpack alluded to above because such variation is commonly observed in mass balance work on the Place Glacier and elsewhere (Østrem and Brugman, 1991), but it is beyond the scope and information base of this study to include it in the modeling approach.

The off-glacier modeled α results also warrant discussion because, in addition to these being smaller than the glacier modeled α values near the end of the old snow cover period, the difference between them increases from ∼0.01 on day 239 to ∼0.03 on day 243, the last full day of old snow (Fig. 3). This invoked the caution stated in Brock et al. (2000) that the structure of Equation 11 allows errors in α to propagate through MS. Therefore, each of steps 4 and 5 in Table 2 was initiated with the α pattern of Figure 3, continued with optimization of Kb according to the RMSE, and completed with inclusion of the off-glacier α model.

Error propagation in the albedo model could, in distributed modeling terms, translate to errors in snowline migration. So the question was explored further, by reiterating the α values obtained at the conclusion of step 5. The α reduction of ∼0.03 for day 239 changed to ∼0.04 after the first iteration, where it remained through nine more iterations. The changes to new snow were greater, with α for day 247 rising to its limit of 0.9 by the fourth iteration, α for day 251 falling from ∼0.8 to ∼0.5. Given the small mass of the fresh snow cover, the effect on Σ(δh) was small, changing its value from −1.196 to −1.225 m. Therefore, with the caveat of greater error sensitivity for thin, fresh snow cover, it seems that Equation 11 is a good choice for distributed modeling, as exemplified in the latest version of the Storglaciären model (Reijmer and Hock, 2008).

More problematical for distributed modeling are the outcomes of steps 4 and 5 because they stem from attempting to fit a continuous modeling structure to what is essentially an intermittent flow regime. Parameterization of Kkat according to air temperature invokes the assumption of turbulent mixing even if the wind has nearly ceased, as is the case for all of day 232 and many hours on other days of the sample displayed in Figure 9. If the assumption is wrong, the error is compounded in Kkat as well as in the air temperature and humidity differences in Equation 7, the result being QH and QE outcomes that are at least double the values listed in Table 2, even with Kb  =  0. The results look more realistic with optimization according to a fixed value for Kb, with Kkat  =  0 (Fig. 9), but even so there is the inherent assumption of continuous flow and therefore little agreement with bulk transfer results. It is possible that the model could be improved by incorporating an “on/off flow switch,” but glacier-based wind speed data would be required, thus defeating the goal of reliance upon off-glacier data.

Furthermore, it may not be the case that the wind regime at one glacier site typifies the whole glacier. The data suggest that Place Glacier may be too small to sustain a katabatic flow regime, but a different picture emerges from wind direction data retrieved from the micrometeorological study conducted further toward the glacier terminus in 1991 (Fig. 9, inset). The data show marked deviations of wind direction from the expected direction of glacier winds, extremely so on days 230 and 236 of that summer. The deviations were often associated with notable reductions in wind speed (not shown), but with none that were below 0.75 m s−1. This observation, and the seemingly glacier-like winds recorded at the off-glacier AWS, which also show extreme directional deviation (Fig. 6), suggest that there is more to consider for the wind regime at the glacier AWS than the failure of a small glacier to maintain a katabatic flow regime.

If the off-glacier AWS is recording katabatic winds, the glacier AWS wind record could be interpreted in terms of insufficient fetch down-glacier for katabatic flow to be established. The contrasting summer wind regimes of the two sites would then be interpreted in terms of insufficient fetch at 2044 m but sufficient fetch at 1840 m, notwithstanding the 500 m of separation between the glacier tongue and the off-glacier AWS. The 1991 results would fit well into such a picture. However, without wind direction data at 2044 m and concordant timing of measurements further down the glacier, the picture is inconclusive.

Also to consider is the unusual shape of the glacier: a small upper glacier area that feeds a larger lower glacier area comprised of two tongues. The glacier AWS appears to be sufficiently close to where airflow streams might diverge to each tongue to suggest that conflicting flow directions could also be a factor to suppress wind speed. It is also possible that the small size of the glacier makes it particularly susceptible to topographic steering of entrained regional airflow, such as to conflict with the establishment of the glacier wind. In this regard, it is noteworthy that the summer and winter wind direction regimes at the off-glacier site appear to be much the same (Fig. 6).

Concluding Remarks

If the results are typical of what to expect of small glacier wind regimes, then they underline what is already well known: that good radiation modeling is central to the task of effective distributed modeling of glacier ablation. The keys to this aspect of the task are accurate digital elevation models to handle the radiation geometry and effective representation of the surface albedo field. It appears that Equation 11, with suitable modification of its optical depth scale, is suitable for the task of surface albedo modeling. Given the small size of the turbulent transfer terms on the Place Glacier, radiation modeling alone might suffice as a first approximation to its surface melt field.

The addition of turbulent transfer to gain a better approximation poses a problem for distributed modeling because it appears that Equation 7 will not work without modification and that, given the constraint of working with off-glacier data, an “on/off switch” is not a tenable modification. Two possibilities do present themselves, depending upon whether glacier size per se, or fetch of cold air drainage is the controlling agent:

  1. If glacier size is the controlling agent, such as to prevent katabatic flow over the glacier, then Kb alone could be used over the whole extent of the glacier surface to achieve a fit with observed melt patterns.

  2. If down-glacier fetch controls katabatic flow, a suitable strategy might be to apply Equation 7 with Kb + Kkat below a threshold elevation that signifies sufficient down-glacier fetch for katabatic flow to be sustained, but with Kb alone at higher elevations. In either case, successive patterns of diminishing glacier snow cover could be used to guide the modeling effort (e.g., Marosz-Wantuch, 2004).

The nature of the base AWS record is also important because it should sample the characteristics of the regional air mass. If in fact it is really sampling a glacier wind modified regime, there is circularity in using katabatically modified air to model katabatic air flow. This can be made to work, but to the extent that modifications needed to make it work are different from those used elsewhere, it will not be clear whether it is because of regional or local influences.


Support for this work has been provided by Environment Canada through its Cryospheric Systems program, the Natural Sciences and Engineering Research Council of Canada, and, currently, the Canada Foundation for Climate and Atmospheric Science. Infrastructure support has been provided through the Canadian Glacier Variations Monitoring and Assessment Network, headed by M. N. Demuth, Natural Resources Canada. The authors have also benefited from collaboration with R. D. Moore, University of British Columbia, as well as from the constructive comments of F. Pellicciotti, ETH, Switzerland, and one other reviewer.

References Cited


E. L. Andreas 2002. Parameterizing scaler transfer over snow and ice: a review. Journal of Hydrometeorology 3:417–432. Google Scholar


Anonymous 1956. Snow Hydrology Summary Report of the Snow Investigations, North Pacific Division of Corps of Engineers, U.S. Army, Portland, Oregon. 437. Google Scholar


B. W. Brock, I. C. Willis, and M. J. Sharpe . 2000. Measurement and parameterization of albedo variations at Haut Glacier d'Arolla, Switzerland. Journal of Glaciology 46/155:675–688. Google Scholar


L. Cheng 2001. Validation of quantitative precipitation forecasts during the intermountain precipitation experiment. M.Sc. Thesis. University of Utah, 105 pp. ⟨⟩. Google Scholar


J. Dozier and J. Frew . 1990. Rapid calculation of terrain parameters for radiation modeling from digital elevation data. IEEE Transactions of Geosciences and Remote Sensing GE28/5:963–969. Google Scholar


B. E. Goodison, H. L. Ferguson, and G. A. McKay . 1981. Measurement and data analysis. In D. M. Gray and D. H. Male . Handbook of Snow. Oxford Pergamon. 191–274. Google Scholar


W. Greuell, W. H. Knap, and P. C. Smeets . 1997. Elevational changes in meteorological variables along a mid-latitudinal glacier during summer. Journal of Geophysical Research 102/D22:25941–25954. Google Scholar


R. W. Hamming 1989. Digital Filters. New York Dover. 284. Google Scholar


D. M. Hannah and A. M. Gurnell . 2001. A conceptual linear reservoir runoff model to investigate melt season changes in cirque glacier hydrology. Journal of Hydrology 246:123–141. Google Scholar


R. Hock and B. Holmgren . 2005. A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden. Journal of Glaciology 51/172:25–36. Google Scholar


E. J. Klok and J. Oerlemans . 2002. Model study of the energy and mass balance of Morteratschgletscher, Switzerland. Journal of Glaciology 48/163:505–518. Google Scholar


T. Konzelmann, R. S. W. van de Wal, J. W. Greuell, R. Bintanja, E. A. C. Henneken, and A. Abe-Ouchi . 1994. Parameterization of global and long-wave incoming radiation for the Greenland ice sheet. Global Planetary Change 9/1–2:143–164. Google Scholar


M. Marosz-Wantuch 2004. Modelling snowline migration and runoff response for Place Glacier Basin. M.Sc. thesis. University of Toronto. 107. Google Scholar


D. S. Munro 1989. Surface roughness and bulk heat transfer on a glacier: comparison with eddy correlation. Journal of Glaciology 35/121:343–348. Google Scholar


D. S. Munro 1990. Comparison of melt energy computations and ablatometer measurements on melting ice and snow. Arctic and Alpine Research 22:153–162. Google Scholar


D. S. Munro 2004. Revisiting bulk heat transfer on the Peyto Glacier in light of the OG parameterization. Journal of Glaciology 50/171:590–600. Google Scholar


D. S. Munro and G. J. Young . 1982. An operational net short-wave radiation model for glacier basins. Water Resources Research 18/2:220–230. Google Scholar


J. Oerlemans 2000. Analysis of a 3 year meteorological record from the ablation zone of Morteratschgletscher, Switzerland: energy and mass balance. Journal of Glaciology 46/155:571–579. Google Scholar


J. Oerlemans and B. Grisogono . 2002. Glacier wind and parameterization of the related surface heat flux. Tellus 54A/5:440–452. Google Scholar


J. Oerlemans and E. J. Klok . 2002. Energy balance of a glacier surface: analysis of automatic weather station data from the Morteratschgletscher, Switzerland. Arctic, Antarctic, and Alpine Research 34/4:477–485. Google Scholar


J. Oerlemans and W. H. Knap . 1998. A 1 year record of global radiation and albedo in the ablation zone of Morteratschgletscher, Switzerland. Journal of Glaciology 44/147:231–238. Google Scholar


T. R. Oke 1987. Boundary-Layer Climates. London Methuen. 435. Google Scholar


G. Østrem and M. Brugman . 1991. Glacier Mass Balance Measurements. National Hydrology Research Institute Report No. 4, Supply and Services Canada. 224. Google Scholar


S. Pohl, P. Marsh, and G. E. Liston . 2006. Spatial-temporal variability in turbulent fluxes during spring snowmelt. Arctic, Antarctic, and Alpine Research 38/1:136–146. Google Scholar


C. H. Reijmer and R. Hock . 2008. Internal accumulation on Storglaciären, Sweden, in a multi-layer snow model coupled to a distributed energy- and mass-balance model. Journal of Glaciology 54/184:61–72. Google Scholar


F. Sauberer 1955. Zur Abschätzung der Globalstrahlung in verschiedenen Höhenstufen der Ostalpen. Wetter Leben 7:22–29. Google Scholar


P. W. Suckling and J. E. Hay . 1977. A cloud layer-sunshine model for estimating direct, diffuse and total solar radiation. Atmosphere 15:194–207. Google Scholar


C. J. Willmott, S. G. Ackleson, R. E. Davis, J. J. Fedderma, K. M. Klink, D. R. Legates, J. O'Donnell, and C. M. Rowe . 1985. Statistics for the evaluation and comparison of models. Journal of Geophysical Research 90/C5:8995–9005. Google Scholar
D. Scott Munro and Marzena Marosz-Wantuch "Modeling Ablation on Place Glacier, British Columbia, from Glacier and Off-glacier Data Sets," Arctic, Antarctic, and Alpine Research 41(2), 246-256, (1 May 2009).
Accepted: 1 January 2009; Published: 1 May 2009

Back to Top