Glacier subsurface heat-flux characterizations for energy-balance modelling in the Donjek Range, southwest Yukon, Canada

Abstract

We apply a point-scale energy-balance model to a small polythermal glacier in the St Elias Mountains of Canada in order to investigate the applicability and limitations of different treatments of the glacier surface temperature and subsurface heat flux. These treatments range in complexity from a multilayer subsurface model that simulates snowpack evolution, to the assumption of a constant glacier surface temperature equal to 0°C. The most sophisticated model includes dry densification of the snowpack, penetration of shortwave radiation into the subsurface, internal melting, refreezing of percolating meltwater and generation of slush layers. Measurements of subsurface temperature and surface lowering are used for model validation, and highlight the importance of including subsurface penetration of shortwave radiation in the model. Using an iterative scheme to solve for the subsurface heat flux as the residual of the energy-balance equation results in an overestimation of total ablation by 18%, while the multilayer subsurface model underestimates ablation by 6%. By comparison, the 0°C surface assumption leads to an overestimation of ablation of 29% in this study where the mean annual air temperature is about −8°C.

Introduction

Snow- and ice-melt modelling in alpine watersheds has a variety of important applications, from predictions of seasonal runoff, streamflow and freshwater-resource availability (e.g. Moore and Demuth, 2001; Lemke and others, 2007) to estimations of sea level associated with climate change (e.g. Raper and Braithwaite, 2006). Since 1961, glaciers in southeast Alaska and the Coast Mountains are thought to have contributed more to sea-level rise than any other glaciological source (Kaser and others, 2006), though volume estimates of this contribution have recently been revised (Berthier and others, 2010). Glaciers in the St Elias Mountains are among the most significant contributors in this region (Arendt and others, 2002). Despite these high rates of mass loss, few detailed studies of glacier energy and mass balance have been undertaken in the St Elias Mountains.

Accurate models of glacier energy balance are essential for understanding the detailed physical mechanisms linking climate and glaciers, and for predicting future climate-driven glacier change. Energy- and mass-balance modelling on cold (polar) glaciers requires consideration of sub-freezing surface temperatures and the heat flux into the glacier (e.g. Greuell and Oerlemans, 1986; Konzelmann and Braithwaite, 1995; Klok and others, 2005). Even on temperate glaciers, however, changes in the temperature of the glacier surface layer can be associated with the storage and release of substantial amounts of energy, which can reduce total ablation, and delay the seasonal and daily onset of both ablation and streamflow (e.g. Kattelmann and Yang, 1992; Hock, 2005; Pellicciotti and others, 2008). Serious questions remain regarding the necessity of taking sub-freezing surface and englacial temperatures into account for glacier melt modelling in regions such as the St Elias Mountains, where both temperate and polythermal glaciers exist (Clarke and Holdsworth, 2002).

The main objective of this work is to evaluate the applicability of glacier surface temperature and subsurface heat-flux characterizations for point-scale energy-balance modelling on glaciers in the Donjek Range, St Elias Mountains. We test four different treatments of the glacier surface/subsurface within an energy-balance model using meteorological and melt data collected in 2008 on an unnamed valley glacier. The treatments range in complexity from one that employs a multilayer subsurface model to one that assumes a constant glacier surface temperature of 0°C. We use cumulative ablation and measured ice-temperature profiles to evaluate model performance. We then compare surface temperatures, hourly ablation and daily and total energy fluxes simulated with the various models. While previous studies have compared the results of energy-balance models that include different treatments of surface temperature and the glacier heat flux (e.g. Reijmer and Hock, 2008; Pellicciotti and others, 2009), the combination of models considered here is unique.

Background

Numerous treatments of glacier surface and subsurface temperatures have been incorporated in energy-balance models, but the applicability of a given method for a wide range of situations is not well established. A constant, 0°C glacier surface temperature is often assumed in energy-and mass-balance studies, and seems to provide satisfactory results for many glaciers (e.g. Hock and Noetzli, 1997; Brock and others, 2000; Oerlemans, 2000; Willis and others, 2002). In other studies, conceptual approaches that compensate for energy deficits before allowing melt (e.g. Braun and Aellen, 1990; Van de Wal and Russell, 1994) or iteratively adjust the glacier surface temperature in order to balance energy deficits (Escher-Vetter, 1985; Braithwaite and others, 1998; Braun and Hock, 2004; Hock and Holmgren, 2005), have been assumed to be sufficient to account for changes in surface temperature, ablation onset time and total ablation on temperate or near-temperate glaciers. In order to simulate englacial temperatures and the conductive heat flux into a glacier, subsurface models have been applied (e.g. Klok and Oerlemans, 2002; Corripio, 2003). An additional level of complexity can account for the transfer of heat related to percolation and refreezing of surface meltwater (e.g. Brun and others, 1989; Jordan, 1991; Greuell and Konzelmann, 1994; Bartelt and Lehning, 2002; Andreas and others, 2004). Huang (1990) and Braithwaite and others (2002) suggest that where mean annual air temperatures are below −6 to −8°C, significant refreezing and percolation of meltwater can occur.

When energy-balance models are extended from the point to the glacier scale (e.g. Arnold and others, 1996; Hock and Noetzli, 1997; Klok and Oerlemans, 2002), many variables cannot be measured and must be parameterized in each model gridcell. Some of these variables, namely the turbulent heat fluxes and the outgoing longwave radiation, are highly sensitive to assumptions made about the glacier surface and subsurface temperature conditions. It is therefore advisable to evaluate the validity of any such assumptions at the point scale before attempting to extend models in space.

Study Area and Field Measurements

The study site is located at 2280 m a.s.l. on a small valley glacier in the Donjek Range (60°50′ N, 139°10′ W) of the St Elias Mountains (Fig. 1a). Separated from the Gulf of Alaska by the highest peaks in the St Elias Mountains, the Donjek Range has a notably continental climate despite being <100 km from the coast. The study glacier is ∼5.3 km² in area and ranges in altitude from 1970 to 2960 m a.s.l. (Fig. 1b). Its present equilibrium-line altitude is ∼2550 m a.s.l. (Wheler, 2009). It has a southerly aspect over most of its length and occupies a valley on the southeast side of the range crest. It has a history of surging and is currently thought to be undergoing a slow surge (De Paoli and Flowers, 2009). The mean annual air temperature at the study site (∼2300 m a.s.l.), in the mid-ablation area of the glacier, is about −8°C. The glacier internal temperature structure appears to be weakly polythermal (De Paoli and Flowers, 2009). It is thus unclear to what extent the glacier heat flux and the variability in glacier surface temperature will impact the surface energy balance and melt regime.

