GIS and field data based modelling of snow water equivalent in shrub tundra

An approach for snow water equivalent (SWE) modelling in tundra environments has been developed for the test area on the Yamal peninsula. Detailed mapping of snow cover is very important for tundra areas under continuous permafrost conditions, because the snow cover affects the active layer thickness (ALT) and the ground temperature, acting as a heat-insulating agent. The information concerning snow cover with specific regime of accumulation can support studies of ground temperature distribution and other permafrost related aspects. Special attention has been given to the presence of shrubs and microtopography, specifically ravines in a modelling approach. The methodology is based on statistical analysis of snow survey data and on GIS (Geographical Information System) analysis of a range of parameters: topography, wind, and shrub vegetation. The topography significantly controls snow cover redistribution. This influence can be expressed as increase of snow depth on concave and decrease on convex surfaces. Specifically, snow depth was related to curvature in the study area with a correlation of R=0.83. An index is used to distinguish windward and leeward slopes in order to explain wind redistribution of snow. It is calculated from aspect data retrieved from a digital elevation model (obtained by field survey). It can be shown that shrub vegetation can serve as a ‘trap’ for wind-blown snow but is not a limiting factor for maximum snow depth, since the snow depth can be higher or lower than shrub height dependent on other factors.


Introduction
The systematic observations of snow cover in Russia have been initiated by A.Voeykov, who noted that snow was an important component of the environment (Voeykov 1889).In permafrost regions, snow cover plays an important role because it di-rectly influences the thermal regime of frozen ground, as it is a natural thermal insulator (Dostovalov & Kudryavtsev 1967).Snow water equivalent (SWE) can be used as a proxy for the integral characteristics of snow depth and density, which define the insulating properties of snow.SWE allows the calculation of the land surface snow stor-FENNIA 193: 1 (2015) Yury Dvornikov et al. age and is also used for the assessment of the water regime of rivers and lakes and the activation of various cryogenic processes (Kuz'min 1960;Gray & Male 1981;Marsh et al. 1995).In Arctic tundra strong winds also play a significant role in snow drift and can be responsible for the increase in snow density (Kuz'min 1957).Zhitkov (1913), based on his field observations in Yamal, Western Siberia, Russia, already noted that on the flat surfaces of tundra, unrelated to hilltops, snow accumulates up to 20-30 cm depth.However, it is documented in Trofimov (1975) that snow with a depth up to 3-4 meter can accumulate in depressions, while in most areas snow depth is usually less than 15 cm.This small depth is caused by low winter precipitation in the Arctic (Kuz'min 1960).However, there are recent observations, which show the increase of maximal snow storage in Russia, including the northern part of Western Siberia (Krenke et al. 2000;Kitaev & Kislov 2008).
According to the "Map of snow depth" (Richter 1948), mean values of snow depth do not exceed 30-50 cm for the Yamal peninsula.These data are consistent with those published by Kopanev (1978), who noted that the average perennial snow depth for Yamal does not exceed 30 cm in the second half of March.More recent studies indicate that the maximal snow storage for Yamal is about 150 mm SWE (Kotlyakov 2004) or 50 cm of snow with the average density of 0.3 g/cm 3 (Trofimov 1975).
Previous studies on snow depth and SWE modelling at the local scale have been carried out for regions with different environmental conditions: alpine (Purves et al. 1998;Winstral et al. 2002;Geddes et al. 2005;Clow et al. 2012), and arctic tundra (König & Sturm 1998;Essery & Pomeroy 2004).Practically the same parameters (independent variables) have been used in these studies for modelling, such as topography, wind, vegetation, differing only in the importance assigned to each variable.
The main purpose of this study is to model the snow accumulation processes for tundra landscape (elevation range 56 meters) where woody vegetation is absent.Three main parameters are considered: topography, wind, vegetation.The degree of each control's weight in the process of snow redistribution is discussed as well.Additionally the use of GIS enables the spatial and landscape consideration of all the parameters.Analysis is characterized by a high level of detail corresponding to micro-level research (Gray & Male 1981).

Study area
The "Vaskiny Dachi" research station is located in Central Yamal (70º20'N, 68º51'E) to the southeast of the Bovanenkovo gas field (Leibman & Kizyakov 2007).The topography is represented by stepped plain, dissected by ravines, lake basins, small rivers and complicated by a complex of cryogenic processes, mainly cryogenic landsliding (Leibman et al. 2012).The field observations of snow cover in this area in the beginning of the 20th century have shown that the strong dissection has some influence on the distribution of snow depth: the deflation from hilltops and accumulation in depressions, such as the ravines and foot of the slopes (Zhitkov 1913), under the influence of strong winds is a very significant component of the study area (Trofimov 1975).According to the Marre-Sale weather station records (www.rp5.ru), prevailing wind direction for the cold season of 2012-2013 was from south-east.During our snow survey in late March 2013 the same prevailing ESE wind direction was observed with a maximum wind speed of 14 m/sec, and an average wind speed of 7.4 m/sec.The total sum of precipitation during the cold season up to the end of the snow survey amounted to 118 mm SWE.
The study area is located in the bioclimatic subzone D (CAVM Team 2003), which is characterized by about 9 °C mean July temperature.Zonal vegetation on gently sloping upland surfaces consists of sedges, prostrate and erect dwarf (<40 cm tall) shrubs and mosses (http://www.a r c t i c a t l a s .o r g / m a p s / t h e m e s / c p / cpbzUnit?queryID=D).
Shrub willows (Salix glauca, S. lanata) and dwarf-birch (Betula nana) are widely distributed here (Rebristaya & Khitun 1998).Plant communities with dense shrub layer mostly occupy valley bottoms and gentle hill slopes.Willows grow up to 2 m high in some places, so the structure of plant communities may affect snow distribution significantly.
Modelling was carried out for the local area represented by a 1.65 km long and 250 m wide transect (Fig. 1).The transect crosses the main geomorphologic units of the study area ("Vaskiny Dachi" research station) and is characterized by different surface properties and vegetation complexes.The transect includes a site for active layer depth monitoring (CALM) with an extent of 100x100 m (Fig. 1, CALM site) (Brown et al. 2000, http://www.gwu.edu/~calm/).

Methodology Snow survey
The field snow survey in the study area was undertaken between 16th-31st March 2013.Snow depth was measured using a metal ruler with 1 mm interval and snow density with the snow sampler VS-43.This device is mostly used in Russian meteorological stations and designed for snow density measurements during snow surveys (Slabovich 1985).It consists of a metal cylinder and scales.At one end of the cylinder is a ring with cutting elements while the other end is closed by a lid.
Snow depth was measured at 233 points, including 121 points on the CALM site, and density was measured at 55 points.

Digital elevation model
The application of a GIS allows the spatial analyses of snow depth (Evans et al. 1989;Purves et al. 1998).Since topography significantly controls the redistribution of snow in Arctic landscapes, a detailed DEM is of high importance for SWE modelling (Litaor et al. 2008).
In summer 2011, a topographic survey was carried out using a tachymeter TopCon GTS-235 (accuracy of angle measurements 5') to produce a detailed DEM with 5x5 m cell size.Based on this raster dataset, terrain derivatives were calculated.Slope aspects within 0-360 degree range were grouped according to the main directions: N-NE-E-SE-S-SW-W-NW, as well as a curvature ranging from -4 to +4.Such a range of curvature corresponds to the dissected topography of the site (Marchand & Killingtveit 2001).Curvature shows the convexity and concavity degree of topographic patterns, i.e. of one raster cell relative to other eight surrounding raster cells (Zeverbergen & Thorne 1987).Negative curvature values correspond to concave and positive to convex patterns.The significance of this index in the studies of snow cover has been noted by many researchers (e.g.Freindlin & Shnyparkov 1985;Golding 1974;Woo et al. 1983;Sexstone & Fassnacht 2014), from which it has been observed that greater snow accumulation occurs on concave slopes rather than on convex hilltops.
A correlation analysis was conducted to reveal the relationship, if any, of snow depth with the above topographic parameters.For each point with measured snow depth, and coverage by the detailed DEM, aspect and curvature values were extracted according to their spatial referencing.

Wind-induced redistribution of snow
It is well documented that snow is blown away from windward slopes and accumulated on leeward slopes (Evans et al. 1989;Winstral et al. 2002;Litaor et al. 2008).The prevailing wind direction is also considered in our model.It is selected as southeastern according to the Marre-Sale weather station records.To consider the influence of the wind, an empirical correction W, which depends on the aspect value of each raster cell, is introduced (Eq.1): where W is a correction for the initial value of simulated snow depth, A the slope aspect, and K a coefficient which is calculated empirically and depends on the amount of snow precipitation as well as on wind speed.

Vegetation-induced snow redistribution
Shrub vegetation on the key site is an important factor for the redistribution of snow masses, being a 'trap' for snow drift as well as a factor decreasing the wind activity (Benson & Sturm 1993) for both windward and leeward slopes (Essery & Pomeroy 2004).However, it is also documented that these trapping processes do not always take place as snow might be partially blown away from shrub vegetation complexes (Pomeroy & Gray 1995)

Snow water equivalent model
A methodology was developed to build a map of SWE, including the consideration of parameters: topography, wind, and vegetation (independent variables).The general scheme of GIS-based modelling is shown in Figure 2. The DEM was converted to a point vector data model with 5 m spacing.
Curvature and aspect values for each point have been added to the attribute table.A shrub height measure was assigned based on the field descriptions to shrub vector polygons retrieved from satellite imagery interpretation.The initial value of snow depth was calculated based on the established dependence between the field measured snow depth and extracted curvature values.To obtain final values of calculated snow depth, corrections based on the wind and vegetation data have been made.The field measurements of snow density were used in transition from snow depth values to SWE.If the density was measured layer-by-layer, the average density for the entire section was used in the model.Comparison of two parameters (snow depth and SWE) showed a linear de-pendence (R2 = 0.99), which has been used in the calculation of SWE from snow depth.After all, the raster surface was interpolated.

Snow survey
The measured snow depth varies, depending on the type of terrain, from 0 to 315 cm.The average snow depth ranges on sub-horizontal surfaces between 15 and 30 cm.In depressions, snow depth exceeds 1 m and can reach several meters.Density also varies depending on the geomorphology: from 0.17 g/cm 3 on the flat hill tops to 0.67 g/cm 3 on concave portions of the slopes.The total results of the snow survey are summarized in Table 1.
The data presented in Table 1 show that snow depth is a very irregular parameter for our site mainly due to relief control.In general, the results of snow depth and snow density measurements match that of previous research (Richter 1948;Trofimov 1975;Kopanev 1978;Kotlyakov 2004).An anomalously high maximum value of snow density measured in the field on the hilltop surface (0.67 g/cm 3 ) points to a strong wind influence that has led to snow cover compression.

Topographic controls
The relation of measured snow depth to curvature appears to be very strong (Fig. 3), and hence the topography has been selected as a main factor for calculation of the initial snow depth value for the transect.
The curvature has been used before as a secondary independent variable in snow modelling (Gold-  FENNIA 193: 1 (2015) Yury Dvornikov et al. ing 1974;Woo et al. 1983;Winstral et al. 2002, Clow et al. 2012).A study of the dependence between snow depth and curvature for the area of Norway showed a very low correlation (Marchand & Killingtveit 2001).The authors analyzed a linear model of the dependence and used a low resolution 100x100 meters DEM, which -as they mentioned -can adversely affect the degree of correlation.In contrast, studies of snow cover in Khibiny mountains (Freindlin & Shnyparkov 1985;Kontsevaya et al. 1989) show the high correlation between these two parameters and the authors used a linear regression model for snow depth calculation.Kasurak et al. (2009) and Sexstone and Fassnacht (2014) also confirm that curvature is well correlated with snowpack properties even for alpine terrain.
In this paper, this variable is used as a main parameter.It appears to be obvious that strong inverse correlation (-0.83) between this parameter and snow depth was revealed only because the detail of the used DEM (5x5 m, which corresponds to 1:1000 mapping scale) allows an accurate description of topography and therefore a more detailed derivative surface.The DEM with lower resolution used for the same area results in the substantially lower correlation between these parameters (Fig. 4).Curvature is very sensitive to the source DEM resolution because the depth values of the eight surrounding cells are taken into account in the calculation procedure for each raster cell.It was also observed in the field that the influence of microtopography is high too.
The measured snow depth varied up to 10-20 cm in an area of less than 4-10 m 2 .However, the acquisition of a DEM with such a detail is a very timeconsuming process.Processed stereo-pairs of airborne and satellite images or LIDAR data may serve as an additional source along with a field topographic survey.

Modelling the wind-induced redistribution of snow
Equation (1) allows to index leeward slopes positively and windward slopes negatively relative to the prevailing wind direction.In Figure 5,   1) were established according to the mean measured snow depth analysis on the slopes of different aspect (Fig. 6).
Inventories of wind impact on snow redistribution using GIS have been carried out by other researchers (Purves et al. 1998;Winstral et al. 2002;Clow et al. 2012) with a similar scheme: 1) the identification of a prevailing wind direction; and 2) indexing a raster surface according to the main established pattern; blowing of snow away from windward slopes and accumulation on leeward slopes.
Following our approach for a more accurate analysis of snow redistribution by wind, the field measurements were subdivided into two groups: 1) collected within the sites with positive curvature values (convex sites), and 2) collected within sites with negative curvature values (concave sites).The data shown in Figure 6 demonstrate that, on average, accumulation is up to 20-30 cm more on the leeward slopes than on windward slopes from which snow is blown away.In order to be consistent with these data, the coefficient K for the correction of W was assumed to be 25.This value allows for the expansion of the value of corrections up to ± 20 cm.This wind correction was entered in the initial calculated snow depth value according to the graph in the attributive table with aspect values.So far, the dependence between snow depth and curvaturewind group is not combined into one index because we are yet unable to account for the prevailing wind direction variations from year to year.It's the main reason why the initial value and wind correction are used separately and this combination remains to be solved.Such combinations have been successfully applied using binary regression trees (Elder et al. 1995;Erxleben et al. 2002;Winstral et al. 2002;Molotch et al. 2005).

