Surface energy-balance model

We use the distributed energy-balance model developed by Reference Hock and HolmgrenHock and Holmgren (2005). The model solves the surface energy balance (Equation (1)) and is briefly described here. For more details, the reader is referred to Reference Hock and HolmgrenHock and Holmgren (2005).

S _in at each gridpoint is based on the observed S _in at a climate station. To obtain S _in at each gridpoint, S _in is divided into a direct _{S in dir} and a diffuse _{S in diff} part. This division is based on an empirical relation between the ratio of observed S _in to incoming shortwave radiation at the top of the atmosphere S _ToA, and the ratio of S _{in diff} to S _in. S _{in dir} at each grid-point is calculated using the ratio of the measured _{S in dir} to potential direct radiation at the climate station multiplied by potential direct radiation at the gridcell. The potential direct radiation is the amount of solar radiation that can reach the surface at a given gridpoint under clear-sky conditions and is approximated from the solar zenith angle, atmospheric transmissivity and surface slope and aspect. S _{in diff} is divided into a diffuse part from the sky and diffuse part from the terrain using a sky-view factor. The sky-view factor is based on the integrated elevation angle of the horizon for each gridpoint.

The parameterization of the surface albedo α is from Reference Oerlemans and KnapOerlemans and Knap (1998) and Reference Zuo and OerlemansZuo and Oerlemans (1996) and describes a decrease of α after a snowfall event that depends on the elapsed time since the event and the total thickness of the snow layer. It includes the effect of water on the surface (Reference Zuo and OerlemansZuo and Oerlemans, 1996) and the diurnal cycle in α due to variations in the solar zenith angle (Reference Segal, Garratt, Pielke and YeSegal and others, 1991). The ice albedo α _ice and the superimposed ice albedo α _sup in the parameterization are constant in space.

Longwave incoming radiation L _in is based on measurements at the climate station. L _in is divided into a sky and terrain irradiance based on the sky-view factor. The sky irradiance is assumed constant over the glacier, while the terrain irradiance is distributed using the sky-view factor. The long-wave outgoing radiation L _out is calculated from T ₀ using the Stefan–Boltzmann relationship with an emissivity of 1.

The turbulent fluxes of heat H and moisture LE are calculated using the Monin–Obukhov similarity theory. The fluxes are calculated between the surface and 2 m, the observation level at the climate station, using air-temperature, humidity and wind-speed data from the climate station. Surface humidity is based on saturation values at the given T ₀. T ₀ is calculated using the snow model described in the following subsection. The snow model also computes the subsurface energy flux Q _G and provides the amount of melt. The energy supplied by the rain Q _R is determined by the temperature difference between the air and the surface as well as the observed rainfall rate.

The snow model is an addition to the Reference Hock and HolmgrenHock and Holmgren (2005) model. In the original model, Q _G is set to 0 Wm⁻² or set a priori. T ₀ is set to 0°C in case of a positive energy balance when excluding melt (i.e. the lefthand side of Equation (1) is positive). Q _M then closes the balance with values larger than 0 W m⁻². The energy is converted to water equivalent melt which immediately contributes to mass loss in the computed mass budget. In case of a negative energy balance when excluding melt, T ₀ is lowered iteratively in small steps until a balance is reached in which Q _M is 0 W m⁻². In effect, H, LE, L _out and Q _R are changed to compensate a zero Q _G.

Multi-layer snow model

We use the SOMARS model developed by Reference Greuell and KonzelmannGreuell and Konzelmann (1994). The snow model calculates T ₀ and provides Q _G and Q _M, including the effects of meltwater percolation and refreezing. Reference Greuell and KonzelmannGreuell and Konzelmann (1994) used it in combination with energy-balance calculations to study the energy and mass balance at the Swiss ETH (Federal Institute of Technology) camp on Greenland. The model was tested extensively by Reference Bougamont, Bamber and GreuellBougamont and others (2005) at further locations in West Greenland.

The model calculates three variables on a grid extending vertically from the surface to an arbitrary depth (in this case, 40 m), namely temperature T _g, the density of the dry part of the snowpack ρ and water content w. The distance between snow model gridpoints increases exponentially from 0.04 m at the surface to a maximum of 5 m. On this grid, the thermodynamic energy equation becomes

