Introduction

The remarkable mass losses that currently affect glaciers in the European Alps may lead to significant future alterations in the hydrology of glaciated catchments. The glacial stream flow regime shows a typical maximum in mid- to late summer, meeting the human requirements for water supply. In addition, glaciers have a damping effect with respect to climatic variations, storing mass during cold and wet periods and releasing it in warm, dry years. This mitigating capability sharply decreases as climatic forcing continues, leading to a strong reduction in glacier area and volume.

Knowledge of the spatial and temporal distribution of glacier mass balance is crucial both for climatic sensitivity studies and for understanding the hydrological behaviour of headwater basins (Reference Fountain, Jansson, Kaser and DyurgerovFountain and others, 1999). For hydrologic investigation purposes, particularly, the glaciers that feed the main stream of a catchment should be regarded as a single system, characterized by complex behaviour and a strongly differentiated response to climatic changes.

Direct mass-balance measurements are only feasible on a few suitable glaciers (Reference Østrem and BrugmanØstrem and Brugman, 1991; Reference Kaser, Fountain and JanssonKaser and others, 2003). Extensive surveys are not practical because measurement zones should meet certain requirements such as appropriate size, altitudinal range, geometry, accessibility and safety. Such constraints often lead to data being only partially representative of the whole glaciated area. The spatial interpolation and extrapolation of point mass-balance data are prone to errors, in particular when data are few and sparse. In this case, computer plotting results are often poor when compared with hand-drawn mass-balance maps (Reference Kaser, Fountain and JanssonKaser and others, 2003); however, manual mass-balance interpolations and extrapolations involve some subjectivity. The mass-balance gradient with altitude can be reliably estimated even with a small network of stakes, provided they are scattered over the elevation range of the glaciated area (Reference Fountain and VecchiaFountain and Vecchia, 1999). However, special care should be taken both in the case of strong lateral mass-balance gradients (Reference Holmlund and JanssonHolmlund and Jansson, 1999) and when extrapolating field data to nearby glaciers with different sizes, geometries, feeding sources, debris cover and surrounding topography. Reference KrimmelKrimmel (1999) attributed the discrepancy between mass-balance estimates by direct and geodetic methods on South Cascade Glacier, Washington, USA, to possible uncertainties in the area-integration procedure, and to sinking stakes.

The mass-balance spatial distribution pattern is quite similar from year to year, since it is mainly controlled by the local topography, which affects both accumulation and ablation processes. Therefore, it is interesting to verify the feasibility of using topographically based attributes to assess the spatial variability of these dominant processes, by means of regression analyses.

The aim of this work is to develop and test a simple method to extrapolate direct mass-balance measurements in unsampled areas of benchmark glaciers or in nearby unsurveyed glaciers. After a description of the study area and field campaigns, we explain how the extrapolation scheme was implemented and applied to the glacial system of a headwater basin of the central Italian Alps. Thereafter, we focus on the results and critically discuss the applicability of the technique, its limits and possible improvements.

Study Area Description and Data Collection

The area we investigated is the upper sector of the Val de La Mare (UVLA), Trentino, Italy (see Fig. 1). Its physical and geographic characteristics make it representative of alpine, high-altitude catchments. The basin has an area of 36.2 km² and an average altitude of 2906 m, ranging from 1931 to 3769 m, at the Cevedale peak. The catchment lies in the southern sector of the Ortles-Cevedale massif, the most glaciated area in the Italian Alps. At present, 25% of its surface is covered by four glaciers (9.2km²) and the annual runoff regime is dominated by snow- and ice melt. In spite of its high average altitude, the area is easily accessible and comprises the glacier with the longest mass-balance data series in Italy (Careser). Moreover, for the Careser Diga station (2605 m) the weather data series extends back to 1930.

Fig. 1. Geographical setting of the study area.

The Val de La Mare is located near the inner dry Alpine zone (Reference SchwarbSchwarb, 2000). Therefore precipitation is less abundant than in the nearby Adamello–Presanella and Brenta massifs and the average equilibrium-line altitude (ELA) is higher. At Careser Diga, the mean annual precipitation is 927 mm and the October–May precipitation averages 556 mm. The mean annual air temperature is –0.7°C, while from June to September it averages 6.0°C. The altitude of the mean annual zero isotherm is 2500 m.

