## Introduction

In the last decade, distributed mass-balance models of varying complexity have been increasingly used on glaciers to study the impact of climate change on glacier mass balance (Reference Arnold, Willis, Sharp, Richards and LawsonArnold and others, 1996; Reference Escher-VetterEscher-Vetter, 2000; Reference Klok and OerlemansKlok and Oerlemans, 2002; Reference Bougamont, Bamber and GreuellBougamont and others, 2005; Reference Hock and HolmgrenHock and Holmgren, 2005). One of the main reasons to start using distributed models was to account for the spatial variability in the shortwave incoming radiation and surface albedo. Some of these models use the energy-balance method to calculate the amount of melt since, compared to empirical temperature index models, they are more capable of describing the large spatial and temporal variability in melt rates encountered over (small) glaciers.

The models using the energy-balance method solve the equation:

where *S*
_{in} is the shortwave incoming radiation, *α* is the surface albedo, *L*
_{in} and *L*
_{out} are the incoming and outgoing longwave radiation, *H* and LE are the sensible and latent heat fluxes, respectively, *Q*
_{R} is the energy supplied by rain, *Q*
_{G} is the subsurface energy flux and *Q*
_{M} is the energy used for melt. Fluxes towards the surface are defined positive. The models calculate the different components of Equation (1) using parameterizations and observations at a climate station on or close to the glacier. More recently, output from global climate models has also been used as input in order to study longer time periods (e.g. Reference Hock, Radic and de WoulHock and others, 2007).

An important parameter necessary to solve Equation (1) is the surface temperature *T*
_{0}, which is used to calculate *L*
_{out}, *H*, LE and *Q*
_{G}. Observations of *T*
_{0} are generally not available, and therefore *T*
_{0} is calculated internally in the models. The simplest solution is to assume *T*
_{0} = 0°C and *Q*
_{G} = 0 W m^{−2}. This is a reasonable assumption on temperate glaciers during the melt season. However, when calculating multi-year energy and mass balances or applying these models on glaciers in (cold) Arctic regions, these assumptions are incorrect. Under these conditions, *T*
_{0} can have values below 0°C, *Q*
_{G} is non-zero and to determine *Q*
_{G} and *T*
_{0} it becomes necessary to calculate the englacial temperature distribution and the subsurface energy flux.

A non-zero englacial temperature distribution affects the energy and mass balance in several ways. Before melting starts, energy is consumed to remove the winter’s cold wave and raise the englacial temperature to the melting point. Thus, this energy is no longer available for melting. Part of the energy used to increase the englacial temperature results from the refreezing of percolating melt- and rainwater. Since this water is retained in the snowpack, runoff is reduced or delayed. The introduction of a snow model can therefore significantly affect the temporal and spatial variability in calculated energy and mass balance (Reference Janssens and HuybrechtsJanssens and Huybrechts, 2000).

Several snow models are available and have been tested on glaciers, sea ice and seasonal snow-covered surfaces (e.g. *Crocus* (Reference Brun, Martin, Simon, Gendre and ColéouBrun and others, 1989), SNOWPACK (Reference Bartelt and LehningBartelt and Lehning, 2002; Reference Obleitner and LehningObleitner and Lehning, 2004), SNTHERM (Reference JordanJordan, 1991; Reference Andreas, Jordan and MakshtasAndreas and others, 2004), Snow-SVAT (Reference Tribbeck, Gurney, Morris and PearsonTribbeck and others, 2004) and SOMARS (Simulation Of glacier surface Mass balance And Related Sub-surface processes; Reference Greuell and KonzelmannGreuell and Konzelmann, 1994)). These models all use similar parameterizations of the subsurface processes and simulate snow properties such as density, temperature and liquid water content. They have been shown to perform reasonably well in comparison to observations. To date, these models have mostly been applied at point locations and have rarely been coupled to distributed models.

Reference Lefebre, Gallée, Ypersele and GreuellLefebre and others (2003, Reference Lefebre2005) and Reference Bougamont, Bamber and GreuellBougamont and others (2005) were among the first to couple a subsurface snow model to a distributed model. Reference Lefebre, Gallée, Ypersele and GreuellLefebre and others (2003, Reference Lefebre2005) coupled a snow model based on *Crocus* to a regional climate model (MAR) and validated the results for individual locations on the Greenland ice sheet. Reference Fettweis, Gallée, Lefebre and van YperseleFettweis and others (2005) then used this model to study the snow water content over the Greenland ice sheet in comparison with satellite observations. Reference Bougamont, Bamber and GreuellBougamont and others (2005) also focus on the snow water content of the Greenland ice sheet using SOMARS.

Here we use the distributed surface energy- and mass-balance model developed by Reference Hock and HolmgrenHock and Holmgren (2005) in combination with the snow model SOMARS developed by Reference Greuell and KonzelmannGreuell and Konzelmann (1994), to study the contribution of internal accumulation to the mass balance of Storglaciären, a small valley glacier in northern Sweden.

Internal accumulation consists of two components: refreezing of meltwater percolating in cold firn in spring, and refreezing of water held by capillary forces when the cold wave penetrates into the firn in winter. Traditional mass-balance measurements only consider density and mass changes down to the summer surface of the previous year. Hence, summer mass balances may be overestimated since any meltwater refreezing when the cold wave penetrates below the summer surface is counted as ablation. Internal accumulation may constitute a significant term in the glacier mass balance, but estimates are scarce.