where c_{p i} is the heat capacity of ice (2009 J kg⁻¹ K⁻¹) and δQ _t /δz is the absorption of energy arriving from the atmosphere. Q _t is the sum of the radiative and turbulent fluxes and only impacts on the uppermost layer. M and F are the melt and refreezing rate, respectively, and L _f is the latent heat of fusion (0.334 × 10⁶ J kg⁻¹). K is the effective conductivity and describes energy exchange through conduction, convection, radiation and vapour diffusion processes. K is a function of snow properties and is generally described as a function of density ρ (e.g. Reference Sturm, Holmgren, König and MorrisSturm and others, 1997). We use the equation of Reference Sturm, Holmgren, König and MorrisSturm and others (1997) (ρ in kg m⁻³):

Figure 1 presents an overview of the procedure followed in the snow model. First, a temperature profile is calculated without considering the effects of melt or refreezing. The water content and the effect of melting and refreezing on the temperature profile are then calculated. Melt occurs when the temperature T _g of a layer exceeds 0°C. T _g is set to 0°C and the remaining heat is used for melt. The meltwater is added to the water content of the layer. Rainwater and water on the surface are added to the water content of the top layer.

Fig. 1.

Flow chart of the multi-layer snow model. Steps 4a–f are calculated from the top down while 4g is calculated from the bottom up: T _g is subsurface (glacier) temperature; ρ is dry snow density; w is snow water content; and K is effective conductivity.

When T _g is below 0°C, refreezing occurs. The amount of refreezing is limited by either the temperature, which cannot be raised above melting point, the available amount of meltwater or the available pore space in the snow. After calculation of refreezing, the remaining water percolates downwards under the influence of gravity. A small amount is held in the layer against gravity by capillary and adhesive forces: the irreducible water content θ _mi. θ _mi is the ratio of mass of irreducible liquid water to the total mass of the layer, and can also be expressed as a fraction of the pore volume θ _pi. We describe θ _mi as a function of density according to an empirical relationship described by Reference Schneider and JanssonSchneider and Jansson (2004), based on measurements on Storglaciären and in a laboratory (Reference Coléou and LesaffreColéou and Lesaffre, 1998):

Here, n is the porosity of the snow layer, i.e. the ratio of pore volume to total volume, and θ_m is the gravimetric liquid water content defined as the ratio of the mass of liquid water m _liq to the total mass of the layer (m _liq + m _sn) where m _sn is the mass of the snow/firn. ρ _tot is the density including the water content, ρ _i is the ice density (900 kg m⁻³) ρ _w is the water density. θ _mi increases with increasing n or decreasing ρ. If the total water content w exceeds the maximum irreducible water content, w is set to this maximum and the rest percolates further downwards into the next layer.

Water may percolate through the successive vertical layers until it reaches impermeable ice. On top of the ice layer, water can accumulate and form a slush layer where all pore spaces are occupied by water. Water is removed from the slush layer as a function of a timescale of runoff t _runoff which in turn depends on the surface slope β (Reference Zuo and OerlemansZuo and Oerlemans, 1996):

The coefficients c ₁, c ₂ and c ₃ are based on the runoff time-scales for surface water on a steep slope τ _steep (β > 5°), on a horizontal surface τ _hor and on a slope of 1° τ _1°. The amount of water leaving the slush layer corresponds to the model ablation. Using this formulation, the fraction of water in the slush layer that is lost per time-step increases with decreasing t _runoff, which decreases with increasing β. The speed at which runoff occurs is assumed to be greater at the surface than within the snowpack, and t _runoff(surface) is therefore multiplied by a constant factor to obtain t _runoff(snow)_. The slush layer can reach the surface at any point of the glacier; on the surface a water layer can form and also refreeze. Table 1 presents values of the different timescales and factors as presented in the literature and used in this study. Note that no measurements are available to verify these timescales.

Table 1.

Timescales of runoff (days) used in Equation (8), as presented in the literature

The densification of the dry snowpack is described by an empirical relation developed by Reference Herron and LangwayHerron and Langway (1980), modified by Reference Li and ZwallyLi and Zwally (2004), and depends on snow temperature T _g and accumulation rate a (in mm w.e.):

where K ₀ is the rate factor and E the activation energy. Both are a function of T _g. R is the universal gas constant (8.3144 J K⁻¹ mol⁻¹) and J is the vapour flux. The factor b is approximately equal to 1. The accumulation rate a represents the change of overburden pressure, which is set to the initial snow layer thickness. Note that densification caused by internal melting or refreezing (through their impact on T _g) will dominate the densification of the total snowpack.