Direct mass-balance measurements have been performed since 1967 on Careser glacier, and extend until the present without interruption (Reference Carturan and SeppiCarturan and Seppi, 2007). As can be seen in Table 1, the Careser glacier geometry is peculiar when compared to the entire glacial system, 95% of its surface being concentrated between 2858 and 3200 m. A long series of strongly negative mass balances has led to rapid surface lowering and fragmentation as the entire glacier lies below the recent years’ average ELA (∽3400 m). In addition, albedo and elevation feedbacks emphasize the climate forcing.

Table 1. Physical characteristics of UVLA glaciers. The elevations are in metres above sea level. ELA65 and ELA50 are the equilibrium-line altitudes corresponding to accumulation-area ratios of 0.65 and 0.50, respectively

The other three glaciers (Vedretta Rossa, Venezia and La Mare) are located on the western slope of the basin and are characterized by higher mean altitudes. Their geometry suggests a possible minor mass change in response to the recent rise of ELA and a lower response time than Careser glacier (Reference Furbish and AndrewsFurbish and Andrews, 1984; Reference Benn and EvansBenn and Evans, 1998; Reference Pelto and HedlundPelto and Hedlund, 2001). This is confirmed by the terminus behaviour in the last several decades, characterized by continuous retreat of Careser and temporary but substantial readvance of the other three glaciers from 1970 to 1985. The presence of crevassed areas and the observed displacement of boulders on the surfaces are evidence of flow on these three glaciers, while repeated global positioning system (GPS) surveys of ablation stakes on Careser glacier denote a complete absence of movement. The different activity of the four glaciers is also evident from the accumulation-area ratio (AAR) values from 2004 to 2006 (Table 2). Vedretta Rossa, Venezia and La Mare appear less in disequilibrium with the present climatic conditions.

Table 2. Accumulation-area ratios for the UVLA glaciers in the years 2004–06 and corresponding anomalies in accumulation season precipitation and ablation season temperature from the 1967 to 2006 averages, at the Careser Diga weather station

In order to set up an adequate experimental basis for our research activity, a series of field surveys was undertaken in the past 5 years. In the summer of 2003, we started direct mass-balance observations over the more representative Vedretta de La Mare, with the aim of completing our knowledge about the mass-balance behaviour of the entire glacier system. Winter and summer balance measurements have been carried out since 2004 in order to understand the climatic causes of mass changes. An automatic weather station (AWS) has been operated from the summer of 2005 on the surface of the glacier, and in 2007 portable albedo measurements were conducted on Careser and La Mare glaciers. Repeated mapping of snow cover on the four glaciers was performed using terrestrial photographs at 2 week intervals during summer.

Vedretta de La Mare can be considered a ‘difficult’ glacier since it presents highly crevassed areas, seracs and steep slopes, and it is not possible to survey the whole surface in one day. Therefore, we focused on the southern sector, which feeds the main tongue (Figs 1 and 2).

Fig. 2. The four glaciers of the UVLA catchment. (a) Vedretta del Careser, (b) Vedretta de La Mare, (c) Vedretta Venezia and (d) Vedretta Rossa. Photographs taken by Reference Carturan and SeppiL. Carturan in 2007.

Winter balance (b _w) measurements were carried out on both glaciers in the last 10 days of May, combining snow-depth soundings with density measurements in snow pits for the calculation of point snow water equivalent (w.e.). On Careser glacier, the entire area was investigated. On La Mare glacier the soundings were limited to safe areas, and in 2004 and 2005 specific points were selected in accordance with the snow-cover patterns observed in previous summers. In 2006, the soundings were performed at regular distances along rectilinear transects. The b _w data were little affected by ablation, which was limited to the lowest area of La Mare glacier.

Summer balance (b _s) measurements were conducted in

the last 10 days of September with 2m aluminium stakes, drilled into the ice or firn with a hand auger; the stakes located in firn were plugged at the lower end with insulating material, in order to avoid sinking. In the

accumulation area, snow-depth soundings were taken and snow pits were dug both for density measurements and for the identification of previous-year layers. The net accumulation was measured on the maximum accumulation areas, where residual snow persisted even in the extremely warm 2003 summer.

