Introduction

Predictive use of watershed runoffmodels and soil–vegetation–atmosphere transfer (SVAT) models often fails in landscapes with a heterogeneous forest cover and a complex topography, because these characteristics strongly influence the micrometeorological surface conditions. This applies especially to winter conditions in subalpine semi-forested areas where both the snow accumulation and the snow ablation are strongly affected by canopy interception of water, snow and radiation, as well as by the complex wind field within such a terrain. In order to improve predictive use of models for such areas, we need to include information on the spatial structure of the landscape at a much lower resolution than the watershed scale, preferably at the scale of the forest mosaic.

In recent years, most of the remote-sensing studies addressing snow-cover distribution had the goal of determining the areal snow coverage within a certain region. The quantitative areal information about the snow coverage is an important feature for mountain hydrology (Reference Rango and ShalabyRango and Shalaby, 1998). While in the beginning the objective was to classify the pixels as snow-covered or snow-free (Reference DozierDozier, 1989), the more recent works tried to map the fractional snow coverage at subpixel resolution (Reference Rosenthal and DozierRosenthal and Dozier, 1996). In combination with information from the ground or simple snowpack models such maps were then used to estimate the snow water equivalent and the runoff during the snowmelt in spring (Reference SchaperSchaper, 2000; Reference Schaper, Seidel and MartinecSchaper and others, 2001).

The success of this strategy is strongly influenced by the image resolution and the overflight frequency of the satellite: the Landsat Thematic Mapper (TM) has a good resolution, but a limited number of overflights per winter season; the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite records an image of the selected area every day, but with a rather poor resolution.

Unfortunately, mapping fractional snow cover using optical remote sensors is very difficult for forested areas (Reference Solberg, Hiltbrunner, Koskinen, Guneriussen, Rautiainen and HallikainenSolberg and others, 1997) because of the highly complicated solar illumination (direct or diffuse) of snow below a canopy. In a current analysis of Landsat TM images, Reference Vikhamar, Solberg, Wunderle and NaglerVikhamar and Solberg (2001) investigate the influence of different tree species and densities on the reflectance of forest with a snow-covered ground in the Jotunheimen mountains, Norway. Their results confirmed the complexity of the reflectance pattern caused by, among other things, heterogeneous illumination of the forest ground, the reflection properties of different tree species or snow interception. A detailed reflectance model including shadowing effects was able to explain no more than about 60% of the variation in reflectance. Another study of the reflectance of snow-covered forest was presented by Reference Klein, Hall and RiggsKlein and others (1998). They improved mapping of snow-covered area for forests from MODIS and Landsat TM images by using two parameters (normalized-difference snow index and normalized-difference vegetation index (NDVI)).

In contrast to the areal extent of the snow cover, it is not feasible to map the snow-depth distribution directly from optical sensor satellite images. Therefore, indirect methods have to be applied to model the snow depth, based on assumptions about the influence of topography and forest density on the snow pattern. Reference ForsytheForsythe (1999) showed that multiple linear regression models using, among other things, Landsat TM channels as explanatory variables of vegetation density were able to reproduce the snow-cover distribution in the Kalkhochalpen region of Austria and Germany fairly well.

In this study we use spatial data at different resolutions from (i) digital elevation maps (gridsize 5 m), (ii) a Landsat TM image (georectified and resampled to a pixel size of 25 m), (iii) aerial photographs (pixel size 0.8 m) and (iv) canopy photographs from the soil surface (support: 60 m²) to derive attributes of topography and vegetation density. We hypothesize that such spatial data allow us to predict the snow-depth distribution in heterogeneous subalpine areas with a high accuracy and a fine spatial resolution. The main objective of this work is to develop a detailed snow-depth distribution map for a subalpine semi-forested catchment based on (a) a parsimonious linear regression model which includes the attributes of topography and vegetation indices as explanatory variables, and (b) geostatistical interpolation techniques.

Site Description And Measurements

The study area, Erlenbach, is a 0.7 km² hydrological catchment situated in the central Swiss subalpine zone. It rises gently from the lowest point at 1100 ma.s.l. up to an altitude of 1600ma.s.l. facing westwards (Fig. 1). Typically for this type of landscape, the forest in the catchment consists of a heterogeneous mosaic of coniferous trees with marshy meadows and open pasture in between (about 40% forest). In contrast to alpine areas, redistribution of snow by wind is not significant in this catchment, apart from the open pasture close to the crest where snow may be relocated by wind drift.