The surface temperature T ₀ is calculated by linear extrapolation of the temperature of the two uppermost grid layers. To obtain the most accurate T ₀, the two uppermost layers must be reasonably thin, of the order of a few centimetres.

Modelled mass balance and internal accumulation

Area-averaged summer balance is derived from the specific summer balances b _s from individual stakes. The stake observations consider density and mass changes down to the summer surface of the previous year (Reference Holmlund, Jansson and PetterssonHolmlund and others, 2005). However, if the cold wave of the winter penetrates below the summer surface of the previous year, water that is still present in the firn will refreeze and is therefore retained. In other words, internal accumulation occurs and, as a result, is overestimated.

Internal accumulation b _i is the sum of refreezing of meltwater percolating in cold firn in spring b _p and refreezing of meltwater held by capillary forces in winter b _c, i.e. b _i = b _p + b _c. b _i can be estimated from firn temperature and density profile observations. For a given temperature profile, the amount of mass that can be refrozen in the firn b _{p max} can be determined by (Reference Schneider and JanssonSchneider and Jansson, 2004):

where ΔT is the temperature increase possible at depth z. H _sf is the depth of the snow–firn interface and H ₀ is the maximum depth of the 0°C isotherm at the end of winter. For a given density profile, the amount of mass that can be refrozen in the firn, b _{c max}, can be determined by:

Both equations give a maximum estimate assuming that the complete cold wave is removed in spring by refreezing water and the firn is filled to its maximum with water held by capillary forces at the end of the summer.

The model we use is capable of distinguishing between the different terms contributing to the mass budget: snowfall, sublimation/condensation, melt and refreezing melt- and rainwater. The mass balance is then calculated as the mass added to the surface minus the water removed from the slush layer (excess melt- and rainwater) and thus includes internal accumulation b _i. The summer balance including b _i will be denoted by b _si.

Equivalent to the observations, which report the mass changes with respect to the previous year’s summer surface, the model computes . The difference between the methods is b _i and includes all water that is still present at the end of the calculations but will run off before refreezing occurs. In addition, the two components of b _i are explicitly calculated. b _p is determined from the actual mass refrozen in spring minus the amount of refrozen mass melted later in the season. b _c is determined from the water content held by capillary forces at the end of the summer. These do not necessarily equal b _{p max} and b _{c max}.

Model input

The model needs two types of input: (1) weather station data that drive the model, and (2) gridded fields that provide boundary and initial conditions. The gridded fields are a digital elevation model (DEM) from which slope, aspect of the slope and the sky-view factor are calculated. As initial conditions, the model needs grids of initial snow cover defined as snow water equivalent on top of either ice or firn, firn layer depth and profiles of snow temperature, density and water content.

The model is forced by hourly data of air temperature T _a, relative humidity RH, wind speed ws, short- and longwave incoming radiation S _in and L _in and precipitation for a location on the glacier. To distribute these over the glacier, T _a is changed using a constant lapse rate. ws and RH are assumed to be constant in space. Precipitation is divided into rain and snow using a linear transition from rain to snow when air temperature ranges between ±1°C of a threshold temperature T _t. Below T _t – 1°C, all precipitation is snow; above T _t + 1°C, 100% rain is assumed. Precipitation is increased by a constant percentage to account for the gauge undercatch error and is assumed to increase with a constant percentage with elevation.

Storglaciären

Our study area is Storglaciären, a small valley glacier in northern Sweden, close to the Swedish research station Tarfala (67°55′ N, 18°35′ E). The 30 × 30 m DEM is based on a map from 1990 (Fig. 2) (Reference HolmlundHolmlund, 1996). The area of the glacier is 3 km², and the elevation spans 1120–1730 m a.s.l. The glacier is polythermal, with a cold surface layer up to 50 m deep overlaying the temperate ice in the ablation area (Reference Pettersson, Jansson and BlatterPettersson and others, 2004). The study period is the 1999 melt season. Different types of observations were carried out on the glacier during this period (e.g. automatic weather station (AWS); snow water content; snow, firn and ice temperatures; ground-penetrating radar surveys; and mass balance using stakes at 53 sites on the glacier (e.g. Reference Jonsell, Hock and HolmgrenJonsell and others, 2003; Reference Schneider and JanssonSchneider and Jansson, 2004)).

Fig. 2.

DEM of Storglaciären. The thick line represents the outline of the glacier. The black dots indicate the stake locations, and the stars indicate locations mentioned in the text.