Fig. 1. Study site. (a) Donjek Range study area within the St Elias Mountains of southwest Yukon with study glacier outlined. (b) Contoured surface elevation of study glacier (20 m contour interval) with locations of the ablation stakes (dots) and the AWS (×). Image in (a) provided through NASA’s Scientific Data Purchase Project and produced under NASA contract by Earth Satellite Corporation.

An automatic weather station (AWS) was installed in the ablation area as part of a field program initiated in summer 2006 at this site. In 2008 an albedometer was added to the AWS. Meteorological instruments, installed at a nominal height of 2 m above the surface, measure variables at 30 s intervals, and a Campbell CR1000 data logger records 5 min averages of these quantities (Table 1). An ultrasonic depth gauge (USDG) is installed on a structure drilled into the ice several metres away from the AWS and provides half-hourly surface-lowering data used for model validation. Englacial temperatures were measured daily at 1 m intervals down to ∼12 m depth at the AWS site in July–October 2008. These data were collected using a custom cable embedded with digital temperature sensors and comprise an additional source of model validation. The cable failed in October 2008 for reasons not yet understood.

Table 1. Instrumentation used in this study. Precision taken from manufacturer documentation. SW: shortwave

We have constructed a map of ice-surface elevation from kinematic GPS surveys conducted in 2006–07, and we maintain a network of 18 ablation stakes on the study glacier to measure mass balance (accumulation and ablation) (Fig. 1b). Profiles of snowpack density and temperature at the AWS location were determined from snow pit measurements made in May 2008. Densities and temperatures were measured at 10 cm depth intervals in the snow pit using, respectively, a 250 cm³ stainless-steel wedge and spring scale, and a 5 in dial-stem snow thermometer. The necessary data for energy-balance modelling are recorded at the AWS and comprise a dataset covering a complete ablation season in 2008 (Fig. 2).

Fig. 2. Hourly meteorological variables at a nominal height of 2 m above the glacier surface at the AWS in 2008. (a) Net shortwave radiation, SW_net. (b) Incoming longwave radiation, LW_in. (c) Air temperature, T _z. (d) Wind speed, u _z. (e) Vapour pressure, e _z. T _z and u _z are measured directly, e _z is calculated from the measured barometric pressure, SW_net is the sum of measured SW_in and SW_out, and LW_in is calculated from Q _N = SW_net+LW_in−LW_out, where Q _N is measured directly and LW_out is calculated by model M1.

Energy-Balance Model Formulation

Conservation of energy requires that incoming energy be balanced by the amount outgoing and the amount consumed. The energy balance over a snow or ice surface can be written as (e.g. Hock, 2005):

(1)

where Q _N is net radiation, Q _R is the sensible heat supplied by rain, Q _H and Q _L are the turbulent fluxes of sensible and latent heat, respectively, Q _G is the glacier or subsurface heat flux and Q _M is the energy available for melt. Heat fluxes toward the surface are considered positive (downward for Q _N, Q _H and Q _L and upward for Q _G). The small quantity of rain observed at the study sites, combined with the low air temperatures, allows us to neglect the heat flux from rain, Q _R.

We evaluate four model formulations, M1–M4, in which the glacier surface temperature, T _s, and glacier heat flux, Q _G, are treated in an increasingly simple manner (Table 2). M1 includes a multilayer subsurface model (MSM) in order to compute the glacier heat flux, Q _G, and the surface temperature, T _s (following Greuell and Konzelmann, 1994). The MSM presented by Greuell and Konzelmann (1994) has been successfully applied to studies of energy and mass balance at point locations, at 30 m spatial resolution across the small valley glacier Storglaciären, Sweden (Reijmer and Hock, 2008), and at 20 km resolution across the Greenland ice sheet (Bougamont and others, 2005). Of the models used in this study, the MSM (included in M1) has the strongest physical basis and is thus used as a reference against which to compare the results obtained with other models (M2, M3, M4). Below we describe our treatment of each component of the energy balance in the various models. Unless otherwise noted, the treatment applies to all models.

Table 2. Model formulations M1–M4, distinguished by their treatments of glacier surface temperature, T _s, and glacier heat flux, Q _G: multilayer subsurface model (MSM), residual Q _G method, iterative temperature scheme (ITS) or T _s = 0°C assumption

Shortwave, longwave and net radiation (Q _SW, Q _LW and Q _N)

Net radiation is directly measured in this study by a Kipp & Zonen NR-LITE net radiometer. The spectral response of this instrument has two peaks (short- and longwave) between 0.3 and 30 μm. Incoming and outgoing shortwave radiation are also measured independently by a Kipp & Zonen CMA6 albedometer (spectral range 310–2800 nm). Both radiation sensors are installed to be horizontal and are relevelled periodically; the albedometer was installed in May 2008, and both sensors were relevelled on 14 July 2008, at which time the albedometer was found to be <1° from horizontal. Measurements of shortwave radiation have been adjusted for the glacier surface slope (southerly aspect, 6° slope). To do this, the potential direct solar radiation is first computed separately for a flat surface and for a surface of 6° slope, determined from the glacier surface digital elevation model. This value of the slope is within the 5–7° range measured with a compass inclinometer at the AWS location. The difference between potential direct solar radiation computed for the AWS location for the horizontal surface and a 6° surface is added to the measured value of incoming shortwave radiation at 5 min intervals. Details of this correction, including the method used to compute potential direct solar radiation, are given by Wheler (2009). This difference averages to 7.8 W m⁻² over the melt season. Note that the incoming shortwave radiation in early morning and late evening for a 6° surface is less than that for a flat surface, due to the northerly component of incoming radiation.