Both glaciers are currently experiencing strong mass losses, which are mainly controlled by b _s and only occasionally heightened by low b _w (Reference ZanonZanon, 1982; Reference Carturan and SeppiCarturan and Seppi, 2007).

For the years 2004, 2005 and 2006, an average number of 125 b _w and 32 b _s measurements for Careser and La Mare glaciers was usable for the implementation and validation of the extrapolation scheme. As reported in Table 2, the accumulation season precipitation was slightly below average in 2004 and 2005 and well below average in 2006. Ablation seasons were warmer and longer than normal. In 2004 the peak of ablation was reached in late

summer; in 2005 and 2006 it occurred between June and July, followed by a cool August with snowfalls that inhibited ablation and stopped it at the higher elevations.

The topographic attributes used as independent variables in this study were derived from a 10 m grid resolution digital terrain model (DTM) generated by processing high-resolution airborne laser altimetry (lidar) data, surveyed in September 2006. In order to assure a wider applicability of the extrapolation method, no meteorological variables were included.

The Extrapolation Method

In the three observation years, residual snow patches were scarce and mostly retreated in remote areas at high altitude. Therefore, the net accumulation data that could be acquired were much fewer than the net ablation data, and not sufficient for the implementation of an extrapolation technique in the accumulation area.

Consequently, we decided to focus on b _w and b _s data rather than net balance, b _n. This approach also has the advantage of separately handling the accumulation and ablation processes, whereas the net quantities result from the combination of both.

The tested methods of extrapolation are single linear regressions and multiple regressions. The regression coefficients of the relationships obtained from the surveyed areas are used to calculate the mass balance on other, unsurveyed, areas.

Winter balance extrapolation

Both wind and avalanches relocate snow, forming a deposition pattern that is primarily governed by topography and is quite similar from year to year. This deposition pattern affects the runoff timing and magnitude during the melt season and the mass-balance distribution on glaciers (Reference Luce, Tarboton and CooleyLuce and others, 1998; Reference HartmanHartman and others, 1999).

The independent variables selected to explain the winter balance spatial variability are reported in Table 3. In addition to the elevation, E, slope, S, and aspect, α, we used a relative elevation attribute, REA _r , calculated as follows:

Table 3. Correlation, R, for b _w and independent variables (E: elevation; S: slope; α: aspect; REA _r : relative elevation attribute). Correlations significant at the 0.05 level are bold type; degrees of freedom are in parentheses

where is the average elevation in a circular area with radius r (in metres) around the gridcell with elevation E _DTM. This topographic attribute expresses the degree of wind exposure of each gridcell. Its values are positive in hollows, depressions and at the foot of rock walls or steep slopes and negative on crests, peaks and generally convex areas. An interesting feature of this attribute lies in the possibility of adjusting the radius and ultimately the extent of the area that affects the local attitude to over- or under-accumulate snow.

As can be seen in Table 3, b _w is directly correlated with E, S and REA _r , but not with α. The direct correlation with E is explicable by the positive precipitation gradient, even if above a certain altitude a reversal of the gradient is expected to exist; moreover, in most glaciers, b _w is likely to decrease in the upper accumulation area because of snowdrift (Reference Oerlemans and HoogendoornOerlemans and Hoogendoorn, 1989). In addition, the inclusion of E among the independent variables allows us to account for possible ablation losses in the lower part of the glaciers. The direct correlation with S does not have a similar physical basis and appears to be coincidental, since snow accumulation is expected to decrease with increasing S, due to avalanche action.

The b _w is directly correlated with REA _r and correlation increases up to a radius of 80–100 m, then decreases. These results differ from a preliminary application of REA _r to a bare-ground area of the same catchment, located between 2600 and 2800 m. In that case, the correlation between the snow w.e. and REA _r was higher (0.689), and optimized by a radius of 50 m. This is probably due to greater terrain roughness, different snow metamorphism and lower wind speeds.

Two types of extrapolation schemes were implemented and compared. The first, W1, was a simple linear regression between E and b _w, while the second, W2, was a multiple linear regression among b _w, E and REA _r .

Summer balance extrapolation