Modelling the vegetation-induced redistribution of snow
Shrub willow distribution highly depends on relief.The tallest shrubs (1-1.5 m) occupy concavities: the valleys of small rivers and slopes of erosion features (Ukraintseva 1997) with the negative values of curvature.Research carried out in the tundra of Alaska also showed an increase of snow density in relation to canopy height, branch diameter (McFadden et al. 2001), as well as shrub density (Sturm et al. 2001).According to the data, shrub communities accumulate 27% more snow than hummocky tundra.
The relationship between snow depth, shrub presence and shrub height was similar to the one with the wind: the measurements were compared depending on surface convexity/concavity. Figure 7 shows the correlation between shrub height and snow depth versus the convexity/concavity index.The data were graded according to this index to study the impact of shrubs on snow specifically at convex sites, as at concave areas this relationship is overridden completely by relief as shown in Fig. 7.In the first case, strong dependence of snow depth on shrub height is shown, in the second case this dependence is not observed.This is due to the biologically caused impossibility of shrub growth over 1.5 m in sites where snow depth exceeds this value (Leibman 2004).
Figure 8 shows sorted values from the dataset for convex slopes.This confirms that the shrub vegetation can serve as a 'trap' for snow (Essery & Pomeroy 2004).Following from this, sites with low shrubs (less than 15 cm) control the snow depth, but for higher shrubs (more than 15 cm) only a part of blown snow is trapped and so they do not control the snow depth.It could be seen through the behavior of curves on Figure 8: at lower shrub height they correlate, but above 15 cm the "shrub height" curve starts to be markedly higher than the "snow depth" curve, which generally is going up simultaneously.This threshold possibly depends on the total amount of fallen snow.
According to the regression equation (Fig. 7), the limitation value for calculated snow depth with the presence of shrubs with certain height is (Eq.2): (2) where X is snow depth, and Y is shrub height.The value for each modelling point, which are defined by the shrub vector polygons, was compared to the initial SWE value and corrected for the wind redistribution impact.When the limitation value was higher, then the calculated value was assumed equal to the limitation value.In this case, the shrub vegetation was accepted as a limit for snow accumulation according to the field data.Results of this study provide the evidence that shrub height can affect the snow depth, but does not fully control this parameter.The relationship is better understood when the topography factor is excluded.The correlation between these parameters is stronger in this case (e.g.Fig. 7).
The assumption that the vegetation acts as a secondary factor contradicts data for other regions, where shrubs -and not topography -have been established as a main control for snow depth (e.g.McFadden et al. 2001).For the area, which is densely dissected by erosion processes, topography plays a more important role for the redistribution of snow.Our analysis shows a linear relationship of the 'trap' effect of shrubs with snow depth.