Meteorological data

The AWS was placed on the glacier at 1370 m a.s.l., reasonably close to the equilibrium line. The station was operational from 9 May to 2 September 1999 and collected all data needed to force the model. In addition, the AWS provided data of reflected shortwave radiation S _ref, outgoing longwave radiation L _out and surface height changes (obtained from a sonic altimeter). Time series of the observations are provided in Figure 3 and time averages in Table 2.

Fig. 3.

Hourly weather station observations of air temperature (a), relative humidity (b), wind speed (c), shortwave incoming radiation (d), albedo (e) and precipitation (f) over the period 9 May–2 September 1999.

Table 2.

AWS observations averaged over the period 9 May–2 September 1999. Except for α, the mean, minimum and maximum values are based on hourly data. For α, they are based on daily averages of S _in and S _ref

The average temperature over the period was 2.7°C which, according to Tarfala station data (1965–2002), was an average summer. This is also reflected in the for 1999 (−1.51 m w.e.) compared with the average during 1946–2001 (−1.69 ± 0.49 m w.e.). Although the average wind speed is relatively low (3.2 m s⁻¹), the near-surface flow over the glacier is dominated by a katabatic wind, indicated by an average wind direction of 267° (downslope of the glacier) and a directional constancy of 0.64. During the summer season, 819 mm w.e. precipitation was recorded, of which about 20% was snowfall. From the end of July the melt season was characterized by frequent snowfall, evident from the albedo record (Fig. 3).

Snow, firn and ice data

The initial snow cover for the model run (equivalent to the winter balance b _w) was determined from roughly 280 snow probings in a 100 m regular grid and density measurements in several snow pits carried out in April and May 1999 (Reference Holmlund, Jansson and PetterssonHolmlund and others, 2005). The firn depth grid was estimated from the ground-penetrating radar observations presented by Reference Schneider and JanssonSchneider and Jansson (2004) (Fig. 4). The firn limit was determined from aerial photographs taken in 1998/99.

Fig. 4.

Firn depth (m firn) based on radar measurements from May 1999 (Reference Schneider and JanssonSchneider and Jansson, 2004). Firn depth is 0 m in the grey area.

At location S29 and the AWS site, firn and snow temperatures were monitored hourly. S29 observations were made throughout the melt season; at the AWS site, observations were made until 22 June. At the AWS site, five sensors were placed at a fixed distance above the snow–ice interface. On 9 May, these depths corresponded to 0.05, 0.6, 1.15, 1.75 and 2.05 m below the snow surface. The results are presented in Figure 5. The snow surface reaches the melting point on 21 May at the AWS site and on 9 June at S29 (not shown). The snowpack starts to warm rapidly due to percolation and refreezing of meltwater. At the AWS site, the whole 2 m snow-pack is at the melting point on 14 June. At S29, firn temperatures were monitored in 1 m intervals down to the firn–ice transition (Reference Schneider and JanssonSchneider and Jansson, 2004), and the entire firn pack is temperate on 4 July (not shown).

Fig. 5.

Observed and modelled (HH05, GK94 and GK94-rd) (a) hourly surface temperature and (b–e) snow temperatures at the AWS site at four different depths below the surface. On 9 May 1999, these depths were (b) ∼0.05 m; (c) ∼0.6 m; (d) ∼1.15 m; and (e) ∼1.75 m. The initial snow thickness was 2.30 m snow, reducing to 2.10 m on 7 June, 1.55 m on 13 June and 1.25 m on 23 June. See Results for a description of the designations HH05, GK94 and GK94-rd.

The observed snow temperatures on 9 May at both locations were used to initialize the model subsurface temperatures T _g. A second-order polynomial was fit through the uppermost observations providing a temperature profile for the first 2.5 m. Below 2.5 m, a linearly increasing temperature was assumed, limited by the melting point. The temperature profiles were distributed over the gridpoints using a constant lapse rate of −0.01 Km⁻¹, and temperature was set to 0°C below 25 m depth at all gridpoints. Examples of the initial temperature and density profiles for three locations and an observed profile (S29) are given in Figure 6. The density profiles are based on observations in snow pits and a firn core in April and May 1999.

Fig. 6.

Initial model (a) temperature and (b) density profiles at three locations on the glacier, including the observed temperature profile on 10 May 1999 at S29 (Reference Schneider and JanssonSchneider and Jansson, 2004). Note that the density profile in snow is the same at all gridpoints.