In general, most of the energy used for melt on glaciers is supplied by radiation, followed by turbulent heat fluxes. Air temperature is a good predictor of melt since it is highly correlated with incoming longwave radiation and turbulent heat fluxes (Reference OhmuraOhmura, 2001; Reference HockHock, 2005). In addition, it controls the fraction of precipitation that falls as snow, which is crucial during the ablation season since it stops the ablation and increases the surface albedo. For this reason, the summer balance on glaciers is strongly altitude-dependent.

However, our measurements at the AWS of La Mare glacier demonstrate that the net shortwave radiation dominates the energy flux during melt, confirming other observations with AWS operated on glaciers (Reference Hock and HolmgrenHock and Holmgren, 1996; Reference OerlemansOerlemans, 2000; Reference Greuell and SmeetsGreuell and Smeets, 2001; Reference Strasser, J., Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004). Therefore, the use of measured or calculated shortwave radiation in addition to E should optimize the capability to account for the high spatial variability of b _s.

In confirmation of this, simple temperature-index ice- and snowmelt models have been appreciably improved by the inclusion of clear-sky direct or global radiation, yielding a more realistic temporal and spatial distribution of melt rates (Reference Cazorzi and Dalla FontanaCazorzi and Dalla Fontana, 1996; Reference HockHock, 1999; Reference Cazorzi, Carturan and Dalla FontanaCazorzi and others, 2005; Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005).

Since the albedo differences among melting substrata (snow, firn and ice) strongly affect the melt rate, b _s data showed a spatial variability which was related to the relative amount of these substrata, and not just a function of topography. In order to spotlight the topographic effect on spatial data variability, they were preliminarily homogenized, using the average albedo differences among snow, firn and ice, as measured by the AWS and by a portable albedometer.

A preliminary application of the grid-distributed melt model developed by Reference Cazorzi, Carturan and Dalla FontanaCazorzi and others (2005) in the UVLA basin yielded encouraging results, but a spatial aggregation of simulation errors was clearly observed. The tuned model underestimated ablation on east-facing slopes and overestimated it on west-facing ones. The comparison of AWS-measured shortwave global radiation versus computed clear-sky radiation revealed a persistent diurnal cycle of cloud cover (mainly cumuliform clouds) during summer, with a minimum around sunrise and a maximum at sunset. This cycle is mainly triggered by the sun through thermal convection, which explains the timing of minima and maxima of cloud cover (Reference Barry and ChorleyBarry and Chorley, 1998).

Three different extrapolation schemes for b _s were implemented and compared. The first, S1, was a simple linear regression with E; the second, S2, was a linear multiple regression that included E and daily-averaged clear-sky radiation (CSR). The third, S3, was similar to the second, but used a normalized cloud-cover factor (NCCF) to reduce the CSR as a function of aspect.

Due to the observed diurnal cycle of cloud cover, the fraction of CSR available for melt is largest for the east-exposed sites and lowest for the west-exposed slopes. This phenomenon can be statistically simulated using a simple harmonic function that sets NCCF to 0 for the gridcells with aspect, α, matching the sunrise azimuth and to 1 for the gridcells with aspect matching the sunset azimuth:

in which α ₀, α ₁ and ξ depend on daily sunrise and sunset azimuth. The reduced clear-sky radiation (RCSR) is then calculated as follows:

where k is a calibration parameter ranging from 0 to 1, which expresses the maximum reduction of CSR by cloud cover for the gridcells with α matching the sunset azimuth. It should be noted that the parameter k does not require any additional meteorological data for its determination; in fact k is selected by subsequent trials to fit the experimental b _s data as well as possible.

Results and Discussion

The determination coefficients of the linear regressions between b _w and independent variables are shown in Figure 3. The statistics denote a sharp and persistent improvement by including REA _r on the linear regression; E alone accounts for less than half of the variance explained by multiple regression. The relationship of each regressor with the dependent variable is significant at the 0.05 level in t tests. The optimal radius for the calculation of REA _r is 120 m, although in 2006 the statistics indicate nearly identical results with 120, 140 and 160 m. Better results, in terms of extrapolations, have to be expected for 2004 and 2005, since in 2006 the correlation is noticeably lower.

Fig. 3. Determination coefficient, R ², of the linear regression between b _w and independent variables.