Creation of SWE map
The SWE map is shown in Figure 9.The resolution of the final raster model is equal to the source DEM with 5x5 m.
Modelling and mapping of snow cover in this paper is done at the local scale.Therefore, results cannot well describe the distribution of snow within different areas, such as neighboring floodplains.Similar conclusions were suggested by other researchers, who dealt with small test sites (Litaor et al. 2008).Environmental features of the area specify the number of independent variables describing the distribution of snow in tundra more accurately (Evans et al. 1989;McFadden et al. 2001;Sturm et al. 2001;Essery & Pomeroy 2004).Even though independent variables, which control SWE, do not correlate linearly with snow distribution (Elder et al. 1995;Molotch et al. 2005), linear regression models have already been used (Golding 1974).Generally, it is obvious that all the factors influencing the redistribution of snow in the site are interdependent and complex.

Validation of results
The reliability of the model for snow depth distribution was estimated through comparison with the model based exclusively on field measurements.For this purpose, the results of snow depth measurements from the entire transect including     CALM site (see Fig. 1) were used, as well as a detailed DEM based on the ground survey.Field data from the CALM site deliberately were not used for the statistical analysis of the snow depth controls.Snow depth was modelled for each point on the transect.The analysis showed that there is a quite high correlation factor between the two value arrays (R=0.82).Field and modelled datasets for CALM site are shown in Figure 10: the volume of snow was calculated for both the in situ and the modelled version (Table 2).This shows that the values of modeled snow depth and accordingly the values of SWE reflect the real snow cover distribution well.