Model sensitivity

Table 3 presents the results of the sensitivity tests of GK94, focusing on the sensitivity of the spatial and temporal variations of the mass balance. In general, the modelled mass balance is very sensitive to the parameterization of the surface albedo α, the chosen roughness length of momentum z _0m, heat z _0h and moisture z _0q, and to all parameters related to snowfall.

The sensitivity to the parameterization of α is not surprising since shortwave radiation provides up to 60% of the melt energy on Storglaciären. A change in α of 0.1 due to fresh snow (α _frsn) results in a 26% change in and ±32% change in b _s at the AWS site. A similar change in α of ice (α _ice) changes by 4% and b _s by ±8% at the AWS site. Both have the largest impact in the ablation area. The representation of the spatial variability of b _s, expressed in the residual standard deviation of b _s at the stakes, also varies considerably and is best when a relatively α _ice is chosen.

Changing the roughness lengths changes the exchange of sensible and latent heat at the surface considerably, resulting in changes in amount of melt. Typical values of z _0m for glacier surfaces are 0.1–5 mm (Reference Brock, Willis and SharpBrock and others, 2006). On Storglaciären over ice, z _0m values of 2.7 mm (Reference Hock and HolmgrenHock and Holmgren, 1996) and 0.1 mm (Reference Grainger and ListerGrainger and Lister, 1966) have been reported. Unfortunately, observations of z _0h and z _0q are scarce. However, Reference AndreasAndreas (2002) shows that the surface renewal model of Reference AndreasAndreas (1987), in which z _0h and z _0q are a function of z _0m and the friction velocity, gives reasonable results over glaciers. Varying z _0m within its uncertainty limits changes by more than 5% and b _s by an even greater amount. Using the model from Reference AndreasAndreas (1987), a change in z _0m also results in a change in z _0h and z _0q. Changing only z _0h and z _0q has similar effects to changing z _0m.

All parameters related to snowfall, such as the threshold temperature, the undercatch correction and the elevational gradient, directly and indirectly affect the mass balance: directly by increasing/decreasing accumulation, indirectly by changing α. Changing the threshold temperature has the largest impact: a change of 1°C changes by as much as 20%.

The modelled mass balance is not very sensitive to most of the parameters in the snow model, which is fortunate since many of these parameters are not well known. Within the range of uncertainty, is not very sensitive to the initial density profile, the densification of the dry snowpack, the conductivity of the snow/ice, or the values of the constants used in calculating the slush layer. The latter is due to the fact that the surface slope over most parts of the glacier is larger than 5° where t _runoff approaches τ_steep. is sensitive to the initial temperature profile and the irreducible water content, and is more sensitive than . Since minus is , is also most sensitive to the initial temperature profile and the irreducible water content. It is noticed that these are also the only parameters, besides the firn density, to determine b _{i max} using Equations (13) and (14).

Model performance

We tuned the model to best represent the spatial and temporal variations in b _s by comparing modelled b _s with stake observations and the surface lowering at the AWS site (Figs 7 and 8; Table 4). For an additional performance check, we used AWS observations of L _out and S _ref and observations of T _g at the AWS site and S29, which are not used as input for the model.

Fig. 7.

Modelled and observed hourly values of surface lowering with respect to the snow–ice interface on 9 May 1999 at the AWS site. Observations are from a sonic altimeter; modelled surface lowering results are from modelled changes in snow thickness and density and melted ice.

Fig. 8.

Modelled vs observed specific summer mass balance b _s at 53 stakes for the model runs GK94, GK94-rd and HH05. For GK94, the summer balance including internal accumulation b _si is also plotted.

Table 4.

Model performance averaged over the modelled period for experiments GK94, GK94-rd and HH05

Figure 7 shows that GK94 represents the temporal variations in the surface lowering at the AWS site reasonably well. The transition from snow to ice is within a day of the observations; the average difference is 0.03 m w.e., mainly caused by a small overestimation of the thickness of the snowpack in the first part of the run. b _s (0.84 m w.e.) is almost equal to the observation (0.81 m w.e.). GK94-rd and HH05 better represent the temporal variations when the surface is snow-covered. The surface lowers slightly faster, which is caused by immediate removal of water in the snowpack, while in GK94 meltwater percolation and refreezing delays the decrease. When the snow has melted, ice melt in GK94-rd and HH05 is severely overestimated.