Fig. 1. Digital elevation model superimposed with an aerial photograph of the Erlenbach catchment.

The spatial snow measurements presented in this paper are from two winters, 1998/99 and 1999/2000, that were rich in snow, with peak snow depths of 2–3m in the upper part of the area. On 4 days of the first winter and 2 days of the second winter, we measured snow depth at 100–200 locations, fairly evenly distributed within the catchment (Fig. 3a). The support of the snow-depth measurements was approximately 4 m². At most of the forest locations, we took a photograph of the canopy with a hand camera and digitized the picture. Using a simple threshold operation in the blue channel, we classified the pixels forming the canopy and used their relative frequency as a measure for the canopy closure, f_c (support: ~50m²; Reference Stähli, Papritz, Waldner and ForsterStähli and others, 2000).

Fig. 3. (a) Spatial distribution of the snow depth measured on 23 December 1999 and 4 April 2000, (b) areal forest fraction determined from the aerial photograph, and (c) solar exposure (Equation (3)) with positive values facing southwards distributed over the Erlenbach catchment.

A Landsat TM image from 15 June 1996 was selected, showing a cloudless early-summer landscape with no snow. The image was georectified using a digital elevation model provided by the Swiss Federal Office of Topography and resampled to a pixel size of 25×25 m2. No atmospherical correction was made. From the Landsat TM channels 1–5 and 7 we calculated vegetation indices (Table 1) that had been reported as feasible indicators of forest density in earlier studies (e.g. Reference Crist, Laurin and CiconeCrist and others, 1986; Reference CohenCohen, 1991; Reference Borry, de Roover, Leysen, Wulf and GoossensBorry and others, 1993; Reference Koch, Förster and MünstererKoch and others, 1993).

Table 1. Vegetation indices used for the multiple regression analysis

Analysis Of The Forest Structure In The Erlenbach Catchment

Since the forest density strongly influences the snow accumulation and ablation in the Erlenbach catchment (see next section), we decided to analyze the spatial structure of the forest canopy. To this end, we calculated sample variograms of selected channels of the aerial photograph and the Landsat TM image of June 1996 that were assumed to be sensitive to the forest density. The semivariance of a variable a(x), which in this study could be either the grey value of a pixel or measured snow depth, may be estimated by:

(1)

where h is the lag distance between two points, and n the number of locations separated by lag distance h. An attribute is said to be spatially dependent if γ(h) increases with h. Details on geostatistics can be found in Reference Webster and OliverWebster and Oliver (1990).

The sample variograms of the green channel of the aerial photograph and the Landsat TM 4 image showed that there are two distinct scales of spatial dependence: one representing the small-scale heterogeneity of the forest canopy (~10m; single-tree scale, Fig. 2a), the other representing the size of the forest patches (~100 m; forest mosaic scale, Fig. 2b). Although the heterogeneity at the single-tree scale obviously cannot be captured by the Landsat TM image, there are methods available to estimate the areal fraction of forest per pixel. One such method is to take a high-resolution aerial photograph, classify the forest (e.g. using a simple threshold algorithm) and then average the classified image to the larger pixel size of the Landsat TM image.

Fig. 2. Semivariograms for pixel values of the aerial photograph (left) and Landsat TM4 (right).

We applied this method to our catchment using the green channel of the aerial photograph. The resulting map of areal forest coverage (Fig. 3b) was related to the vegetation indices determined from the Landsat TM image. A multiple linear regression model with the three explanatory variables TM4, TM7 and brightness (TC_br) explained 72% of the variation in forest coverage. The single vegetation index that represented the forest coverage best was TC_br with an R₂ value of 0.66, followed by TM4 with an R₂ value of 0.55. On the other hand, the forest coverage was only weakly correlated with NDVIs (R₂<0.29) and the simple ratios (R₂<0.34).

From this analysis we concluded that TM4, TM7 and TC_br offered the most promising way of representing the forest density in this type of subalpine landscape. We therefore expected that these vegetation indices would be best suited to model the spatial distribution of the snow depth in the Erlenbach catchment.

Multiple Linear Regression Analysis Of The Snow Depth