Conclusions
Our observations of the snow cover in the study area showed that the topography of the Central Yamal tundra caused the extremely uneven distribution of snow.They also confirmed the importance of wind drifting from the hilltops and accumulation of snow in the concavities, as observed by other researchers.
Snow redistribution is due to strong winds and the lack of significant natural barriers.The data analysis confirmed a pattern specified by many authors previously: blowing away of snow from the windward slopes and accumulation on the leeward slopes in accordance with the prevailing wind direction.To account for this pattern in the general model proposed here we introduce a correction based on the slope aspect.This correction allows for an indexing of the slopes with respect to the prevailing wind direction and facilitates automation.
The impact of shrubs is estimated as 'traps' for the snow or limiting control, depending on the combination of shrub height and surface shape.The relationship between shrub height and snow depth is well expressed when concave landforms do not prejudge the accumulation of snow exceeding the height of shrubs.As high shrubs in the key area are linked to depressions, the limited impact of shrub communities on the depth of the snow is observed.However, shrubs retain part of the snow cover, providing a limiting effect which is used in the model.
Based on field data and identified relationships between various controls and snow cover depth, a technique for GIS-based modelling of the snow water equivalent was developed, which is applicable to a particular type of landscape.These studies can be applicable in the detailed permafrost mapping and modeling of the number of parameters, such as active layer or ground temperature on the local scale, as well as for the hydrological studies of the tundra environments with the link to the various cryogenic processes' activation.

Fig. 2 .
Fig. 2. Workflow for the calculation of Snow Water Equivalent (SWE) using GIS and remote sensing (RS) data.

Fig. 3 .
Fig. 3.The relation of field measured snow depth (cm) to surface curvature.

Fig. 4 .
Fig. 4. The relation of field measured snow depth (cm) to surface curvature derived from a 25 meters resolution DEM based on 1:25,000 topographic map.

Fig. 5 .
Fig. 5. Correction values to the modelled snow depth according various slope aspects in relation to prevailing wind direction (see Eq. 1).

Fig. 7 .
Fig. 7.The relation between snow depth and shrub height.The dataset is subdivided into two groups: points located on convex places (curvature > 0), and points located on concave places (curvature < 0).

Fig. 8 .
Fig. 8.The dependence between snow depth and shrub height (samples located on convex places).

Fig. 9 .
Fig. 9. Snow water equivalent map for central part of the transect.

Fig. 10 .
Fig. 10.Snow depth on the CALM site (a -interpolation of field data, b -calculated values).

Table 2 .
The comparison of field and calculated snow depths for the transect.