A statistical approach to modelling permafrost distribution in the European Alps or similar mountain ranges

Estimates of permafrost distribution in mountain regions are important for the assessment of climate change effects on natural and human systems. In order to make permafrost analyses and the establishment of guidelines for e.g. construction or hazard assessment comparable and compatible between regions, one consistent and traceable model for the entire Alpine domain is required. For the calibration of statistical models, the scarcity of suitable and reliable information about the presence or absence of permafrost makes the use of large areas attractive due to the larger data base available. We present a strategy and method for modelling permafrost distribution of entire mountain regions and provide the results of statistical analyses and model calibration for the European Alps. Starting from an integrated model framework, two statistical sub-models are developed, one for debris-covered areas (debris model) and one for steep bedrock (rock model). They are calibrated using rock glacier inventories and rock surface temperatures. To support the later generalization to surface characteristics other than those available for calibration, so-called offset terms have been introduced into the model that allow doing this in a transparent and traceable manner. For the debris model a generalized linear mixed-effect model (GLMM) is used to predict the probability of a rock glacier being intact as opposed to relict. It is based on the explanatory variables mean annual air temperature (MAAT), potential incoming solar radiation (PISR) and the mean annual sum of precipitation (PRECIP), and achieves an excellent discrimination (area under the receiver-operating characteristic, AUROC = 0.91). Surprisingly, the probability of a rock glacier being intact is positively associated with increasing PRECIP for given MAAT and PISR conditions. The rock model is based on a linear regression and was calibrated with mean annual rock surface temperatures (MARST). The explanatory variables are MAAT and PISR. The linear regression achieves a root mean square error (RMSE) of 1.6 C The final model combines the two sub-models and accounts for the different scales used for model calibration. The modelling approach provides a theoretical basis for estimating mountain permafrost distribution over larger mountain ranges and can be expanded to more surface types and sub-models than considered, here. The analyses performed with the Alpine data set further provide quantitative insight into larger-area patterns as well as the model coefficients for a later spatial application. The transfer into a mapbased product, however, requires further steps such as the definition of offset terms that usually contain a degree of subjectivity.

We present a strategy and method for modelling permafrost distribution of entire mountain regions and provide the results of statistical analyses and model calibration for the European Alps.Starting from an integrated model framework, two statistical sub-models are developed, one for debris-covered areas (debris model) and one for steep bedrock (rock model).They are calibrated using rock glacier inventories and rock surface temperatures.To support the later generalization to surface characteristics other than those available for calibration, so-called offset terms have been introduced into the model that allow doing this in a transparent and traceable manner.
For the debris model a generalized linear mixed-effect model (GLMM) is used to predict the probability of a rock glacier being intact as opposed to relict.It is based on the explanatory variables mean annual air temperature (MAAT), potential incoming solar radiation (PISR) and the mean annual sum of precipitation (PRECIP), and achieves an excellent discrimination (area under the receiver-operating characteristic, AUROC = 0.91).Surprisingly, the probability of a rock glacier being intact is positively associated with increasing PRECIP for given MAAT and PISR conditions.The rock model is based on a linear regression and was calibrated with mean annual rock surface temperatures (MARST).The explanatory variables are MAAT and PISR.The linear regression achieves a root mean square error (RMSE) of 1.6 • C. The final model combines the two sub-models and accounts for the different scales used for model calibration.
The modelling approach provides a theoretical basis for estimating mountain permafrost distribution over larger mountain ranges and can be expanded to more surface types and sub-models than considered, here.The analyses performed with the Alpine data set further provide quantitative insight into larger-area patterns as well as the model coefficients for a later spatial application.The transfer into a mapbased product, however, requires further steps such as the definition of offset terms that usually contain a degree of subjectivity.

