Permafrost Distribution Modeling in the Semi-Arid Chilean Andes

Mountain permafrost and rock glaciers in the dry Andes are of growing interest due to the increase in human activities in this remote area. Empirical models of mountain permafrost distribution based on the spatial analysis of intact and relict rock glaciers and mean annual air temperature (MAAT) have been established as a tool for regional-scale assessments of permafrost favorability across entire mountain ranges; however, this kind of model approach has never been applied for a 10 large portion of the Andes. In the present study, this methodology is applied to map permafrost favorability throughout the semi-arid Andes of central Chile (29° S-32° S), excluding areas of exposed bedrock. After spatially modeling MAAT distribution from scarce temperature records (116 station-years) using a linear mixed-effects model (LMEM), a generalized additive model (GAM) was built to model the activity status of 3524 rock glaciers. A Permafrost Favorability Index (PFI) was obtained by adjusting model predictions for conceptual differences between permafrost and rock glacier distribution. The 15 results indicate that the model has an acceptable performance (median AUROC: 0.76). Conditions highly favorable to permafrost presence (PFI  0.75) are predicted for 1051 km2 of mountain terrain, or 2.7 % of the total area of the watersheds studied. Favorable conditions are expected to occur in 2636 km2, or 6.8% of the area. Especially in the Elqui and Huasco watersheds in the northern half of the study area, where a substantial surface portion (11.8 % each) was considered to be favorable for permafrost presence, while predicted favorable areas in the southern Limarí and Choapa watersheds are mostly 20 limited to specific sub-watersheds. In the future, local ground-truth observations will be required to confirm permafrost presence in favorable areas, and to monitor permafrost evolution under the influence of climate change.


Introduction
Mountain permafrost is widely recognized as a phenomenon that may influence slope stability (Kääb et al., 2005;Harris et al., 2009) and hydrological systems (Bommer et al., 2010;Caine, 2010;Haeberli, 2013), which poses a challenge to economic development in high mountains (Taillant, 2015).Its extent and characteristics and its response to climate change are of major concern, but the present knowledge is still very limited especially in remote environments such as the high Andes.Some cases of rock glacier destabilization have been observed in the Andes (Bodin et al., 2012), raising the question of sensitivity of ice-rich permafrost to climate warming.In addition, the cryosphere of the central Andes is a growing interest zone due to the increase in human activities related to mining and border infrastructure as well as societal concern about possible environmental damages and new environmental acts (Brenning, 2008;Brenning and Azócar, 2010a).
Empirical models describing the distribution of mountain permafrost based on geomorphological permafrost indicators and topographic and climatic predictors are a simple yet effective approach toward a first assessment of its distribution at a regional scale (Lewkowicz and Ednie, 2004;Janke, 2005;Boeckli et al., 2012a, b;Sattler et al., 2016).In the Andes, the boundaries of permafrost have been determined based on global climate models at coarse resolution scales (Gruber, 2012;Saito et al., 2015); however, regionally calibrated permafrost mapping based on geomorphological evidence has not been attempted previously in the Andes.The recent availability of accurate global digital elevation models (DEMs), the compilation of new rock glacier inventories, and the improved access to meteorological data provide a unique opportunity to model permafrost distribution at high resolution for a large portion of the Andes.
The goal of this research was therefore to create an index map of potential permafrost distribution in the semi-arid Chilean Andes between ∼ 29 and 32 • S based on statistical modelling of topographic and climatic conditions and rock glacier activity status.For this purpose, a rock glacier inventory was compiled in order to obtain a variable indicating the presence and absence of permafrost conditions according to rock glacier activity status.Considering the potential incoming solar radiation (PISR) and the mean annual air temperature (MAAT) as potential predictors, they are modelled as predictor variables for the favourability for permafrost occurrence.A generalized additive model (GAM) was then used to map a permafrost favourability index (PFI) in debris surfaces within the study area.