The spatial variability in b _s is less well represented (Fig. 8). Comparing modelled b _s for all three runs with stake observations, the model overestimates melt in the higher parts of the firn area and underestimates melt in the lower parts of the ablation area. To illustrate the effect of b _i, b _si (GK94) is also plotted in Figure 8. The average b _s over the stake locations is reasonably well reproduced in GK94, but this is due to these compensating errors. Melt in GK94-rd and HH05 is on average overestimated, especially in the firn area. Since the same over- and underestimation pattern is found in all runs, it is probably not associated with the treatment of T ₀, meltwater percolation or refreezing. A possible explanation is the assumption of a constant wind speed over the whole glacier. Assuming an increase in wind speed from the upper parts to the tongue, wind speeds are lower in the accumulation area, resulting in smaller turbulent fluxes and less melt. On the lower parts, wind speeds increase, causing the turbulent fluxes and amount of melt to increase. However, there was no clear relation found between altitude and wind speed in observations carried out at different locations on the glacier (Reference Konya, Hock and KlingbjerKonya and Hock, 2004; Reference Hock and HolmgrenHock and Holmgren, 2005). Other possible explanations are errors in the initial snow-cover field, spatial variations in snowfall, snowdrift or roughness length.

The seasonal mean S _ref is reasonably well reproduced in each run. From the observed S _ref, α is derived. In GK94-rd and HH05, α is on average too high when the surface is snow-covered and too low when ice is exposed (Fig. 9). Since the albedo parameterization is the same in all runs, the underestimation of α can be partly explained by the fact that in GK94-rd and HH05 superimposed ice does not form: α _super = 0.5 and α _ice = 0.3. Furthermore, fresh snowfall increases α, but measurements of snowfall are difficult. Comparing the sonic altimeter to the precipitation record shows that observed snowfall is often underestimated. A constant factor is applied to correct this undercatch problem. However, problems remain in the cases where no snowfall is recorded but the sonic altimeter data show a surface height increase. Finally, in GK94-rd and HH05, water percolation and refreezing are neglected; in GK94, these processes delay the melting of fresh snow and therefore also delay the decrease in albedo.

Fig. 9.

Modelled and observed surface albedo based on daily averaged S _ref and S _in at the AWS site.

The discrepancies in the temporal and spatial variations of the mass balance b _s described above are partly explained by the problems in modelling α. An α which is too low impacts on b _s by increasing melt. The overestimation can be improved by increasing α _ice in the parameterization of α. However, although an α which is too low results in the overestimation of melt, the underlying problem is mainly the input precipitation data and the neglect of superimposed ice in GK94-rd and HH05. It is therefore not desirable to use α _ice to compensate for these problems. Note that when using observed α, which is independent of modelled snowfall and superimposed ice formation, modelled melt at the AWS site will improve.

Surface temperature T ₀ is derived from observed L _out. In all runs, modelled T ₀ is on average too high, except in the first 10–15 days when air temperatures are still below freezing point. In this period, the daily cycle in T ₀ is overestimated, with largest differences during the night in HH05. During cold spells, Q _G is positive; heat is transported from lower layers to the surface. When set to 0, as in the iterative method in HH05, this heat transport is compensated for by reducing the other surface fluxes by lowering T ₀, which increases the turbulent fluxes and decreases L _out. As a result, T ₀ is underestimated, especially in night-time stable conditions. Introducing the snow model and especially melt and refreezing decreases the variability and increases T ₀ during the night, thereby increasing the average T ₀ (Table 4; Fig. 5).

The cold bias in T ₀ in the first 10–15 days of the model is transported down into the snow. This is seen in the temperature profiles at the AWS site (Fig. 5) and at S29 (not shown). The model runs show a daily cycle which is too large and daily average snow temperatures T _g which are too low, similar to T ₀. The T _g are transported downwards and reach 1.15 m depth before the onset of the melt season. Since the bias originates from the surface, it is likely caused by inadequacies in the calculations of the surface energy fluxes.

Figure 5 also illustrates the importance of including percolation and refreezing of water for correct representation of T _g, as was shown in other studies (e.g. Reference Greuell and KonzelmannGreuell and Konzelmann, 1994). For the first ∼15 days, the snow is dry, and runs GK94 and GK94-rd show similar T _g. The percolation and refreezing of meltwater from 23 May onwards raises T _g to the melting point fairly quickly. The observations show a more gradual increase. Without percolation and refreezing, the snow below the upper ∼0.1 m does not reach the melting point, clearly emphasizing the dominant role of refreezing of percolating water in removing the cold content of the snowpack.