We do not measure longwave radiation directly. For M1, the net longwave radiation, Q _LW, is determined by subtracting the measured net shortwave radiation, Q _SW, from the measured net all-wavelength radiation, Q _N. Incoming longwave radiation, for all models, M1–M4, is then computed as the difference between Q _LW and the outgoing longwave radiation simulated by M1. Outgoing longwave radiation is calculated individually for each model as a function of the modelled glacier surface temperature, according to the Stefan–Boltzmann law: , where ∊ is the emissivity of the surface which we take as 1.0 (Arnold and others, 1996; Hock and Holmgren, 2005), σ = 5.6704 × 10⁻⁸ J s⁻¹ m⁻² K⁻⁴ is the Stefan–Boltzmann constant and T _s is the surface temperature.

Turbulent heat fluxes, Q _H and Q _L

We compute the turbulent fluxes of sensible (Q _H) and latent (Q _L) heat using the bulk aerodynamic method (Munro, 1989). The glacier surface temperature and vapour pressure are either computed or assumed constant (0°C and 611 Pa), depending on the surface/subsurface model applied. The turbulent fluxes are thus calculated based on Monin–Obukhov similarity theory using a single level of wind speed, temperature and relative humidity measurements by the equations:

(2)

(3)

where ρ _a is the density of air computed from the measured air pressure, P, c_p = 1010 J kg⁻¹ K⁻¹ is the specific heat capacity of air, L _v = 2.514 × 10⁶ J kg⁻¹ is the latent heat of vaporization, and C _H and C _E are the turbulent exchange coefficients for heat and vapour pressure, respectively. The wind speed, u _z, temperature, T _z, and air pressure, P, are measured at height z (nominally 2 m). T _s and e _s are the temperature and vapour pressure at the surface and e _z is the 2 m vapour pressure computed from measured relative humidity.

The turbulent exchange coefficients C _H,E are computed as (e.g. Hock, 2005):

(4)

where k = 0.4 is the von Kármán constant, L is the Monin–Obukhov length, Ψ_M,H,E are the integrated stability constants and z _M,H,E are the roughness lengths for momentum, heat and water vapour, respectively. We adopt a value of z _M = 1 mm for both snow and ice, and set z _H = z _E = z _M/100 (e.g. Hock and Holmgren, 2005). Values of z _M in the range 0.01–1.1 mm for ice (mean of 0.65 mm) and 0.015–1.8 mm for snow (mean of 0.38 mm) have been measured at the study site (MacDougall and Flowers, in press). Order-of-magnitude changes in the value of z _M have little effect on our results. For example, varying the value of z _M for ice from 0.5 mm to 10 mm yields only a 0.06 m (or 3%) difference in calculated total ablation over the melt season. This suggests a low model sensitivity to z _M; however, a more robust sensitivity analysis is required to confirm this. The turbulent heat fluxes are small contributors to the energy balance in our study area, but their importance in the modelled energy balance could increase, for example, if stability corrections were not applied.

C _H,E are computed for stable atmospheric conditions using the functions presented by Beljaars and Holtslag (1991) and for unstable conditions using the nonlinear functions of Paulson (1970) and Dyer (1974). In order to compute the stability corrections, an iterative procedure is employed following Munro (1990). To initiate the iteration, neutral atmospheric conditions are assumed and preliminary values of the turbulent exchange coefficients, C _H,E, and the sensible and latent heat fluxes, Q _H and Q _L, are computed. Nonzero values of the stability correction functions are then estimated and used to correct C _H,E for atmospheric stability and compute corrected values of Q _H and Q _L. Full correction is achieved when the value of Q _H stabilizes, which we define as a mean change from the previous step of <0.1 W m⁻² and a maximum change of <0.5 W m⁻².

Melt energy, Q _M

The melt energy, Q _M, is computed as the residual of the energy balance in Equation (1). The simulated water-equivalent (w.e.) ablation rate, M _s (m w.e. s⁻¹), is determined by:

(5)

where ρ _w = 1000 kg m⁻³ is the density of water and L _f = 3.34 × 10⁵ J kg⁻¹ is the latent heat of fusion. Only positive values of Q _M are considered in the ablation-rate calculation, but both negative and positive values of Q _L are considered. Both sublimation/deposition and evaporation/condensation are expected to occur, so the application of the latent heat of vaporization is an arbitrary choice (following e.g. Oerlemans, 2000). Depending on the partitioning of the negative latent heat flux between evaporation and sublimation, this assumption may result in ablation differences of a few millimetres to centimetres water equivalent for the dataset under consideration.

Glacier surface temperature, T _s, and heat flux, Q _G

M1: multilayer subsurface model

There are a variety of ways that heat is transferred from the surface into the glacier. In addition to refreezing meltwater, shortwave radiation penetrates below the surface (Paterson, 1994) and conduction through the snow/ice occurs wherever a temperature gradient exists. We follow the approach used in the SOMARS (Simuation Of glacier surface Mass balance And Related Sub-surface processes) subsurface model of Greuell and Konzelmann (1994), which takes all of these heat transfer mechanisms into account.

For the one-dimensional MSM, the following thermodynamic equation is solved in each layer of the subsurface grid (Greuell and Konzelmann, 1994):

(6)

with temperature T, material density ρ and specific heat capacity c = 152.5 + 7.122(T + 273.15) (Paterson, 1994). The term on the right-hand side of the equation is the conductive heat flux, Q _C, where K _e is the effective conductivity and ∂T/∂z is the subsurface temperature gradient. Quantity ∂Q _a /∂z is the absorption of energy from the atmosphere, where Q _a = Q _N + Q _H + Q _L is the sum of the radiative and turbulent fluxes. M and F are the melt and refreezing rates, respectively, and L _f is the latent heat of fusion.

To implement the MSM we establish a grid comprising 25 subsurface layers, with the upper boundary defined as the snow or ice surface and the lower boundary set to a depth of ∼40 m (i.e. well below the depth of ice affected by seasonal temperature fluctuations). The thickness of the subsurface layers varies from 0.04 m near the surface to ∼5 m for the deepest layers. Precise layer thicknesses are adjusted so the depth of the snow–ice interface at the beginning of the melt season (i.e. the previous year’s summer surface) coincides with a boundary between two subsurface layers (following Greuell and Konzelmann, 1994).

Glacier surface conditions at the beginning of the melt season must be prescribed to initialize the model. The simplified MSM applied here requires only initial snow depth and profiles of density and temperature, which we obtain from snow pit measurements made adjacent to the AWS in May 2008. Following Greuell and Konzelmann (1994), a constant initial density equal to the integrated snow pit density is assigned to all snow layers, and a constant ice density, ρ _i, of 900 kg m⁻³ is assigned to the ice. Measured snow pit temperatures are linearly interpolated onto the subsurface model grid.