Study area
The study area comprises a large portion of the semi-arid Chilean Andes, covering from north to south the upper sections of the Huasco (9766 km 2 ), Elqui (9407 km 2 ), Limarí (11 683 km 2 ), and Choapa (7795 km 2 ) river basins between ∼ 28.5 and 32.2 • S. Regarding the altitudinal distribution of the studied watersheds, about 35 % of the terrain is located above 3000 m a.s.l., with the Huasco and Elqui basins bearing the highest elevations (i.e.Cerro de Las Tórtolas is 6160 m a.s.l.).Their median elevations are 2995 and 2536 m a.s.l., respectively, compared to median values of approximately 1300 m a.s.l. for the Limarí and Choapa basins (Fig. 1), limiting the presence of glaciers, snow, and permafrost in the latter two basins.
Population is scarce, but lowland populations rely on water resources from the high-elevation headwater areas (Gascoin et al., 2011).Climatically, it is located in a transition zone between arid and semi-arid climates where the presence of the South Pacific anticyclone inhibits precipitation, favouring clear skies and high solar radiation.Precipitation at high elevations almost exclusively occurs as snow as it is concentrated between May and August, i.e. austral winter (Gascoin et al., 2011).Regarding the recent climatic trends in this part of the Andes, atmospheric warming reaches 0.2 to 0.4 • C decade −1 (period 1979-2006;Falvey and Garreaud, 2009), whereas precipitation exhibits no trend (∼ 1986-2005;Favier et al., 2009); these trends are, however, subject to large uncertainties due to limited data availability.
Vegetation is nearly absent above 3500 m a.s.l.(Squeo et al., 1993).Glaciers are present where summit elevations are higher than the modern equilibrium line altitude (ELA) of glaciers, which is located at about 5000 m a.s.l. at ∼ 30 • S and also depends on local climatic anomalies (Kull et al., 2002;Azócar and Brenning, 2010).Most of the surface ice bodies in the study region correspond to small glaciers or "glacierets" (< 0.1 km 2 ).In contrast, 2 % of the glaciers in the study region are greater than 1 km 2 , mostly located in the Elqui watershed at and near Cerro Tapado and in the upper Huasco watershed (Nicholson et al., 2009;Rabatel et al., 2011;UGP UC, 2010).Active and inactive rock glaciers as ice-debris landforms in the study area are present above ∼ 3500 m a.s.l.near the southern limit of the study region and above ∼ 4200 m a.s.l.towards its northern limit (Brenning, 2005a;Azócar and Brenning, 2010; UGP UC, 2010).Rock glaciers with ground ice and current movement have been detected in various cases (UGP UC, 2010 ;Monnier, et al., 2011Monnier, et al., , 2013;;Monnier andKinnard, 2013, 2016;Janke et al., 2015).While this pattern broadly follows the latitudinal trend of the 0 • C isotherm altitude (ZIA) of present-day MAAT, in the southern part of the study region rock glaciers also exist at positive regional-scale MAAT levels (Brenning, 2005b;Azócar and Brenning, 2010).Rock glacier presence is furthermore constrained by topographic suitability for rock glacier development, and the scarce availability of quaternary glacial sediments may further limit their development in the southern Limarí and northern Choapa watersheds (Azócar and Brenning, 2010).Recent geophysical investigations on rock glaciers in the semi-arid Andes revealed complex internal structures and landform evolution as well as highly variable ice contents, while confirming the presence of ice-rich permafrost in intact rock glaciers and the general relationship between landform appearance and presence of ground ice in this region (Monnier and Kinnard, 2015a, b;Janke et al., 2015).

Methods
In adopting the empirical, regional-scale approach to permafrost favourability modelling of Boeckli et al. (2012a), several preprocessing steps were necessary to compile the necessary data (Fig. 2).To create a geomorphological indicator variable representing likely permafrost presence and absence, rock glaciers were mapped and then classified according to their likely activity status (Sect.3.1.1).Further, predictor variables such as MAAT and PISR were calculated (Sect.3.1.2),a linear mixed-effects model (LMEM) was used to regionalize point measurements of MAAT at weather stations to the landscape scale using a DEM.Finally, the GAM was chosen to create a PFI for the study area (Sect.3.2.),extending the generalized linear modelling approach chosen by Boeckli et al. (2012a).(Birkeland, 1973;Haeberli et al., 2003).They have a tongue or lobe shape, with ridges and furrows on their surface that are indicative of their present or past deformation.Rock glaciers were identified following the criteria of classification proposed previously by Barsch (1996) and Roer and Nyenhuis (2007) and slightly adapted by Azócar (2013).Frequently, active and inactive rock glaciers (grouped together as intact forms in this study) have a steep front with visual unstable rocks (see Fig. S1 in the Supplement); in contrast, an irregular and collapsed surface due to thawing of ice commonly indicates that rock glacier is in its relict form (Putman and David, 2009).Because of uncertainties in this classification, the "intact" category was considered as a whole, including active and inactive landforms.Rock glaciers were digitized as points located at the rock glacier toe using an on-screen map scale of 1 : 7000.Due to the lack of ground observation data about the dynamic status of rock glaciers, a validation process cannot be accomplished in this research.Nevertheless, the authors have performed measurements of rock glacier dynamics and active layer status within and near the study region (UGP UC, 2010) and have knowledge of additional geophysical evidence of several rock glaciers in the dry Andes (Janke et al., 2015).Although these observations cannot be used for independent validation since they were previously known to the authors, their direct evidence of rock glacier dynamic status and ground ice presence is generally consistent with our assessment based on remote-sensing imagery.