Introduction
Many models already exist for estimating the spatial distribution of mountain permafrost in regions of the European Alps (Hoelzle, 1992;Keller, 1992;Imhof, 1996;Gruber and Hoelzle, 2001;Lambiel and Reynard, 2001;BAFU, 2005).These models are difficult to compare or combine because they have different empirical or statistical approaches and differing types of indices as output.Their extrapolation is difficult as they are usually calibrated for specific regions and, as a consequence, no estimation of permafrost distribution exists in large regions of the European Alps to date.
First estimations of permafrost occurrence in the Alps were based on the "rules of thumb" (Haeberli, 1975), which use basic relations of permafrost occurrence with L. Boeckli et al.: Statistical permafrost distribution model topographic attributes and surface characteristics.These relationships were first implemented within a GIS environment by Keller (1992) and later incorporated in further studies to predict spatial permafrost occurrence (Imhof, 1996;Frauenfelder et al., 1998;BAFU, 2005;Ebohon and Schrott, 2008).Besides topographic variables, climatic information and other direct proxy variables of the surface energy balance (such as MAAT and PISR) are often used in statistical or empirical permafrost models.The basal temperature of snow (BTS), introduced by Haeberli (1973) as an indicator of permafrost occurrence, has been widely used for model calibration (Hoelzle, 1992;Keller et al., 1998;Riedlinger and Kneisel, 2000;Gruber and Hoelzle, 2001;Stocker-Mittaz et al., 2002;Lewkowicz and Ednie, 2004;Brenning et al., 2005).Measurements in boreholes and near the ground surface have been used for model evaluation (Gruber et al., 2004;Heggem et al., 2005;Etzelmüller et al., 2006Etzelmüller et al., , 2007;;Allen et al., 2009).Other studies used rock glacier inventories to infer the occurrence of permafrost (Janke, 2004), to identify the lower boundary of discontinuous permafrost (Nyenhuis et al., 2005), or for model assessment (Imhof, 1996;Gruber and Hoelzle, 2001).
Existing permafrost distribution models typically do not distinguish between different surface characteristics such as fine or coarse-grained debris slopes or steep bedrock, even though their differences in snow-cover and material characteristics can cause strongly differing ground temperatures.Because the amount and type of available permafrost data differs between the various surface types this often results in an unknown degree of extrapolation.Models based on BTS for example are also used to predict permafrost occurrence in bedrock without snow or on ground that is nearly always blown free of snow.Similarly, models based on rock glacier occurrence are used to predict permafrost also in fine-grained substrate.The output of most established models consists of permafrost zonation classes such as "probable permafrost" or an index that is then related to a qualitative probability.These quantities are usually defined beforehand (e.g.rules of thumb (Haeberli, 1975); the BTS relationship (Haeberli, 1973); or based on the concept of permafrost zonation limits (Lambiel and Reynard, 2001)) and cannot be validated quantitatively later on.A notable exception is the work by Lewkowicz and Ednie (2004) but applied to more diverse terrain, extrapolation to differing surface conditions may again be problematic.
With this paper we aim to base a statistical model on reliable permafrost data for a larger area and to make the extrapolation to surface types for which no suitable calibration data is available transparent by introducing dedicated offset terms.While this does not solve the problem of difficult validation due to strong heterogeneity and little data, it provides a basis for better model calibration with the data available and for better separation of quantitative statistical analysis and of subjective adjustment.The objectives of the present study are thus to (a) introduce a suitable strategy and method for statistic-based modelling the permafrost distribution for large mountain regions, (b) to discuss generic difficulties of such models regarding model application to an entire landscape, and (c) to provide the results of our statistical analysis (i.e.model coefficients) for the European Alps.