As a next step towards a model for the snow-cover distribution, we determined a parsimonious multiple linear regression model for the snow depth, z_snow, within the catchment. We started from a full model including all the attributes of topography and vegetation indices derived from the Landsat TM image (Table 1) as explanatory variables, and reduced the set of explanatory variables by a stepwise procedure to a minimum set that best explained the variation of snow depth in space. We fitted a common model to the data of all the six measurement dates. The change of the average snow depth in time was modelled by a time-dependent constant, leading to the following model:

(2)

where z_alt is altitude a.s.l., I_i(t) are time indicators for the respective measurement dates, ε is the error term, and e_t is the solar exposure of the terrain assuming a midday-sun zenith angle, α_sz, of 60° (which corresponds approximately to an average winter-season sun zenith at the latitude of the Erlenbach catchment, 47° N; Fig. 3c):

(3)

where α_t is the aspect deviation from south and s_t is the slope angle of the terrain.

Obviously, the parsimonious model found by the stepwise reduction included only some of the vegetation indices that were found to represent forest density best (see previous section): TM4 and TC_br. On the other hand, SR₄₃ and TC_wet had a larger influence on the snow depth than TM7. The coefficient of determination, R₂, of model Equation (2) was 0.704 (Table 2). A linear regression model with the three vegetation indices representing the forest density best gave only a slightly smaller coefficient of determination than the model Equation (2).

More important than the difference between vegetation indices was that the topography, i.e. z_alt and e_t, and the time-dependent constants alone resulted in an _R2 value only marginally smaller than the coefficient of determination of model Equation (2). Thus, we concluded that the vegetation indices of the Landsat TM image explained hardly any of the variation in snow depth. This applies not only to locations that were mapped as open land during the snow measurements, but also to the forested or semi-forested measurement locations.

For single measurement dates (exemplified in Table 2 for 23 December 1999, an early-winter situation, and 4 April 2000, a late-winter situation), we found best-fitting linear regression models that differed slightly from the parsimonious model for all measurements (Equation (2)). However, the gain in precision using the best-fitting models for given dates against model Equation (2) was only minor. We also observed that the vegetation indices were more influential, and the topography less dominant, when the model was fitted to single dates.

To sum up, the Landsat TM image contained hardly any information on the forest density that could be used to model variation of snow depth at the single-tree scale. To see if the information on the forest density at the single-tree scale would improve the explanation, we went on to fit a linear regression model that included the canopy closure, f_c, determined from the photographs taken from below the canopy as the only vegetation-related explanatory variable. Clearly, this variable was much more important than the Landsat TM vegetation indices. This was found for the pooled dataset of all six measurement dates, but also when fitting the model to single measurement dates. A particularly impressive gain of precision was obtained for the data of 23 December 1999 where the model with the canopy closure explained as much as 76% of the snow-depth variation, whereas the R₂ of the best regression models with Landsat TM images was not higher than 38%.

So far, we have only discussed regression models where z_snow was untransformed. Unfortunately, we noted for the pooled dataset of all six measurement dates that the residuals of the fitted models were not normally distributed. A square-root transformation of z_snow resulted in residuals that were approximately normally distributed. Obviously, this transformation did not significantly impair the coefficient of determination (Table 2). Fitting the models to the untransformed z_snow of the single measurement dates resulted in approximately normally distributed residuals.

Interpolation Of The Snow Depth Using Universal Kriging

To analyze the spatial structure of the error term σ, we computed the residuals, r, of the square-root-transformed snow depth for the model in Equation (2), fitted to the full dataset (n = 853), shown in Figure 4 (top) for 23 December 1999 and 4 April 2000. Considering only the locations on the 75 675 m² grid, the patterns imply some continuity in the spatial distribution of the residuals. Over shorter distances, however, there were pronounced differences in the magnitude of the residuals: on 4 April, for example, large differences between adjacent locations were found at (697.5, 211.1) and (697.3, 211.3).

Fig. 4. Spatial distribution of the residuals of the regression model (Equation (2)) on 23 December 1999 and 4 April 2000 (top). The area of the circular symbols denoting the position of the measurement locations is proportional to the residuals shifted by a positive constant. The bottom graph shows the sample variograms for the corresponding dates. Here the area of the circles is proportional to the number of data pairs for a given lag class.