Regionalization of MAAT
Air temperature in mountain areas is mainly controlled by latitude, altitude, and topography (Barry, 1992;Whiteman, 2000), considering in particular the effects of global atmospheric circulation patterns, global as well as local differences in PISR, and adiabatic temperature lapse rates.In order to regionalize (or interpolate) weather station data, most studies therefore utilize regression or hybrid regressioninterpolation approaches with a combination of predictors representing elevation, geographic position, and local climatic phenomena (e.g.cold air pools), depending on data availability and size of the study region (Lee and Hogsett, 2001;Hiebl et al., 2009;Lo et al., 2011).
Since MAAT strongly varies latitudinally and altitudinally throughout the study region, a statistical regression model with latitude and elevation as predictors was used for the regionalization of this variable.Additional regression predicwww.the-cryosphere.net/11/877/2017/The Cryosphere, 11, 877-890, 2017 tors were not used since the limited sample size would not allow for an accurate estimation of their coefficients.Considering the small number of and large distances between weather stations, no attempt was made to spatially model or interpolate the residuals.The heterogeneous temporal coverage of different weather stations made it necessary to choose a model that is able to account for year-to-year temperature variations instead of regionalizing observations of MAAT directly.We therefore chose a LMEM in order to relate annual average temperature (AAT) data to latitude and elevation as so-called fixed-effects variables while accounting for yearto-year variability by including a random-effect term for the year of observation.

Response and predictor variables
To build a response variable representing air temperature conditions, AAT was calculated using data from eleven weather stations for a time period of between 1 and 30 years .Temperature data for eight weather stations were provided by Chile's water administration, Dirección General de Aguas (DGA), and additional data were obtained from mining projects distributed throughout the study area.AAT for a particular year was calculated as the arithmetic average of that year's mean monthly temperatures.Since the consistency of elevation references (e.g.above sea level) of available weather station elevations was in doubt, consistent elevation values were extracted from ASTER GDEM.Only weather stations located above 2000 m a.s.l. and at least 100 km inland from the coast were selected in order to reduce oceanic influence (Hiebl et al., 2009) and focus on mountain climate.Moreover, stations located slightly north and south of the study area were also included to better represent latitudinal changes.Based on these criteria, 116 station years from 11 weather stations were available (Table 1).Homogenization of temperature data was not possible due to the short length and partial overlap of most available time series.
The predictors' elevation (metres) and northing (coordinate in metres) were used, taking elevation data from the ASTER Global Digital Elevation Model (GDEM) version 2 with a 30 m × 30 m resolution and a vertical precision around 15 m (Tachikawa et al., 2011).Easting was not considered due to the limited extent of the study area in west-east direction and possible confounding with the dominant altitudinal MAAT gradient.
LMEMs are appropriate models when data are organized in hierarchical levels and observations are therefore grouped (Pinheiro and Bates, 2000).Unlike ordinary linear regression models, LMEMs are able to account for the dependency among observations that arises from grouping (Twisk, 2006).In a LMEM, the predictors can contain random and fixed effects.Random effects can be thought of as additional error terms or variation in coefficients that are tied to different grouping levels.In a climatic context, this may include a year's overall departure from long-term mean temperature, for example.Fixed effects are comparable to ordinary linear regression predictors whose coefficients do not vary by group.
In this study, fixed and random effects were estimated by maximum likelihood (ML).AAT data were considered as being grouped by year (30 groups, one for each year).Thus, the following model specification of the statistical temperature distribution model for AAT at a weather station i in year j was used: (2) In this model, AAT ij is the AAT at a particular weather station i in a particular year j , b 0 represents the overall mean intercept of AAT, which may vary from year to year by a random amount u 0j , whose variance is estimated during model fitting.ε ij denotes the residual error for each weather station and year.Thus, the model error is split into two com-ponents: the residual variation between years (u 0j ) and the residual variation among weather stations within a particular year (ε ij ).
In order to predict spatially varying MAAT, year-to-year and station-specific residual variation were disregarded since their statistical expected value is zero, and therefore only the fixed-effects portion is used to predict the expected value of MAAT at any location x with known elevation and latitude within the study region: (3) The overall fit of the LMEM was evaluated by examining the residual standard error (RSE) in • C. The LMEM implementation of the "nlme" package in R was used (Pinheiro and Bates, 2000;R Core Team, 2012).
The PISR, defined as the annual sum of direct and diffuse incoming solar radiation, was calculated using SAGA GIS version 2.1.0(Conrad et al., 2015) and ASTER GDEM data.PISR was calculated at 10-day intervals with a half-hour time resolution between 04:00 and 22:00; latitudinal effects were accounted for.To simulate the extremely clear and dry skies, a lumped atmospheric transmittance of 0.9 was used (Gates, 1980).

Statistical permafrost favourability model
In order to create a response variable Y indicating likely permafrost presence/absence, intact rock glaciers were used as indicators of permafrost presence (Y = 1) and relict rock glaciers as indicators of permafrost absence (Y = 0).A GAM was applied to model this response variable (Sect.3.2.1),and several adjustments were implemented a posteriori in order to account for biases related to the use of rock glaciers as indicators (Sect.3.2.3).
MAAT predicted by the LMEM and DEM-derived PISR were used as predictor variables.Due to the lack of snow cover information, annual PISR was preferred over PISR because seasonal PISR is highly correlated with snow cover periods.Nevertheless, the statistical interaction (Hosmer and Lemeshow, 2000) of PISR and MAAT was included in the model; this interaction term is capable of representing temperature-dependent (altitudinal) variation in the influence of PISR on permafrost occurrence, thus capturing the effects of snow cover duration on absorbed solar radiation (Brenning and Azócar, 2010b).