Conceptual background
The lack of sufficient and reliable data for calibration and validation probably is one of the most important limitations for permafrost modeling and it is important to devise strategies for the efficient use of existing data.While much progress has been made in therms of physics-based modelling of mountain permafrost in recent years (cf.Riseborough et al., 2008;Harris et al., 2009), also those methods are challenged in terms of their validation for diverse conditions and we therefore decided to use a statistical approach.
To benefit from a large data basis, the presented model is based on an Alpine-wide collection of permafrost evidence (Cremonese et al., 2011).Overall, only MARST measurements and the rock glacier inventories offer enough data to support the fitting of statistical models.The other permafrost observations described by Cremonese et al. (2011) were not used for model calibration because, (a) they are not sufficient in number to allow consistent statistical analysis; (b) the integration and homogenization of heterogeneous permafrost observations are subject to large uncertainty and subjectivity, and (c) observations are strongly biased towards permafrost existence and less observation in non-permafrost conditions are available.The evidence collection is subject to some degree of homogenization, especially the rock glacier inventories in it may be subject to differences due to slightly differing conventions and data used that must be accounted for in statistical models.To utilize rock glacier inventories and MARST data, a combination of two sub-models is required: one for debris-covered areas and one for steep bedrock.With the term "steep bedrock" we refer to terrain that (a) is not or only marginally affected by a snow cover in wintertime, (b) does not contain large amounts of blocks and/or debris, and (c) is without vegetation coverage.The statistical models are calibrated based on rock glacier activity status as a binary variable, and rock temperatures as a continuous variable, respectively.While this model combination results in probabilities sensu strictu, their application to areas that are not rock glaciers or steep bedrock is difficult.Because we only have limited evidence for permafrost occurrence or absence, this extrapolation is an integral part of permafrost modeling.To make this transparent, we use offset terms that we embed in the model to allow later subjective adjustment.Because the measurements of MARST used for model calibration are representative on a scale of few tens of centimeters only, the application of results in a model based on a DEM of several tens of meters is difficult and requires dedicated attention.Because MAAT has regional trends besides its local dependance on elevation, we prefer to use MAAT instead of elevation as an explanatory variable.The modelling approach we present is suitable for large mountain regions and is here demonstrated for application to the European Alps (43

Debris model
In debris slopes, intact (active and inactive) rock glaciers are a diagnostic and well visible geomorphological feature to detect the presence of permafrost whereas relict forms indicate its absence (e.g.Haeberli, 1985).Active and inactive rock glaciers are grouped as "intact" rock glaciers because of their existing ice content (cf.Haeberli, 1985;Ikeda and Matsuoka, 2002;Roer and Nyenhuis, 2007;Lilleoren and Etzelmüller, 2011) and the reliable indication of permafrost they offer.Relict rock glaciers do no longer contain ice, show a collapsed surface due to melting of the ice, and they often have a vegetation cover (Roer and Nyenhuis, 2007).
The possibility of mapping rock glaciers from e.g.aerial photographs or field observations makes this an attractive and unique data source.The existence of rock glaciers depends on suitable debris production and transport mechanisms, implying that permafrost can also exist in areas where rock glaciers are absent (Imhof, 1996).Because here, no visible features indicate permafrost, a model-based estimation is especially valuable.Due to a cooling effect of the coarse block surface (Harris and Pedersen, 1998) and the creeping of rock glaciers, an estimation of permafrost distribution based on rock glacier activity generally results in an overestimation of the amount of permafrost below surfaces that are not rock glaciers.This effect may be compensated by offset terms to describe the permafrost status of those surfaces relative to the status of rock glaciers.
The debris model is calibrated using status information of rock glaciers resulting in a binary response (permafrost yes/no).While generalized linear models (GLMs) are commonly used to model binary response variables such as presence/absence of permafrost (Lewkowicz and Bonnaventure, 2008), an extension of this model that is able to account for random inventory effects is required here.Such random effects may be related to different observation techniques and interpretation criteria being applied in the compilation of inventories by different research groups, which may result in an inventory-specific bias.The generalized linear mixed model (GLMM; Venables and Ripley, 2002) is able to account for such random inventory effects and is therefore applied in this study.It is implemented as "glmmPQL" in the R package "MASS" (Venables and Ripley, 2002).GLMs and GLMMs for binary response variables can be specified using either the probit or, more commonly, the logistic link function (Hosmer and Lemeshow, 2000;Gelman and Hill, 2007).While both are nearly identical (Aldrich and Nelson, 1984;Gelman and Hill, 2007), a probit link function is preferred in this study due to mathematical advantages in the combined framework, as outlined in Sect.3.1.

Rock model
Near-surface temperatures in steep bedrock have been observed in the Alps since 2002 (Gruber et al., 2003;PERMOS, 2010).Based on this data and other analyses (Gruber et al., 2004;Allen et al., 2009), it was shown that short-wave radiation is the major controlling factor for the lateral variability of MARST in steep and homogeneous rock making it suitable for a linear statistical model.Because this is based on near vertical slopes that have no blocky layer and no snow, the extrapolation into the most common type of rock slope that is heterogeneous, fractured and partly snow-covered requires special attention.Differing mechanisms and effects have been postulated (Gruber and Haeberli, 2007;Noetzli et al., 2008) andmeasured (Hasler et al., 2011), and require an offset term for their inclusion in a permafrost distribution model.Measurements of ground surface temperatures (GST) or BTS in less steep terrain were not included in our analysis because of their large inter-annual variability caused by the strong influence of the snow cover (Hoelzle et al., 2003;Brenning et al., 2005).

Model combination
In our case of Alpine-wide permafrost distribution modelling, we wish to integrate two models that are fitted separately in two different model domains: debris surfaces and steep bedrock.While the debris model is based on an Alpinewide digital elevation model (DEM) with coarse grid spacing, the rock model is calibrated using locally measured terrain attributes, which refer to fine-scale topographic information.When combining these two models we have to consider scale effects with particular emphasis on the situation where an empirical model developed using fine-scale in situ measurements is applied at a coarser resolution for regional-scale application.

Statistical method
The introduction of our statistical approach for permafrost modeling starts by using the probit model formulation to show the formal equivalence between permafrost probabilities derived from temperature as a continuous random variable, and presence/absence as a binary one (Sect.3.1).This lays the foundations of the proposed approach, because this allows us to establish a unified framework for the rock model and debris model, which are based on these two different types of response variables.A simple approach for combining permafrost probabilities from rock and debris models is then proposed (Sect.3.2).This may result in scale issues, which require further attention (Sect.3.3).Two practical aspects of the application of our approach are then www.the-cryosphere.net/6/125/2012/The Cryosphere, 6, 125-140, 2012 discussed, the discrimination between rock and debris surfaces (Sect.3.4), and the assessment of the models' goodness of fit (Sect.3.5).

Model formulation
Permafrost is defined thermally by the permanent presence of zero or negative ground temperatures ( • C) over two entire years (van Everdingen, 1998).Because we are interested in depths where the variability of annual ground temperatures can be neglected, we assume maximum and mean ground temperatures to be equal.We therefore express the probability p of permafrost occurrence at a given location by the probability of mean ground temperature ϑ being ≤ 0 • C: If the ground temperature ϑ is modeled linearly (as we will later do in the rock model), where ε is a normally distributed residual error term with mean 0 and variance σ 2 , the x i and β i are the model's k explanatory variables and their coefficients, and α + represents an intercept term that is explained in detail in Sect.3.3.Throughout this work, model coefficients with a tilde refer to the temperature scale as in Eq. ( 2), while model coefficients at the probit scale will carry no tilde.In a predictive situation, this model will allow us to predict ϑ with a variance σ 2 pred ≥ σ 2 , which can be estimated from the model.In this situation, the permafrost probability p is therefore predicted to be where is the cumulative standard normal distribution function.The negative sign is due to the fact that we are interested in the probability of negative rather than positive temperatures.
On the other hand, direct evidence of permafrost presence or absence (debris model) allows us to model the permafrost probability p directly using generalized linear models.The probit link function is used in this study, because of its relation to the cumulative normal distribution function (Aldrich and Nelson, 1984;Gelman and Hill, 2007).Here, the probability of permafrost presence is modeled linearly not at the probability scale but at the probit scale, which is obtained from an inverse cumulative distribution function of the standard normal distribution: Thus, and if we introduce an additional (thermal) offset term (Sect.3.2) into the probit model, we write the debris model as where the x i and β i are the model's k explanatory variables and their corresponding coefficients, and α is the model intercept.
From Eqs. (3), ( 4) and ( 5) it becomes evident that a linear regression of ϑ is equivalent to a probit regression of p with scaled coefficients: This relationship between the temperature-based model and the model based on rock glacier presence/absence allows us to convert the temperature-based model in Eq. ( 2) into a probit-based probability model (Eq.5).It will later be shown that in this specific context, this is relatively insensitive to the estimation of the prediction variance, and we will therefore use a conservative variance estimator σ 2 pred , which will later be specified.The above equivalence allows us to integrate continuous-and binary-response permafrost distribution models within the formal framework of a linear model with comparable model coefficients.

Integration of continuous-and binary-response models
The coefficients of model M d (debris surface) are derived from the debris model, resulting in the coefficients α d ,β d,1 ,β d,2 adopted from Eq. ( 5).d is introduced into this model as a fixed offset value that can be used for adjusting effects such as rock glacier movement; this value is not estimated from the data but represents the possibility to later introduce an expert-defined adjustment term.
Model M r (rock surface), by contrast, is derived from the rock model (Eq.2) and partly uses the same explanatory variables as model M d , with the exception of a difference in spatial scale (discussed in Sect.3.3).It is important to note that in this model formulation, the adjustment offset r can be directly interpreted as a thermal offset of the nearsurface ground temperature (MARST) minus the temperature at the top of permafrost (TTOP).Given this model's prediction variance σ 2 pred,r , we estimate the probit-scale coefficients of M r from Eq. ( 6), i.e. by dividing all temperaturescale model coefficients by − σ pred,r .
In practice, the spatial distribution of different surface types is usually not well known and may exhibit transitions such as spatially varying debris or snow cover thicknesses.We represent this in a simple way through a (spatially varying) degree of membership in a land cover class, m τ , with values between 0 and 1 that sum up to 1 at each location.The integrated model is then defined to be which has an obvious generalization to more than two land cover classes.Probabilities of permafrost occurrence can be obtained from this integrated probit value by applying the inverse probit transformation, and probit-scale prediction variances are integrated in a similar way as the weighted sum of each model's prediction variances.

Scaling issues
With scale effects we refer to the fact that model coefficients may change at different scales or levels of aggregation as coarser-scale explanatory variables tend to show a smaller range of values and less scatter.This situation is related to the change of support problem (e.g.Gotway and Young, 2002), but instead of a geostatistical interpolation setting we need a solution that is tailored to the situation of integrating two linear models.
We start by looking at the scaling problem encountered in the situation where the rock model is fitted at a fine scale (parameters with index "F") and applied at a coarser scale (index "C") and consider initially only a linear model with one explanatory variable (k = 1) and no offset term F = 0.
Thus, from Eq. ( 2), where the residual variance is var ε F = σ 2 F .In a predictive situation, we have to approximate the finescale x F with its coarse-scale equivalent x C .We therefore predict x F using a scaling model, where the residuals shall be assumed to be independent and identically distributed according to a normal distribution with mean 0 and variance varε C = σ 2 C .The function f represents an arbitrary predictive model, such as a linear regression in x C .More generally, we could approximate x F using a model built on multiple variables other than x C .Thus, where the residuals are Since the spatial predicitons are to be made at the coarse scale, where one grid cell is composed of N fine-scale grid cells, we have where we make use of the fact that x C does not vary within a coarse-resolution grid cell.We refer to the last two terms, which involve the fine-and coarse-scale residuals, as the residual of ϑ C .
The estimation of the residual variance of ϑ C is not an easy task because the within-cell residuals ε F,i and ε C,i , respectively, can certainly not be considered to be independent due to the likely presence of (positive) spatial autocorrelation over these short distances.The variances of the partial residual terms would be expected to decrease proportionally to N −1 under the assumption of independent within-cell replication.The estimation variance of the mean value of positively autocorrelated random variables, by contrast, decreases more slowly with increasing sample size.In the extreme case of perfect within-cell dependence, averaging over N identical pseudo-replications would not reduce estimation variance compared to using only one replication.We adopt this conservative approach by assuming that averaging over N finer-scale grid cells does not reduce the uncertainty varε in the statistical model ground temperature at the aggregate scale (Hurlbert, 1984).In addition, we replace β 2 F with the square of a one-sided (upper) 95-% confidence limit of |β F |, β 2 F,cl .Thus, as a conservative estimator for varε , we use Consequently, the residual variance of the scaling model adds to the residual variance of the scaled model, using the regression coefficient for variance weighting, and the equation would be expanded by additional β 2 i,F,cl σ 2 i,C for each additional explanatory variable to be scaled.Estimates of β F and σ 2 F can be obtained from the fine-scale rock temperature model, and an estimate of σ C from the scaling model.
In a predictive situation, σ 2 F can be replaced with the corresponding prediction variance of the rock temperature model, which is generally slightly greater than σ 2 F .The prediction variance varies, however, slightly between samples.In the present study, the prediction variance is inflated only by 6 % on average, with a maximum of 11 %, and we therefore increase σ 2 F by 6 % in general in this study as a first-order approximation.

Surface types
To distinguish between the two model domains (debris vs. bedrock) one of the following approaches can be applied: (1) an index describing the degree of membership in the exposed bedrock rock surface class, (2) a statistical model of land cover such as a logistic regression or generalized additive model (Hastie and Tibshirani, 1990) or (3) remotelysensed or map-based land cover products.

Model evaluation
To assess the accuracy of the debris model, the area under the receiver-operating characteristics curve, which is known as AUROC, was calculated.This value ranges between 0.5 (random model behavior) and 1.0 (perfect model; Hosmer and Lemeshow, 2000).AUROC values reported in this study  2011) are based on model predictions that include the inventory random effect.A 10-fold cross-validation was performed to assess how transferrable to independent test data sets the model is (Hand, 1997).The original data set was randomly partitioned into 10 sub-samples.Of these 10 sub-samples, a single sub-sample was retained for testing the model, and the remaining 9 subsamples were used as training data.This process was repeated 10 times using each of the 10 sub-samples exactly once as the validation data.The 10 results from the folds were combined to produce a single estimation which then was used to measure the AUROC.
The goodness-of-fit for the rock model was obtained by calculating the R 2 and the root mean square error (RMSE).Furthermore, the RMSE resulting from a 10-fold crossvalidation was calculated.