We then computed the sample variograms (Equation (1)) of the residuals for the individual measurement dates. Figure 4 (bottom) shows these statistics for 23 December 1999 and 4 April 2000. In general, the semivariance depended hardly at all on the lag distance when analyzed for the six dates individually: for example, the average squared difference over the shortest lag distance was as large as (23 December) or even larger than (4 April) the differences found for longer lags. Since there was little to no spatial dependence discernible in the variograms of individual dates, we computed the variogram of the pooled dataset, and this statistic implied that the residuals were spatially correlated up to a distance of 0.5 km with a nugget/sill ratio approximately equal to 0.5.

We checked this model in a cross-validation exercise using universal (or external drift) kriging (Reference CressieCressie, 1993), with the pooled variogram and altitude, solar exposure and the vegetation indices as explanatory variables for the mean function. In such a leave-one-out scheme we eliminateda given datum from the dataset and predicted its snow depth from the remaining data. Then we compared the predicted and the measured values and computed their difference (prediction error). We repeated this procedure for all the observations of the set, and computed various statistics (coefficient of determination, mean squared error of prediction, mean absolute error of prediction, etc.). We used only the observations that had been recorded on the same date as the target datum when computing the cross-validation predictions. Surprisingly, when we used the variogram estimated from the pooled dataset, the universal kriging predictions were no better than the predictions obtained from the linear regression model, i.e. when the correlation of the error term was ignored. The coefficients of determination were: 0.687 (universal kriging) and 0.693 (linear regression), and the respective mean square errors of prediction were 3.48 and 3.39, respectively. Kriging was marginally better with respect to the mean absolute prediction errors, but the differences were negligible. The results of the cross-validation exercise are in agreement with the lack of spatial dependence exhibited by the sample variograms of individual dates. It therefore appears that the spatial correlation signalled by the pooled variogram is an artifact mainly caused by the superposition of unrelated patterns.

Since the cross-validation analysis confirmed that there was no spatial dependence, we used the coefficients determined by the regression analysis to map the snow depth in the Erlenbach catchment. Figure 5 shows the spatial distribution on 23 December 1999, a precipitation-free day with a mean air temperature close to 0°C. On that day, the snow depth varied in the catchment between 20 and 110 cm. The model predicted increasing snow depth with increasing altitude, but at the smaller scale the structure of the prediction surface was largely determined by the vegetation. In particular, east of 697 the snow depth is closely related to the forest coverage (cf. Fig. 3b). Avery deep snow cover was predicted around location (698.5, 211). Inspection of the aerial photograph revealed a patch of bare ground at this location which obviously influenced the spectral signature recorded by Landsat TM.

Fig. 5. Spatial distribution of the snow depth on 23 December 1999, predicted from the regression model Equation (2) fitted to the full dataset.

Discussion

One main conclusion resulting from this study is that indices of the forest density derived from the Landsat TM image were poor explanators of snow-depth variation measured within the Erlenbach catchment. On the other hand, local estimates of the canopy closure from photographs taken from below the canopy significantly improved our spatial snow-depth simulation. The key problem related to this result already showed up in the initial analysis of the forest structure: measurements and image information are recorded at partly contrasting scales. Whereas the snow-depth measurements, aerial photographs and canopy photos basically reproduce the pattern of single trees and glades, the Landsat TM vegetation indices may only capture the general structure of the forest mosaic. Altitude and solar exposure were determined both locally and from a 25 m DEM, which gave very similar results. In the regression analysis we could not strictly adhere to one scale, which, at least partly, explains the modest R₂ values of our regression model Equation (2) for the snow-depth variationin the catchment for single days. Using canopy density information (f_c) recorded at the same scale as the snow-depth measurements resulted in a much better fit of the model. An excellent discussion of such scale issues is provided in Reference BlöschlBlöschl (1999).

The range of goodness of fit obtained with the model that includes f_c is in agreement with the R² values of Reference ForsytheForsythe (1999) who applied a similar multiple regression model to snow-course measurements in the Kalkhochalpen region. As in our study, he found that the goodness of fit of such a model varied significantly within the season and between winters.

In spite of the moderate success in the present study, we are still convinced that one has to take into account vegetation indices derived from satellite images when modelling the spatial distribution of the snow cover. The main limitation is the spatial resolution of the remotely sensed spectral data from which we have to derive the vegetation indices.

A final conclusion from this study is that geostatistical approaches do not improve the prediction of the snow-depth variation in this type of landscape, so we have to be content with applying the multiple linear regression model.