Generalized additive model
A GAM was utilized as the statistical model of permafrost favourability, and its spatial predictions are referred to as the PFI.This type of model has been applied in environmental sciences, including ecology (Guisan et al., 2002), forestry (Janet, 1998), periglacial geomorphology (Brenning and Azócar, 2010b), and landslide research (Goetz et al., 2011).A GAM can be described as a generalized linear model (GLM) in which some or all of the linear predictors are specified in terms of smooth functions of predictors (Hastie and Tibshirani, 1990;Wood, 2006).Like the GLM, the GAM can be applied to model binary categorical response variables such as the presence (Y = 1) versus absence (Y = 0) of permafrost.Specifically, the probability p = P (Y = 1|MAAT, PISR) of permafrost occurrence given known values of the predictors MAAT and PISR can be modelled in a GAM with a logistic link function as where β 0 is the intercept and f 1 and f 2 are smoothing functions for the predictors MAAT and PISR.The ratio p/(1−p) is referred to as the odds of permafrost occurrence, and its logarithm as the logit.
When an interaction term MAAT and PISR is included in the above equation, the model becomes where f is a bivariate smoother.In this model, the relationship between PISR and the response may depend on MAAT, and conversely the relationship between MAAT and the response may vary with PISR.
The GAM has the advantage over the GLM that it increases model flexibility by fitting nonlinear smoother functions to the predictors (Wood, 2006).In this study, the smoothing terms were represented using a local regression or "loess" smoother with 2 degrees of freedom.In this method local linear regressions are fitted using a moving window that describes a smoothly varying relationship between predictor and response variables.One of the advantages of this method is that assumptions about the form of the relationship are avoided, allowing the form to be discovered from the data.

Model evaluation
The overall performance of the GAM was evaluated using the area under the receiver operating characteristic (ROC) curve (AUROC).In the present context the curve shows the probability of detecting observed permafrost occurrences (sensitivity) and absences (specificity) for the whole range of possible decision thresholds that could be used to dichotomize predicted odds into permafrost presence/absence (Hosmer and Lemeshow, 2000).The AUROC can range from 0.5 (no separation) to 1 (complete separation of presence and absence by the model).
The performance of the predictive model was evaluated using spatial cross validation where testing and training data sets are spatially separated (Brenning, 2012).k-means clustering was used to partition the subsets randomly into k = 10 equally sized subsamples for k-fold spatial cross validation.Cross validation was repeated 100 times.

Model adjustments
In order to distinguish between two surfaces with different thermal properties, bedrock and debris-covered areas were identified.This was necessary since the present model was based on rock glaciers as evidence of permafrost conditions in areas with rock debris as the surface material.Therefore, the permafrost model cannot be extrapolated to non-debris surfaces such as steep bedrock slopes, which have different thermophysical properties and limited winter snow cover.A slope angle ≥ 35 • in the ASTER GDEM was used to identify probable bedrock surfaces and to exclude these areas from the predictive modelling following Boeckli et al. (2012b).Although bedrock exposures may be found on gently inclined terrain, we disregarded these and tentatively classified slopes less steep than 35 • as debris zones.
While rock glaciers are commonly used as indicators of permafrost conditions in mountain regions, they are also known to overestimate permafrost distribution in adjacent debris-covered areas for several reasons, including (1) a cooling effect in coarse blocky material that is often present on the surface of rock glaciers, (2) the long-term creep of rock glaciers downslope into non-permafrost terrain, and (3) a delayed response of ice-rich permafrost to climatic forcings (Boeckli et al., 2012b).Boeckli et al. (2012b) suggest that the last two effects can be compensated by the use of a temperature offset term.However, the first effect cannot be easily accounted for due to the lack of information about surface characteristics of rock glaciers.Nevertheless, qualitatively speaking, in our experience surface grain size of rock glaciers in the semi-arid Andes is not necessarily much different from surrounding areas, and the complexity of relationships between the structure of the surface layer and thermal characteristics makes it difficult to estimating the size of this potential bias (Brenning et al., 2012).
In this research, a temperature offset term was applied to estimate the above-mentioned effects (2) and (3).It was approximated by the mean altitudinal extent of rock glaciers, which was derived from the mean slope angle and mean length of rock glaciers inventoried by Azócar (2013) and UGP UC (2010).The calculated mean altitudinal extent was 89 m, which corresponds to a temperature offset of −0.63 • C, using the lapse rate of −0.0071 • C m −1 obtained in this work.This temperature offset was chosen and added to MAAT values (renamed as "adjusted MAAT") prior to permafrost model fitting.In other words, the adjusted MAAT used for model fitting was made 0.63 • C cooler, while us-ing the warmer original MAAT for prediction; this effectively shifts the high PFI values towards higher elevations.