Initial ice temperatures are linearly interpolated onto the grid between the measured 12 m ice temperature (recorded in September 2008) and the measured temperature at the snow–ice interface at the beginning of the ablation season (−10°C in early May 2008). With the initial subsurface temperature and density profiles defined, the surface temperature, T _s, is calculated by linear extrapolation of the temperature of the two uppermost model layers. This extrapolation is preferable to setting T _s equal to the temperature of the uppermost layer, because of the large temperature gradients just below the surface (Greuell and Konzelmann, 1994). According to Reijmer and Hock (2008), in order to obtain the most accurate surface temperature, the two uppermost layers must be of the order of a few centimetres thick; we use thicknesses of 4 and 5 cm. The surface temperature obtained from the energy balance at a given time-step is then used to calculate the outgoing longwave radiation, the turbulent heat fluxes (Equations (2) and (3)) and the surface energy balance (Equation (1)) for the following time-step. We found a negligible difference in results between using a linear versus shape-preserving spline function to compute the surface temperature.

The penetration of shortwave radiation beneath the surface is taken into account following Greuell and Oerlemans (1986) by assuming 36% is absorbed entirely at the surface (represented by the uppermost model layer). The radiation absorbed at the surface corresponds approximately to wavelengths >0.8 μm. The remaining 64% is assumed to penetrate below the surface and is absorbed by all layers (including the uppermost model layer) according to the Beer–Lambert law:

(7)

where S(z) is the shortwave radiation flux at depth z, Q _SW is the net shortwave radiation incident at the surface and κ. is the extinction coefficient. Snow/ice absorption of shortwave radiation is predicted to increase with material density up to ∼450 kg m⁻³ and then decrease with density from ∼450 to 900 kg m⁻³ (Bohren and Barkstrom, 1974). In the datasets presented here, the minimum snow density is ∼250 kg m⁻³ and we apply a constant extinction coefficient of κ. = 20 m⁻¹ for all snow densities up to 450 kg m⁻³ (for comparison, the maximum value of κ. given by Greull and Konzelmann (1994) was applied to the summer snow density of 300 kg m⁻³). Following Greuell and Konzelmann (1994), κ. is assumed to be a linear function of density from ∼450 to 900 kg m⁻³ such that κ. = 2.5 m⁻¹ at 900 kg m⁻³. Because the extinction coefficient is taken from the literature, we have not explicitly accounted for any site-specific effects related to contaminants or dust in the study glacier snowpack.

The major steps taken to implement the MSM, after the grid is established, are described below. The conductive heat flux depends on the effective conductivity, K _e, which governs conduction, convection, radiation and vapour diffusion processes. We compute K _e as a function of snow or ice density, ρ (kgm⁻³), following Sturm and others (1997):

(8)

Conduction between adjacent subsurface layers is then determined from K _e and the temperature gradient, ∂T/∂z.

If the computed temperature of any layer exceeds 0°C, the temperature is set to 0°C and the excess energy is applied to melt. Melting can occur in subsurface layers due to the penetration of shortwave radiation. When snow is present at the glacier surface, meltwater percolates until it reaches a layer with a subzero temperature and refreezes. The amount of refreezing in any given model layer is limited by the temperature, which cannot exceed 0°C, and the available pore space in the snow. An amount of meltwater equal to the irreducible water content is retained even if the refreezing limit has been reached, and the remaining liquid meltwater percolates down to the underlying layer. Following Reijmer and Hock (2008), we describe θ _mi, the ratio of the mass of irreducible liquid water to the total mass of the layer, based on an empirical relationship developed by Schneider and Jansson (2004) using data from Storglaciären and laboratory experiments (Coléou and Lesaffre, 1998):

(9)

(10)

(11)

where n is the porosity of the snow layer, ρ _tot is the total density including both solid and liquid, m _liq is the mass of liquid water, m _sn is the mass of snow and θ _m is referred to as the gravimetric liquid water content. Increases in the total density correspond to increases in porosity and to increases in the irreducible water content according to Equations (9) and (10). If the sum of the available meltwater (produced within a given layer or percolating from above) and the liquid water content present in the layer exceeds the irreducible water content, the excess percolates to the underlying layer.

Any meltwater that does not refreeze and is not retained as irreducible water content percolates to the base of the snowpack, where it runs off along the impermeable ice surface. Runoff is computed using a runoff coefficient for snow, c _s = 0.8, which is the ratio of runoff to available water (e.g. Seidel and Martinec, 2004). This means that 80% of the water runs off immediately, while 20% is retained in the snowpack at each time-step. This runoff is equivalent to the model ablation. Meltwater that does not immediately run off is used to build a saturated slush layer that extends upward from the base of the snowpack. In the dataset considered, the results are not particularly sensitive to the value of the runoff coefficient (values from 0.5 to 1.0 were examined), in part because snow makes a much smaller contribution than ice to the total ablation. Furthermore, the glacier surface slopes in the study area are >5°, producing relatively rapid runoff and little opportunity for the formation of slush layers. For situations in which the contribution of snow is more significant, or slopes more conducive to retaining water, the runoff coefficient may warrant more careful investigation.

In addition to changes in the density of snow layers due to melting and refreezing, which will dominate in a melting snowpack, small changes in layer densities can also occur due to dry densification. The dry densification rate, dρ/dt, is computed following the approach of Li and Zwally (2004), which is a modification of the empirical relationship developed by Herron and Langway (1980):

(12)

where K(T) = K ₀(T) exp(−E(T)/RT) is the temperature-dependent rate factor with universal gas constant R = 8.3144 J K⁻¹ mol⁻¹, rate constant K ₀ and activation energy E. E and K ₀ are functions of temperature, as described by Zwally and Li (2002). A is the mean accumulation rate representing the change in overburden pressure and is set to the initial snow layer thickness (Reijmer and Hock, 2008), α is ∼1 and J is the vapour flux. Changes in density due to the vapour flux are computed following Sturm and Benson (1997). Calculated surface lowering in the model thus includes contributions from dry densification.

