The authors propose a model of glacial mass balance based on correlations with meteorological observations and data from climate re-analysis. The minimum input data required include the following: average monthly temperature on the glacier and in its vicinity during summertime for a reference time period, average monthly air temperature, and average precipitation total at the nearest weather station or from re-analysis. This model was used to hindcast the mass balance and its components at Werenskioldbreen (southern Svalbard) over the period 1912–2005. The hindcast specific mass balance was then used to estimate the change in the thickness of the snout of Werenskioldbreen over the period 1958–1990. These results were compared with results obtained using a cartographic method. Comparing the topographic maps, the glacier front lowered 28.7 m on average over 32 years. The average difference in the calculation of the change in glacier thickness between these two methods amounted to 3.7 m (based on meteorological data) and 0.2 m (using ERA-40). The discrepancy of less than 13% confirmed that the method is a reasonably accurate way of predicting past glacier mass balance. The proposed method can find a broad application in hindcasting the mass balances of small Svalbard glaciers where observation data are scarce or nonexistent.
Introduction
The study of changes in glacier geometry, which are an important control on mass balance and movement, constitutes one of the most important research topics in ice-covered areas. The research provides valuable information about changes in glacier volume, trends in climate, and glacier contributions to sea level rise. In addition to helping scientists understand and model geomorphological and hydrological processes, glacier research has many practical applications. Understanding glacier mass balance is crucial to monitoring and forecasting flood risks and to the hydropower potential of glacierized basins.
Accurate mass balance estimates require high-resolution measurements of changes in ice surface elevation and glacier motion through melt seasons. Because these measurements are expensive and time-consuming to make, mass balance estimates are only available for a limited number of glaciers, and the period of record for most of these glaciers is quite short. As a result, models that can estimate the past, present, and future mass balance of unmonitored glaciers are needed.
This study tests a simple approach to modeling and hindcasting the mass balance of Werenskioldbreen, a small valley glacier in southern Svalbard (Fig. 1). This glacier belongs to the largest typological group of glaciers found in Svalbard, which are particularly numerous in Nordenskiöld Land and Andree Land. The project objective was to obtain reliable mass balance results for a glacier with very limited glacial and meteorological data using correlations and spatial relations with data recorded at Svalbard's weather stations, or ERA-40 re-analysis data. Popular models of glacier mass balance normally require a large amount of input data and when this is not available the applicability of such models tends to be limited. A simple mass balance modeling method with minimal, cheaply collected input-data requirements and accuracy of the results comparable to the ones obtained through direct mass balance measurements could be widely applicable to Svalbard glaciers of similar type and size.
Glacier Mass Balance Estimation Methods
There is a variety of ways to determine the mass balance of a glacier. Which method is used is typically determined by the data available and the desired accuracy of the result. Kaser et al. (2003) listed the following as the most popular methods for determining the mass balance of a glacier: (1) the glaciological method, (2) the hydrological method, (3) the cartographic method (geodetic, topographical, photogrammetric), (4) indirect methods based on relationships between the glacier mass balance and the altitude of the glacier equilibrium line (ELA), or its accumulation area ratio (AAR), (5) the flux method, and (6) the flux divergence method.
Mass balance models that use weather data simulate ablation and accumulation processes. Accumulation is mostly computed from precipitation using a threshold temperature to distinguish between rainfall and snowfall (e.g. Hock et al., 2007; Schuler et al., 2007). There are two popular methods of ablation process modeling:
(1) Energy balance can be used to estimate ablation (Hock, 1998; Oerlemans, 1993; Hock and Holmgren, 1996, 2005; Hock and Noetzli, 1997). Energy balance models calculate the amount of snow or ice that can be melted for specified time intervals based on measurements of the amount of energy supplied to and emitted from the glacier surface. Even the simplest energy balance models require a large amount of data inputs, including air temperature, relative humidity, net radiation, cloudiness, and precipitation.
(2) Temperature index methods (Braithwaite and Olesen, 1989; Laumann and Reeh, 1993; Johannesson et al., 1995, Jóhannesson, 1997; Oerlemans et al., 1998; Hock, 1998, 2003) use simple relationships between air temperature and the rate of ablation of ice/snow at various time intervals, which usually range from hourly to seasonal. An advantage of using this method is that it requires minimum data collection for model inputs, as the model only requires air temperature records to calculate the amount of snow or ice melt. The temperature index, however, is firmly based on physical relationships with long-wave atmospheric radiation, and the sensible heat flux, the primary controls on ablation (Ohmura, 2001; Hock, 2003). This method was used in the study to hindcast Werenskioldbreen mass balance.
Study Area and Data Source
Werenskioldbreen is located in Wedel Jarlsberg Land, north of the Hornsund (Fig. 1). It is a land-terminating polythermal valley glacier that has a well-defined basin boundary. Its accumulation field consists of four basins: the northern basin contains Skilryggbreen and Slyngfjellbreen tongues, the central basin consists of Werenskioldbreen proper, and the southern basin has the Angellkroken tongue. The tongues are separated by medial moraines, the most distinctive of which is located between Skilryggbreen and the main stream. The glacier runs longitudinally with its front turning to the north. In 1990, the glacier surface was ca. 28 km2 (Jania et al., 2002). The elevation of the glacier ranges between 30 and 780 m a.s.l., whereas in 1990, 90% of the surface area was below 500 m a.s.l.. Average glacier velocities are a few centimeters a day (Kosiba, 1960; Baranowski, 1977). A folding of the medial moraine may suggest a glacial surge in the past.
MASS BALANCE DATA
For the mass balance study it was necessary to collect appropriate meteorological data and direct information about the elements determining the mass balance. The existing glacier mass balance measurements are incomplete. The earliest available data on ablation and accumulation were recorded in the 1957/1958 and 1958/1959 seasons (Kosiba, 1960). The measurements were carried out at four sites in a vertical profile from the altitude of 65 m a.s.l. to 493 m a.s.l. (two stakes in ablation zone and remaining two in accumulation area). Accumulation measurements are underestimated because they were conducted after the winter accumulation maximum and Kosiba (1960) suggested treating them as minimum estimates. Possibilities to evaluate the mass balance based on these data are rather limited. The next time mass balance measurements were made on Werenskioldbreen was during the 1979/1980 season (Haeberli, 1985; Pulina, 1986). Net mass balance in this season was determined basing on sounding of snow cover thickness at the end of the accumulation season and hydrological observations in summer. During the 1993/1994 season, mass balance was calculated from measurements of changing snow and ice thickness at 8 stakes located on a vertical profile along the centerline of the glacier (Jania, 1994). None of the stakes were located on the northern tributary glacier of Skilryggbreen. The observations were carried out at the end of the ablation and accumulation seasons. Additional mass balance monitoring using ablation stakes was carried out during the 1999–2002 seasons. These observations were conducted using 11 stakes systematically located on the main stream of the Werenskioldbreen and its tributary glaciers (Fig. 1). Because several stakes were installed at similar altitudes, the resolution of the mass balance gradient in elevation is relatively poor. In spite of this, significant differences can still be observed between stakes located on a centerline of the main stream of Werenskioldbreen and Skilryggbreen. Smaller accumulation, larger ablation, and lower net balance were recorded at the tributary glacier compared to the main stream. It was assumed that measurements of snow and ice thickness at the stakes had an error of ±10 mm and the error for estimates of snow/firn/ice density was about 14% (±0.05 g cm-3).
METEOROLOGICAL DATA
The direct air temperature at Werenskioldbreen forefield are very poor. The following data were used: 10-day and monthly average temperature measurements for the 1970–1974 ablation seasons (Pereyma, 1988) and 1979 (Pulina et al., 1984), 5-day and monthly average temperature measurements for the 1983 (Pereyma and Piasecki, 1988) and 1985 (Brázdil et al., 1988) ablation seasons, and mean daily temperature data for the 1986 summer season (unpublished data from the University of Wrocław and University of Silesia expedition). Matching records from the Hornsund station were also used to complement these data (Pereyma, 1988; Pulina et. al., 1984). From 1982 onwards, the average daily meteorological data for Hornsund come from the meteorological data set of the Institute of Geophysics of the Polish Academy of Sciences. Air temperature data, including both monthly average and daily records, also came from some other Svalbard weather stations: the Isfjord Radio, Svalbard Lufthavn, and Ny-Ålesund (the Norwegian Meteorological Institute database).
For the mass balance model for wintertime (October–May), the study used precipitation totals recorded at the Isfjord Radio (1912–1975), Svalbard Lufthavn (1976–1978 and 1982), and Hornsund (1979–2005) stations (the Norwegian Meteorological Institute and the Institute of Geophysics of the Polish Academy of Sciences [ http://www.glacio-topoclim.org] databases).
Homogeneity sequences of both precipitation and temperature applied in this work were tested by Førland et al. (1997) for data from the following stations: Svalbard Lufthavn, Ny-Ålesund, Isfiord Radio, and by Marsz and Styszyńska (2007) for data from Hornsund station.
An alternative application of data from the European Centre for Medium-Range Weather Forecasts re-analysis (ERA-40) for a gridcell 2.5° × 2.5° including the glacier was proposed. The data from the ERA-40 re-analysis were successfully applied to mass balance models of Storglaciären in northern Sweden (Radić and Hock, 2006; Hock et al., 2007), as well as Midre Lovènbreen on NW Svalbard (Rye et al., 2010). The data of the ERA-40 re-analysis consist of a three-dimensional database which was produced from a digital climatic model, meteorological observations, and satellite data (Uppala et al., 2005). The database includes the period from mid-1957 to mid-2002. From this database, we retrieved daily air temperature from 2 m above the surface (June–September) to model summer mass balance, and daily precipitation totals in winter seasons (October–May) to model winter balance. These data were also used to calculate monthly temperature and monthly precipitation totals. The ERA-40 temperature and precipitation are related to mean altitude of a gridcell, i.e. 143 m a.s.l. (based on GLOBE, 1999). In order to adjust the re-analysis data to local conditions, we scaled temperature based on a lapse rate of 0.8 K/100 m. The lapse rate factor comprises temperature variations depending on elevation, horizontal distance, and model bias (Hock et al., 2007).
The local temperature gradient was calculated from temperature data provided by Pereyma (1983) and Pereyma and Piasecki (1988) from stations located at 20 m a.s.l. (direct forefield of Werenskioldbreen) and on the glacier surface (380 m a.s.l. in August 1979, 250 m a.s.l. in August–September 1983). The observed temperature lapse rate was in the range of 0.67–1.26 K/100 m (10-day mean) with an average of 0.85 K/100 m. A relatively high summer temperature gradient is a result of locating the meteorological sites on different surfaces (moraine, rock, ice). According to the temperature data from Werenskioldbreen area (Pereyma, 1983; Pereyma and Piasecki, 1988; Nasiolkowski and Pereyma, 2007), the temperature shift between the places located at the same altitude on unglaciated and icy surfaces is c. -1 K. This constant shift factor was adopted to calculate the temperature data over the glacier. To extrapolate the temperature across the glacier the lapse rate of 0.66 K/100 m, based on the temperature data from unglaciated Bratteggdalen valley near the glacier (Nasiółkowski and Pereyma, 2007), was used. The precipitation was scaled according to the applied precipitation gradient (60% increase of precipitation recorded at the sea level at every 100 m of altitude; Grabiec, 2005). The applied coefficients show changes of meteorological parameters according to altitude, distance, and model biases.
Figure 2 compares mean monthly summer temperatures and winter precipitation totals in the period of 1959–2002, downscaled to sea level (typical for the Werenskioldbreen forefield), calculated from the meteorological data and ERA-40 re-analysis. These temperature records are highly correlated (r = 0.93), which is statistically significant at the level of α = 0.05. Based on the data from the meteorological stations, the mean monthly summer temperature was 3.6 °C (standard deviation SD = 1.7) comparing well with the average of 3.5 °C (SD = 1.5) calculated from ERA-40 re-analysis data. Precipitation records are less well-correlated (r = 0.41), also statistically significant at the level of α = 0.05. Average winter precipitation was respectively 231.2 mm w.e.; SD = 75.4 (data from meteorological stations) and 227.7 mm w.e; SD = 61.1 (ERA-40).
Methods
MASS BALANCE MODEL
The net mass balance at a point of known altitude can be derived from the following formula:
where bn (meters of water equivalent [m w.e.]) is net mass balance, bw (m w.e.) is winter balance, and bs (m w.e.) is summer balance.WINTER BALANCE
In this study, we derived the winter balance using the accumulation algorithm (Grabiec, 2005):
where bw(l,hi) (m w.e.) is winter balance on a glacier located at a distance l (km) from open water at a point situated at an elevation of hi (m a.s.l.); Bh0 (m w.e.) is winter accumulation at a reference station at an elevation h0 (m a.s.l.); Lc is the location coefficient; (hi - h0) is the parameter calculated as the difference between the elevation of a point on the glacier for which the accumulation is determined (hi) and the altitude of the reference meteorological station (h0); and τ the accumulation gradient factor—the rate of accumulation increase per 100 m of elevation expressed as fraction of total winter accumulation on sea level.The winter accumulation at the selected reference station (Bh0) can be calculated from total winter precipitation (Grabiec, 2005)
where Ph0 (m w.e.) is winter (October to May) precipitation at a reference station at elevation h0; and K is the accumulation efficiency correction coefficient.The accumulation efficiency correction coefficient K derives accumulation rates from the sum of winter precipitation. The K factor is determined mainly by precipitation measurements error, which increases K, and the influence of liquid precipitation, which decreases K (Grabiec, 2005).
If only winter precipitation totals are known, but there is no way to distinguish between rain and snow, K is assumed to be equal to 1.1. This value of the K coefficient is justified when the solid precipitation consists of c. 60% of the total precipitation recorded (Førland and Hanssen-Bauer, 2000; Grabiec, 2005).
The Lc coefficient explains the statistical relationship between the winter precipitation totals at a glacier ELA located at d (km) from open water and precipitation at a meteorological station located at w (km) from open water (Grabiec, 2005). These relation-ships were determined by comparing the long-term average winter precipitation records from weather stations in Svalbard with their distance from the open sea.
For Werenskioldbreen, d was taken to be 6 km (approximate distance between the glacier ELA and the open sea).
The amount of sea-level accumulation achieved by Equation (3) is then converted into accumulation at a measurement altitude by applying the accumulation lapse rate factor τ. Following Grabiec (2005), τ = 0.6 for the purpose of this study.
The study uses winter precipitation totals from the Isfjord Radio (1912–1975), Svalbard Lufthavn (1976–1978 and 1982) and Hornsund (1979–2005) stations. The 1931 and 1941–1947 records from Isfjord Radio were incomplete and the gaps were filled by using the Werenskioldbreen (bww) winter balance formula based on the statistical relationship with altitude (A) in the observation seasons 1994 and 1999–2002:
Winter mass balance estimation may be much simpler using ERA-40 precipitation data. Re-analysis data represent average conditions of the gridcell. In consequence, sum of precipitation refers to mean grid elevation (143 m a.s.l.) and the data require scaling to local elevation using the precipitation gradient. As indicated above, applied precipitation gradient factor equals accumulation gradient factor (i.e. τP = 0.6). The gridcell including the glacier consists of ∼50% sea surface and ∼50% land. In such a case, the location coefficient (Lc) may be neglected. Taking into consideration these issues, winter balance may by expressed as: where bw(hi) (m w.e.) is the winter balance on a glacier at a point located at an altitude of hi (m); PERA-40 h0 (m w.e.) is ERA-40 precipitation at average grid elevation h0 (m); and τP is precipitation gradient factor—rate of precipitation increase per 100 m of elevation expressed as fraction of total winter precipitation at sea level.SUMMER BALANCE
The summer balance (bs) was determined from its relationship to the sum of positive degree-day (ΣPDD) in the summer season:
where ΣPDD (°C) is the sum of positive degree-day on particular elevation range in the summer months between June and September; and a (mm w.e./°C) is an ablation coefficient per PDD (degree-day factor).The set of ΣPDD for Werenskioldbreen was derived for the 1912–2005 seasons using monthly average temperatures that were calculated from correlations of air temperatures recorded in Werenskioldbreen forefield with the data from the weather stations at Hornsund, Isfjord Radio, and Svalbard Lufthavn (Fig. 3). ERA-40 temperature data from the gridcell including Werenskioldbreen for the period 1959–2002 were also used as an alternative method of calculating summer balance.
When average monthly temperature was >2.5 °C (Fig. 4, Part A), the monthly sum of PDD (ΣPDDm>25°c) was derived using the following formula:
where tm (°C) is the average monthly temperature, and d is the number of days in a month.If average monthly air temperatures are in the range of 2.5 to - 3.5 °C, this procedure may generate errors linked to the possible occurrence of negative daily temperatures. As a result, when temperatures fall in this range, we employ the following empirically derived formula (Fig. 4, Part B):
This formula was obtained by comparing average monthly air temperatures to the sum of PDD computed using average daily temperatures. Spitsbergen weather station data from the 1978–2005 seasons were used (Hornsund, Svalbard Lufthavn, Ny-Ålesund).For months with average air temperature below -3.5 °C the PDD sum equals 0.
The degree-day factor a calculates the amount of snow and/or ice that is melted per PDD. The a factor varies in time and space and is linked to the physical parameters of ice and snow. Example values of the a factor were summarized by Hock (1998, 2003). In this model we used different melt factors for snow and ice: 6 mm w.e./1PDD and 8.3 mm w.e./1PDD, respectively. The applied melt factors were calculated on the basis of the temperature data and ablation rate recorded by sonic rangers in the 2008 summer season on neighboring Hansbreeen (180 m a.s.l.). Those melt factors are consistent with the degree day factor at the same glacier obtained by Szafraniec (2002), which ranged from 2 to 8 mm w.e.
The refreezing of meltwater was reconstructed by a simple approach. The adopted Pmax coefficient (Reeh, 1989) determines the rate of winter accumulation that is retained during the ablation season. When snowmelt exceeds the local refreezing potential, melted water starts to contribute to ablation. Such a method of refreezing calculations was successfully employed in the Austfonna (NE Svalbard) mass balance model by Schuler et al. (2007). As the reference mass balance data from Werenskioldbreen did not include the internal accumulation, in the model development we tested two summer balance outputs: with neglecting refreezing Pmax = 0, and assuming value Pmax = 0.6 (Reeh, 1989; Schuler et al., 2007).
Specific net balance of Werenskioldbreen was calculated based on temporal variability of surface area and area-elevation distribution between 1936 (1:100,000; NPI), 1957–1959 (1:5000; Lipert, 1961), and 1990 (1:25,000; Jania et al., 2002). During 1936–1990, glacier surface area diminished by 4.3 km2 (13%). In the calculation of the glacier balance, a constant rate of change of the relation of areas in individual altitude zones from the state in 1936 to the state in 1990 was assumed. For the earlier and later periods the values were linearly extrapolated.
CHANGES OF GLACIER GEOMETRY AND VELOCITY FIELD
We verified net mass balance estimates from the 1990 front reach to the Angellfjellet-Wernerknatten cross section (Fig. 1) by comparing the glacier thickness change between 1958 and 1990 as determined from: (a) comparison of relevant digital elevation models (DEM); and (b) net balance estimations. The following maps were used: Werenskiold Glacier, lower and middle sections, scale 1:5000, drawn in 1957–1959, based on field images and with contours spaced at 2.5 m (Lipert, 1961); and Werenskioldbreen, Orthophotomap 1:25,000, based on aerial photos from 1990 with contours spaced at 25 m (Jania et al., 2002). DEMs based on these maps had a 100 m grid and were calibrated to UTM WGS 84. Based on DEM differencing, average thickness change of the glacier terminus between 1958 and 1990 was 28.7 m, or 0.9 m per year.
The change of thickness during a net balance year at any given point of the glacier (h) is expressed by the continuity equation (Cogley et al., 2011):
where (m of ice equivalent/year) is the surface mass balance rate (rate of surface accumulation + surface ablation); (m of ice equivalent/year) is the flux divergence.The negative flux divergence in the frontal part of the glacier is called emergence velocity EV (vertical vector of movement pointing upwards in the ablation zone) (Cogley et al., 2011). The average emergence velocity (EVm) of the frontal part of the glacier can be obtained from the analysis of ice volume which annually flows across the Angellfjellet-Wernerknatten cross section (EVm):
where Qi (m3/year) is the annual volume of ice flowing across the Angellfjellet-Wernerknatten cross section; and S (m2) is the glacier surface downstream from the cross section.The annual volume of ice flowing across the section (Qi) is:
where Vi (m/year) is the average annual ice flow velocity at the cross section; and Sp (m2) is the cross-section area.The ice flow velocity across the section roughly equals the average surface flow velocity (Paterson, 1994). Glacier flow velocity was assumed to be constant because there was no evidence of a glacial surge during the study period. The average ice flow velocity at the Angellfjellet-Wernerknatten cross section was determined to be 7.3 m/year based on a study carried out during 1956–1959 (Kosiba, 1960).
The area of the Angellfjellet-Wernerknatten cross section (Sp) was determined using the results of a ground-penetrating radar survey performed along that section in 2008, as well as the glacier surface profiles from 1958 and 1990 (Fig. 5). During that period, the glacier surface level was lowered by 6.2 m on average. Average lowering included a maximum reduction of 18.6 m in the middle part of the cross section, which was partially offset by a maximum increase of 10.1 m near the Angellfjellet slopes. As a consequence, there was little difference in glacier cross-sectional area between 1958 and 1990 (0.33 and 0.31 km2, respectively).
The average emergence velocity of the lower section of Werenskioldbreen (EVm) ranged between 0.41 m/year in 1958 and 0.59 m/year in 1990.
Calculating the vertical component at every point of the glacier requires detailed data on glacier movement and slope inclination. Possible glacier movement vectors in the ablation zone were discussed by Arnold (1981) and Jania (1988). An emergence velocity EV (m/year) may be calculated based on the following formula (Jania, 1988):
where Vh (m/year) is the horizontal velocity; α is the inclination angle of the glacier slope; and β is the angle between a horizontal plane and vector of movement (positive when the vector of movement is directed above and negative when it is directed below the horizontal plane). DEMs from 1958 and 1990 were used to generate the value of inclination angle of the glacier front surface.In order to obtain horizontal velocity in individual grid cells, we used glacier velocity along Angellfjellet-Wernerknatten transversal profile (Fig. 1) from the period 1956–1959 (Kosiba, 1960), 1970–1971 (Baranowski, 1977), 1981–1982 (Jania, 1988), and the velocity maps of the glacier front movement in the summer season of 1970 (Baranowski, 1975) and winter season of 1982 (Jania, 1988). The comparison of data from these periods suggests there were no significant changes in velocity across the Angellfjellet-Wernerknatten profile. As a result, we assume the glacier velocity was stable during 1957–1990. The model of glacier velocity at the terminus in the years 1958 and 1990 is shown in Figure 6.
Vertical angle (β) of the vector and the movement velocity of the frontal part of Werenskioldbreen was calculated by Jania (1988). The results of measurements indicate that in the lower part of the glacier the movement vector directed above the horizontal surface predominates, however, other cases should not be excluded. The mean value of β angle for the frontal part of the glacier was calculated from Equation (13) based on horizontal velocity, slope inclination, and average emergence velocity (EVm). The values of β were determined to range from 0.6° for 1958 to 0.8° for 1990. The estimated values of EV for 1958 and 1990 are shown in Figure 7. Intermediate values were adopted for other years. Maximum EV values occur in the central and northeastern part of the glacier front (near Angellfjellet-Wernerknatten profile line), reaching 0.9 m/year in 1958 and 1.3 m/year in 1990. The differences of values and changes in spatial distribution of vertical component of the movement between 1958 and 1990 are the effect of changes in glacier geometry which, in turn, results in changes of front inclination and changes of movement velocity. The glacier front below the Angellfjellet-Wernerknatten cross section increased its average inclination in the analyzed period from 4.8° to 5.5°. During the same time period, the glacier surface area decreased by almost 2 km2 (34% of 1990 surface area of the front), which corresponded to a decrease in volume of ice supply of 6%.
Results
MODELED VERSUS OBSERVED MASS BALANCE DATA
Figure 8 shows observed and estimated net balance data in selected measurement seasons. As the field records of mass balance neglected refreezing of melted water, they were compared with the model assuming refreezing potential, Pmax = 0. The results obtained from these methods are not substantially different from each other. Using the input data from the meteorological observations, standard deviation (SD) of the net balance modeled values from all the individual values observed on the stakes was estimated to be 0.34 m w.e., which compares to 0.32 m w.e. calculated from ERA-40 data. In the case of winter balance, the values of SD are 0.15 and 0.16 m w.e., respectively, and for summer balance, 0.28 and 0.30 m w.e. The model best fits the net balance in the season 1993/1994. SD of the net balance values from the 1993/1994 data is 0.07 m w.e., according to meteorological data, and 0.18 m w.e. based on the ERA-40 re-analysis data. The average correlation coefficient of the modeled and observed mass balance in the 5 seasons (1994, 1999–2002) was 0.67 (meteorological data) and 0.70 (ERA-40). The average correlation coefficients of winter balance are considerably smaller (r = 0.36 for the series based on meteorological data, and r = 0.33 for the ERA-40 data series) than the correlation coefficient of summer balance (r = 0.72 and r = 0.66, respectively).
Based on the observation data of mass balance from the years 1994 and 1999–2002, a specific net mass balance for the whole glacier was estimated together with its winter (bw) and summer components (bs). These results, together with the data from the season 1979/1980 (Haeberli, 1985), were used to compare a specific net mass balance with the modeled results (Fig. 9). The absolute error of the estimated specific mass balance is listed in Table 1. The results obtained from the meteorological data show lower values of the average absolute error than the ones based on ERA-40 re-analysis.
LONG-TERM MASS BALANCE RECONSTRUCTION
Based on the results of the model, the net mass balance for the whole glacier in the period 1912–2005 was prepared. Both the net mass balance and its components (summer and winter balance) were determined. Summer and net balance were calculated both with respect to refreezing potential (Pmax = 0.6) and neglecting this component (Pmax = 0). All the calculations of a mass balance were referred to a balance year in “stratigraphic system” of mass balance study (Østrem and Brugman, 1991). The glacier mass balance determined on the model is constant at the same altitude. Changes of the Werenskioldbreen's area and altitude from 1912 to 2005 were taken into account when estimating net balance of the whole glacier.
TABLE 1
Absolute error of estimated specific mass balance for selected seasons, m w.e. = meters of water equivalent, bs = summer balance (m w.e.), bw = winter balance (m. w.e.), bn = net mass balance (m w.e.).
Two mass balance data were calculated based on two different sets of input data. For the period 1912–2005, mass balance was estimated based on meteorological data from selected Spitsbergen stations. The second series of mass balance calculations included the period 1959–2002 and was based on the data from the ERA-40 re-analysis. We then statistically compared both techniques. Summer balance estimates correlate similarly as winter balance ones (r = ∼0.7). Correlation coefficients of winter and summer balances are statistically significant at the level α = 0.05. The results of net balance determined on different input data are correlated with each other at the level r = ∼0.6 and they are not statistically significant. The statistics showing the comparison of the results of mass balance modeling in the years 1959–2002 are shown in Table 2.
TABLE 2
Statistics of the results of glacier mass balance modeling in the years 1959–2002 basing on the input data from meteorological observations and ERA-40 re-analysis.
The results of Werenskioldbreen mass balance assessment (specific net balance, winter balance, summer balance neglecting refreezing) over the period 1912–2005 based on data recorded on meteorological stations and ERA-40 re-analysis are presented in Figure 10. A negative mass balance predominated over the period examined. The average specific net balance of Werenskioldbreen over the study period was -0.43 m w.e. (based on meteorological data). The glacier mass budget was only positive during a small number of seasons, mostly in the 1960s. This was mainly due to considerable positive values of the winter balance during that period. However the greatest positive specific net balance (0.67 m w.e.) was calculated for 1927/1928.
The cumulative balance neglecting refreezing over the 70 seasons shows nearly 40 m w.e. loss of mass, whereas taking into consideration Pmax = 0.6 the cumulative balance is approximately -20 m w.e. (Fig. 11). The results of cumulative balance are c.4 m higher when employing re-analysis input data.
The specific winter balance of Werenskioldbreen was calculated between 1.7 m w.e. (1927/1928 season) and 0.33 m w.e. (1947/1948 season), with the average of 0.8 m w.e. The hindcast specific summer balance was consistently a bit less variable from season to season than the specific winter balance and ranged from -2.1 m w.e in 2001/2002 season (-1.93 m w.e. in 1989/1990 season using ERA-40 data) to -0.65 m w.e. in 1916/1917 season (-0.54 m w.e. applying ERA-40 data in 1961/1962 season), and the average specific summer balance was -1.23 m w.e. No significant trend of the specific net balance was identified over the period 1912–2005. Indeed, the low positive specific winter balance trend and the negative specific summer balance trend cancel each other over a long-term period and neither of these trends is statistically significant.
ESTIMATION OF EQUILIBRIUM LINE ALTITUDE (ELA)
According to the model, the equilibrium line altitude (ELA) ranged between 203 m and 600 m (Fig. 11). As a result, the accumulation area ratio (AAR) ranged widely from 0.77 to 0.02. At the average ELA of 412 m the AAR was c. 0.35.
There is very little direct ELA data available (Fig. 11). Kosiba (1960) determined its value at 320 m in the season 1957/1958 and at 330 m in 1958/1959. The model estimated the corresponding seasonal values of 450 and 252 m, respectively. The ELA was also calculated for the seasons when the components of the mass balance were observed directly, i.e. 1994, 1999–2002. In both cases, the altitudes produced by the model may be quite different from estimates based on observation data (Fig. 11). ELA, however, is generally known to be difficult to determine and it varies in a wide elevation range (e.g. Holmlund et al., 2005; Engeset et al., 2002)
ANALYSIS OF THE MODEL SENSITIVITY TO CHANGES OF MODEL PARAMETERS
Here we analyze the sensitivity of our model of summer and winter balance to individual input parameters (Table 3). The model sensitivity was calculated as average change (%) of Werenskiold-breen mass balance elements (bw; bs) based on the data obtained from meteorological stations as well as re-analysis data in 1959–2002. The following parameters were analyzed: temperature lapse rate, precipitation gradient (τ), precipitation correction factor (K), and degree day factors (asnow, aice).
The summer balance model is determined by more parameters than the winter one. The temperature lapse rate, melt factor, as well as amount of accumulated snow (conditioned by τ and K) influence the ablation rate. The summer balance model is more sensitive to changes of asnow than aice. The sensitivity of bs model is lower when neglecting refreezing (Pmax = 0). Generally, the range of mass balance changes was similar irrespective of the applied input data set (meteorological, ERA-40).
CHANGES IN THE GEOMETRY OF THE GLACIER SNOUT DURING THE PERIOD 1958–1990 BASED ON CARTOGRAPHIC DATA AND ON ESTIMATES
The average change in glacier thickness within its 1990 outline was computed using Equation (10) at 32.4 m (1 m/year) based on meteorological data and 29 m (0.91 m/year) using ERA-40 data. A comparison of the maps yielded a change of 28.7 m (0.9 m/year). The difference between results based on mass balance modeling and geodetic is 3.7 m and 0.2 m, respectively, over 32 years (0.12–0.01 m/year).
While there is only a minor difference between the changes in average glacier thickness at its snout downstream of the Angellfjellet-Wernerknatten cross section (within the 1990 outline), the spatial error distribution ranges between 47.7 m and - 17.8 m (Fig. 12). Differences of ± 10 m are found on over 58% of the surface of the glacier snout. The glacier geometry model that uses estimated net specific balance data tends to overestimate the glacier thickness in its terminus and to underestimate it in the southeastern part (near the Angellfjellet-Wernerknatten cross-section profile; Fig. 12).
Discussion
Research presented here offers a method for hindcasting the mass balance of small Spitsbergen glaciers. We used two series of input data (observation data from meteorological stations and data from the ERA-40 re-analysis) to drive mass balance models. In both techniques, the results of the summer balance are more consistent than those of the winter balance. This results from higher correlation coefficient of both series of air temperature in summer season than winter precipitation (compare Fig. 2). The correlation of precipitation totals determined from meteorological data and ERA-40 re-analysis is weaker than the results obtained by Radić and Hock (2006). This may derive from the differences in the size of the re-analysis gridcells. Radić and Hock (2006) applied 0.5° × 0.5° grid, whereas in this work the grid was 2.5° × 2.5°. However it is difficult to estimate unequivocally the quality and usability of both series of winter precipitation totals used in the work because also data based on meteorological observations use only statistical relations between the sum of precipitation and the distance of the measurement site from the open sea. Additionally, spatial variability of precipitation may also complicate correlations.
TABLE 3
Sensitivity of the mass balance model on changes of basic parameters: temperature lapse rate, melt factor (a), precipitation/accumulation gradient factor (τ) and precipitation correction factor (K).
Correlation coefficients between the observed and modeled ablation stake balances on Midre Lovènbreen by ERA-40 derived mass balance model (Rye et al., 2010) are better than those obtained for Werenskioldbreen. Correlation coefficient for net, winter, and summer balance for Midre Lovènbreen are as follows: 0.83, 0.84, and 0.77, whereas for Werenskioldbreen they are 0.67, 0.36, and 0.72, respectively. However, Rye et al. (2010) developed a more detailed surface energy balance model than the temperature indexprecipitation model applied in this work. Rye et al. (2010) were only comparing their model to stakes located along the glacier centerline, whereas we are comparing our model to more spatially variable stake networks on Werenskioldbreen (Fig. 1).
While our model provides a satisfactory hindcasting of the glacier specific mass balance, its weakness comes from using only very rudimentary and average meteorological data from weather stations that are up to 160 km away as well as meteorological reanalysis averaged for broad area. Despite significant correlations between weather data recorded at the various Spitsbergen stations and re-analysis data, this simplification may potentially generate errors in the estimation of mass balance. Particularly serious doubts are associated with the maximum winter specific mass balance in the season of 1927/1928, especially based on the precipitation total at the Isfjord Radio.
A potential source of errors in estimating winter balance and, as a consequence, net mass balance, is the assumption of a temporally constant accumulation efficiency correction coefficient (K). As mentioned before, K integrates both the errors of precipitation measurement and the fraction of solid precipitation in the total of winter precipitation. The real precipitation in Nordic areas represents the total of measured precipitation, precipitation loss by wetting and evaporation, multiplied by correction factor due to aerodynamics effects (Førland et al., 1996). Precipitation correction factor is a function of the precipitation type, wind velocity, air temperature (in the case of solid precipitation), and the intensity of liquid precipitation, all of which change in time. Førland and Hanssen-Bauer (2000) assumed a correction coefficient of 1.15 for liquid precipitation in Svalbard, and 1.85 for solid precipitation. According to Førland and Hanssen-Bauer (2003) the fraction of solid precipitation in Svalbard has diminished in the last several decades. During the accumulation season, solid precipitation represents 44% of precipitation in Hornsund station (1978–2002) (Łupikasza, 2003). Mixed precipitation constitutes 35%, and the rest is liquid precipitation. During the same time period, the increasing trend of the fraction of liquid precipitation was recorded in the accumulation period by 40 mm w.e. per decade with a simultaneous decreasing trend of the fraction of solid precipitation by 24 mm w.e. per decade. This trend overlaps with a general increase in precipitation in Spitsbergen. Although this trend is less significant in the accumulation period than in the ablation period, the observed increase of spring precipitation in Svalbard Lufthavn station (1912–2001) amounted to 2.2% per decade (Førland and Hanssen-Bauer, 2003).
The method assumes the net mass balance to be constant at the same altitudes, which does not take into account local conditions that can affect mass balance components. The modeling results are the most satisfactory along the centerline of the glacier, which may differ considerably from conditions found on the glacier's lateral sections where the local conditions, especially the local topography, can considerably modify the accumulation and ablation processes on the glacier. This effect can be observed by comparing the modeled data with the results of net mass balance in the period 1998–2002 (Fig. 8). The best agreements between modeled and observed values were recorded in 1993/1994, when stakes were only located along the centerline of the main stream of Werenskioldbreen. The problem of spatial variability may be solved to some extent by using a distributed temperature-index model that includes potential solar radiation (e.g. Hock, 1999; Schuler et al., 2005; de Woul et al., 2006; Hock et al., 2007) and considers topographic effect on ablation distribution. That model incorporates spatiotemporal variability in degree day factors that result from variable direct solar radiation.
Local differences in mass balance are also observed in comparisons of glacier front thicknesses change calculated by the carto-graphic method and the balance-dynamic method involving Equation (10) (Fig. 12). Especially intriguing is that models estimate lowering of the southern section of the glacier snout whereas the 1958 and 1990 maps suggest that it was in fact thickening (Fig. 5). Thickening in this area may be explained either by a faster rate of emergence velocity linked to the overdeepening and a larger volume of ice or, alternatively, by ice flow instability over time. Other local factors influencing mass balance components may also have played a role in the buildup. Higher wintertime accumulation rates were observed in the area to result from snow drifts in the northern Skilryggbreen accumulation field (Grabiec et al., 2006), but also lower rates of ablation may be expected due to lower solar radiation in the shadow of the Angellfjellet. Neither can a cartographic error be ruled out, especially in the map printed in 1961. The altitude of this section in 1990 corresponds with a GPS survey of 2008. Indeed, a full model of the glacier geometry change can only be achieved when the mass balance estimation method is supplemented with an accurate model of the glacier movement dynamics and with a glacial energy balance model.
Other errors involved in the method are generated by the inaccuracy of the PDD sum that was computed on the basis of average monthly air temperatures in summer (June–September). Errors result from negative temperatures occurring throughout the year on Svalbard, including summer. This error margin was largely reduced by using an empirical formula [Equation (9)] to calculate PDD sums during months with average temperatures between 2.5 °C and - 3.5 °C. Finally, the PDD sum was converted into ablation using a degree day factor (a). The model takes into consideration changes of this factor in time and space while applying a different value for snow and ice.
Our mass balance model assumes bs without summer accumulation. This assumption is not always valid because snow can occur during the ablation season, thereby increasing albedo and slowing down melt rate. Additional energy necessary to melt out summer accumulation is not captured in our model. This effect should be small, however, because the mean sum of summer (VI–IX) solid precipitation (1978–2002) observed in Hornsund is only 6.6 mm (Lupikasza, 2003), which only requires c. 1 PDD to melt. In winter, short-term thawing may cause surface melting, but the majority of this meltwater percolates through the snow cover and refreezes, and therefore does not contribute to mass loss.
In the mass balance model, the amount of water refreezing within snow cover has been tested with the application of a local refreezing potential coefficient (Reeh, 1989). When applying Pmax = 0.6, the net mass balance of the glacier increases by 0.22 m w.e. on average, i.e. 27% of average winter balance. Application of the Pmax coefficient to the model reduced the negative cumulative balance of Werenskioldbreen by c. 50% over the whole studied period (Fig. 11). The calculation is in conformity with other findings on internal accumulation and superimposed ice on Arctic glaciers. Internal accumulation on Storglaciären (N Scandinavia) is estimated by different authors to be as little as 3–5% of annual accumulation (Schneider and Jansson, 2004) and as much as 20% of glacier winter balance (Reijmer and Hock, 2008). In an even wider range, Trabant and Mayo (1985) estimated internal accumulation on Alaska glaciers to be 7–64% of net accumulation. Super-imposed ice on Midre Lovènbreen (NW Spitsbergen), which type and size are similar to Werenskioldbreen, was estimated between 16 and 37% of net accumulation (Wadham and Nuttall, 2002; Wright et al., 2005). In the case of tidewater Kongsvegen, Brandt et al. (2008) calculated the superimposed ice rate in the range 5–10% of the winter balance.
While there will always be uncertainty in any model, we believe our model can reconstruct observations within an acceptable degree of error.
Conclusions
The method of hindcasting glacial mass balance and its components presented in this paper were able to reconstruct observations for Werenskioldbreen reasonably well. The minimum data necessary to run the model include the average monthly summer temperature at the glacier (or in its vicinity) and a precipitation total at the nearest weather station. If weather stations are not available, it is possible to use data provided by the climatic re-analysis projects, such as European Centre for Medium-Range Weather Forecasts Re-analysis (ERA-40). Our hindcasting results were verified by comparing limited observation data, which included comparing the thicknesses of the glacier snout derived from the model with those from cartographic comparisons. The model's reasonable accuracy combined with low data requirements means that the method could be widely used to estimate specific mass balances for Spitsbergen glaciers of a similar type and size.
Acknowledgments
The authors would like to thank the Institute of Geophysics of the Polish Academy of Sciences, and the Norwegian Meteorological Institute for granting access to meteorological data and mass balance data as well as ECMWF for access to the ERA-40 data set. We also would like to express our thanks to Dr. Leszek Kolondra for making the DEM of the glacier available and Dr. Jan Leszkiewicz for unique meteorological data. This research has been financially supported for scientific research in 2007–2010 by research grants from the Polish Ministry of Science and Higher Education (IPY/269/2006). The final stage was supported from the Polish-Norwegian Research Found on the project “Arctic Climate and Environment of the Nordic Seas and the Svalbard-Greenland Area—AWAKE” (PNRF-22-AI-1/07). The authors wish to thank two anonymous reviewers and the Associate Editor, whose comments improved this paper. Jason Gulley helped in smoothing the text.