Response variables
Most of the rock glacier inventories used to fit the debris model were provided by the permafrost observation collection of the PermaNET project (Cremonese et al., 2011).This collection was complemented by inventories from Switzerland published at the Seventh International Conference on Permafrost ("Yellowknife inventories"; ICP Yellowknife, Canada, 23-27 June 1998; Delaloye et al., 1998;Frauenfelder, 1998;Hoelzle, 1998;Imhof, 1998;Reynard and Morand, 1998;Schoeneich et al., 1998) and an inventory from the Upper Engadine, Switzerland (Frauenfelder et al., 2001;Frauenfelder, 2005).The final data set used as basis for the model development includes 2184 intact and 4218 relict rock glaciers from Austria, France, Italy and Switzerland (Table 1, Fig. 1).
For each rock glacier, information concerning its activity is available.The activity information from the different inventories was reclassified into the two classes (1) intact and (2) relict.The inventories from the PermaNET data contain polygon information for each rock glacier.From this inventories a stratified random sample was selected that resulted in one random point within the polygon for each rock glacier.For the Yellowknife inventories the centroids of the rock glaciers were used instead, because polygon information was unavailable.Finally, from each of the inventories an equal number of intact and relict rock glacier samples was drawn randomly in order to obtain balanced samples.
MARST data from France, Germany, Italy and Switzerland contained in the PermaNET inventory were used (Pogliotti, 2006;Pogliotti et al., 2008;PERMOS, 2010;Cremonese et al., 2011) and complemented with additional measurements from Switzerland (Hasler et al., 2011, Table 1, Fig. 1).All 57 sensors were located in rock walls > 55 • steep and several meters above flat ground to ensure snowfree conditions.The data originate from eight areas (Fig. 1) within which a wide range of aspects and elevations has been sampled.Measurement depths are on the order of 10 cm.The rock types sampled vary between areas and include limestone, granite and gneiss.The attribute data for each logger contains elevation, slope angle and aspect (measured in the field) as well as the observation period (logger years) taken for the calculation of MARST values.With this information MARST measurements of single or few years duration were adjusted to longer-term temperature trends according to Allen et al. (2009): longer-term MAAT from Piz Corvatsch (Upper Engadina, MeteoSchweiz, 2010) for the period 1961-1990 (MAAT = −6 • C) were compared with MAAT of the period corresponding to the specific logger years.The difference between these temperatures was used to correct the MARST values.The underlying assumption is that the difference of MAAT to its longer-term mean and the difference of MARST to its longer-term mean are equal (cf.Fig. 3.1 of PERMOS, 2010).By using the period of 1961-1990 as reference, the air temperature warming especially in the past decade due to climate change is neglected.This leads to an optimistic estimation (biased towards an overestimation of permafrost distribution) of MARST, but is in line with the debris model that also follows an optimistic approach.In comparison with the high number of rock glaciers available, 57 measurement points are few.They are, however, used for describing a system that is much less complicated than rock glaciers because the influence of snow, phase change, a mixed-media active layer and the downslope displacement of ice-rich material is minimal or non-existent.