Mass balance

Figure 10 and Table 5 present the observed winter b _w, modelled summer b _s and net mass balance b _n for 1998/99. Since the mass-balance measurements do not include internal accumulation, we use the GK94 model results without internal accumulation b _i taken into account for comparison with specific and area-averaged mass balances and area-averaged mass balances derived from the measurements.

Fig. 10.

(a) Winter balance b _w 1998/99 based on observations carried out in April and May 1999; (b) modelled (GK94) summer mass balance b _s until 2 September 1999; and (c) net mass balance b _n.

Table 5.

Observed and modelled mean mass balance and mass balance at S29 (m w.e.)

b _w is based on observations and shows an increasing accumulation with elevation. On the lower part of the glacier, small spots with higher amounts of accumulation are observed which are likely due to convergence of snowdrift. The glacier-averaged is 1.34 m w.e.

The modelled b _s mainly reflects the altitudinal variations of the glacier. The spatial average for GK94 is −1.36 m w.e., which is smaller than that derived from stake observations (−1.42 m w.e.; see Table 5). This can be explained by the underestimation of melt in the lower parts of the ablation area dominating the overestimation of melt in the higher parts of the accumulation area. In HH05 and GK94-rd, is overestimated, which is mainly due to a general overestimation of melt over the whole glacier with the exception of the lowest parts. Note that in both the observed and modelled estimate of , percolation and refreezing of water (which may result in b _i) is not taken into account. It is assumed that all mass lost is removed from the glacier. Taking b _i into account, (GK94) is smaller than (−1.08 m w.e.). At the end of the model run, some water still present in the firn will drain before the cold wave refreezes it. This will add a small amount to at a later date, resulting in and an estimated of 0.25 m w.e.

Figure 10b also shows that b _s is not negative over the whole glacier. In the highest parts of the glacier, b _s is slightly positive. However, observations do not indicate the presence of an area with positive b _s for 1999.

Figure 10c presents b _n for 1999 as the sum of b _w (observed) and b _s (GK94). In 1999, the equilibrium line is at about 1500 m which is above the firn line shown in Figure 4. Furthermore, b _n is positive in several parts of the area where firn depth was defined to be 0 m which is mainly due to reasonably high amounts of accumulation in winter. Combining the observed (1.34 m w.e.) with the modeled (−1.36 m w.e.) results in a of −0.02 m w.e. From observations of b _s, is estimated to be −0.08 m w.e. In both estimates, b _i is not taken into account.

Internal accumulation

Table 5 presents the modelled (GK94) area-averaged internal accumulation and its compenents and ; Figure 11 presents their spatial variability. The figure shows that b _i increases with distance from the firn line, which is due to the increasing thickness of the firn layer. is +0.25 m w.e., constitutes a contribution of about 20% to and results in a positive balance for 1999 instead of the negative balance found without taking into account. and contribute 20% and 80% to , respectively.

Fig. 11.

Modelled internal accumulation (GK94): (a) total amount of water refrozen in spring b _p; (b) irreducible water content at the end of the model run b _c; and (c) b _i, the sum of (a) and (b). Internal accumulation is zero in the grey area.

When calculating b _p,max and b _c,max (Equations (13) and (14)) from the initial model temperature profiles and the final model density profiles, the resulting values for are similar to that presented in Table 5. Thus, the total winter cold wave is removed by refrozen water, and the maximum capacity of the firn to hold water against gravity is filled at the end of the melting season. At the end of the melting season, the water content is not at its maximum only at the higher parts of the glacier. This is also the area where b _s is positive.

Compared with the values presented by Reference Schneider and JanssonSchneider and Jansson (2004), the value for b _i presented here is much larger, mainly due to a much larger value of b _c. From observed density and temperature profiles, Schneider and Jansson estimate b _i to be 0.10–0.22 m w.e. at S29. Modelled values indicate values of ∼0.39 m w.e. at the same location, which is ∼18% of the accumulation at this site.

The sensitivity tests showed that is sensitive to the initial temperature profile, especially , and the irreducible water content θ _pi, especially . The reasonable agreement between modelled and observed b _p at S29 is explained by the fact that the initial temperature profile is based on observations at S29. However, θ _pi is also partly based on observations at S29 and therefore cannot explain the difference between model and observed b _c at S29.