The scatter plots in Figure 4 show comparisons between measured and calculated b _w, for the two extrapolation methods. The multiple regressions better match the observed data, the higher values in particular. A point-by-point comparison of measured vs calculated values revealed that the altitude alone underestimates b _w on snowdrifts and overestimates it for wind-scoured areas; whereas REA _r tends to overestimate b _w at lower altitudes and to underestimate it at higher altitudes. Combining the two attributes with a multiple linear regression considerably overcomes these drawbacks, even if two limitations are implied in the simplicity of this extrapolation method: the first is a scale dependency of snow redistribution, which varies with altitude and is not captured by a constant radius of 120m in the calculation of REA _r ; the second is the inability to catch the drifts near to abrupt slope changes.

Fig. 4. Measured versus calculated b _w using the W1 method (b _w and E) and the W2 method (b _w, E and REA _r ). Empty squares: La Mare; filled squares: Careser. The straight line indicates a 1 : 1 relationship.

The high correlation between altitude and b _s ensures better estimates on unsampled areas compared to b _w. However, if b _s data come from an unrepresentative area of a glacier system, care should be taken in the extrapolation since systematic errors are possible.

Figure 5 shows the R ² values for the three different extrapolation schemes of b _s. For S3 four RCSR _k were used, in order to optimize the value of the calibration parameter, k. The elevation alone explains a large proportion of the variability in b _s. The effect of CSR on b _s is statistically significant in 2004 (t-test significance <0.05), while in 2005 and 2006 this improvement is not so clear (t-test significance = 0.128 for 2005 and 0.054 for 2006). Better results for these two years are achievable by the simplified simulation of the diurnal cycle of cloud cover to reduce CSR. For 2005 and 2006 a k value of 0.2 proves to be optimal and the trends of R ² with the various combinations of independent variables are quite similar. The 2004 results show different behaviour, and the best correlation is achieved with a k value of 0.1. These differences are likely to be due to the climatic conditions of the three ablation seasons, as described above. The t tests for calibrated RCSR _k provided a significance level always <0.05.

Fig. 5. Determination coefficient, R ², of the linear regression between b _s and independent variables.

The scatter plots in Figure 6 compare the observed and calculated b _s using the three extrapolation methods. The plots of measured vs observed values show a closer alignment with the 1 : 1 line when RCSR _k is used. The scatter of points also diminishes, and with few exceptions the improvement from using CSR and RCSR _k is evident. As observed in the regression analysis, the inclusion of the cloud-cover daily cycle yields better results in 2005 and 2006 than in 2004.

Fig. 6. Measured versus calculated b _s with the S1 (b _s and E), S2 (b _s, E and CSR) and S3 (b _s, E and RCSR _k ) methods. Empty squares: La Mare; filled squares: Careser. The straight line indicates a 1 : 1 relationship.

The results of the validation tests are shown in Tables 4 and 5, where the accuracy of each extrapolation method is expressed in terms of R ² and root-mean-square error (rmse). For winter balance (Table 4), the tests confirm the findings of the regression analysis, since the multiple regressions yield higher R ², lower rmse and b _w that ranges closer to the observed values, in almost all cases. For summer balance (Table 5), the method which combines E and RCSR _k provides consistently better results. In fact, the CSR-based extrapolation is sometimes worse than the simple linear regression of summer balance vs altitude. Specifically, this happens when Careser data are used to infer La Mare data, or vice versa, leading to systematic under- or overestimations of b _s. The different prevalent aspect of the two glaciers makes them differently affected by the cloud-cover daily cycle. The reduction of CSR by clouds is more effective on Careser glacier due to a wider presence of southern and western slopes.

Table 4. Accuracy of the two b _w extrapolation methods, W1 (b _w and E) and W2 (b _w, E and REA _r ). Bold type indicates a determination coefficient, R ², significant at the 0.05 level and n is the sample size

Table 5. Accuracy of the three b _s extrapolation methods S1 (b _s and E), S2 (b _s , E and CSR) and S3 (b _s , E and RCSR _k ). Determination coefficients, R ², are significant at the 0.05 level and n is the sample size