In the final step, the subsurface grid is redefined so it is effectively shifted down to compensate for the change in thickness of the uppermost layer due to ablation. Summer snow accumulation is not accounted for in this version of the MSM. Mass (solid and liquid) and energy (i.e. cold content) are conserved when layers are shifted. The thickness of the bottom snow layer is incrementally reduced as snowmelt proceeds, until it is <1 cm, at which time the layer is merged with the overlying snow layer.

In summary, the MSM computes the glacier heat flux, Q _G, as the integrated change in cold content in all model layers at each time-step and facilitates calculation of outgoing longwave radiation, Q _LW, and the sensible (Q _H) and latent (Q _L) heat fluxes by estimating the surface temperature. Melting is computed directly from the energy balance in each model layer.

M2: solving Q_G as the residual energy flux

In this formulation, the glacier heat flux, Q _G, is solved as the residual of the energy-balance equation when the surface energy flux is negative (e.g. Klok and Oerlemans, 2002). The resultant changes in surface temperature, and related changes in the turbulent fluxes and the surface energy balance are then recomputed. For each 5 min time-step, the surface energy balance is initially calculated assuming T _s = 0°C. If the surface energy flux is negative, the glacier heat flux, Q _G, is computed as the residual required to achieve energy balance (Equation (1)) with Q _M. The temperature change in the 5 cm thick surface layer is then computed from Q _G as

(13)

where Δt is the time-step in seconds, and ρ _s, c _s and d _s are the density, specific heat capacity and thickness of the surface layer, respectively. Varying the layer thickness, d _s, affects the temperature, T _s, but does not affect the amount of energy associated with warming/cooling the surface layer, since the amount of energy is determined by Q _G. The surface energy balance is then recalculated with the surface temperature, T _s, assumed equal to the new temperature of the 5 cm surface layer. The iteration continues until the value of Q _G stabilizes, which usually occurs after two or three iterations.

This simple approach to estimating the glacier heat flux and the surface temperature is a compromise between the MSM (M1) and the iterative temperature scheme (M3). M2 extends the iterative temperature scheme (Braun and Hock, 2004) by taking into account the energy associated with temperature changes in a 5 cm thick surface layer and has some similarities to the two-layer subsurface model of Klok and Oerlemans (2002). M2 is expected to be most applicable on temperate or near-temperate glaciers where the conductive component of the glacier heat flux is small.

M3: iterative adjustment of surface temperature; M4

The third method we use to determine the glacier surface temperature, T _s, referred to as the iterative temperature scheme (ITS), is applied following Braun and Hock (2004). In this method we assume that T _s = 0°C when the surface energy flux (Equation (1)) is positive. If the surface energy flux is negative, Q _M is set to zero and the surface is assumed to cool (Escher-Vetter, 1985; Braithwaite and others, 1998). The surface temperature is then lowered iteratively in 0.25°C steps and the turbulent heat fluxes recomputed at each step. A lowered surface temperature increases the turbulent heat fluxes and decreases the outgoing longwave radiation in order to compensate for the negative surface energy flux, and the iteration continues until balance is achieved (Braun and Hock, 2004; Hock and Holmgren, 2005). The glacier heat flux, Q _G, is neglected in this method. Reijmer and Hock (2008) compare this method to the MSM of Greuell and Konzelmann (1994) and find that the latter reproduces observed surface temperatures with significantly greater skill. This is expected, since the MSM accounts for the relevant surface and subsurface processes in a physical manner. However, the potential advantage of the iterative scheme is that it does not require knowledge of the temperature or density structure of the ice or snowpack. M3 is differentiated from M2 in that it assumes an instantaneous adjustment of surface temperature, whereas M2 allows the temporary storage of heat in a shallow subsurface layer.

In the extreme case where wind speed u _z = 0, the turbulent heat fluxes will be zero and the surface temperature iteration will never converge. Surface temperatures are therefore capped at a minimum of −30°C for M3. We assume the surface energy is balanced when temperatures reach this threshold. When the melt energy is positive, M3 (ITS) is identical to M4, in which the surface temperature is assumed constant at 0°C. Melt rates simulated with these two methods are thus identical.

Results and Discussion

We present results for the 2008 melt season, beginning with a comparison of simulated and measured ice temperatures and cumulative ablation. This is followed by an intercomparison of model results, including modelled surface temperatures, hourly ablation and daily and total energy fluxes.

Subsurface temperature structure

The MSM applied in M1 computes temperature and density changes down to ∼40 m below the glacier surface. In this study, the primary purpose of the MSM is to calculate the glacier heat flux, as outlined above. The MSM provides a sound physical basis for modelling subsurface temperatures and densities, and yields better agreement with the observed englacial temperature profiles when shortwave penetration into the subsurface is included (Fig. 3).

Fig. 3. Subsurface temperature profiles measured (diamonds) and modelled with M1 (crosses) for 18 July (top row), 15 August (middle row) and 11 September (bottom row). Left column: M1 initialized on 5 May with initial temperature profile linearly interpolated between temperature measured at snow–ice interface on 5 May at 12 m ice depth. Centre column: M1 initialized on 18 July with measured ice-temperature profile. Right column: as in centre column, but with shortwave penetration neglected (i.e. all radiative energy absorbed entirely by uppermost model layer).

The temperature profile is captured within 1°C for the model initialized on 5 May with a simple temperature profile linearly interpolated between the measured temperature at the base of the snowpack and at 12 m ice depth (left column, Fig. 3). If the initial ice-temperature profile is instead assumed constant and equal to the measured temperature at the ice–snow interface at the beginning of the melt season, the modelled temperature profile is substantially colder. For this dataset, where the initial snowpack was very thin, the difference was up to ∼ 6°C from July to September; this difference is expected to decrease with an increasing thickness of the initial snowpack.

Based on visual inspection, good agreement with the observations is achieved when the measured temperature profile is used as an initial condition for a simulation beginning on 18 July (middle column, Fig. 3). It could be argued that the measured near-surface temperatures are better reproduced by the simulation shown in the left column of Figure 3. Neglecting shortwave penetration into the subsurface leads to modelled near-surface temperatures that depart substantially from the measured values late in the season (right column, Fig. 3). In this simulation, these departures extend to 5 m depth, with modelled temperatures being up to ∼2.5°C colder near the surface by 11 September 2008. The model neglecting shortwave penetration would presumably also under-predict near-surface temperatures in the early melt season.

Cumulative ablation