Reference Trabant and MayoTrabant and Mayo (1985) estimated that internal accumulation was 7–64% of the net accumulation on Alaskan glaciers. Reference Schneider and JanssonSchneider and Jansson (2004) estimated the contribution of internal accumulation to the mass balance of Storglaciären to be 3–5% of the annual accumulation in an average year. However, their results are based on observations of snow water content, temperature and density at a single location. Using a distributed model will enable us to make a better estimate of the area-averaged contribution of internal accumulation to the mass balance and also study the spatial variability.

We first introduce both models and our study area. We then present the results of the model sensitivity study and model performance. Finally, we present the mass balance of Storglaciären including the contribution of internal accumulation.

## Model

### 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*
_{0} 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*
_{0}. *T*
_{0} 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^{−2} or set a priori. *T*
_{0} 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^{−2}. 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*
_{0} is lowered iteratively in small steps until a balance is reached in which *Q*
_{M} is 0 W m^{−2}. 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*
_{0} 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

^{−1}K

^{−1}) 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

^{6}J kg

^{−1}).

*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

^{−3}):

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.

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

^{−3})

*ρ*

_{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*
_{1}, *c*
_{2} and *c*
_{3} 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.

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*
_{0} 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^{−1} mol^{−1}) 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*
_{0} is calculated by linear extrapolation of the temperature of the two uppermost grid layers. To obtain the most accurate *T*
_{0}, 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*
_{0} 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.

## Field Observations and Model Initialization

### 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^{2}, 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)).

### 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.

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^{−1}), 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.

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).

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^{−1}, 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.

## Results

We present results of three model runs, GK94, GK94-rd and HH05, including a sensitivity and performance analysis of GK94. GK94 uses the Reference Hock and HolmgrenHock and Holmgren (2005) model including the snow model of Reference Greuell and KonzelmannGreuell and Konzelmann (1994). The model parameters are tuned within their uncertainty limits such that the observed spatial and temporal variations of the mass balance *b*
_{s} are best represented. Table 3 presents the values of key parameters in the model and marks the parameters mainly used for tuning. The GK94-rd run is identical to GK94 except that percolation and refreezing of meltwater is excluded. HH05 uses the original iterative scheme to calculate *T*
_{0} and *b*
_{s} (Reference Hock and HolmgrenHock and Holmgren, 2005). Runs GK94-rd and HH05 are not tuned to obtain the best spatial and temporal *b*
_{s}, but the best-fit parameter set from GK94 was assumed, since we aim to illustrate the effect on the results of including meltwater percolation and refreezing.

### 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.

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*
_{0}, 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.

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*
_{0} is derived from observed *L*
_{out}. In all runs, modelled *T*
_{0} 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*
_{0} 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*
_{0}, which increases the turbulent fluxes and decreases *L*
_{out}. As a result, *T*
_{0} is underestimated, especially in night-time stable conditions. Introducing the snow model and especially melt and refreezing decreases the variability and increases *T*
_{0} during the night, thereby increasing the average *T*
_{0} (Table 4; Fig. 5).

The cold bias in *T*
_{0} 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*
_{0}. 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.

*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.

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.

## Summary and Concluding Remarks

To investigate the mass balance including internal accumulation on Storglaciären, we coupled a snow model (Reference Greuell and KonzelmannGreuell and Konzelmann, 1994) to a distributed energy- and mass-balance model (Reference Hock and HolmgrenHock and Holmgren, 2005). The snow model describes the temperature, density and water-content evolution of the snow/ice pack and includes the processes of percolation and refreezing of melt- and rainwater. Including percolation and refreezing of water is necessary to correctly represent the subsurface temperature in the model and makes it possible to take internal accumulation into account in the estimation of the mass balance of the glacier.

Sensitivity tests show that the model results are mostly sensitive to the albedo parameterization, the chosen momentum and scalar roughness lengths and all parameters related to snowfall. After tuning, the temporal variations in mass balance at the AWS site are well reproduced. The average mass balance based on stake observations is also reasonably well reproduced. However, this is the result of compensating errors. Melt is underestimated in the lower ablation area, while it is overestimated in the upper accumulation area. This is not an artefact of the snow model; the same is seen in the model run without the snow model. At the moment we do not have an explanation.

The area-averaged modelled internal accumulation is +0.25 m w.e., which amounts to about 20% of (+1.34 m w.e.) and would result in a positive net balance for 1999 of +0.23 m w.e. This is much larger than the amount of internal accumulation estimated by Reference Schneider and JanssonSchneider and Jansson (2004), which was 3–5%. The two components of internal accumulation, refreezing of percolating meltwater in spring and refreezing of capillary water in winter, amount to about 20% and 80% of the internal accumulation, respectively. Compared with estimates from observations, the modelled amount of internal accumulation seems high, which is mainly due to a very high value of refreezing of capillary water in winter. The modelled amount of internal accumulation is sensitive tothe initial snowtemperature profile and the irreducible water content, which are based on observations.

In conclusion, the addition of the snow model provides additional information on the mass budget of the glacier since it makes it possible to estimate the contribution of internal accumulation. However, results are relatively sensitive to input parameters and the model does not include all processes contributing to the specific mass balance on the glacier. For example, snowdrift and the horizontal movement of water from one gridcell to another have not yet been taken into account.

## Acknowledgements

We are grateful to the Department of Physical Geography and Quaternary Geology of Stockholm University for their hospitality. T. Schneider is gratefully acknowledged for providing the firn temperature data and W. Greuell for valuable information regarding SOMARS. Two anonymous reviewers are acknowledged for their useful comments. R. Hock is a Royal Swedish Academy of Science Research Fellow supported by a grant from the Knut Wallenberg Foundation. This work was financed by a Spinoza award from the Netherlands Organization for Scientific Research (NWO) to J. Oerlemans.