Small-scale variation of snow in a regional permafrost model

The strong winds prevalent in high altitude and arctic environments heavily redistribute the snow cover, causing a small-scale pattern of highly variable snow depths. This has profound implications for the ground thermal regime, resulting in highly variable near-surface ground temperatures on the metre scale. Due to asymmetric snow distributions combined with the nonlinear insulating effect of snow, the spatial average ground temperature in a 1 km2 area cannot be determined based on the average snow cover for that area. Land surface or permafrost models employing a coarsely classified average snow depth will therefore not yield a realistic representation of ground temperatures. In this study we employ statistically derived snow distributions within 1 km2 grid cells as input to a regional permafrost model in order to represent sub-grid variability of ground temperatures. This improves the representation of both the average and the total range of ground temperatures. The model reproduces observed sub-grid ground temperature variations of up to 6 C, and 98 % of borehole observations match the modelled temperature range. The mean modelled temperature of the grid cell reproduces the observations with an accuracy of 1.5 C or better. The observed sub-grid variations in ground surface temperatures from two field sites are very well reproduced, with estimated fractions of sub-zero mean annual ground surface temperatures within±10 %. We also find that snow distributions within areas of 1 km2 in Norwegian mountain environments are closer to a gamma than to a lognormal theoretical distribution. The modelled permafrost distribution seems to be more sensitive to the choice of distribution function than to the fine-tuning of the coefficient of variation. When incorporating the small-scale variation of snow, the modelled total permafrost area of mainland Norway is nearly twice as large compared to the area obtained with grid-cell average snow depths without a sub-grid approach.


Introduction
High altitude and arctic environments are exposed to strong winds and drifting snow 20 can create a small-scale pattern of highly variable snow depths. Seasonal snow cover is a crucial factor for the ground thermal regime in these areas (e.g. Goodrich, 1982;Zhang et al., 2001). This small-scale pattern of varying snow depths results in highly variable ground temperatures on the meter scale; up to 6 • C within areas of less than Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a large number of snow surveys in the Northern Hemisphere shows a large spread of CV sd values, in particular within this land use class, ranging from 0.1 to 0.9 (Clark et al., 2011). This illustrates the need for improved representation of snow distribution within this land use class. An accurate representation of the small scale snow variation highly influences the 5 timing and magnitude of runoff in hydrological models, and a detailed picture of the sub-grid variability is of great value for the hydropower industry and in flood forecasting. Adequate representations of the snow covered fraction in land surface schemes are important for enhanced realism of simulated near surface air temperatures, ground temperatures and evaporation due to the considerable influence of snow cover on the 10 duration of melt season and the surface albedo.
In this study we derive functional dependencies between distributions of snow depth within 1 km × 1 km grid cells and CV sd , based on an extensive in-situ data set from Norwegian alpine areas. In a second step, we employ the resulting snow distributions as input to the permafrost model CryoGRID1, a spatially distributed, equilibrium per- 15 mafrost model . From the sub-grid representation of ground temperatures, permafrost probabilities are derived, hence enabling a more realistic, fuzzy permafrost boundary instead of a binary, sharp transition. With this approach, we aim to improve permafrost distribution modelling in inhomogeneous terrains.
2 Setting 20 The model is implemented for the Norwegian mainland, extending from 58 to 71 • N. Both the topography and climate in Norway is dominated by the Scandes, the mountain range stretching south-north through Norway, separating the coastal western part with steep mountains and deep fjords from the eastern part where the mountains gradually decrease in height. The maritime climate of the west coast is dominated by low- 25 pressure systems from the Atlantic Ocean resulting in heavy precipitation, while the eastern parts of the Scandes have a more continental and drier climate. Mountain permafrost is present all the way to the southern parts of the Scandes, with a gradient in the lower limit of permafrost from ∼ 1400 to 1700 m from east to west in central southern Norway, and from ∼ 700 to 1200 m from east to west in northern Norway . While permafrost is also found in mires at lower elevations both in southern and northern Norway, most of the permafrost is located in exposed terrain 5 above the tree line. This environment is dominated by strong winds resulting in heavy redistribution of snow.
In-situ records of snow depth data used to establish the snow distribution scheme were collected at the Hardangervidda mountain plateau in the southern part of the Scandes (Fig. 1). It is the largest mountain plateau in northern Europe, located at ele- 10 vations from 1000 to above 1700 m a.s.l., with occurrences of permafrost in the highest mountain peaks. The terrain is open and slightly undulating in the east, while in the west it is more complex with steep mountains divided by valleys and fjords. The mountain range represents a significant orographic barrier for the prevailing westerly winds from the Atlantic Ocean, giving rise to large variations in precipitation and strong winds, two 15 agents promoting a considerably wind-affected snow distribution. Mean annual precipitation varies from 500 to more than 3000 mm over distances of a few tens of kilometres, and maximum snow depths can vary from zero to more than 10 m over short distances (Melvold and Skaugen, 2013). 20