Figure 4 shows observed and simulated cumulative ablation from 5 May to 11 September 2008. The USDG ablation record is constructed from differences in daily average surface lowering using the integrated snow pit density for the period of snowmelt, and a density of 900 kg m⁻³ for the period of ice melt. Note that ablation of summer snowfall is included in the USDG ablation record, while snowfall is not explicitly included in the model. However, because the model uses measured albedo, the influence of summer snowfall on surface albedo is implicitly included. Figure 4 also shows several ablation measurements at the stake location (a few metres from the AWS) and their attendant uncertainties. These uncertainties reflect those estimated for individual stake measurements based on the complexity of the glacier surface in the immediate vicinity of the stake; they do not reflect the standard deviation between measurements or include uncertainty in the representativeness of the stake for a larger region.

Fig. 4. Cumulative ablation from 5 May to 11 September 2008 simulated by M1–M4, compared to ablation-stake measurements and daily average values of ultrasonic depth gauge (USDG) surface lowering. Initial snow depth is 0.41 m (0.103 m w.e.) and total ablation measured at the stake is 1.80 m w.e. while that measured by the USDG is 1.78 m w.e. M1 underestimates ablation by 6%, while M2, M3 and M4 overestimate ablation by 18%, 27% and 29%, respectively.

Significant differences in cumulative ablation as simulated by the models arise during the period of snowmelt prior to 10 June, with M1 underestimating ablation near the beginning of this period. Models M2–4 perform better than M1 from roughly 15 to 25 May; M1 performs significantly better than the other models immediately after this, from 25 May to 10 June. All models overestimate the USDG measured ablation rate similarly and significantly during the period 10–25 June. We interpret this mismatch as being due to the evolution of the glacier surface during this phase of the melt season: a heavily weathered and uneven cryoconite surface emerges from beneath the snow cover and later collapses. This results in depressed rates of surface lowering (and hence ablation) as recorded by the USDG prior to the cryoconite collapse. From about 25 May onward, M2–M4 generally over-predict ablation, resulting in modelled ablation that exceeds the measured value by 18%, 27% and 29% for M2, M3 and M4, respectively, by the end of the season (Table 3). M1 underestimates cumulative ablation by 6%, increasing to 10% if a uniform initial ice-temperature profile is assumed.

Table 3. Averages of simulated hourly energy fluxes (W m⁻²) for the entire 2008 ablation season, and for periods of melting only, with M1–M4: outgoing longwave radiation, L _out, sensible heat flux, Q _H, latent heat flux, Q _L, and glacier subsurface heat flux, Q _G. Also reported is the duration of surface melting in days, t _M; the ablation due to sublimation, M _QL (mw.e.); and the total ablation, M _TOT (m w.e.), which includes mass lost due both to melt and sublimation. If the melting period for M1 is calculated based on simulation of internal (subsurface) melting, rather than surface melting, L _out = −310, Q _H = 6, Q _L = −5, Q _G = −26 and t _M = 60.8

The subsurface model in M1 produces a reasonably accurate simulation of the timing of the snow-to-ice surface transition. According to the USDG record, the ice surface became temporarily snow-free on 7 June for a few hours, was covered by snow later that night, and again became snow-free on 10 June, after which ice remained exposed at the surface, aside from a few small snow events later in the melt season. The M1-simulated surface transition occurs on 11 June. The actual surface transition manifests itself in all the modelled ablation records as an inflection point, due to the increase in net radiation as the surface albedo drops substantially in the transition from snow to ice. The surface transition date is important for simulating changes in surface albedo, subsurface heat-transfer processes, supraglacial hydrology and associated impacts on the melt and runoff. These changes have a particularly large impact on spatially distributed melt modelling, where albedo, outgoing longwave radiation and glacier surface conditions are simulated rather than measured or observed (e.g. Arnold and others, 1996; Klok and Oerlemans, 2002; Hock and Holmgren, 2005).

It should be noted that the snowpack was relatively thin (0.41 m) in 2008 at the USDG location, compared with some other recent years where we have measured >1.5 m snow depth at the same site. For the thin snowpack, including dry densification in the model did not result in significant surface lowering. However, percolation and refreezing were significant and resulted in 0.1 m w.e. of internal accumulation; this is comparable to the winter accumulation for 2008. A deeper snowpack will provide a more informative test of the model in some respects; however, we have chosen to focus on the year for which we have subsurface temperature measurements. A future study will examine the performance of the MSM over several years, encompassing a range of initial snowpack conditions.

Surface temperature, T _s

Glacier surface temperatures simulated with M1–M3 are shown in Figure 5; M4 assumes T _s = 0°C. The mean difference between surface temperatures simulated by M1 initialized with a constant as opposed to a linearly varying initial temperature profile was 0.04°C; the maximum difference was 2.32°C. The modelled temperature records are distinct for M1–M3, with M1 consistently producing the highest daily temperature minima, often greater than −7°C during the mid-ablation season. M1 also produces multi-day periods of subfreezing surface temperatures, particularly late in the melt season. Sub-freezing temperatures rarely persist for more than a day in M2 or M3 after 15 May. Surface temperatures simulated with M3 are highly variable and exhibit the lowest daily minima. The arbitrary threshold of −30°C is rarely reached in M3 after 15 May, though the temperature frequently dips below −20°C.

Fig. 5. Modelled surface temperature, T _s, at the AWS in 2008. (a) M1. (b) M2. (c) M3.

While M1 and M2 both account for the thermal inertia of the snow or ice layers, hence limiting the magnitude of surface temperature variability, there is still a significant difference between including a full subsurface model (M1) and only a thin subsurface layer (M2). The penetration of shortwave radiation below the surface in M1 contributes to the higher modelled surface temperature minima: warmer subsurface layers are able to transfer heat toward the surface as it cools, thus limiting the surface temperature drop. The glacier heat flux, Q _G, is generally positive (upward), transferring heat from depth toward the surface, during subfreezing conditions at the surface.