The best-performing regression methods were combined to extrapolate the field data and calculate the net annual mass balance, b _n, on the surface of the entire UVLA glacier system. A multiple regression of E and REA₁₂₀ was used to extrapolate b _w. For b _s a multiple regression of E and RCSR _k was employed, with k ranging from 0.1 (for 2004) to 0.2 (for 2005 and 2006). The difference between b _w and b _s (b _n) calculated for the entire glacier system was –450, –1064 and –1011mmw.e. in 2004, 2005 and 2006, respectively.

The reliability of the extrapolation method was assessed against b _n measured at the ablation stakes and the snow-cover maps surveyed during the summer. This was done with the assumption that the isolines of annual mass balance run roughly parallel to the transient snowline (Reference Kamniansky and PertzigerKamniansky and Pertziger, 1996).

In Figure 7, the 2004 snow-cover distribution on four dates is reported, alongside the calculated b _n. As can be seen, the multiple-regression method better matches the lower boundary of the snow cover at the end of the ablation season, in particular on north-facing slopes. Above the equilibrium line, however, the calculated residual snow cover is more continuous than in reality. Furthermore, b _n appears to be overestimated at the highest altitudes and on north-facing steep slopes. In the ablation area, the distribution of b _n calculated by linear regression with E appears to be excessively simplified with respect to the observed snow-cover pattern. The use of multiple regression yields more realistic results, and some specific features are better represented. Among these is the reduction (increase) of ablation on north-facing (south-facing) slopes and in high-accumulation (low-accumulation) areas.

Fig. 7. Snow-cover maps for four dates in 2004 (left) and net-balance maps for the same year, obtained by single linear regression of altitude (centre) and by multiple regression (right). The black dots are the ablation stakes. Careser glacier has been translated to fit in the figure.

The scatter plots in Figure 8 show the observed and calculated b _n at the ablation stakes in 2004. With few exceptions, the results of the multiple-regression method are in closer agreement with measurements. The improvement is greater on Careser glacier. On La Mare glacier, good results are obtained both in the ablation and accumulation areas, even if differences up to 1200mm persist at single points. Similar results are obtained for 2005 and 2006. A steady overestimation of net ablation can be seen in proximity to the observed ELA, which is probably related to the strong mass-balance gradient in that area.

Fig. 8. Measured versus calculated b _n at the ablation stakes of Careser glacier (filled squares) and La Mare glacier (empty squares) in 2004, on the left by simple linear regression with E, on the right by multiple regression.

The discrepancies observed in the final b _n output (local overestimation of accumulation areas and of mass balance) are primarily due to incorrect b _w estimates in areas of avalanche redistribution, at high altitude and near abrupt slope changes. The strong mass-balance gradient observed in proximity to the ELA is probably due to fresh summer snow, which tends to persist longer on snow-covered than on bare-ice areas. To capture the complexity of these processes is beyond the aims of this work.

Concluding Remarks

Topographic attributes easily derived from a DTM (E, α, CSR, NCCF, REA _r ) allow for a satisfactory interpretation of the dominant processes governing the mass exchange on glaciers. A regression analysis of accumulation and ablation data surveyed during three years was conducted separately for b _w and b _s, and was more problematic for the former. The amount of explainable variability in the observed data was nearly double for b _s with respect to b _w. Moreover, whereas use of the altitude alone allows for reasonable extrapolation of ablation data, the same cannot be said for accumulation.

The use of an REAk that accounts for the wind exposure of a site improves the ability to explain the spatial variability of the accumulation processes. Indeed, snow redistribution on glaciers is mainly controlled by wind action. The use of clear-sky radiation in the summer balance regression equation is greatly improved by the introduction of a simple function that relates aspect to cloud cover (NCCF), allowing for a correction of CSR according to the daily cycle of cloud cover.

The final regression equations were tested and used to extrapolate the net balance to the entire glacier system of the basin. Some aspects remain unsolved. Among these is the snow redistribution by avalanches, the formation of snowdrifts in the proximity of abrupt slope changes or on the lee side of obstacles, and the decrease of winter balance above a certain altitude due to increased effectiveness of wind erosion.

An ultimate approach probably requires a fully distributed energy- and mass-balance model, but it is worth noting that in the method proposed here no other information besides topography is required.