Rock glacier inventory
Information on 3575 rock glaciers was compiled based on existing inventories and the identification of additional rock glaciers in the study area.Of these, 1075 were classified as active, 493 as inactive, 343 as intact, and 1664 as relict forms.Rock glaciers are most abundant in the Elqui (n = 681), Limarí (n = 486) and Huasco (n = 424) watersheds; in contrast, they are less abundant in Choapa watershed (n = 324).The majority of rock glaciers (∼ 60-80 %) are located below the present ZIA obtained in this study from the LMEM of temperature.However, 37 % of active, 21 % of inactive, 26 % intact, and 15% of relict rock glaciers are located above the ZIA.These proportions vary considerably between watersheds.In the Huasco and Elqui watersheds, nearly 50 % of all active rock glaciers are located at negative MAAT; in contrast, in the Limarí and Choapa watersheds in the southern part of the study contain less than 20 % (Fig. 3).
Of the 3575 rock glaciers in the inventory, 51 were removed from the data set due to their isolated or unusual location in order to prevent these from becoming influential points.Of the remaining 3524 rock glacier observations, 1909 were used as indicators of permafrost and 1615 as indicators of the absence of permafrost.In general, MAAT and PISR values at sites with permafrost were lower than the sites without permafrost.Exploratory analysis showed a clear trend towards a greater presence of intact rock glaciers at lower MAAT and lower PISR.

Mean annual air temperature model
Model coefficients indicate an average temperature lapse rate of −0.71 • C per 100 m while accounting for latitude and year-to-year variation (95 % confidence interval: −0.68 to −0.74 • C per 100 m; Table 2).With every 200 km of northward distance, the AAT increased on average by 1.6 • C (95 % confidence interval: 1.4-1.8• C).Therefore, there is a 4 • C MAAT difference between the northern and southern border of the study area at equal elevations.Based on the fitted model, regional-scale MAAT was predicted using the equation  the ZIA was situated at ∼ 4350 m a.s.l. in the northern (29 • S) section and dropped southward to ∼ 4000 m a.s.l. at 32 • S (Fig. 4).

Statistical permafrost favourability model
In order to understand the effect of the interaction between PISR and MAAT in the permafrost favourability model, it is necessary to examine the modeled effect of MAAT at sites with different levels of PISR.We will compare and contrast sites with regional-mean PISR with extremely sunny sites with PISR 2 standard deviations above regional-mean PISR, looking at how much less favourable sites with a MAAT of +1 • C are relative to a MAAT of 0 • C. At regional-mean PISR, locations with a MAAT of +1 • C are, on average, associated with ∼ 33% lower predicted odds (or relative probabilities) of permafrost occurrence.At extremely sunny sites, in contrast, the model predicts ∼ 73 % lower odds at +1 • C. Similarly, a large PISR had a greater estimated effect at higher MAAT levels than at lower MAAT.At −1 • C MAAT, a PISR value 1 standard deviation above the average PISR was associated with approximately 27 % lower odds of permafrost occurrence compared to average PISR.However, the same PISR difference at +1 • C MAAT was associated with an estimated 57 % decrease in the odds of permafrost occurrence.
Based on spatial cross-validation estimation of the GAM's overall accuracy using a decision threshold of 0.5, 64 % of observations indicating permafrost conditions were correctly classified by the model (sensitivity: 58 %; specificity: 78 %).Spatially cross-validated AUROC confirmed this acceptable performance (median AUROC: 0.76), and the comparison with non-spatial cross-validation results (median AUROC: 0.75) suggests that the GAM generalized well from the training data.

Spatial distribution of PFI
Considering a PFI ≥ 0.5 and excluding steep bedrock and glacier surfaces, favourable conditions for permafrost occurrence were inferred for ∼ 6.8 % of the study area, or 2636 km 2 .Considering only highly favourable conditions with a PFI ≥ 0.75, the potential permafrost area would be limited to 2.7 % of the area, or 1051 km 2 (Table 3).The largest potential permafrost areas were concentrated in the Huasco and Elqui watersheds where the PFI ≥ 0.5 covers more than 10 % of each watershed (1150 km 2 in the Huasco; 1104 km 2 in the Elqui).In the Limarí and Choapa watersheds, areas with PFI ≥ 0.5 represented less than 3 % of each watershed (Limarí: 217 km 2 ; Choapa: 192 km 2 ).
The spatial distribution of the predicted PFI in the study area is depicted in Fig. 5.In general, the potential permafrost areas tend to decrease southward.More favourable conditions were concentrated in the highest areas in the central part of the study area (e.g.Pascua-Lama area, Cerro El Toro 6168 m a.s.l., Las Tórtolas 6160 m a.s.l., Olivares 6216 m a.s.l.).In contrast, lower scores (< 0.5) were associated with lower hillslopes and valley bottoms.