Without accounting for thermal inertia and the energy released due to cooling of the surface/subsurface at night, M3 requires that the nightly surface temperature be adjusted to a relatively low value in order to compensate for a negative energy flux at the surface. The negative energy flux is the result of the negative net radiative flux, which occurs on most nights in this dataset. The smaller the energy deficit, the smaller the surface temperature adjustments required to achieve energy balance in M2 and M3, and the smaller the differences in simulated surface temperature minima between models. On a few nights, when the net longwave energy loss is small and winds are relatively high, surface temperature minima simulated with M1, M2 and M3 are similar (e.g. on 30 August). When the surface energy balance remains positive at night, the simulated surface temperature remains at 0°C with all models (e.g. on 5 July). The assumption of a constant surface temperature of 0°C does not necessarily have a large impact on simulated melt rates, since melting does not occur when T _s < 0°C. However, simulation of the energy fluxes during times when T _s < 0°C does have an impact on the timing of melt onset the following day.

Hourly ablation

To understand the detailed differences in simulated ablation between models, we examine hourly ablation rates and surface temperatures over 5 days in August (Fig. 6). The highest hourly mean nightly surface temperatures are obtained with M1 and the lowest temperatures with M3 (Figs 5 and 6b). The instantaneous response of M3 to changes in the surface energy balance results in large and rapid surface temperature fluctuations. Following cold night-time temperatures, melting temperatures at the glacier surface are often achieved earlier in the morning with M3, compared with M1 and M2 (Fig. 6a); the thermal inertia of the snow/ice layers modelled with M1 and M2 generally reduces the rate and magnitude of surface temperature variability. The conductive heat flux simulated with M1 is positive at night, which causes warmer surface temperatures and allows ablation to begin earlier on some days compared with M2 (e.g. 10–11 August).

Fig. 6. Hourly variables simulated for 5 days in August by M1, M2 and M3. (a) Ablation rates. Values are identical for M3 and M4. (b) Surface temperatures, T _s.

Although the timing of daily ablation onset can differ by a few hours between models, there is little difference in simulated peak daily ablation rates. When the surface is melting, simulated ablation rates are equal with all models except M1. Ablation rates are frequently less with M1 than the other models due to the negative glacier heat flux, and due to the refreezing and retention of meltwater during periods of snowmelt. This difference is larger during ice melt than snowmelt, but the difference between peak rates of ice melt simulated by M1 and the other models is generally small (Fig. 6a). However, these differences can accumulate to substantial amounts over time (Fig. 4).

Daily energy fluxes

Differences in modelled surface temperature are important, because they affect the calculation of the outgoing longwave radiation and the turbulent heat fluxes. In Figure 7 we compare the daily energy fluxes simulated with M1–M4 to examine the implications of varying treatments of T _s and Q _G on the energy-balance components. To facilitate comparison between models, in Figure 7b–d we plot the components as differences from those simulated by M1. The most notable features of Figure 7 are: (1) the large difference between M4 and the other models; (2) the similarity of most energy fluxes simulated by M2 and M3; and (3) the minimal discrepancy between M1 and M2 in the second half of the melt season.

Fig. 7. Daily mean energy fluxes: full components for M1 (a), as well as components expressed as differences from those calculated by M1 (b–d) (net longwave radiation, Q _LW, glacier heat flux, Q _G, sensible heat flux, Q _H, and latent heat flux, Q _L). Values are stacked rather than superimposed (i.e. the total length of each bar is the sum of the components). Only values >5 W m² are shown, as smaller values cannot be resolved in this figure. (a) M1. (b) M2 − M1. (c) M3 − M1. (d) M4 − M1.

The large difference in daily energy fluxes between M4 and the other models, especially at the beginning and end of the melt season, is due entirely to the assumption of a 0°C glacier surface temperature. Prior to 16 May the simulated surface temperature (Fig. 5a) rarely reaches 0°C, so M4 simulates much more outgoing longwave radiation than the other models. In contrast to M2 and M3, M4 generally underestimates net longwave radiation relative to M1 throughout the melt season. Differences in the simulated values of net longwave radiation are entirely a product of differences in calculated outgoing longwave, as incoming longwave has been prescribed identically for all models. Q _LW is generally higher for M2 and M3 relative to M1, as M2 and M3 predict lower surface temperatures (Fig. 5) and hence less outgoing longwave. Aside from during the early melt season, when temperatures are notably colder for M2, compared with M1, Q _LW is fairly similar between M1 and M2. Although minimum modelled surface temperatures are significantly different between M1 and M2, the daily mean modelled temperatures are often similar and thus produce similar daily mean values of Q _LW.

The turbulent heat fluxes, Q _H and Q _L, also distinguish M4 from the other models. The sensible heat flux, Q _H, is second in importance to the net radiative flux as an energy source, providing 15% of the daily energy on average. This percentage was calculated by separating the positive and negative contributions of each daily flux and dividing their magnitudes by the respective total positive and negative daily fluxes. The sensible heat flux provides most or all of the daily energy on several days early and late in the melt season. After mid-May, there is little difference in the sensible heat flux (which generally acts as a small source of energy) as simulated by M1–M3. In contrast, Q _H simulated by M4 is a significant energy sink (large and negative), both at the beginning and end of the melt season. Simulations of the latent heat flux with M1, M2 and M3 have very similar patterns of overall variability. All models yield negative fluxes of latent heat, but M4 frequently yields values 5–25 W m⁻² lower than the other models.

The daily mean glacier heat flux, Q _G (only computed in M1 and M2), is the most substantial energy sink (aside from melting, which is not shown in Fig. 7). On average, Q _G is twice the magnitude of the other consistent energy sink, Q _L. With M1, the glacier heat flux is a significant source of energy on only nine days in the melt season. Early in the melt season, the snowpack is rapidly warmed due to percolation and refreezing of meltwater. This warming can result in a positive glacier heat flux when the surface cools at night, but the net glacier heat flux is predominantly negative when daily averages are computed. Other than M1, the models presented here do not consider heat exchange between the surface and subsurface, aside from the thin, 5 cm surface layer considered in M2. Similar discrepancies in Q _G with respect to M1 between M2 and M3 (Fig. 7) suggest the treatment of the glacier heat flux in M2 is inadequate for the near-surface temperature conditions at the study site.

Although Q _G varies diurnally throughout the melt season in M1, there is a shift from values near zero (during snowmelt) to negative values (during ice melt) at the simulated snow-to-ice surface transition (not shown). Two factors contribute to this shift. The most important is the effect of meltwater percolation and refreezing during snowmelt. Refreezing meltwater rapidly warms the upper snowpack, resulting in an isothermal near-surface temperature profile. In the absence of a near-surface temperature gradient, the conductive heat flux is zero. Once the surface is snow-free (i.e. after the transition date), percolation no longer occurs in the model, and without the energy released by refreezing meltwater, subsurface temperatures remain below zero when the surface is melting. This establishes a near-surface temperature gradient in which heat is conducted away from the surface, into the glacier, and Q _G is thus negative.