A statistical model for snow depth variation
The Winstral terrain-based approach (Winstral et al., 2002)  The terrain-based exposure parameter (Sx), described in detail in Winstral et al. (2002), quantifies the extent of shelter or exposure of the grid-cell considered. Sx is determined by the slope between the grid-cell and the cells of greatest upward slope in the upwind terrain. The upwind terrain is defined as a sector towards the prevailing wind direction d constrained by the maximum search distance (d max = 100 m) 5 and a chosen width (A) of 30 • with the two azimuths extending 15 • to each side of d (see Fig. 2). The cell of the maximum upward slope is identified for each search vector, separated by 5 • increments. Sx for the given grid-cell is finally calculated as the average of the maximum upward slope gradient of all seven search vectors: where d is the prevailing wind direction, (x i , y i ) are the coordinates of the considered grid-cell, and (x v , y v ) are the sets of all cell coordinates located along the search vector defined by (x i , y i ), A and d max . This gives the degree of exposure or shelter in the range −1 to 1, where negative values indicate exposure.
To estimate a realistic degree of exposure based on the observed wind pattern 15 at a local site, Sx was computed for each of the eight prevailing wind directions d = [0, 45, 90, 135, 180, 225, 270, 315 • ], and weighted based on the wind fraction (wf d ).
wf d accounts for the amount of different exposures in the terrain at various wind directions, and represents the fraction of hourly wind direction observations over the accumulation season for the eight wind directions. The accumulation season is here chosen 20 as January to March. Wind speeds below a threshold of 7 ms −1 are excluded, as this threshold is considered a lower limit required for wind drifting of dry snow (Lehning and Fierz, 2008;Li and Pomeroy, 1997). The calculated Sx parameter values are used as predictors in different regression analyses to describe the CV sd within 1 km × 1 km derived from the ALS. The coefficient 25 of variation of exposure degrees (CV Sx ) within each 1 km × 1 km grid cell is computed by aggregating the Sx map from 10 m to 1 km resolution according to: Introduction Sx values below the 2.5th and above 97.5th percentiles of the Sx distributions are excluded, giving Sx ≈ [−0.2, 0.2]. Three regression analyses were performed to reduce the RMSE between CV Sx and observed CV sd , where additional predictors such as elevation above treeline (z) and maximum snow depth (µ) successively have been in-5 cluded (Table 1). Elevation above treeline is chosen as predictor to account for the increased wind exposure with elevation. Ideally, wind speed should be included as predictor. However, the NORA10 dataset (Sect. 4.1) does not sufficiently reproduce the local variations in wind speeds over land, especially not at higher elevations and for terrain with increased roughness. Because of the strong gradient in treeline and gen-10 eral elevation of mountain peaks from high mountains in the south to lower topography in the north of Norway, applying only elevation as predictor would result in an underestimation of redistribution in the north.