Rock glacier inventory
The rock glacier inventory prepared for this work expands our previous knowledge of rock glacier distribution in the semi-arid Andes as it adds rock glaciers to previous compilations that were mostly based on lower-quality imagery or statistical sample surveys (Azócar and Brenning, 2010a;Nicholson et al., 2009;UGP UC, 2010).In comparison to the recent inventory of rock glaciers realized by UGP UC (2010) in the Elqui, Limarí, and Choapa watersheds and Huasco by Nicholson et al. (2009), the present inventory increases the number of known active rock glaciers from 621 to 933 (increase 60 %), inactive rock glaciers from 151 to 415 (increase 275 %), and intact rock glaciers from 135 to 249 (increase 184 %) within these watersheds.This has been possible since to rock glaciers are recognized using high-resolution satellite imagery in comparison to previous studies.Moreover, this work added relict rock glaciers (n = 1664), which were missing from previous inventories (Nicholson et al., 2009; UGP UC, 2010) and were mapped for the first time systematically throughout this study area.This has been possible because in the current research, rock glaciers are recognized using images with better resolution (finer than 2 m × 2 m) than in previous studies, allowing the identification of small landforms (below 0.1 km 2 ).The considerable proportion of intact rock glaciers located at (present, regional-scale) MAAT above 0 • C is confirmed by this study (Brenning, 2005;Azocar and Brenning, 2010).It has previously been attributed to the delayed response of ice-rich permafrost to MAAT increases since the Little Ice Age, rock glacier advance towards lower elevations, and the possible preservation of ice-rich mountain permafrost in favourable local topoclimatic conditions (Azócar and Brenning, 2010;Boeckli et al., 2012b).The interaction of permafrost with buried glacier ice from Holocene/Little Ice Age glacier advances may play an additional important role in the development of rock glaciers in these areas (Monnier and Kinnard, 2015a).

Temperature distribution model
The temperature distribution model indicated a present (1981-2010 mean) ZIA at ∼ 4350 m a.s.l. at 29 • S, which drops southward to ∼ 3900 m a.s.l. at 32 • S.Although the result cannot be directly compared with other studies that may refer to different time periods or provide only general statements, our ZIA is somewhat higher than previously thought (Brenning, 2005a; ZIA ∼ 4000 m at 29 • S, ∼ 3750 at 32 • S).
The temperature lapse rate obtained in this study (−0.71 • C per 100 m) is reasonably close to the average temperature decrease in the free atmosphere (∼ −0.6 • C per 100 m; Barry, 1992).
The precision (RSE) of the MAAT distribution model of 0.8-1.08 • C (with 95 % confidence) highlights important uncertainties in this key predictor of permafrost favourability in this region.Sensitivity analyses using leave-one-stationout cross validation confirmed the validity of this value.As example of the change in predicted MAAT in leave-onestation-out cross validation compared to full data set ( • C), using 4200 m a.s.l. as reference elevation point indicate that largest individual changes in predicted MAAT were obtained when leaving out the temperature data from Los Bronces (−0.73 • C change) and Embalse La Laguna (−0.56 • C); however, all these values are smaller than the model's RSE.
This level of uncertainty is roughly comparable to MAAT products used in permafrost distribution models for the European Alps (RSE < 1 • C in Hiebl et al., 2009) and at a global scale (RSE ≈ 1 • C in Gruber, 2012), although such comparisons may only provide very general guidance.Even though our research approach focused on incorporating all available high-elevation temperature data into a hierarchic regression model, the limited availability of temperature data can affect permafrost prediction.
Considering the overrepresentation of station years from valley locations, which tended to be associated with positive residuals, overall the MAAT regionalization may have a warm bias at high elevations and on the upper slopes.Such bias would, however, be substantially smaller than the residual standard error of 0.93 • C due to averaging effects.As a consequence, permafrost favourability may be underestimated at the highest elevations or on the upper slopes.