Explanatory variables
As potential explanatory variables we consider PISR, MAAT, PRECIP, and a seasonal precipitation index (SEASONAL).PISR was derived from the Global Digital Elevation Model (GDEM; Hayakawa et al., 2008) with a grid spacing of 1 (approximately 30 m) using RSAGA (Brenning, 2008) and the algorithm of Wilson and Gallant (2000).PISR was calculated for one year with an hourly temporal resolution and clear sky conditions (100 % atmospheric transmittance) and is calculated for a latitudinal extent of 1 • (6 bands according to the total latitudinal extent of our study area).ASTER GDEM covers the entire Alpine arc and shows an overall vertical accuracy on a global basis of approximately 20 m at 95 % confidence (USGS et al., 2009).Alpine-wide MAAT data for the period 1961-1990 (Hiebl et al., 2009) was provided by the Central Institute for Meteorology and Geodynamics (ZAMG, Austria).MAAT is based on the GTOPO30 elevation model (Center, 1997) with an approximate resolution of 1000 m and shows a monthly standard error of less then 1 • C (Hiebl et al., 2009).A constant lapse rate of 0.65 • C 100 m −1 was used to interpolate the coarse MAAT based on more precise elevation information from the ASTER GDEM.
Alpine-wide monthly precipitation data (Efthymiadis et al., 2006) is available for 1800-2003, gridded at 10 resolution (approximately 15 km, available from ALP-IMP, http: //www.cru.uea.ac.uk/cru/data/alpine/).Based on this data, PRECIP for the period 1961-1990 was calculated.As potential explanatory variable for the debris model, PRECIP was centered (cPRECIP) by subtracting its mean value of 1271 mm.Centering PRECIP allows to directly compare the coefficients of the different models including and excluding PRECIP as explanatory variable.Additionally, an index describing the seasonality of precipitation (SEASONAL) was computed by dividing the mean sum of summer precipitation (May-October) by the mean sum of winter precipitation (November-April).
For the locations of the MARST loggers, the usage of locally measured terrain parameters is necessary for the characterization of MARST, because they strongly depend on micro-topographic radiation effects such as sun exposure or terrain shading (the resolution of the ASTER GDEM is too coarse for this purpose.).For increased accuracy, PISR was therefore calculated following Corripio (2003)   fish eye lens (Gruber et al., 2003) and considered in the PISR calculations.Further, the MAAT provided by ZAMG was adjusted for the logger locations using local elevation information measured in the field.

Alpine-wide permafrost model
This section contains the model calibration and interpretation for the debris and the rock model.For each sub-model, three different sets of explanatory variables were used to compute regression models.Afterwards, the final model was chosen based on goodness-of-fit-statistics.In the last subsection of this chapter, the combination of the sub-models is presented.

Debris model
A total of 3580 rock glacier points (Table 2) were used for model calibration.While MAAT and PISR in the debris model show a clear relation to the activity of rock glaciers, the correlation with precipitation is less obvious in a univariate analysis (Fig. 2).Nevertheless, cPRECIP was included in the final model based on the high significance of the Wald test (Table 3).The seasonality of precipitation (SEASONAL) shows no significant contribution within the debris model and was therefore omitted for the final model.The chosen GLMM includes MAAT, PISR and cPRECIP as fixed effects and the membership of each point in the different inventories as random effects (Table 3).All explanatory variables show a high significance (p-value).When considering random effects, the debris model achieves an AU-ROC of 0.91, respective 0.91 for the 10-fold cross-validation (AUROC cv ), which both are "outstanding" discriminations according to Hosmer and Lemeshow (2000).
The coefficients of the final model indicate: a difference in cPRECIP of 400 mm is identical with a change of 0.52 on the probit scale.A difference in MAAT of 1 • C is equivalent to a probit-change of 0.91.Thus, a change in cPRECIP of 400 mm is identical to a difference in MAAT of 0.57 • C and leads to a dislocation of the limit between intact and relict rock glaciers of 88 m (assuming a constant lapse rate of 0.65 • C 100 m −1 ).An increase of 240 W m −2 (approximate difference in PISR of a south vs. north exposed slope with an angle of 30 • ) is associated with a decrease of 1.78 on the probit scale.This change is equivalent to an increase in MAAT by 1.96 • C or approximately 300 m in elevation.

Rock model
For all 57 locations, MARST are higher than MAAT (Table 4, Fig. 3, top left) and the difference between MARST and MAAT increases with higher PISR (Fig. 3, top right).PRECIP was not included in the rock model, because the variable showed no high significance and it deteriorates the Akaike Information Criterion (AIC), which measures model fit while penalizing for model size (Table 5; Gelman and Hill, 2007).SEASONAL was omitted from the final model because its range of values on the present training sample was too narrow (from 0.76 to 1.66) to allow for an Alpine-wide application of this empirical relationship (SEASONAL between 0.50 and 2.47).The coefficients of the chosen model indicate that MARST are generally warmer than the corresponding MAAT.An increase in PISR of 240 W m −2 is associated with a decrease in MARST of 4.6 • C and is equivalent to a change in MAAT of 4.2 • C. Thus, a change in slope aspect from south to north has a similar influence on MARST as a change in elevation of approximately 650 m.

Scaling model and model combination
A LIDAR DEM covering South Tyrol (data provided by Autonomous Province of Bolzano -South Tyrol, Italy) with a resolution of 2.5 m was used to estimate the prediction variance of the scaling model.The other variance component was estimated from the rock model.PISR derived from the LIDAR DEM refers to local, "real world" estimates and can be compared with PISR values calculated for the rock logger locations.
The following linear regression was fitted to a random sample of 28 640 points within South Tyrol above 2000 m and relates finer-scale (2.5 m, LIDAR DEM) PISR to coarsescale values calculated from a reduced-resolution and reduced quality (30 m, ASTER GDEM) equivalent: This model resulted in an R 2 = 0.72 and a residual standard error of 46 W m −2 .
For the two other explanatory variables (MAAT and cPRE-CIP), no scaling correction was necessary because both variables show negligible spatial variation within ASTER GDEM grid cells.
The conservative estimation of σ (Eq.14) obtained a value of 1.95 • C and was used for converting predicted MARST into the corresponding probit-based values (Eq.6).The individual variance components are displayed in Table 6.The adjustment parameters r (Eq.2) and the d (Eq.5) were set to zero.Both models (debris and rock model) were then combined using Eq. ( 7).Probabilities of a rock glacier being intact as opposed to relict, respectively probabilities of MARST ≤ 0 • C in steep bedrock, were obtained by applying the inverse probit transformation (Fig. 4).A sample application of the model showing a map-based output product is presented in Fig. 5.

Use and limitations of the model
The presented model approach is based on statistical relations and thus limited in the ability to represent physical processes such as snow redistribution by avalanche and wind that is known to have an impact on mountain permafrost occurrence (Haeberli, 1975;Hoelzle et al., 2001).To account for the different thermal responses related to surface conditions in two domains (debris and rock cover) an adjustmentoffset can be applied in our model for each sub-domain model individually (Sect.3.1).Identifying suitable average adjustment parameters for each domain is challenging because of the large spatial variation of the offsets between different locations (Hoelzle and Gruber, 2008).
The explanatory variables MAAT and cPRECIP are derived from existing data sources (Sect.4.2).PISR estimates are based on a DEM.For an Alpine-wide model application, the ASTER GDEM can be used to calculate PISR values.For regional model application (e.g.South Tyrol, Italy), where more precise DEM data is available, this could be used to derive the PISR values.Functions similar to Eq. ( 15) are then needed to address the scaling from fine to coarse resolution for the debris model.The prediction is also possible based on two different DEMs: a coarse elevation model (e.g.ASTER GDEM) for the debris model representing the mesoscale characteristics of rock glaciers, and a more precise DEM for the rock model because MARST values more strongly depend on accurate PISR estimates.The ASTER GDEM, which is used in this study to calculate PISR and to rescale MAAT for the debris model, shows limitations in the Alps when compared to a more reliable DEM (Frey and Paul, 2011).However, no better DEM is available at the moment for the entire Alps.The spatial distribution of the rock glaciers used for model calibration (debris model) nearly covers the entire Alps.In contrast, only 57 MARST measurements mostly from the central part of the Alps were available.This inevitably requires a strong generalization of the rock model, especially regarding the precipitation (Sect.6.2).The temporal extrapolation of MARST values to the period 1961-1990 was addressed by using long-term MAAT measurements.However, the corrected data is sensible to inter-annual variability, because some of the measurement series were only one year long.

The
The transition zone between debris covered slopes and steep bedrock requires further investigation.Some ground surface temperature (GST) measurements exist in this zone but as mention in Sect. 2 the large inter-annual variability makes this data unsuitable for statistical modelling.

Influence of precipitation
The precipitation variable in the debris model can be seen as a simple proxy for the amount of snow in a regional context or the reduction of short wave insolation by cloud cover.The positive coefficient of precipitation in the regression model (Sect.5.1) implies that in areas with higher precipitation, rock glaciers are more likely to be intact, or equivalently, the limit between intact and relict rock glaciers tends to shift towards lower elevations.According to our model, this means that for given MAAT (or elevation in a local context) and PISR conditions, the boundaries of permafrost occurrence in debris-covered and wet areas of the Alps are on average approximately 220 m lower than in relatively dry areas with 1000 mm lower PRECIP.This contrasts with several studies that state that permafrost boundaries are lower in dry or continental areas (e.g.Barsch, 1978;King, 1986), but it is consistent with regional-scale trends in the lower limit of intact rock glacier distribution in the Andes of Central Chile (Brenning, 2005;Azócar and Brenning, 2010).Debris rock glaciers, referring to rock glaciers developed in strong relation to a glacier (e.g.Barsch, 1996;Hughes et al., 2003), may offer an explanation for the positive coefficient of PRECIP in the debris model.However, their precise definition is difficult and as a consequence their influence on the debris model cannot be assessed.
The positive influence of precipitation regarding the intactness of a rock glacier is also shown in Fig. 6, where three different models without precipitation as explanatory variable for drier, normal and relatively wet areas are compared.The three models were calibrated using three subsamples of the entire data set representing drier, normal and relatively wet inventories (drier: mean PRECIP = 1105 mm, normal: mean PRECIP = 1291 mm, and relatively wet: mean PRECIP = 1679 mm).The variable precipitation is not just significant, but also relevant regarding possible model prediction as shown in Fig. 7.The predicted values modelled with cPRECIP as explanatory variable differ with a maximum of 1.5 • C (or approximately 200 m of elevation) from the model prediction without cPRECIP included as explanatory variable.
To further investigate possible relationships between precipitation and the spatial density of rock glaciers, we compared data from two different rock glacier inventories, for which the inventory perimeters were manually digitized (IGUL, Tecino, Switzerland and GEOL,Trentino, Italy); inventory boundaries are currently not available for the other inventories.The results show that rock glacier density in the Alps tends to be higher in areas with less precipitation (Fig. 8).This could explain the widespread notion that also permafrost boundaries occur at lower elevation in dry areas.
The correlation of MARST and precipitation is weak (Fig. 3, bottom right) and PRECIP shows no significance in the rock model (Table 5).However, the observed significance and magnitude of the influence of SEASONAL suggests that further research on the physical relationship of precipitation seasonality on rock temperatures would be desirable, and that a larger rock temperature data basis would allow us to incorporate an additional relevant predictor variable into the model.According to Gruber et al. (2004) 3 calculated for a randomly selected probability range of 0.475-0.525for intact rock glacier occurrence.Black crosses: debris model without cPRECIP as explanatory variable, red bubbles: debris model including cPRECIP as explanatory variable.

Conclusions
We have presented an approach and statistical model designed to cater for the specific needs of permafrost distribution estimation for entire mountain regions.Based on this, rock glacier inventories and MARST measurements were used to calibrate model coefficients for the European Alps.By using intact and relict rock glaciers as calibration data, the prediction of the debris model is biased towards and overestimation of the permafrost distribution, while the rock model generally underestimates the current permafrost distribution.Current data does not permit extending the statistical analyses to other types of surface cover with the same statistical rigor, as a high number of observations would be required.This is a fundamental challenge to all statistical permafrost distribution modelling, and equally to the validation of physically-based numerical models.However, the quantitative statistical model presented already contains offset terms to allow later subjective adjustment for the extension to other surface types.
By allowing analyses (i.e.model calibration) in larger areas, more robust and new insight can be derived because of more available data and due to larger environmental gradients covered by one analysis.In this way, a shift of the limit between intact and relict rock glaciers towards lower elevation has been detected to coincide with increasing precipitation.This influence of precipitation needs further investigation because it conflicts with previous but less quantitative studies.However, this is the first investigation known to the authors, which systematical analyses the spatial rock glacier distribution in relation to precipitation patterns with a large data sample in the Alps.
The model presented is build on reliable data and subjective adjustment is clearly separated from statistical analysis.In this sense, this paper is a step towards an Alpine permafrost map, but the model and coefficients presented cannot be directly applied to an entire landscape and they require subjective adjustments of the offsets d and r .The following steps are needed to use the presented approach for Alpine-wide model application and to provide a map-based output product: (a) the definition of offsets terms for surfaces other than steep bedrock and rock glaciers; (b) the derivation of scaling functions to correct PISR estimates derived from different DEMs; (c) the preparation of gridded land cover maps quantifying the membership in the two classes debris and bedrock slopes as well as any other surface type considered with an offset term; and (d) the establishment of a legend and interpretation guidelines for map users.After the inclusion of offset terms, model results are no longer probabilities but should be considered to be an index, which describes the permafrost occurrence per grid cell and than can be shown in a map (cf.Boeckli et al., 2012).

Fig. 1 .
Fig. 1.Spatial distribution of intact/relict rock glaciers (blue crosses) and the locations of the rock surface temperature loggers (red crosses).

Fig. 2 .
Fig. 2. Frequencies of intact as opposed to relict rock glaciers conditional on potential explanatory variables.Bar widths in these spinograms are proportional to the empirical frequency of the given interval of values of explanatory variable.The figure does not account for random inventory effects.

Fig. 3 .
Fig. 3. Scatterplots illustrating the relation of MARST to MAAT (top left) and of the difference between MARST and MAAT to PISR (top right) for the 57 MARST values.In the lower panels, the model residuals (mean = 0 • C, standard deviation = 1.52 • C) of the rock model are plotted against PRECIP and SEASONAL to visualize possible relations of these two variables.

Fig. 7 .
Fig. 7. Prediction values for the two first models fromTable 3 calculated for a randomly selected probability range of 0.475-0.525for intact rock glacier occurrence.Black crosses: debris model without cPRECIP as explanatory variable, red bubbles: debris model including cPRECIP as explanatory variable.

Fig. 8 .
Fig. 8. Spatial density of rock glacier occurrence in inventories with high (IGUL Tecino, Switzerland; mean PRECIP = 1900 mm) and low (GEOL, Trentino, Italy; mean PRECIP = 1072 mm) precipitation.The relevant area is calculated for areas above 2000 m of elevation excluding glaciered areas (data provided by Paul et al., 2009) and steep slopes (slope angle < 50 • ).

Table 1 .
Overview of data used for model calibration (RG Rock glacier; AT Austria, CH Switzerland, DE Germany, FR France, IT Italy).

Table 2 .
Summary statistics of the randomly sampled points representing potential explanatory variables for the debris model (Q25 lower quantile, Q75 upper quantile).

Table 3 .
Model coefficients (and standard errors in parentheses) of debris models using different sets of explanatory variables, and the corresponding goodness-of-fit statistics.Debris model 1 was chosen as final model.(The units for the explanatory variables are given in square brackets for each coefficient.)< 0.05, * * < 0.01, * * * < 0.001 *

Table 4 .
Summary statistics and Pearson correlations between MARST and potential explanatory variables for all MARST locations for the rock model (Q25 lower quantile, Q75 upper quantile).

Table 5 .
Model coefficients (and standard errors in parentheses) of rock models using different sets of explanatory variables, and the corresponding goodness-of-fit statistics.Rock model 2 was chosen as final model.(The units for the explanatory variables are given in square brackets for each coefficient.)< 0.05, * * < 0.01, * * * < 0.001 *

Table 6 .
Variance components used for combining the rock and debris models.