CryoGRID 1 with an integrated sub-grid scheme for snow variation
The equilibrium permafrost model CryoGRID 1 Westermann et al., 15 2015) provides an estimate for the MAGST (Mean annual ground surface temperature) and MAGT (Mean Annual Ground Temperature) from freezing (FDD a ) and thawing (TDD a ) degree days in the air according to and 20 MAGT = where P is the period that FDD a and TDD a are integrated over, r k is the ratio of thermal conductivities of the ground in thawed and frozen states, while nT and nF are semiempirical transfer-functions including a variety of processes in one single variable. The winter nF factor relates the freezing degree days at the surface to the air and thus accounts for the effect of the winter snow cover, and likewise the nT factor relates the 5 thawing degree days at the surface to the air and accounts for the surface vegetation cover: FDD s = nF · FDD a and TDD s = nT · TDD a . nT = −0.13 · µ + 1.1. We assume that the distribution of maximum snow depths within a grid cell with a given CV sd and average maximum snow depth (µ) follows a gamma distribution with a probability density function (PDF) given by: with a shape parameter α = CV −2 sd and a rate parameter β = µ×CV 2 sd (e.g. Kolberg and 5 Gottschalk, 2006;Skaugen et al., 2004). Corresponding n factors are computed for all snow depths (x) based on Eqs. (6) and (7), and related to the PDF (Eq. 8). The model is run for each nF from 0 to 1 with 0.01 spacing, giving 100 model realizations.
Each realization corresponds to a unique snow depth, represented with a set of nF and nT factors. Based on the 100 realizations a distribution of MAGST and MAGT are calculated for each grid cell, where the potential permafrost fraction is derived as the percentage of sub-zero MAGT. To assess the sensitivity of the choice of the theoretical distribution function, the model was also run with PDFs following a lognormal distribution, given by (e.g. Liston, 2004):

Model evaluation
The CV sd was derived for 0.5 km × 1 km areas based on the ALS snow depth data (Sect. 4.1) resampled to 10 m × 10 m resolution. Each 0.5 km × 1 km area includes 500 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of fit evaluations for the theoretical lognormal and gamma distributions applying the Anderson-Darling test in MATLB (adtest.m, Stephens, 1974) were conducted for each distribution. Parameters for gamma (shape and rate) and lognormal (mu, sigma) distributions were estimated by maximum likelihood as implemented in the MATLAB functions gamfit.m and lognfit.m.