Permafrost favourability model assumptions
The use of rock glaciers as the empirical foundation for permafrost favourability models requires researchers to make several assumptions, not all of which are equally known in the literature.On the one hand, altitudinal biases related to the movement and characteristics of rock glaciers have been known (Boeckli et al., 2012b), and adjustments are available and have been used in this work, albeit not in all studies relying on rock glaciers (Janke, 2005;Deluigi and Lambiel, 2011;Sattler et al., 2016).However, additional regional calibration will be necessary in the future in order to obtain more precise adjustments.Geophysical soundings or direct borehole evidence would be suitable for this, at least if placed representatively according to a meaningful sampling design.
On the other hand, the transfer of relationships between permafrost presence and predictor variables (e.g.MAAT, PISR) from rock glaciers to non-rock glacier areas also requires the additional assumption that these relationships (e.g.model coefficients) remain more or less the same in debris areas as within rock glaciers.We will refer to this assumption as the transferability assumption.This assumption has previously not been made explicit in the literature despite the frequent application of rock-glacier-based permafrost distribution models.
In combining permafrost models based on rock temperatures and rock glacier activity status, Boeckli et al. (2012a) pointed to a mathematical relationship between (probit) presence/absence models on the one hand and linear regression models on the other (see Sect. 3.1 of Boeckli et al., 2012a).This mathematical relationship may shed some light on the transferability assumption.Based on this relationship between probit and linear regression models, a sufficient condition for the transferability assumption is that (1) ground temperatures in both model domains (i.e.rock glaciers and other debris surfaces) show similar relationships with MAAT and PISR and (2) such linear regressions of ground temperature would have similar residual standard deviations, or precisions.Evidence for or against this is, unfortunately, scarce, since sufficiently replicated ground temperatures have only been measured at shallow depths, and previous studies have paid little or no attention to such differences between rock glaciers and debris surfaces.In the semi-arid Andes, Apaloo et al. (2012) andCEAZA (2012) examined regression relationships between near-surface ground temperatures (NGST) and topoclimatic predictors including elevation (as a proxy for MAAT) and PISR both within and outside of rock glaciers.These studies showed no convincing evidence of differences in ground temperature between ice-debris landforms and other debris surfaces under otherwise equal conditions.While interaction terms of ice-debris landforms with elevation or PISR were not examined in these studies, a reanalysis of data from Apaloo et al. (2012) showed no evidence of relationships between NGST and air temperature or NGST and PISR varying between rock glaciers and other debris surfaces.Thus, while the transferability assumption merits further evaluation in future studies, we have no concrete evidence against this assumption in this climatic setting.

Permafrost favourability model
Debris areas with a permafrost favourability score ≥ 0.5 in our study cover a spatial extension of 6.8 % (2636 km 2 ) of the study area.Although the model includes the main factors that control the regional distribution of permafrost in the semi-arid Chilean Andes, such as the temperature and the potential amount of solar radiation in relation to the elevation and latitude (Azócar and Brenning, 2010), the permafrost model does not account for the effect of specific local environmental factors in debris areas, such as spatially variable soil properties and the effects of long-lasting snow patches that can influence ground thermal regimes locally (Hoelzle et al., 2001;Apaloo et al., 2012).Therefore, all these local factors must be considered when interpreting PFI values (Boeckli et al., 2012b).In the Andes of Santiago directly south of the study region, for instance, 30 days of additional snow cover in spring/summer were with 0.1-0.6 • C lower mean ground surface temperatures (MAGST), while openwork boulder surfaces were about 0.6-0.8• C cooler than finer debris (Apaloo et al., 2012).
In summary, our model suggests that in areas with PFI ≥ 0.75, permafrost will occur in almost all environmental conditions; in contrast, in areas where PFI ranges between 0.5 and 0.75, permafrost will be present only in the favourable cold zones described before.In areas with PFI < 0.5, permafrost may be present in exceptional local environmental circumstances.
Comparisons of modelled permafrost distribution with independent ground-truth observations is naturally difficult due to the bias toward permafrost presence sites and sampling design (Boeckli et al., 2012b).Such comparisons were therefore not performed in this study.Direct permafrost observations outside of rock glaciers are still scarce in the semi-arid Andes and mostly limited to unpublished data from a mining context.Additional systematic ground-truthing is therefore required for a quantitative assessment of permafrost extent and to reduce model uncertainties (Lewkowicz and Ednie, 2004).Model-model comparisons are therefore currently the only means for assessing uncertainties in permafrost distribution.
The permafrost zonation index (PZI) of Gruber ( 2012) is an independent global-scale empirical modelling effort at a 1 km × 1 km resolution based on downscaled reanalysis data and within-pixel relief.For comparison, PFI was resampled to PZI resolution (1 km) by averaging all PFI pixels that fall within a given PZI pixel (Azócar, 2013).While we do not suggest that equivalent PZI and PFI values should be interpreted the same way, it should be pointed out that potential permafrost areas with a PZI ≥ 0.75 are substantially smaller (209 km 2 , including exposed bedrock areas) than the areas with PFI ≥ 0.75 predicted by our approach (1051 km 2 ), even though the latter only considers debris surfaces.Even PZI ≥ 0.50 covers only 653 km 2 .Since vast areas with acrock glaciers present PZI values below 0.50, we conclude that the PZI is a very conservative measure of permafrost favourability in this region.However, more research is needed to confirm the appropriateness of the bias adjustment used in this and earlier rock glacier based studies (Boeckli et al., 2012) and thus the adequate calibration of such models.
Overall, compared to previous empirically based studies (Lewkowicz and Ednie, 2004;Janke, 2005;Boeckli et al, 2012a, b;Sattler et al., 2016), the permafrost modelling approach used in this work integrates the regionalization of MAAT distribution from scarce data with more flexible, nonlinear models of potential permafrost distribution in comparison with previous research approaches.
In this study, the discrimination of rock glacier classes based on PFI values results in an AUROC of 0.76 which is similar to the AUROC value obtained in recent permafrost modelling studies in the Alps (Boeckli et al., 2012b).Other studies have reached AUROC values ≥ 0.9 (Deluigi and Lambiel, 2011;Sattler et al., 2016), which means nearly perfect separation.It is worth to mentioning that excessively high AUROC values may result in numerical difficulties in the estimation of logistic regression coefficients (Hosmer and Lemeshow, 2000), and in the case of Sattler et al. (2016) the high value may be explained by the omission of inactive rock glaciers as a permafrost landform located in topographic conditions between active and relict ones.