The second factor contributing to the shift in Q _G from the period of snow- to ice melt is the conductivity itself. As defined in Equation (8), effective conductivity depends on density and is thus significantly greater in ice than in snow. For example, for typical snow and ice densities of 500 and 900 kg m⁻³, respectively, the conductivity of ice is more than four times that of snow. The increased conductivity facilitates more efficient energy exchange between the surface and subsurface layers. The combined results of changes in refreezing and conductivity described above are: relatively cold temperatures in the upper subsurface layers during ice melt (due to the lack of percolation and refreezing) and enhanced energy exchange between the cold subsurface layers and the melting surface (due to the increased conductivity).

The major difference between the glacier heat flux in M1 versus M2 is that M2 does not include conduction in the subsurface, but simply calculates changes in the temperature of the 5 cm surface layer. The glacier heat flux, Q _G, is solved as the residual in the surface-energy-balance equation in M2 when the sum of the radiative and turbulent fluxes is negative. This effectively lumps the energy that is physically associated with temperature changes in both the surface and subsurface layers into a single surface layer. While M2 captures the diurnal reversal of Q _G reasonably well (not shown), contributions of opposite sign tend to cancel one another over the course of a day in this model. This explains the significant discrepancy in daily values of Q _G between M1 and M2 in Figure 7.

Total energy fluxes

The mean contributions of each energy flux to the surface energy balance are given in Table 3. Considering the full period, the sensible (Q _H) and latent (Q _L) heat fluxes are similar for M1–M3, with values of 6–7 and −3 W m⁻², respectively. The glacier heat flux, Q _G, simulated with M1 is −6 W m⁻² when averaged over the full period of record, indicating a net flux of heat from the surface to depth. During the period of melting, Q _G averaged −30 W m⁻².

Energy fluxes computed with M3 and M4 are identical during melting conditions, as the two models are identical when the surface energy balance is positive. The total ablation simulated by M3 and M4 (2.26 and 2.30 m w.e., respectively) differs only by the contribution of sublimation. The negative latent heat fluxes simulated with M4 result in the largest contribution of all the models to ablation by sublimation. The T _s = 0 assumption results in an assumed 611 Pa surface vapour pressure, which is typically larger than the vapour pressure in the atmosphere (see Fig. 2). This increases simulated sublimation (Equation (3)), especially when air temperatures are low (e.g. at night when the surface is not melting). However, even the relatively large 0.05 m w.e. of sublimation simulated with M4 is only a few per cent of the total ablation. Models that incorporate the 0°C assumption, such as M4, generally only compute ablation when the surface energy balance is positive (e.g. Arnold and others, 1996), and often ignore sublimation because of its generally small contribution to the total ablation (e.g. Willis and others, 2002).

Values of L _out, Q _H and Q _L are similar between M2 and M3 for the full period and for the period of melting. The difference in total ablation predicted by these two models is therefore primarily a function of the modelled duration of melting, 40.3 and 50.9 days for M2 and M3, respectively. Values of L _out and Q _H are also similar between M1 and M2, but the subsurface heat flux, Q _G, comprises a significant heat sink for M1. In M2, Q _G acts as a short-term (i.e. daily to a few days) energy storage term and thus has a value of zero when averaged over longer timescales. In contrast, the glacier heat flux simulated with M1 is an estimate of the true heat exchange between the glacier surface and subsurface. The result is a total ablation of 1.68 m w.e. simulated by M1 versus 2.10 m w.e. simulated by M2. While the duration of surface melting predicted by M1 is only 26.4 days, the duration of surface and subsurface melting combined is 60.8 days.

Summary and Conclusions

We quantify the differences between various treatments of the glacier surface temperature and subsurface heat flux in a point-scale energy-balance model, aiming to identify an optimal balance of model simplicity and accuracy. The most physically defensible model, M1, is a MSM that includes dry densification of snow, penetration of shortwave radiation and subsurface melting, percolation and refreezing of meltwater in the snowpack and generation of slush layers. This model simulates the measured ice-temperature profile well, and provides a reasonably robust simulation of total ablation with different initial temperature profiles. The subsurface heat flux, Q _G, plays a non-negligible role in reducing the ablation simulated by M1 compared with the other models. While the other models significantly overestimate the total ablation measured at the AWS site in 2008, the ablation modelled with M1 was within 6% of the measured value. Given the advantages of this model, an effort should be made to collect at least the snow temperature and density data (i.e. vertical profiles) required to initialize the model. Such data can be readily collected with minimal additional effort as part of most traditional mass-balance studies.

Model M2 solves for the glacier heat flux, Q _G, as the residual in the energy-balance equation when the surface energy flux is negative, and computes the temperature of a 5 cm surface layer based on the value of the glacier heat flux. Because of the significance of Q _G as a heat sink in this study, and the limitation of M2 furnishing only short-term heat storage, cumulative total melt was overestimated by 18% with this model. Surface temperatures simulated by M2 were systematically lower than those simulated by M1. M3 employs an iterative temperature scheme whereby the surface temperature is adjusted until the energy fluxes are balanced. This model produces highly variable and significantly lower surface temperatures than M2 or M1, and ignores any delay in the onset of melting after a period of subfreezing temperatures. This results in a 27% overestimation of total ablation by the end of the melt season, primarily as a result of overestimating the duration of melt. M3 performs only marginally better than the simplest model, M4, which assumes a constant surface temperature of 0°C.

While M1 is the only model that arguably performs well in this study, M2 has an appeal in that it does not require any knowledge of the subsurface snow or ice properties. M2 produces results intermediate between the physically based subsurface model in M1 and the commonly made 0°C assumption in M4. In particular, ablation and surface temperature are simulated with greater accuracy by applying M2 than by applying M3 or M4. For studies that require treatment of sub-freezing surface temperatures and the glacier heat flux, but that aim to avoid the complexity of a full subsurface model, a shallow multilayer hybrid model may be a productive compromise (e.g. MacDougall and Flowers, in press).