5
The results of the permafrost model are evaluated with respect to the average MAGST and MAGT within each grid cell, as well as the fraction of sub-zero MAGST. The model runs are forced with climatic data for the hydrological year corresponding to the observations. The performance in representing fractional permafrost distribution is evaluated at two field sites where arrays of 26 (Juvvasshøe) and 41 (Finse) 10 data loggers have measured the distribution of ground surface temperatures within 500 m × 500 m areas for the hydrological year 2013 (Gisnås et al., 2014). The general lower limits of permafrost are compared to permafrost probabilities derived from BTS (basal temperature of snow) -surveys (Haeberli, 1973;Lewkowicz and Ednie, 2004), conducted at Juvvasshøe and Dovrefjell (Isaksen et al., 2002). The model performance 15 of MAGST is evaluated with data from 128 temperature data loggers located a few cm below the ground surface in the period 1999-2009 (Farbrot et al., , 2011Isaksen et al., 2008Isaksen et al., , 2011Ødegaard et al., 2008). The loggers represent all vegetation classes used in the model, and spatially large parts of Norway (Fig. 2). Four years of data from 25 boreholes (Isaksen et al., 2007(Isaksen et al., , 2011Farbrot et al., 2011Farbrot et al., , 2013 are used 20 to evaluate modelled MAGT (Fig. 2).

Forcing and evaluation of the snow distribution scheme
Wind speeds and directions during the snow accumulation season are calculated from the boundary layer wind speed and direction at 10 m in the Norwegian Reanalysis TCD 9, 2015 Small-scale variation of snow in a regional permafrost model ERA-40 to a spatial resolution of 10-11 km, with hourly resolution of wind speed and direction (Reistad et al., 2011). The dataset is originally produced for wind fields over sea, and underestimates the wind speeds at higher elevation over land (Haakenstad et al., 2012). Comparison with weather station data revealed that wind speeds above the tree line are underestimated by about 60 % (Haakenstad et al., 2012). For these 5 areas the forcing dataset has been scaled accordingly. The snow distribution scheme is derived from an Airborne Laser Scanning (ALS) snow depth over the Hardangervidda mountain plateau in southern Norway (Melvold and Skaugen, 2013). The ALS scan is made along six transects, each covering a 0.5 km × 80 km area. The survey was first conducted between 3 and 21 April 2008, 10 and repeated in the period 21-24 April 2009. The snow cover was at a maximum during both surveys. A baseline scan was performed 21 September 2008 to obtain the elevation at minimum snow cover. The ASL data are presented in detail in Melvold and Skaugen (2013). Distributions of snow depth, represented as CV sd , are calculated for each 0.5 km×1 km area, based on the snow depth data resampled to 10 m×10 m reso-15 lution. About 400 cells of 0.5 km×1 km exist for each year, when lakes and areas below treeline are excluded.
The snow distribution scheme is validated with snow depth data obtained by ground penetrating radar (GPR) at Finse ( in open, non-vegetated alpine landscapes with major wind re-distribution of snow. However, they differ with respect to elevation (1300/1450 m a.s.l.), mean maximum snow depth (∼ 2 m/∼ 1 m), average winter wind speeds (7-8/10-14 m s −1 ) and topography (very rugged at Finse, while steep, but less rugged at Juvvasshøe). The timing of the snow surveys were late March to April (2009April ( , 2012April ( -2014 around maximum snow 25 depth, but when the snow pack was still dry. The GPR surveys at Finse are constrained to an area of 1 km × 1 km, while at Juvvasshøe they cover several square kilometres, but with lower observation density. The GPR data from the end of the accumulation season in 2013 are presented in Gisnås et al. (2014), and the data series from the other years are obtained and processed following the same procedures, described in detail in Dunse et al. (2009). The propagation speed of the radar signal in dry snow was derived from the permittivity and the speed of light in vacuum, with the permittivity obtained from snow density using an empirical relation (Kovacs et al., 1995). The snow depths were determined from the two-way travel time of the reflection from the 5 ground surface and the wave-speed. Observations were averaged over 10 m×10 m grid cells, where grid cells containing less than three samples were excluded. The CV sd for 1 km × 1 km areas are computed based on the 10 m resolution data.

Permafrost model setup
The climatic forcing of the permafrost model is daily gridded air temperature and snow 10 depth data for the period 1961-2013, called the seNorge dataset, provided by the Norwegian Meteorological Institute (Mohr and Tveito, 2008;Mohr, 2009) and the Norwegian Water and Energy Directorate (Engeset et al., 2004;Saloranta, 2012). The dataset is based on air temperature and precipitation data collected at the official meteorological stations in Norway, interpolated to 1 km × 1 km resolution. Snow depths are derived 15 from the air temperature and precipitation data, using a snow algorithm accounting for snow accumulation and melt, temperature during snow fall and compaction. Freezing (FDD a ) and thawing (TDD a ) degree days in the air are calculated as annual accumulated negative (FDD) and positive (TDD) daily mean air temperatures, and maximum annual snow depths (µ) are derived directly from the daily gridded snow depth data. 20 The CryoGRID 1 model is implemented at 1 km × 1 km resolution over the same grid as the seNorge dataset. Soil properties and surface cover is kept as in Gisnås et al. (2013)

Observed snow distributions in mountain areas of Norway
CV sd within 1 km × 1 km areas in the ALS snow survey at Hardangervidda ranged from 0.15 to 1.14, with mean and median of respectively 0.58 and 0.59. According to the Anderson-Darling goodness of fit evaluations 70 out of 932 areas had a snow distri-5 bution within the 5 % significance interval of a gamma distribution, while only 1 area was within the 5 % significance interval of a lognormal distribution. Although the null hypothesis rejected more than 90 % of the sample distributions, the Anderson-Darling Test Score was all over lower for the gamma distribution, indicating that the observed snow distributions are closer to a gamma than to a lognormal theoretical distribution 10 ( Fig. 4). For lower lying areas with less varying topography and shallower snow depths, in particular in the eastern parts of Hardangervidda, the observed snow distributions were similarly close to a lognormal as to a gamma distribution. In higher elevated parts with more snow to the west of the plateau the snow distributions were much closer to a gamma distribution. Based on these findings a gamma distribution was used in the 15 main model runs, while a model run with lognormal distributions of snow was made to evaluate the sensitivity towards the choice of the distribution function (Sect. 3.2).

Evaluation of the snow distribution scheme
Three regression models for CV sd as a function of the terrain-based parameter Sx, elevation (z) and mean maximum snow depth (µ) were calibrated with the snow dis-20 tribution data from the ALS snow survey over the Hardangervidda mountain plateau (Table 1). Model 1 results in a root mean square error (RMSE) of only 0.14, however, the correlations of the distributions are significantly improved by including elevation as predictor (Model 2; R 2 = 0.52). By including maximum snow depth as additional predictor (Model 3) the model improves slightly to R 2 = 0.55 (Fig. 5). The distribution of rougher topography (western side of Norway) and higher elevations (central part following the Scandes), with maximum CV sd up to 1.2 in the Lyngen Alps and at peaks around Juvvasshøe (Fig. 1, site 2 and 4). The lowest values of 0.2-0.3 are modelled in larger valleys in south eastern Norway, where elevations are lower and topography gentler.

5
The regression models for CV sd are validated with data from GPR snow surveys at Juvvasshøe and Finse (Table 1). The correlation for Model 1 is poor, with R 2 = 0.04 and Nash-Sutcliff model efficiency (ME) = −0.7 (Table 1). Model 2 improves the correlation significantly, while the best fit is obtained with Model 3 (Fig. 5, RMSE = 0.094, R 2 = 0.62 and ME = 0.61). The improvement in Model 3 compared to Model 2 10 is more pronounced in the validation than in the fit of the regression models, and is mainly a result of better representation of the highest CV sd values. The validation area at Juvvasshøe is located at higher elevations than what is represented in the ALS snow survey data set and undergoes extreme redistribution by wind. The representation of extreme values therefore has a high impact in the validation run.

Modelled ground temperatures for mainland Norway
The main results presented in this section are based on the model run with 100 realizations per grid cell, applying gamma distributions over the CV sd from Model 3. According to the model run, in total 25 400 km 2 (7.8 %) of the Norwegian mainland is underlain by permafrost in an equilibrium situation with the climate over the 30-year period 1981-20 2010 (Fig. 1). 12 % of the land area features sub-zero ground temperatures in more than 10 % of a 1 km grid cell, and is classified as sporadic (4.4 %), discontinuous (3.2 %) or continuous (4.3 %) permafrost (Fig. 1). In comparison, the model run without a subgrid variation results in a permafrost area of only 13 460 km 2 , corresponding to 4.1 % of the model domain ( Table 2). The difference is illustrated for Juvvasshøe (Fig. 6a) and 25 Dovrefjell (Fig. 6c), where the sub-grid model very well reproduces the observed lower limit of permafrost based on borehole temperatures and BTS-surveys. In contrast, the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | model without sub-grid variability indicates a hard line for the permafrost limit at much higher elevations ( Fig. 6b and d). At Juvvasshøe, the model without sub-grid distribution still reproduces the permafrost limit to some extent because of the large elevation gradient. At Dovrefjell, where the topography is much gentler, the difference between the models is much larger and the approach without sub-grid distribution is not capable 5 of reproducing the observed permafrost distribution. The modelled permafrost area for model runs applying the other models for CV sd and theoretical distribution functions are summarized in Table 2. The standard deviations of the modelled sub-grid distribution of MAGT range from 0 to 2.5 • C (Fig. 7, right panel). The highest standard deviation values are found in the Jotunheimen area, where modelled sub-grid variability of MAGT is up to 5 • C. Also at lower elevations in south eastern parts of Finnmark standard deviations exceed 1.5 • C.
Here, the CV sd values are below 0.4, but because of cold (FDD a < −2450 • C) and dry (max SD < 0.5 m) winters even small variations in the snow cover result in large effects on the ground temperatures. 15 Close to 70 % of the modelled permafrost is situated within open, non-vegetated areas above treeline, classified as mountain permafrost according to Gruber and Haeberli (2009). This is the major part of the permafrost extent both in northern and southern Norway. In northern Norway the model results indicate that the lower limit of continuous/sporadic mountain permafrost decreases eastwards from 1200/700 m a.s.l., re-20 spectively, in the west to 500/200 m in the east. In southern Norway, the southernmost location of continuous mountain permafrost is in the mountain massif of Gaustatoppen at 59.8 • N, with continuous permafrost above 1700 m a.s.l. and discontinuous permafrost down to 1200 m a.s.l. In more central southern Norway the continuous mountain permafrost reaches down to 1600 m a.s.l. in the western Jotunheimen and  all boreholes except of one, this being 0.2 • C outside the range. All the average modelled MAGT are within ±1.6 • C of observations, while 90 % are within 1 • C. The RMSE between the observed and modelled average MAGT is 0.6 • C (Fig. 8, right panel).
The evaluation of the model runs with all three CV sd -models, as well as lognormal instead of gamma distribution functions are summarized in Table 2. The highest corre-20 lation between observed and mean MAGST and MAGT was obtained by Model 3, but Model 2 yielded similar correlations. All three model runs capture 58 % of the observed MAGST and more than 98 % of the observed MAGT within the temperature range of the corresponding grid cell. The total area of modelled permafrost is 9 % less when applying the simplest snow distribution model (Model 1) compared to the reference model

The effect of a statistical representation of sub-grid variability in a regional permafrost model
The total distribution of modelled permafrost with the sub-grid snow scheme corresponds to 7.8 % of the Norwegian land area, while the modelled permafrost area with-5 out a sub-grid representation of snow is ∼ 4 %. This large difference in total modelled permafrost area stems exclusively from differences in the amount of modelled permafrost in mountains above the treeline. In these areas the snow distribution is highly asymmetric with a majority of the area having below average snow depths. Because of the non-linearity in the insulating effect of snow cover the mean ground temperature of 10 a grid cell is not, and is often far from, the same as the ground temperature below the average snow depth. Often, the majority of the area in high, wind exposed mountains is nearly bare blown with most of the snow blown into terrain hollows. Consequently, most of the area experiences significantly lower average ground temperatures than with an evenly distributed, average depth snow cover. In mountain areas with a more gentle to- 15 pography and relatively small spatial temperature variations, an evenly distributed snow depth will result in large biases in modelled permafrost area, as illustrated at Dovrefjell in Fig. 6. This study is clear evidence that the sub-grid variability of snow depths should be accounted for in model approaches targeting the ground thermal regime and permafrost distribution.

20
The model reproduces the large range of variation in sub-grid ground temperatures, with standard deviations up to 2.5 • C. This is in accordance with the observed smallscale variability of up to 6 • C within a single grid cell (Gisnås et al., 2014;Gubler et al., 2011). Inclusion of sub-grid variability of snow depths in model approaches allows for a more adequate representation of the gradual transition from permafrost 25 to permafrost-free areas in alpine environments, and thus a better estimation of permafrost area. With a warming of the climate, a model without such a sub-grid representation would respond with an abrupt decrease in permafrost extent. In reality, bare 6677 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | blown areas with mean annual ground temperatures of −6 • C need a large temperature increase to thaw. Increased precipitation as snow would also warm the ground; however, bare blown areas may still be bare blown with increased snow accumulation during winter. A statistical snow distribution reproduces this effect, also with an increase in mean snow depth.

5
CryoGRID1 is a simple modelling scheme delivering a mean annual ground temperature at the top of the permanently frozen ground based on near-surface meteorological variables, under the assumption that the ground thermal regime is in equilibrium with the applied surface forcing. This is a simplification, and the model cannot reproduce the transient evolution of ground temperatures. However, it has proven to capture the re-10 gional patterns of permafrost reasonably well Westermann et al., 2013). Because of the simplicity it is computationally efficient, and suitable for doing test-studies like the one presented in this paper and in similar studies (Westermann et al., 2015).
For the model evaluation with measured ground temperatures in boreholes 15 (Sect. 5.4), the modelled temperatures are forced with data for the hydrological year corresponding to the observations. Because of the assumption of an equilibrium situation in the model approach, such a comparison can be problematic as many of the boreholes have undergone warming during the past decades. However, with the majority of the boreholes located in bedrock or coarse moraine material with relatively high 20 conductivity, the lag in the climate signal is relatively small at the depth of the top of permafrost.
The large amount of field observations used for calibration and evaluation in this study is mainly conducted in alpine mountain areas. The large spatial variation in winter snow depths is a major controlling factor also of the ground temperatures in peat 25 plateaus and palsa mires, and is a driving factor in palsa formation (e.g. Seppälä, 2011). The sub-grid effect of snow should therefore also be implemented for mire areas, where comparable data sets are lacking.

Model sensitivity
The sensitivity of the model for CV sd to the modelled ground temperatures is relatively low, with only 9 % variation in permafrost area, although the performance of the snow distribution scheme varies significantly between the models when evaluated with GPR snow surveys (Table 1). In comparison, a lognormal instead of a gamma distribu-5 tion function reduces the permafrost area by 18 % ( Table 2). The choice of distribution function therefore seems to be of greater importance than the fine tuning of a model for CV sd . This result contradicts the conclusions by Luce and Tarboton (2004), suggesting that the parameterization of the distribution function is more important than the choice of distribution model. With a focus on hydrology and snow cover depletion curves, equal 10 importance was given to both the deeper and shallower snow depths in the mentioned study. In contrast, an accurate representation of the shallowest snow depths is crucial for modelling the ground thermal regime. The low thermal conductivity of snow results in a disconnection of ground surface and air temperatures at snow packs thicker than 0.5-1 m, depending on the physical properties of the snow pack (e.g. Haeberli, 1973). 15 In wind exposed areas prone to heavy redistribution, large fractions of the area will be entirely bare blown (Gisnås et al., 2014). These are the areas of greatest importance for permafrost modelling. In order to reproduce the gradual transition in the discontinuous permafrost zone, where permafrost is often only present at bare blown ridges, shallow snow covers must be satisfactorily represented. Compared to a gamma func-20 tion, a lognormal distribution function to a larger degree underestimates the fraction of shallow snow depths, resulting in a less accurate representation of this transition. Several studies include statistical representations of the sub-grid variability of snow in hydrological models, most commonly applying a two-or three-parameter lognormal distribution (e.g. Pomeroy et al., 2004;Liston, 2004;Nitta et al., 2014;Donald et al., 25 1995). Observed snow distributions within 1 km × 1 km in the ALS snow survey presented in this paper are closer to a gamma than to a lognormal distribution, supporting the findings by Skaugen (2007) and Winstral and Marks (2014)  in non-forested alpine environments. However, the difference is not substantial in all areas; the two distributions can provide near-equal fit in eastern parts of the mountain plateau where the terrain is gentler and the wind speeds lower. We suggest that the choice of distribution function of snow is important in model applications for the ground thermal regime, and recommend the use of gamma distribution for non-vegetated high 5 alpine areas prone to heavy redistribution of snow. While a gamma distribution offers improvements over a lognormal distribution, the bare blown areas are still not sufficiently represented. One attempt to solve this is to include a third parameter for the "snow free fraction" (e.g. Kolberg and Gottschalk, 2010;. We made an attempt to calibrate such a parameter for this study, however, no correlations to any of the predictors were found. It is also difficult to determine a threshold depth for "snow free" areas in ALS data resampled to 10 m resolution, where the uncertainty of the snow depth observations are in the order of ten centimetres (Melvold and Skaugen, 2013).
In this study a high number of realizations could be run per grid cell because of the 15 low computational cost of the model. To evaluate the sensitivity of sampling density, the number of realizations was reduced from 100 to 10 per grid cell. This resulted is a 2.6 % increase in total modelled permafrost area relative to the reference model run. This demonstrates that a statistical downscaling of ground temperatures as demonstrated in this study is robust and highly improves the model results with only a few additional 20 model realizations per grid cell.

Conclusions
We present a modelling approach to reproduce the variability of ground temperatures within the scale of 1 km 2 grid cells based on probability distribution functions over corresponding seasonal maximum snow depths. The snow distributions are derived from 25 climatic parameters and terrain parameterizations at 10 m resolution, and are calibrated with a large scale data set of snow depths obtained from laser scanning. The model re- -The same permafrost model without a sub-grid representation of snow produces almost 50 % less permafrost. Because of the non-linearity in the insulating effect of snow cover in combination with the highly asymmetric snow distribution within 10 each grid cell, sub-grid variability of snow depths must be accounted for in models representing the ground thermal regime.
-Observed variations in ground surface temperatures from two logger arrays with 26 and 41 loggers, respectively, are very well reproduced, with estimated fractions of sub-zero MAGST within ±10 %. 94 % of the observed mean annual temperature 15 at top of permafrost in the boreholes are within the modelled ground temperature range for the corresponding grid cell, and mean modelled temperature of the grid cell reproduces the observations with an accuracy of 1.5 • C or better.
-The sensitivity of the model to the coefficient of variation of snow (CV sd ) is relatively low, compared to the choice of theoretical snow distribution function. How-20 ever, both are minor effects compared to the effect of running the model without a sub-grid distribution.
-The observed CV sd of snow within 1 km 2 grid cells in the Hardangervidda mountain plateau varies from 0.15 to 1.15, with an average CV sd of 0.6. The distributions are generally closer to a theoretical gamma distribution than to a lognormal 25 distribution, in particular in areas of very rough topography, thicker snow cover In areas subject to snow redistribution, the average ground temperature of a 1 km 2 grid cell must be determined based on the distribution, and not the overall average of snow depths within the grid cell. Furthermore, modelling the full range of ground tempera-5 tures present over small distances enables representation of the gradual transition from permafrost to non-permafrost areas and most likely a more accurate response to climate warming. This study is clear evidence that the sub-grid variability of snow depths should be accounted for in model approaches targeting the ground thermal regime and permafrost distribution.

10
Acknowledgements. This study is part of the CryoMet project (project number 214465; funded by the Norwegian Research Council). The field campaigns at Finse were partly founded by the hydropower companies Statkraft and ECO, while the field work at Juvvasshøe was done in collaboration with Ketil Isaksen (Norwegian Meteorological Institute). The Norwegian Meteorological Institute provided the NORA10 wind data and the seNorge gridded temperature data. 15 The Norwegian Water and Energy Directorate provided the seNorge gridded snow depth data and the ALS snow survey at Hardangervidda. Kolbjørn Engeland gave valuable comments to the statistical analysis presented in the manuscript. We gratefully acknowledge the support of all mentioned individuals and institutions. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Gubler, S., Fiddes, J., Gruber, S., and Keller, M.: Scale-dependent measurement and analysis of ground surface temperature variability in alpine terrain, The Cryosphere Discuss., 5, 307-338, doi: 10.5194/tcd-5-307-2011, 2011. Haakenstad, H., Reistad, M., Haugen, J. E., andBreivik, Ø