Permafrost distribution and effects of climate change
According to the PFI model, mountain permafrost in the semi-arid Chilean Andes between 29 and 32 • S is widespread above ∼ 4500 m a.s.l. and more scattered between ∼ 3900 and 4500 m a.s.l.Permafrost areas near the lower limit of permafrost distribution can be more sensitive to degradation processes due to the possible effects of climate change (Haeberli et al., 1993) considering also its location at or near the present ZIA.A rise in air temperature can potentially lead to permafrost thaw and long-term degradation of icerich frozen ground (e.g.intact rock glaciers) as well as the acceleration and sudden collapse of rock glaciers (Schoeneich et al., 2015).In addition, this warming could lead to geotechnical problems related to high-altitude road or infrastructure (Brenning and Azócar, 2010a;Bommer et al., 2010).
In this context, PFI maps can serve as a first resource to assess permafrost conditions and uncertainties in mountain research and practical applications such as infrastructure planning (Boeckli et al., 2012b).

Conclusions
The statistical permafrost distribution model proposed here provided more detailed, locally adjusted insights into mountain permafrost distribution in the semi-arid Chilean Andes compared to previous coarser-resolution results from a global permafrost distribution model such as PZI (Gruber, 2012).General climatic and topographic patterns proved to be useful for mapping broad permafrost distribution patterns while local environmental factors (such as substrate properties or snow cover duration), which are not included in the model, could determine permafrost presence locally.Data from rock glacier inventories combined with topographic and topoclimatic attributes were used to model the probability of permafrost occurrences in the semi-arid Chilean Andes.The GAM is particularly suitable for modelling these relationships due to its ability to incorporate nonlinear relationships between predictor and response variables.
Using rock glaciers as indicators of permafrost conditions in areas with debris as surface material, the result of the permafrost model cannot be extended to other types of surface covers.Therefore, future studies should address this limitation in order to determine potential permafrost areas in steep bedrock.Furthermore, the effect of a delayed response of rock glaciers with high ice content to climate forcings should be considered in future analyses along with systematic direct ground-truthing of permafrost occurrence outside of rock glaciers.
Rock glacier landforms can be used as a variable indicative of permafrost distribution in mountain areas; however, the lack of permafrost observation outside of the boundaries of rock glaciers must be addressed in future studies.Nevertheless, outside of the rock glaciers boundaries, there is not a systematic way to infer permafrost presence or absence in large areas; therefore, for general studies of mountain permafrost distribution, rock glaciers are possibly one of best proxy variable to infer permafrost outsides of boundaries of rock glacier areas.
However, the results of the MAAT model show that LMEMs can be appropriate to determine temperature distribution with scarce and heterogeneous temperature records from weather stations.Moreover, the results of the MAAT model can provide valuable data to study other climatically driven Earth surface features such as glacier and vegetation patterns.
Finally, the findings of this research contribute mainly to improving the general knowledge about permafrost distribution in the Andes, providing valuable information to government and economic sectors as a starting point for additional local site investigations.Considering the uncertainties inherent in regional-scale modelling, local site conditions such as surface material, snow cover duration, or topoclimatic conditions should be taken into account when interpreting the present permafrost favourability index.Additional research, in particular ground-truthing using a meaningful sampling design, is necessary in order to refine the present model and eliminate possible biases.
Data availability.The PFI and MAAT prediction maps are available for downloading and visualizing at www.andespermafrost.com.Moreover, the rock glaciers inventory is accessible for downloading.Alternatively, the data can be downloaded through Pangaea server (doi:10.1594/PANGAEA.859332).
The Supplement related to this article is available online at doi:10.5194/tc-11-877-2017-supplement.

Figure 1 .
Figure 1.Overview map of the study area in the semi-arid Andes.Watersheds studied are indicated in light green.

Figure 2 .
Figure 2. General methodological framework of permafrost favourability and MAAT models.

Figure 4 .
Figure 4. Altitudinal and latitudinal distribution of modelled AAT (lines) and AAT records per weather station.

Figure 5 .
Figure 5. Permafrost favourability index (PFI) map of the semiarid Chilean Andes based on the permafrost distribution model for debris areas.

Table 1 .
Location of the weather stations and number of years between 1981 and 2010 with observations.

Table 2 .
Model coefficients and goodness-of-fit statistics for the linear mixed-effects model for annual air temperature.
P value of the Wald test * < 0.001.

Table 3 .
Distribution of areas favourable for permafrost in the semi-arid Chilean Andes by watershed.