1. Introduction
In recent years there has been increasing interest in predictive tools for spatially distributed estimates of melt rates (KirnbauerReference Kirnbauer, Bloschl and Gutknecht and others, 1994), further promoted by increased availability of remote-sensing data. Although, for many purposes, melt rates are obtained by simple empirical temperature index methods (Reference HockHock, 2003), more physically based energy-balance models generally better take into account the large spatial and temporal variability in melt rates typically encountered in a valley glacier environment. However, few studies have attempted to extend the surface energy balance at the point scale (for summary see Reference Ohmura, Kasser and FunkOhmura and others, 1992; Reference HockHock, 2005) to a larger spatial unit(e.g. on a catchment scale). Distributed models for snow-covered basins have been presented by, for example, Reference Bloschl, Kirnbauer and GutknechtBloschl and others (1991), Reference OhtaOhta (1994), Reference Ujihashi, Takase, Ishida and HibobeUjihashi and others (1994) and Reference Fierz, Pluss and MartinFierz and others (1997), while models developed for valley glaciers have been suggested by Reference Escher-VetterEscher-Vetter (1985), Reference FunkFunk (1985), Reference Arnold, Willis, Sharp, Richards and LawsonArnold and others (1996), Reference Hock and NoetzliHock and Noetzli (1997), Reference Brock, Willis and SharpBrock and others (2000b), Reference Escher-VetterEscher-Vetter (2000) and Reference Klok and OerlemansKlok and Oerlemans (2002). However, these show one or a combination of the following shortcomings: albedo is prescribed rather than simulated; longwave incoming and outgoing radiation are assumed constant in space; the direct and diffuse components of solar radiation are not calculated individually; or the surface temperature is assumed to be constantly at melting point. Hence, part of the spatial and temporal variability of energy-balance components is neglected.
This paper presents a fully distributed surface energy- balance model and its application to Storglaciären, Sweden, for two melt seasons. The model calculates the components of the surface energy balance for every hour and each gridcell, based on meteorological data collected at one point on the glacier. Direct and diffuse radiation are treated individually, considering topographic effects. Sky and terrain radiation are also considered separately for computation of both the diffuse and longwave irradiance. Albedo and surface temperatures are generated internally by the model. The bulk method is applied for calculating the turbulent fluxes. Results are validated by means of observations of discharge, melt rates and snow-line retreat.
2. Physical Setting and Data Input
Storglaciären is located in northern Sweden (67°55′ N, 18°35′ E), has an area of 3 km2 and ranges in elevation from ~1120 to ~1730 m (Fig. 1). The drainage basin of the glacier is 70% glacierized. The glacier holds the longest continuous and detailed mass-balance record in the world (Reference Holmlund and JanssonHolmlund and Jansson, 1999). The mass balance was positive in 1993 (1.00 m), and negative in 1994 (-0.37 m), although summer balances were relatively similar (-1.25 m in 1993 and -1.43 m in 1994).
Drainage basin of Storglaciären with 25m contour lines. N and S refer to the sites of the water-stage recording stations at the glacial streams Nordjåkk and Sydjåkk. The black dots mark the positions of the ablation stakes in 1994, only slightly differing from those in 1993. The circles labelled A-E denote the sites of five micrometeorological stations on the glacier.
During the 1993 and 1994 melt seasons, up to five automatic weather stations were operated on the glacier (Reference Hock and HolmgrenHock and Holmgren, 1996). The air-temperature, humidity, wind-speed, global radiation, reflected shortwave radiation and net all-wave radiation data collected at station B (Fig. 1), located at 1370ma.s.l. close to the equilibrium line, were used to force the model. However, the data of all meteorological stations on the glacier were analyzed to investigate their spatial variability, as air temperature, relative humidity, wind speed and precipitation have to be extrapolated to each gridcell. As an example, scattergrams of air-temperature, humidity and wind-speed data at stations A and B are given in Figure 2.
Air temperature (°C), relative humidity (%) and wind speed (m s-1) of station A (1250ma.s.l.) on the glacier tongue vs station B (1370ma.s.l.) close to the equilibrium line (Fig. 1), using hourly means of the 1994 melt season.
Air temperature correlated well between stations, with generally higher temperatures at lower elevations. An average lapse rate of 0.55 K(100 m)-1 was obtained from all the data available and applied in the model, although lapse rates on glaciers have been found to vary in time (Reference Braun and HockBraun and Hock, 2004). We investigated this simplification by running the model using a variable lapse rate as measured every hour during the sub-periods of simultaneous temperature measurements on the glacier. Model results differed only negligibly. We also explored the possibility of expressing a temporally variable lapse rate as a simple function of air temperature, but we did not detect any significant relationship. Thus, we adopt a constant lapse rate, allowing the model to be forced by data from only one weather station. Relative humidity varied relatively little among stations and was therefore assumed to be equal over the entire grid. As expected, wind speed displayed the largest scatter. Westerly winds dominated, coinciding with the down-glacier direction. During high wind-speed conditions, wind speed at the station in the accumulation area (station C) tended to be lower than at the meteorological stations at lower altitudes, probably due to the proximity of a precipitous rock wall and resulting shielding effects. For low wind speeds (<3ms—1), especially combined with high air temperatures, an increase in wind speed towards the southern margin (station E) and in particular towards the glacier tongue (station A) was observed (Fig. 2). Despite these apparent systematic variations and considering the large scatter among stations, the data available are considered too sparse to derive a general scheme to quantitatively assess the spatial distribution of wind speed. Thus, wind speed measured at station B was assumed representative for the entire area calculated.
Measured precipitation was increased by 35% to account for the gauge undercatch error and then extrapolated assuming a 10% increase with elevation for the catchment (Reference HieltalaHieltala, 1989). Besides the meteorological data, the model requires a digital elevation model (DEM) sufficiently large to allow for the computation of topographic shading of the area of interest and raster maps containing slope and aspect angles of each gridcell. The DEM was based on the 1 : 20 000 map by Reference Holmlund and SchyttHolmlund and Schytt (1987). In order to distinguish between firn and ice albedo after the retreat of the snow-line, the model also needs a raster map indicating for each glaciated gridcell whether the surface below the winter snow cover is firn or ice. The boundary is approximated from the previous years snow-line.
3. Surface Energy-Balance Model
For every hour and each gridcell the energy at the surface available for melt Q M is calculated from
where I is direct solar radiation, expressed in terms of the beam radiation received on the horizontal, D
s and D
t are diffuse sky and diffuse terrain radiation, respectively, α is the albedo, 
and 
are longwave sky and terrain irradiance, respectively, L ↑ is longwave outgoing radiation, Q
H is the sensible-heat flux, Q
E is the latent-heat flux, Q
G is the ground heat flux in the ice or snow and Q
R is the energy supplied by rain. Energy fluxes directed towards the surface are taken positive, those away from the surface negative. The effects of subsurface melting are not considered. The energy available for melt, Q
M (W m—2), is converted into meltwater equivalent M (mm h—1). The effects of Q
L on the mass balance and on the water available for runoff are included: M provides the input for the runoff modelling. M corrected for the mass transfer by sublimation or resublimation, henceforth referred to as ablation, will be compared to mass-balance data.
3.1. Global radiation
Highly sophisticated models have been proposed for spectral solar radiation in complex terrain (Reference DozierDozier, 1980), but their use is restricted to clear-sky conditions and they require various atmospheric parameters as input. We apply a simpler method using global radiation data at one point and applicable under all weather conditions. Measured global radiation is separated into the direct and diffuse components, which are then extrapolated individually to each gridcell considering terrain effects. The separation is based on the ratio of measured global radiation G to top-of-atmosphere radiation I ToA. With increasing cloudiness, G/I ToA tends to decrease, whilst the diffuse portion D/G increases. In the case of a completely overcast sky with low, dark clouds, the ratio G/I ToA is smallest and all global radiation is diffuse. Empirical relationships between D/G and the G/I ToA have been suggested by, for example, Reference Liu and JordanLiu and Jordan (1960) and Reference Collares-Pereira and RablCollares-Pereira and Rabl (1979). These may vary with space and time as a function of atmospheric conditions and site characteristics such as elevation, surrounding topography and its reflecting characteristics. Therefore, in the present study, a new empirical relationship was established using data collected by the meteorological station of the Swedish Meteorological and Hydrological Institute in Kiruna, 60 km east of Storglaciären (67°50′ N, 18°26′ E;408 m a.s.l.).
Hourly data of diffuse and global radiation for May- September 1987-96 were analyzed. The ratios D/G are plotted against G/I ToA in Figure 3a. Top-of-atmosphere radiation I ToA is calculated using standard methods (Reference SellersSellers, 1965):
where I 0 = 1368 Wm-2 is the solar constant (Reference Frohlich, McBean and HantelFrohlich, 1993), R and R m are the instantaneous and the mean sun- Earth distances, respectively, ϕ is the latitude, δ the solar declination and ω the solar hour angle.
Ratio of diffuse radiation to global radiation D/G vs the ratio of global radiation to top-of-atmosphere radiation G/I ToA at Kiruna: (a) hourly data; (b) daily means and functions to describe the relationship; (c) least-squares fits for May_September based on daily means.
The ratios considered are expected to lose accuracy with decreasing radiation amounts. Therefore, global radiation data or calculated top-of-atmosphere radiation less than 20Wm-2, or data when zenith angles exceeded 85°, were excluded from the dataset, leaving a total of about 22 000 data points. Figure 3a displays large scatter as the relationship is expected to vary with cloud type and zenith angle. Multiple regression techniques indicated that the effect of zenith angle was negligibly small. Therefore, the zenith angle was neglected, and daily means were applied in further analysis, thus also reducing the susceptibility of the involved ratios to errors, and hence the scatter in the correlation. A cubic polynomial function was fitted by the least-squares method to all daily data and to the data of each month individually (Fig. 3b and c), the latter in order to investigate the seasonal variation in the correlation. Whereas the functions for June-September hardly differ from each other, the function obtained for May yields ratios of D/G up to 0.1 higher for a given ratio of G/I ToA. The difference is attributed to the presence of a snow cover in May, resulting in enhanced diffuse radiation due to reflection from adjacent slopes and backscattering of reflected global radiation. The function derived from all daily data reads:
with D/G the ratio of total diffuse to global radiation and x the ratio of global to top-of-atmosphere radiation G/I ToA. The limits of 0.15 and 0.8 were arbitrarily chosen from the plot (Fig. 3b).
The plausibility of the least-squares fit was evaluated by comparing it to the curve obtained by averaging the data over defined classes of ratios G/I ToA (empirical function in Fig. 3b). For each class the ratios G/I ToA and the corresponding ratios D/G were averaged and plotted against each other. The class width was variable and adjusted to include 100 values for each class. The function thus obtained agreed well with the fitted function, indicating that the least-squares fit provides a reasonable description of the correlation (Fig. 3b). The least-squares fit is also in surprisingly good correspondence with the relationship derived from daily data of five different stations in the USA by Reference Collares-Pereira and RablCollares-Pereira and Rabl (1979).
For every hour between sunrise and sunset, Equation (3) together with measured global radiation is used to obtain the diffuse radiation at climate station B, which is then subtracted from global radiation to yield the direct solar radiation at this site, I s. Topographic shading is calculated for every hour and for each gridcell from the path of the sun and the effective horizon. If the climate station is shaded by surrounding topography, any measured global radiation is assumed diffuse.
3.1.1. Direct solar radiation
Following Reference OhtaOhta (1994) and Reference Fierz, Pluss and MartinFierz and others (1997) for every hour and each gridcell, direct radiation I is obtained from
where the subscript s refers to the location of the climate station and c denotes clear-sky conditions. I s is obtained from Equation (3), and I c = I sc at the location of station B. The ratio I s/I sc accounts for the deviation from clear-sky conditions, expressing the reduction of potential clear-sky direct solar radiation mainly due to clouds. It approaches unity under clear-sky conditions and tends to decrease with increasing cloudiness as direct radiation I s decreases. The ratio I s/I sc obtained at the climate station is assumed to be spatially constant. Potential direct solar radiation I c is calculated for each gridcell taking into account the combined effects of slope and aspect angles of each gridcell (Reference Garnier and OhmuraGarnier and Ohmura, 1968):
where ψa is the atmospheric clear-sky transmissivity, P is the atmospheric pressure, P 0 is the mean atmospheric pressure at sea level, Z is the zenith angle, β is the gridcell slope angle, φsun is the solar azimuth angle and φslope is the slope azimuth angle. An average value of 0.75 was derived for clear-sky transmissivity from the measurements of global radiation and an assumed fraction of diffuse radiation of global radiation of 15% (Reference Konzelmann and BraithwaiteKonzelmann and Ohmura, 1995). Variations in transmissivity ψa are not considered since values of I c are subsequently scaled by measurements, and hence the resulting error is considered negligible. The P/P 0 ratio accounts for the effect of altitude on direct radiation, and cos Z in the exponent expresses the variation of the path length with sun altitude. If the gridcell calculated is topographically shaded, direct radiation is set to zero. Potential direct radiation is computed at 10 min intervals, and hourly means are applied in Equation (4).
If the climate station is shaded by surrounding topography, Equation (4) is not applicable since the ratio I s/I sc is no longer defined, as potential direct solar radiation at the climate station, I sc, is zero. This case must be considered since the climate station is located on a valley glacier, and thus potentially subject to topographic shading, while other parts of the glacier may still be in the sunshine. For these gridcells the last ratio that could be obtained before the climate-station gridcell became shaded is applied. Thus, cloud conditions are assumed to remain constant until a new ratio I s/I sc can be determined at the climate station, usually the next morning. The error introduced by this assumption is considered small, as topographic shading of the climate station will only occur during times of low sun altitude, when direct radiation tends to be small.
3.1.2. Diffuse radiation
Diffuse radiation is assumed to be isotropic, although clear- sky diffuse radiation is significantly anisotropic (Reference KondratevKondratev, 1969). However, the error is considered negligible, because under clear-sky conditions the diffuse fraction of global radiation is small (10-20%) and under overcast conditions, when the diffuse component is dominant, the diffuse radiation is almost isotropic with the topographic screening by surrounding topography. Total diffuse radiation D is computed by means of a sky-view-factor relationship taking into account the fraction of the hemisphere that is obstructed by surrounding topography and by considering additional diffuse radiation reflected from adjacent slopes:
where the first term refers to sky radiation and the second term to terrain radiation. D 0 is diffuse radiation from an unobstructed sky, αm is the mean albedo of the surrounding terrain, G is global radiation and F is the sky-view factor defined by (Reference OkeOke, 1987)
where γ is the elevation angle of the horizon, which is integrated over the azimuth φ here in steps of 15°. Values range between zero, indicating a completely obscured sky, and unity, referring to an unobscured sky. The expression includes the variation of radiant intensity with the direction of radiation. Equation (7) assumes a horizontal surface, although strictly speaking the angle between the normal to the inclined surface and the direction of radiation should be used. However, the error due to this assumption is considered small. As diffuse radiation at the climate station is already affected by surrounding topography, the diffuse radiation of an unobscured sky, D 0, is calculated for the location of the climate station from Equation (6) and then assumed constant in space. The mean albedo αm will vary with time, mainly due to the progressive melt of the snow cover on the surrounding rock walls and on the glacier. Therefore, αm is obtained for every hour as the arithmetic mean of the modelled albedo of all gridcells of the entire drainage basin.
3.2. Albedo
Snow and ice albedo are treated individually. Among the large variety of albedo parameterizations (Reference Brock, Willis and SharpBrock and others, 2000a), a simple but widely used approach to account for the tendency of snow albedo to decrease with time is the formulation suggested by USACE (1956). Snow albedo is related to the number of days since the last significant snowfall. In the present study, this approach yielded large discrepancies between measured and modelled albedo, in particular during prolonged periods without snowfall (Fig. 4). In the present study, a new approach is suggested deriving hourly snow albedo from the albedo of the preceding hour. During snowfall, the preceding albedo is increased as a function of snowfall amounts. In the absence of snowfall, snow albedo is lowered as a function of the number of days after snowfall and air temperature. The former is taken as a surrogate to account for the effects of grain-size growth due to snow ripening and increasing impurity content with time. The lowering is assumed to be highest immediately after the snowfall event, decreasing exponentially with time. The inclusion of air temperature accounts for the acceleration of these processes as meltwater production is enhanced. Based on these considerations, the following empirical parameterization was derived:
where T is hourly air temperature (°C), P s is snow precipitation (mm h—1), n d is the number of days since snowfall and a 1-a 4 are coefficients obtained by calibration.
Measured and simulated daily snow albedo at stations B and C using the parameterization defined by Equations (8) and (9) and the formulation by USACE (1956): α = a 0 + a 1 exp (a 2 n d), where a 0 = 0:45, a 1 = 0:44, a 2 = -0:1 and n d is the number of days after snowfall.
To account for cloud effects, snow albedo obtained by Equation (8) is further modified using the ratio of measured global radiation to top-of-atmosphere radiation G/I ToA as an index for cloudiness. We incorporate the empirical relationship between the percentage change in albedo between two successive hourly measurements and the corresponding change in the ratio of G/I ToA as derived from an investigation of detailed albedo measurements on Storglaciären (Reference Jonsell, Hock and HolmgrenJonsell and others, 2003), and snow albedo is obtained by
where b = —20.7 and radiation is in W m-2. In case of snowfall, albedo is delimited to a maximum of 0.9. Snow albedo turned out to drop too rapidly, when hours of snowfall and no snowfall alternated, as the terms subtracted in Equation (8) are largest directly after snowfall events. Therefore, if calculated albedo drops below the value calculated for the time-step before the last snowfall, albedo and n d are set to their values immediately before the snowfall. Albedo will then continue to decrease at the same rate as if no snowfall had occurred. Hence, it is assumed that the new snow has melted.
A reduction in snow albedo when the snowpack becomes shallow is not included in the model. This common assumption was found to be unsubstantiated when tested with extensive field measurements on Haut Glacier dArolla, Switzerland (Reference Brock, Willis and SharpBrock and others, 2000a). In the firn area, when the winter snow has melted, the albedo is lowered once by 0.03 to account for the generally lower albedo of the previous years summer surface. Firn albedo is then treated as snow albedo. The value of 0.03 was obtained on average from a few instantaneous albedo measurements with portable equipment over snow and firn surfaces in the firn area. The coefficients a 1 – a 4 were fitted from the albedo measurements at station C in 1994 (a 1 = 0.005, a 2 = —1.05, a 3 = 0.0005, a 4 = 0.02 h mm—1). To avoid erroneous calibration due to possible effects of the zenith angle (Reference ArendtArendt, 1999) which are not accounted for in our parameterization, or due to measurement inaccuracies caused by surface slope effects (Reference Jonsell, Hock and HolmgrenJonsell and others, 2003), modelled hourly albedos were compiled into daily values from daily averaged global and shortwave reflected radiation prior to comparison with measured albedo. Measured and modelled daily snow albedo for station C in 1994 and station B in 1993 are compared in Figure 4. The model using Equations (8) and (9) performed significantly better than the routine by USACE (1956).
The albedos of ice at stations A and B varied little and did not show a significant decrease with time, as also found on other glaciers by Reference Cutler and MunroCutler and Munro (1996) and Reference Brock, Willis and SharpBrock and others (2000a). Any temporal variation in ice albedo is therefore neglected, and an albedo of 0.355 is assumed at the climate station, as obtained on average from the data available (Reference Jonsell, Hock and HolmgrenJonsell and others, 2003). To account for the tendency of debris to accumulate towards the glacier tongue, this value is extrapolated using an assumed increase of 3%(100m)-1 with increasing elevation, resulting in modelled ice albedo on Storglaciären of 0.33-0.38.
3.3. Longwave radiation
3.3.1. Longwave outgoing radiation
Longwave outgoing radiation L ↑ is approximated by
where ε is the emissivity of the surface (assumed to be 1), σ is the Stefan-Boltzmann constant and T is the absolute temperature (K) of the emitting surface. The surface temperature is assumed to be zero if the energy available for melt is positive. If melt turns negative, the surface temperature is lowered iteratively by steps of 0.25 K until the computed melt for the time-step and the gridcell considered becomes zero (Reference Escher-VetterEscher-Vetter, 1985; Reference Braithwaite, Konzelmann, Marty and OlesenBraithwaite and others, 1998). An iterative loop is needed because a lower surface temperature also affects outgoing longwave radiation, the turbulent heat fluxes and the heat transport supplied by rain.
3.3.2. Longwave incoming radiation
In complex topography, longwave sky radiation is reduced because part of the sky is obstructed, while additional radiation is received from surrounding slopes. Sky and terrain radiation are calculated individually. Longwave terrain irradiance is calculated using the parameterization suggested by Reference Pluss and OhmuraPluss and Ohmura (1997):
where F is the sky view factor (Equation (7)), T
a is the nearsurface air temperature (°C), T
s is the temperature of the emitting surface (°C) and a and b are constants (a = 0.77 W m—2 sr—1, b = 0.54 W m—2 sr—1). The constant L
b = 100.2 W m—2 sr—1 is the emitted radiance of a 0°C black body. The parameterization was derived in a snow- covered alpine environment using a narrowband radiation model (LOWTRAN) and includes the effects of the air between the emitting surface and the receiving surface. Longwave sky irradiance is obtained from
with L
0 the sky irradiance from an unobstructed sky. L
0 is obtained for the location of the climate station involving two steps. First, longwave incoming radiation is determined as a residual in the radiation balance at the climate station using the measurements of net radiation, global radiation and reflected shortwave radiation and the computed value of outgoing longwave radiation. Radiation obtained in this way is the sum of sky and terrain radiation 
Therefore, secondly, sky radiation from an unobscured sky, L
0, is calculated from Equations (11) and (12) for the location of the climate station. This value is assumed invariant in space and used to compute sky radiation of the entire grid. Total longwave irradiance at each grid element is the sum of and .
Deriving longwave incoming radiation from the radiation measurements merits some comment. Reference Halldin and LindrothHalldin and Lindroth (1992) showed that a large variation can be expected by comparing different types of net radiometers. In the present study, a Funk-type net radiometer (Reference KondratevKondratev, 1969) was used. Reference Müller and OhmuraMüller and Ohmura (1993) and Reference Konzelmann and BraithwaiteKonzelmann and Ohmura (1995) showed that this type of instrument is subject to a substantial underestimation of net radiation at an increasing rate with rising global radiation. They suggested a correction formula relating the error to global radiation. Net radiation corrected in this way for sites in Greenland and Switzerland agreed well with net radiation derived from measurements of the individual radiation components. Therefore, the net radiation data after correction are considered sufficiently accurate to allow for the procedure adopted here.
3.4. Turbulent heat fluxes
The turbulent heat fluxes are calculated from the bulk aerodynamic method, assuming the sensible- and latent- heat flux to be proportional to air temperature T z, wind speed u z and vapour pressure e z at height z above the surface:
where p is the air density, the subscript 0 refers to the mean atmospheric pressure at sea level, P 0, c p is the specific heat capacity of air, k is von Karmans constant (0.4), T 0 is the surface temperature, e 0 is the vapour pressure of the surface, z 0w , z 0T and z 0e are the roughness lengths for logarithmic profiles of wind speed, temperature and water vapour, respectively, ΨM, ΨH and ΨE are the stability functions, L is the Monin-Obukhov length, and L v is the latent heat of evaporation (2.514 × 106J kg-1) or sublimation (2.849 × 106J kg-1) as appropriate. Which one is applied is controlled by the direction of the latent-heat flux and the surface temperature. If the latent-heat flux is positive, the surface is expected to experience condensation in case of a melting surface and resublimation (a phase change from vapour to ice) in case of subfreezing surface temperatures. For negative latent-heat fluxes, sublimation is assumed, no matter what the surface temperatures.
Based on Reference Forrer and RotachForrer and Rotach (1997), for the stable case we apply the non-linear stability functions by Reference Beljaars and HoltslagBeljaars and Holtslag (1991):
We assume ΨE = ΨH. For the far less frequent unstable case we apply the commonly used Businger-Dyer expressions (see Reference PaulsonPaulson, 1970). The Monin-Obukhov length L is computed from
where u * is the friction velocity and T z is in Kelvin. Since the calculation of L requires a priori knowledge of Q H and u *, L is determined by iteration for the location of station B and for every time-step following the procedure outlined by Reference MunroMunro (1990). Stability functions determined in this way are then assumed spatially invariant across the glacier.
Detailed profile measurements of temperature, wind speed and water vapour were not available to determine the roughness length z 0w. Reported roughness lengths of ice and snow surfaces cover a range of orders of magnitude (for summary see Reference MooreMoore, 1983; Reference BraithwaiteBraithwaite, 1995) and are expected to vary in space and time (Reference HolmgrenHolmgren, 1971; Reference Greuell and KonzelmannGreuell and Konzelmann, 1994; Reference Pluss and MazzoniPluss and Mazzoni, 1994). A general pattern as to whether or not z0w values over ice are higher than over snow does not emerge in the literature. Hence, obtaining roughness lengths from the literature is problematic. On Storglaciären over ice, z 0w 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. Energy-balance studies have shown a large sensitivity of results to assumed roughness lengths (Reference MorrisMorris, 1982; Reference Harding and OerlemansHarding and others, 1989; Reference Hock and HolmgrenHock and Holmgren, 1996). Based on these considerations, in our model the roughness lengths are treated as model parameters, and hence chosen to yield optimal agreement between simulations and observations considering both melt and discharge. The roughness lengths for ice and snow were assumed equal and also constant in time during both melt seasons (z 0w = 10 mm). The roughness lengths of temperature and vapour pressure were assumed to be two orders of magnitude smaller than that of wind speed according to Reference HolmgrenHolmgren (1971) and Reference AmbachAmbach (1986). When interpreting relative contributions of energy- balance components, it must be borne in mind that the procedure of tuning roughness lengths entails that any errors concerning components other than the turbulent heat fluxes are lumped into the roughness lengths, thus affecting modelled energy partitioning. Nevertheless, considering the uncertainties with respect to the computation of the turbulent fluxes in the glacier environment (Reference HockHock, 2005), and specifically to the specification of roughness lengths and the as yet infant stage of available attempts at their temporal and spatial extrapolation, the method adopted here is considered to suffice in place of a more sophisticated scheme, but should be subject to future model improvement as more datasets become available and new generally applicable parameterizations evolve. Determination of the roughness lengths (or exchange coefficients) from closure of the seasonal energy budget has been adopted by various authors (e.g. Reference OerlemansOerlemans, 2000; Reference Konya, Matsumoto and NaruseKonya and others, 2004).
3.5. Ground heat flux and rain energy
A heat flux into the snow cover and the modelling of the formation of superimposed ice are neglected because snow temperature measurements at station B indicated that the winter cold content of the snow cover had been eliminated by the time the simulation began. On Storglaciären a heat transport into the ice is to be expected, because, although mostly temperate, a cold surface layer exists in the ablation area (Reference Pettersson, Jansson and HolmlundPettersson and others, 2003). Temperature-depth profiles in 1994 and during a micrometeorological experiment in the 1970s (Holmgren, unpublished information) indicated that the ice heat flux was approximately 0-5 W m—2 after the ice had been exposed (Reference Hock and HolmgrenHock and Holmgren, 1996). The amount is small but, because the direction of the flux is constant, neglect of the flux would introduce a systematic error. Therefore, for any ice-exposed gridcell, the ground heat flux is crudely approximated by assuming that it drops from 5 Witt 2 to 0Wm—2 between 1 July and 1 September, with linear interpolation in between. The model neglects the heat energy required to warm the surface up to melting point following periods of negative energy balance. However, the cooling process is not restricted to the surface but rather extends over a volume, whereas during the reverse process it is sufficient to raise the temperature of the surface to 0°C for melt to occur. Subsurface layers are then effectively heated by the latent heat released from the refreezing of percolating surface meltwater. Consequently, the processes of cooling and warming-up of the surface are asymmetrical, and only part of the integrated negative energy balance needs to be compensated for before melt will occur.
The energy supplied by rain Q R is calculated by
where p w is the density of water, c w is the specific heat of water (4.18 kJ kg—1 K—1), R is the rainfall rate and T r and T s are the temperatures of rain (assumed to be equal to air temperature) and the surface, respectively. T r becomes equal to T s before rain leaves the glacier surface.
3.6. Accumulation
Any snow precipitation is accumulated as a snow cover on the glacier. The aggregational state of each precipitation is determined according to a fixed air-temperature divider of 1.5°C (Reference Rohrer and SevrukRohrer, 1989). A mixture of rain and snow is assumed for a transition zone ranging from 1 K above and 1 K below the threshold temperature. Within this temperature range, the snow and rain percentages of total precipitation are obtained from linear interpolation. Redistribution of the original snowfall by wind transport or avalanches is not considered.
3.7. Initial snow cover
The initial snow water equivalent on the glacier was determined from the winter balance survey and ablation stake readings. The winter balance is determined as part of the annual mass-balance programme (Reference Holmlund and JanssonHolmlund and Jansson, 1999) by manual probing in a 100 × 100 m grid across the glacier and several density pits. The amount of melt between the survey and the starting day of simulation was determined by means of approximately 50 ablation stakes and was then subtracted from the winter balance. The snow water equivalent was interpolated to each glacier gridcell using ordinary kriging interpolation (Reference Hock and JensenHock and Jensen, 1999). Initial snow cover on the surrounding valley slopes is difficult to quantify because measurements are inhibited by the steepness of the terrain. In 1994, the surrounding slopes and the glacier forefield were mostly snow-free at the beginning of the simulation, so any snow cover outside the glacier was neglected. In 1993 the areal snow distribution in the forefield was estimated from snow surveys. On the valley slopes snow was distributed to the slope gridcells as a function of elevation, slope and curvature of the terrain following a procedure discussed in Reference Bloschl, Kirnbauer and GutknechtBloschl and others (1991).
4. Results and Discussion
The energy balance was calculated for each gridcell of the catchment for the periods 7 June-17 September 1993 and 5 July-6 September 1994. Simulated results were compared to cumulative ablation at the ablation stakes on the glacier and to observed patterns of snow-line retreat. In order to allow for additional verification by means of observed discharge records, the melt model was coupled to a discharge routine described in detail by Reference Hock and NoetzliHock and Noetzli (1997). The model routes areal melt and precipitation through the glacier by means of three parallel linear reservoirs corresponding to the different storage properties of firn, snow and ice using storage constants of 350, 30 and 16 hours, respectively.
4.1. Discharge
Discharge data were available for the periods 9 July-17 August 1993 and 11 July-6 September 1994. The discharge pattern in 1993 was dominated by rain-induced discharge peaks, whereas melt-induced discharge cycles prevailed during the 1994 melt season. Despite the different weather character, simulated and observed discharge hydrographs for the two years are in very good agreement (Figs 5 and 6). The model performs well with respect to both the diurnal and the seasonal discharge fluctuations. Expressing the agreement in terms of the R2 efficiency criterion by Reference Nash and SutcliffeNash and Sutcliffe (1970), values of 0.84 and 0.86 are obtained for the periods considered in 1993 and 1994, respectively.
Simulated melt averaged over the glacier, M (mm h-1), hourly data of air temperature T(°C) and precipitation P (mm h-1) at station B, and simulated and measured hourly discharge Q (m3 s-1) (black lines) of Storglaciären, 7 June_27 August 1993; and cumulated discharge Q cum (m3) (grey lines) over the period of discharge observations, 9 July_17 August.
Simulated melt averaged over the glacier, M (mm h-1), hourly data of air temperature T(°C) and precipitation P (mm h-1) at station B, and simulated and measured hourly discharge Q (m3 s-1) of Storglaciären, 5 July_6 September 1994; and cumulated discharge Q cum (m3) (grey lines) from 11 July, when the measured discharge record started.
4.2. Cumulative ablation
Simulated ablation was cumulated for the locations of the ablation stakes over the periods 7 June-17 September 1993 and 5 July-25 August 1994, and compared to measured values (Fig. 7). The results for 1994 show a tendency for melt to be overestimated by the model. In 1993 the model tends to slightly under-predict high ablation amounts and overestimates low ablation.
Measured vs simulated cumulative ablation (mw.e.) at the ablation stakes on Storglaciären for the periods 7 June_17 September 1993 and 5 July-25 August 1994.
4.3. Snow-line retreat
Accurate predictions of snow-line retreat are a decisive factor in melt modelling, because the snow-line separates two areas of appreciably different albedo. The observed and modelled patterns of snow-line retreat for the 1994 melt season are illustrated in Figure 8. Results indicate good agreement between observations and simulations, although in the first part of the melt season the model tends to underestimate the areas where the ice is exposed, whereas later in the season this tendency is reversed. The mismatch at the glacier margins later in the season may be associated with the uncertainties involved in determining the initial snow water equivalent. Snow probings showed a tendency for snow depths to increase abruptly at the margins of the glacier tongue. This increase, however, cannot be represented by the applied interpolation scheme, as data were not available in these areas, so the initial snow cover was probably underestimated.
Observed and simulated patterns of snow-line retreat on Storglaciären for the 1994 melt season. The black areas are those where the ice or firn is exposed.
4.4. Energy-balance components
The quality of the radiation model was evaluated by comparing modelled and measured global and net radiation at stations A and C. Results show fairly good agreement between measurements and simulations (Fig. 9). Discrepancies in global radiation are mainly associated with errors in the modelling of the timing of topographic shading. If the climate station is in sunlight (shade), but the model assumes shading (sunlight), global radiation will be underestimated (overestimated). The absence of systematic deviations in the net radiation plots promotes confidence in the accuracy of the new albedo parameterization, although albedo underestimation at one time might be compensated by overestimation at another time. On average, modelled net shortwave radiation fluxes agree well with those measured at the weather stations.
Measured vs simulated hourly global radiation and net radiation at stations A and C (Fig. 1) for the data in 1994.
On average, net radiation only contributed roughly 40% and 60% of the melt energy in 1993 and 1994, respectively, the remaining energy mainly being supplied by the sensible- heat flux (Table 1). As expected, direct radiation exhibited the largest spatial variability, reflecting the effects of topographic shading, slope and aspect (Fig. 10). The spatial variability of diffuse radiation was small, with highest values occurring on steep slopes. The differences between melt and ablation were small (Table 1), indicating that the mass changes by (re)sublimation were negligible compared to melt. On average, in 1993 the glacier experienced condensation, while in 1994 the latent-heat flux was negligible. In contrast to 1993, in 1994 evaporation prevailed in the higher parts of the glacier, whereas condensation was dominant at lower elevations. The large relative contribution of the latent-heat flux in 1993 appears high compared to the results obtained from most single location studies on valley glaciers in the Alps (e.g. Reference HockHock, 2005). However, in particular in more maritime areas, considerable condensation has been found on average (e.g. on Worthington Glacier, Alaska, USA (47Wm-2; Reference Streten and WendlerStreten and Wendler, 1967), Ivory Glacier, New Zealand (23 Witt 2; Reference Hay and FitzharrisHay and Fitzharris, 1988), and Peyto Glacier, Canada (15 Wm-2; Reference FohnFohn, 1973)). The difference in the relative contributions on Storglaciären between the two years is attributed to the difference in prevailing weather character. In contrast to 1994, the melt season in 1993 was subject to frequent rain and fog. In addition, the simulations extend over different time periods, and the relative importance of the energy- balance components may change considerably during the melt season (e.g. Reference OrvigOrvig, 1954; Reference HolmgrenHolmgren, 1971; Reference Wendler and WellerWendler and Weller, 1974).
Direct, diffuse, global and net radiation (Wm-2) and cumulative simulated melt (cm) of Storglacia_ren averaged over the period 7 June_17 September 1993.
Energy-balance components (W m-2) averaged over the glacier and the entire period of computation in 1993 and 1994
4.5. Parameter sensitivity
The sensitivity of results to the choice of various parameters and parameterization schemes was analyzed in order to assess the robustness of the model.
4.5.1. Turbulent heat fluxes
As the sensible-heat flux contributed a considerable amount to the energy available for melt on Storglaciären, the sensitivity of the results to assumptions about roughness lengths was investigated. Figure 11 shows simulated discharge for 1994 for a lowering of z 0w from the optimized value of 10 mm to 0.1 mm as reported for ice on Storglaciären by Reference Grainger and ListerGrainger and Lister (1966). The average sensible-heat flux is reduced from 36 to 15 Wm-2 and discharge is considerably underestimated, yielding R 2 = 0.5. Furthermore Figure 11 displays a model run neglecting a correction for atmospheric stability (Ψ = 0). Although, strictly speaking, a stability correction should be included since stable stratification is common over melting snow and ice, several studies have concluded that the turbulent fluxes tend to be underestimated when correction formulas in terms of the bulk Richardson number or the Monin-Obukhov length are applied over sloping glacier surfaces, in particular under strong stability (e.g. Reference Harding and OerlemansHarding and others, 1989; Reference Konzelmann and BraithwaiteKonzelmann and Braithwaite, 1995; Reference Van den BroekeVan den Broeke, 1996; Reference Forrer and RotachForrer and Rotach, 1997; Reference Martin and LejeuneMartin and Lejeune, 1998). Hence, various studies have neglected stability correction schemes in their energy-balance calculations (e.g. Reference Hock and HolmgrenHock and Holmgren, 1996; Reference OerlemansOerlemans, 2000; Reference Konya, Matsumoto and NaruseKonya and others, 2004). Our results based on neglect of stability correction yield an improved simulation of glacier runoff (R 2 = 0.92) mainly during the period of prolonged high air temperatures when runoff tended to be underestimated in the run including the stability correction. This finding may provide additional evidence that existing stability correction schemes underestimate the turbulent fluxes in the glacier environment. An assumption that the roughness length of ice exceeds that of snow did not lead to improved results, probably because the surface roughness of snow will rapidly increase as melt proceeds. Computing the roughness lengths z 0T and z 0e using the surface renewal model of Reference AndreasAndreas (1987), as, for example, done by Reference MunroMunro (1990) and Reference Van den BroekeVan den Broeke (1996), considerably underestimated discharge, especially during the warm, calm period in mid-July 1994. Generally speaking, results indicate that in regions where the turbulent fluxes play a significant role in the energy balance, the roughness lengths are sensitive parameters, emphasizing the importance of accurate estimates of the roughness lengths for melt modelling.
Hourly measured and simulated discharge Q (m3 s-1) of Storglaciären, 11 July-6 September, 1994, for two model sensitivity runs: reducing z 0w from 10mm to 0.1 mm; and neglecting a correction for atmospheric stability and applying optimized z 0w = 8mm and z 0T = z 0e = 0:008 mm.
4.5.2. Ice albedo
In various studies, albedo has been recognized as a key factor for the energy balance of a glacier (e.g. Reference Van de Wal, Oerlemans and van der HageVan de Wal and others, 1992). In our model, temporal variations in the albedo of ice-exposed gridcells are neglected, and spatial variations are only considered by an elevation dependency. Two sensitivity runs were performed, increasing and decreasing ice albedo by 0.05, thus assigning it 0.305 and 0.405. The results indicate little impact within the range considered, yielding almost identical hydrographs. Concerning the spatially averaged shortwave radiation balance, the overall effect is dampened, as larger parts of the glacier are snow-covered and not affected by variations in ice albedo. The sensitivity to albedo variations is also dampened by the large relative contribution of the turbulent fluxes to melt energy.
4.5.3. Separation of direct and diffuse radiation
The function splitting global radiation into the direct and diffuse components obtained from all available data (Equation (3)) differed significantly from that derived only for the May data (Fig. 3c). Therefore, the latter function was applied in a sensitivity run for comparison. Discharge simulations were almost identical, yielding the R 2 value as before. Spatially averaged melt energy remained the same, although the fractions of direct and diffuse radiation were altered by roughly 10%.
5. Concluding Remarks
A distributed surface energy-balance model has been developed, combining existing and new approaches to the spatial modelling of ice- and snowmelt in rugged terrain on an hourly time-scale. The model is forced by near-surface data of air temperature, humidity, wind speed and global radiation collected atone meteorological station considered representative of the area calculated. In addition, data on reflected shortwave radiation and net radiation were used to derive longwave incoming radiation, which could alternatively be obtained from measurements or from a parameterization based on cloud data. The model was applied to Storglaciären to simulate hourly melt for each grid element of a DEM for two melt seasons. Modelled melt rates varied substantially across the glacier, mainly due to the effects of surrounding topography on direct solar radiation. The contribution of the turbulent fluxes to the energy available for melt was found to be remarkably high, with roughly 60% and 40% in 1993 and 1994, respectively, resembling the energy partitioning generally found in maritime climates.
In general, simulations of melt rates, patterns of snow-line retreat and proglacial discharge corresponded well with observations. In particular, the model performed remarkably well with respect to the simulation of the melt-induced diurnal discharge fluctuations typical of glacial regimes. Such diurnal variability is generally not captured by models based on simpler temperature-index methods (Reference HockHock, 2003). The energy-balance model may be used to provide the meltwater input for runoff modelling in any snow- or ice- covered terrain, particularly when the daily discharge cycles are to be simulated, or for any other purposes requiring the prediction of melt rates with high spatial or temporal resolution.
Sensitivity studies indicated the robustness of the model. However, results were sensitive to the choice of roughness lengths within the range of values reported for ice and snow surfaces in the literature, emphasizing the importance of accurate estimates of the roughness parameters in case the turbulent heat fluxes are substantial contributors to melt. Roughness lengths are difficult to obtain and their temporal and spatial variability poses an additional problem, thus rendering the roughness lengths a tuning parameter. Further research should be directed towards the parameterization of roughness lengths over snow and ice surfaces and their spatial and temporal variability, as well as to more realistic extrapolation of wind speed across glaciers. In addition, the new snow albedo scheme needs further scrutiny using other datasets.
The model has successfully been tested in other climate settings (e.g. in Antarctica (Reference Braun and HockBraun and Hock, 2004) and on a tropical glacier (Reference SicartSicart, 2002)). Our model intends to provide a tool for sensitivity tests and to quantify the response of glacier melt and discharge to predicted climate changes. Although temperature-index models have been used for such purposes (e.g. Reference Braithwaite and ZhangBraithwaite and Zhang, 2000) and are more practical due to lower data requirements, the stability of calibration factors under a different climate remains uncertain, so more physically based methods are preferable if the necessary data are available.
Acknowledgements
The work was financed by the Swiss National Foundation for Scientific Research (grant No. 2100-037839.93) and the Swedish Research Council. Further financial support was given by the Kommission fur Reisestipendien der Schwei- zerischen Akademie der Naturwissenschaften and an exchange grant of the European Ice Sheet Modelling Initiative (EISMINT) within the European Science Foundation. H. Blatter, T. Konzelmann and P. Calanca gave constructive comments on the manuscript. K. Schroff prepared the climate instrumentation. C. Schneider established and provided the DEM. Special thanks to C. Tijm- Reijmer for pointing out a crucial programming error. An anonymous reviewer and the scientific editor M. van den Broeke are acknowledged for valuable remarks and suggestions.
