Future permafrost conditions along environmental gradients in Zackenberg, Greenland

The future development of ground temperatures in permafrost areas is determined by a number of factors varying on di ﬀ erent spatial and temporal scales. For sound projections of impacts of permafrost thaw, scaling procedures are of paramount importance. We present numerical simulations of present and future ground temperatures at 10 m 5 resolution for a 4 km long transect across the lower Zackenberg valley in NE Greenland. The results are based on stepwise downscaling of General Circulation Model-derived future projections using observational data, snow redistribution modeling, remote sensing data and a ground thermal model. Comparison to in-situ measurements of thaw depths at two CALM sites and 10 m ground temperatures in two boreholes suggest 10 agreement within 0.10 m for the maximum thaw depth and 1 ◦ C for annual average ground temperature. Until 2100, modeled ground temperatures at 10 m depth warm by about 5 ◦ C and the active layer thickness increases by about 30 %, in conjunction with a warming of average near-surface summer soil temperatures by 2 ◦ C. While permafrost remains thermally stable until 2100 in most model grid cells, the thaw threshold is ex-15 ceeded for a few model years and grid cells at the end of this century. The ensemble of all 10 m model grid cells highlights the signiﬁcant spatial variability of the ground thermal regime which is not accessible in traditional coarse-scale modeling approaches.


Introduction
The stability and degradation of permafrost areas are extensively discussed in the fate 20 of future climate changes as potentially important source of greenhouse gases (Schuur et al., 2008(Schuur et al., , 2009Elberling et al., 2010Elberling et al., , 2013, as well as in the light of infrastructure stability (Wang et al., 2003(Wang et al., , 2006 and farming potential (Mick and Johnson, 1954;Merzlaya et al., 2008). Depending on the emission scenario, future projections based on coarse-scale General Circulation Models (GCM) suggest a loss of 30 % to 70 % of 25 the current near-surface permafrost extent until 2100, in conjunction with a significant 3908 jections of future greenhouse gas emissions from permafrost areas.
Regional Climate Models (RCMs) facilitate downscaling of GCM output to scales of several kilometers so that e.g. regional precipitation patterns and topography-induced temperature gradients are much better reproduced. Based on RCM output, projections of the future ground thermal regime have been performed for a number of permafrost 10 regions, e.g. NE Siberia (50 km resolution, Stendel et al., 2007), Greenland (25 km resolution, Daanen et al., 2011), and Alaska (2 km resolution, Jafarov et al., 2012). While this constitutes a major improvement, many processes governing the ground thermal regime vary strongly at even smaller spatial scales so that the connection between model results and ground observations is questionable. In high-Arctic and mountain 15 permafrost areas exposed to strong winds, redistribution of blowing snow can create a pattern of strongly different snow depths on distances of a few meters. Since snow is an effective insulator between ground and atmosphere (Goodrich, 1982), a distribution of ground temperatures with a range between average maximum and minimum temperatures of 5 • C and more is created (e.g. Gisnås et al., 2014), which is of sim-forcing crosses the thawing threshold, while near-surface permafrost is conserved for more than two decades in areas with high organic and ground ice contents and/or a dry insulating surface layer. In addition, the soil carbon content in Arctic landscapes is unevenly distributed (Hugelius et al., 2013), and GHG emissions from localized carbonrich hotspots can contribute a significant part to the landscape signal (e.g. Walter et al., 5 2006;Mastepanov et al., 2008). Therefore, both the carbon stocks and the physical processes governing permafrost evolution must be understood at the appropriate spatial scales to facilitate improved predictions on the permafrost-carbon feedback.
In recent years, modeling schemes capable of computing the ground thermal regime at significantly higher spatial resolutions of 10 to 30 m have been developed and ap-10 plied in complex permafrost landscapes (e.g. Zhang, 2013;Zhang et al., 2013;Gruber, 2012, 2014;Fiddes et al., 2013). These approaches can capture smallscale differences in altitude, aspect and exposition, as well as in surface and subsurface properties, but the redistribution of snow through wind drift is only included in a simplified way through precipitation correction factors (Fiddes et al., 2013;Zhang et al., 15 2012). On the other hand, dedicated snow redistribution models of various levels of complexity exist (e.g. Winstral et al., 2002;Lehning et al., 2006) with which the pattern and evolution of snow depths can be simulated.
In this study, we make use of such an approach, the deterministic snow modeling system MicroMet/SnowModel (Liston and Elder, 2006a, b), to achieve high-resolution 20 simulations of the ground thermal regime at the Zackenberg permafrost observatory in NE Greenland (Meltofte et al., 2008) until 2100. MicroMet/SnowModel is employed as part of a sequential downscaling procedure, including the RCM HIRHAM5 (Christensen et al., 1996) and the ground thermal model CryoGrid 2 (Westermann et al., 2013). With a spatial resolution of 10 m, the effect of snow distribution patterns and different subsur- 25 face and surface properties on ground temperatures can be accounted for. The study aims to fill the gap between the coarse-and the point-scale modeling studies on the future ground thermal regime which are available for the Zackenberg valley so far. The 25 km-scale, Greenland-wide assessment of Daanen et al. (2011) puts Zackenberg in 3910 the zone of "high thaw potential" until the end century, with modeled ground temperatures of −5 to −2.5 • C and an active layer thickness of 0.5 to 0.75 m for the period 2065-2075. On the other hand, the detailed point-scale study by Hollesen et al. (2011) suggests a future active layer thickness of 0.8 to 1.05 m for a site with average soil moisture conditions which is not representative for many other sites, e.g. the wetlands, 5 found in the Zackenberg valley. Extending this earlier work, we present simulations for a 4 km-transect cutting across typical vegetation zones in the lower parts of Zackenberg valley which allow estimating the range of ground thermal conditions that could be encountered until the end of the century.
2 The Zackenberg site 10 Zackenberg valley (74 • 30 N, 20 • 30 W) is a wide lowland valley dominated by Quaternary non-calcareous sediments with significant periglacial activity and continuous permafrost (Elberling et al., 2004(Elberling et al., , 2008, with a mean annual air temperature of −9.5 • C (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) according to Elberling et al. (2010). Maximum active layer thickness varies from 40 cm to more than 2 m and has increased significantly by 0.8 cm to 1.5 cm 15 per year between 1996 and 2012 (Elberling et al., 2013), which has been determined at two sites (denoted ZeroCalm 1 and 2) of the "Circumpolar Active Layer Monitoring" (CALM) program (Brown et al., 2000). This increasing thickness of the active layer represents only a fraction of the permafrost degradation taking the high ice contents (40-80 %) into account. 20 From the hilltops towards the depressions, an increase in soil water content is seen from dry to wet conditions at the foot of the slopes due to snow melt water being released during large parts of the summer. Roughly one third of the lowland area in Zackenberg is poorly drained. Given the low summer precipitation, water availability during the growing season is mainly controlled by the location of large snow patches 25 melting during the growing season, resulting in the distinct vegetation zonation around these. The topography, landscape forms and wind direction are main factors controlling both water drainage and snow distribution. These patterns are found on both a landscape scale and a small scale (100-200 m) and can therefore be illustrated conceptually as a transect across typical landscape forms in the valley from hilltops to depressions. The top of the hills are windblown and exposed throughout the year with little or no 5 accumulation of snow. From the hilltops towards the depressions there is an increase in soil water content from dry conditions (even arid conditions and salt accumulation at the soil surface) at the hilltops to wet conditions in the bottom of the depression. The dominant wind pattern during winter leaves large snow-patches on the south facing slopes ensuring high surface and soil water contents during a large part of the growing 10 season. Bay (1998) described and classified the plant communities in the central part of the Zackenberg valley and mapped their distribution. The vegetation zones range from fens in the depressions to fell-fields and boulder areas towards the hill-tops. East of the river Zackenbergelven the lowland is dominated by Cassiope tetragona heaths 15 mixed with Salix arctica snow-beds, grasslands and fens; the latter occurring in the wet, low-lying depressions, often surrounded by grassland. On the transition from the lowland to the slopes of Aucellabjerg (50-100 m a.s.l.), the vegetation is dominated by grassland. Between 150 and 300 m a.s.l., open heaths of mountain avens, Dryas sp., are dominating, and gradually the vegetation becomes more open with increasing 20 altitude towards the fell-fields with a sparse plant cover of Salix arctica and Dryas sp. Grassland, rich in vascular plant species and mosses, occurs along the wet stripes from the snow-patches in the highland (250-600 m a.s.l.).
For monitoring purposes, a transect cutting across the main ecological zones of the Zackenberg valley from sea level to 1040 m a.s.l. at the summit of Aucellabjerg 25 has been established. Along this so-called ZERO ("Zackenberg Ecological Research Operations") line (Fredskild and Mogensen, 1997;Meltofte et al., 2008), changes in species composition and distribution of plant communities are investigated regularly. The most recent analysis in 2005 showed no distinct signs of floristic changes since 3912 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2000, although the major snow-patches in the area have been diminishing during the last five years.

Modeling tools
In order to sample the spatial variability of ground temperatures in the Zackenberg valley, simulations from 1960 to 2100 are performed for grid cells of 10 m resolution for the 5 lower 4 km of ZERO-line from the coast to an elevation of 200 m a.s.l. (in total 437 grid cells). In addition, the 100 by 100 m large CALM sites ZeroCalm 1 and 2 are simulated (in total 200 grid cells). To compile forcing data sets at such high resolution, a multistep downscaling procedure is employed, which is schematically depicted in Fig. 1. It is designed to account for the spatial variability of snow depths, differences in summer 10 surface temperature (due to e.g. different evapotranspiration rates caused by surface soil moisture and land cover), as well as spatially variable ground thermal properties and water/ice contents. Differences in insolation due to exposition and aspect are not accounted for which is acceptable for the gentle topography (average slope 2.8 • ) in the modeled part of ZERO-line. In the following, the different parts of the scheme and their 15 interplay are described.

The permafrost model CryoGrid 2
CryoGrid 2 is a physically-based thermal subsurface model driven by time series of near-surface air temperature and snow depth which has been recently employed to assess the evolution of permafrost extent and temperatures in Southern Norway (West-20 ermann et al., 2013). The physical basis and operational details of CryoGrid 2 are documented in Westermann et al. (2013), so that only a brief overview over the model properties is given here. CryoGrid 2 numerically solves Fourier's Law of conductive heat transfer in the ground to determine the evolution of ground temperature with the thermal conductivity k [W m −1 K −1 ] being a function of the volumetric fractions and thermal conductivities of the constituents water, ice, air, mineral and organic (Westermann et al., 2013) following the formulation of Cosenza et al. (2003). For the thermal 5 conductivity of the mineral fraction of the soil, we assume 3.0 W m −1 K −1 which is a typical value for sedimentary and metamorphic rock with low quartz content (Clauser and Huenges, 1995), as dominant in most parts of the Zackenberg valley (Koch and Haller, 1965). For the organic soil fraction, the standard value of 0.25 W m −1 K −1 (e.g. Côté and Konrad, 2005) for peat is employed. 10 The latent heat from freezing soil water or melting ice is accounted for in terms of an effective heat capacity c eff [J m −3 K −1 ], which increases strongly in the temperature range in which latent heat effects occur. Its shape is determined by the soil freezing characteristic, i.e. the function linking the soil water content to temperature, which is related to the hydraulic properties of the soil in CryoGrid 2 (Dall'Amico et al., 2011) for 15 three soil classes, sand, silt and clay. To account for the build-up and disappearance of the snow cover, the position of the upper boundary is allowed to change dynamically by adding or removing grid cells. Movement of soil water is not accounted for, so that the sum of the soil water and ice contents are constant in CryoGrid 2. For spatially distributed modeling, the target domain is decomposed in independent grid cells each 20 featuring a set of model parameters.

Model initialization:
The initial temperature profile for each grid cell is obtained by a multi-step initialization procedure, which allows to approximate steady-state conditions in equilibrium with the climate forcing for the first model decade (September 1958-August 1968) in a computationally efficient way. The method which is described in more 25 detail in (Westermann et al., 2013) accounts for the insulating effect of the seasonal snow cover as well as the thermal offset (Osterkamp and Romanovsky, 1999

Driving data sets:
To account for small-scale differences in surface soil moisture, which give rise to spatially different surface temperatures, we employ the empirical concept of n factors, which relate average air temperature T air to surface temperature T s by T s = n t T air : This rough treatment of summer surface temperatures (which has been applied in previous modeling studies, e.g. Hipp et al., 2012) is focused on seasonal averages and can not reproduce surface temperatures on shorter timescales, e.g. the daily cycle. As a result, a comparison of temperatures in upper soil layers is less meaningful than for deeper layers, which are only influenced by seasonal or even multi-annual average temperatures. However, the n factor based approach precludes the need to compute the surface energy balance, and allows employing measured historic time series of air temperatures (such as the one from Daneborg, Sect. 3.4) for ground thermal modeling. The summer-time n factor n t is computed according to the "Normalized Difference

15
Vegetation Index" NDVI of each grid cell using The relationship is compiled with n t as the ratio of degree-day sums at the soil surface to that in the air over the summer season at both Zackenberg (74.5 • N) and Kobbefjord (65.6 • N) close to Nuuk in W Greenland. As n t is known to be strongly affected by 20 the cover fraction and the amount of vegetation the widely used NDVI at the maximum of the growing season was used as the controlling factor (Fig. 2). n t values near unity indicate little difference between the degree-day sums at the soil surface and the air and n t = 1 is given by NDVI = 0.2, which equalize a biomass production of 315 kg ha −1 (Jacobsen and Hansen, 1999), a cover fraction of 5 % (Jacobsen and Hansen, 1999)  the soil-surface temperatures are warmer that the air temperatures and this mostly occur on nearly barren mineral soils. The minimum n t values of approx. 0.65 are found in moist fen areas indicating a strong cooling effect during the summer on the mineral soils of these sites. Fig. 2 also shows a strong correlation between n t values (Klene et al., 2001) and NDVI values (Walker et al., 2003) from the Kuparuk River Basin, 5 Alaska, USA.
For each 10 m model grid cell, an NDVI value was accumulated from a 2.5 m multispectral Quickbird 2 image of the Zackenberg area acquired around noon on 7 July 2011. While the acquisition date is close to the annual maximum of NDVI values, it represents a single point in time. Therefore, it must be put in perspective to the 10 strong seasonal and interannual variability in plant growth and consequent evolution of NDVI values (Tamstorf et al., 2007). While this constitutes hard to quantify error source, the general agreement in the coverage of the different vegetation classes (see Sect. 3.2) with field observations suggests that the satellite image is an adequate basis to capture the pattern of surface soil moisture and summer surface temperatures along 15

ZERO-line.
Ground properties: Based on a NDVI-classification, six ecosystems were identified in Zackenberg valley (Bay, 1998;Tamstorf et al., 2007;Ellebjerg et al., 2008). Areas with NVDI < 0.2 are dominated by fell-field with a sparse vegetation. In the high mountains such areas are found on solifluction soils, patterned ground and rocky ravines. Dryas 20 heath dominates areas with NDVI between 0.2 and 0.3. Fell-field and Dryas heath are both situated at exposed plateaus, where snow often blows off during the winter months and hence have thinner snow cover. Here, plant species experience an early snowmelt and hence, an early start of the growing season. Cassiope heath (NDVI between 0.3 and 0.4) depends on a protective snow cover during in winter occurs mainly in the 25 lowland on gentle slopes facing south and leeward from the northerly winds which dominate the winter period . Salix snow-bed features NDVI values are between 0.4 and 0.5. This ecosystem, which is unique for E Greenland, occurs mostly on sloping terrain often below the Cassiope heath belt on the slopes, where the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | snow cover is long-lasting, so that the soil moisture in the Salix snow-bed areas are higher. In the wetland areas with NDVI higher than 0.5, grassland and fen areas are distinguished. Grassland occurs mostly on slightly sloping terrain with adequate supply of water early in the season, while the soil water regime can change from wet to moist later in the season. The fen areas occur on flat terrain in the lowland, where the soil is 5 permanently water-saturated throughout the growing season. In August 2013, a classification of ecosystem classes according to the dominating plant species and qualitative surface moisture conditions was conducted along the modeled part of ZERO-line at spatial resolution of 10 m which resulted in 5 % fell, 20 % Dryas, 35 % Cassiope, 15 % Salix snow-bed and 25 % wetland (fen and grassland areas were not distinguished).

10
Using satellite-derived NDVI values (see previous section), these ecosystem fractions could be well reproduced for fell (9 %), Dryas (22 %) and Cassiope (39 %), while a strong discrepancy was encountered for the Salix and wetland classes. Therefore, Salix snow-bed was merged with wetland, yielding a wetland fraction of 30 %. The "true" Salix class is hereby split between Cassiope and wetland which is reflected in 15 the strong concentration of grid cells with NDVI values around 0.4. This suggests a significant overlap of the NDVI values from the different classes in this region for the particular satellite acquisition date, so that the classes can not be separated by their NDVI value. While the NDVI-derived ecosystem classification constitutes a potentially important source of uncertainty in the modeling chain, it provides the possibility to use 20 satellite images and thus apply the classification procedure for larger regions, e.g. the entire Zackenberg valley, at high spatial resolutions, which can hardly be achieved by manual mapping.
For the remaining four classes fell, Dryas, Cassiope and wetland, typical soil stratigraphies were assigned based on and guided by in-situ measurements in soil samples 0.2 and 0.3 for the active layer at a Cassiope site (Hollesen et al., 2011). For the Dryas and fell classes, large changes in soil moisture were encountered after rain falls which made the values strongly dependent on the timing of the sampling. With 5 % and less, the volumetric contents of organic material are low in all classes and have negligible influence on the thermal properties of the soil. Following measurements of soil cores to 10 2 m depth (Elberling et al., 2010), saturated conditions are assumed below the current active layer for all classes (Table 1), except for fell for which no in-situ data are available. Furthermore, bedrock is assumed below 10 m which is a pure estimate but has limited influence on the outcome of the simulations. Snow properties: In CryoGrid 2, constant thermal properties in space and time are 15 assumed for the snow cover (see Westermann et al., 2013, for details). Following in-situ measurements, a snow density of 300 kgm −3 is employed, which results in a volumetric heat capacity of c snow = 0.65 MJ m −3 K −1 . In the absence of in-situ measurements of the thermal conductivity of the snow cover, we use the empirical relationship between density and thermal conductivity from Yen (1981), which is also employed in 20 the detailed snowpack scheme CROCUS (Vionnet et al., 2012). The resulting value is k snow = 0.25 W m −1 K −1 , slightly lower than those employed in CryoGrid 2 simulations for the mountain environments of Southern Norway where average winter temperatures are higher than in Zackenberg, but predominantly wind-packed snow is encountered as well.

Future climate scenario with HIRHAM
There are several types of uncertainties related to climate projections. Apart from "external" uncertainties such as the future evolution of greenhouse gas emissions, there 3918 Introduction  (Hazeleger et al., 2010(Hazeleger et al., , 2012. The IFS in the current EC-EARTH model is based 10 on ECMWF cycle 31r1 with some improvements from later cycles implemented, including a new convection scheme and a new land surface scheme (H-TESSEL) as well as a new snow scheme (Hazeleger et al., 2012). The atmospheric part of EC-EARTH is configured with a horizontal spectral truncation of T159, which is approximately 125 km × 125 km in latitude and longitude. The vertical resolution is 62 layers.

15
The ocean and sea ice components have 42 vertical layers and a roughly 1 degree horizontal resolution with refinement to 1/3 degree around the equator. EC-Earth is one of the models of CMIP5 (Coupled Model Intercomparison Project) and has been used for the experiments for the IPCC AR5 report.
To resolve the topography of Greenland adequately, a horizontal resolution of 5 km 20 or finer is required . The output of EC-Earth is therefore downscaled to the RCM grid. The RCM used here is HIRHAM5 in its newest version, which includes calculation of the surface mass balance of the Greenland Ice Sheet. A surface snow scheme has been implemented over glaciers. The model setup is described in Rae et al. (2012) except that the resolution here is 0.05 degrees ( namely 1991-2010, 2031-2050 and 2081-2100. The scenario used was RCP 4.5 (Thomson et al., 2011;Clarke et al., 2007;Smith and Wigley, 2006;Wise et al., 2009), which gives an additional radiative forcing in 2100 with respect to preindustrial values of 4.5 W m −2 . In this rather conservative scenario, CO 2 emissions peak around 2040 and decline thereafter, 5 resulting in a CO 2 concentration of 550 ppm in 2100, which is just below a doubling with respect to preindustrial values.

Modeling snow distribution by MicroMet/SnowModel
SnowModel is a spatially distributed snow-evolution modeling system (Liston and Elder, 2006a) which was applied in the Zackenberg study area (14 km by 12 km) to describe 10 the snow distribution through a seven year period covering August 2003 to September 2010. SnowModel consists of three interconnected submodels: Enbal, SnowPack, and SnowTran-3D. Enbal calculates surface energy exchanges and snowmelt (Liston, 1995(Liston, , 1999, SnowPack models the evolution of the snow depth and snow-water equivalent in time and space (Liston and Hall, 1995;Liston and Mernild, 2012), and the trans-15 port of blowing snow is generated by SnowTran-3D (Liston and Sturm, 1998;Liston et al., 2007). SnowModel was coupled with a high-resolution atmospheric model, Mi-croMet (Liston and Elder, 2006b), which spatially distributed the micrometeorological input parameters over the simulation domain. MicroMet requires meteorological station and/or atmospheric (re)analysis inputs of air temperature, relative humidity, pre-20 cipitation, wind speed, and wind direction. Furthermore, available observed incoming shortwave and longwave radiation were included. All meteorological parameters except precipitation were measured by five automatic weather stations distributed in the valley and on mountains contained within the simulation domain ( Table 2). Because of missing data and uncertainties associated with in situ winter precipita-25 tion measurements, MicroMet precipitation inputs were provided by the North American Regional Reanalysis (NARR) (Mesinger et al., 2006). These NARR precipitation fluxes were adjusted using the SnowAssim (Liston and Hiemstra, 2008) data assimila-3920 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tion scheme under the constraint that modeled snow-water-equivalent depth matched observed pre-melt snow depth and snow density at locations where those observations were made. Additionally, a digital elevation model (DEM) and a land-cover map were required for the MicroMet/SnowModel simulations. These distributions were provided over the simulation domain at a 10 m by 10 m spatial resolution. The DEM was based on an August 2000 aerial survey, and the land-cover map was based on the Elberling et al. (2008) vegetation classification. In the land-cover map, vegetated areas were assigned a canopy height which defined the snow-holding depth (shd), i.e. the depth to which the vegetation is able to hold the snow and prevent snow transport by wind (snow over this depth is available for wind redistribution). Non-vegetated areas were assigned a shd of 0.01 m. The modeled snow depths and the timing of snowmelt are validated against automated and manual snow depth observations conducted in two areas of 100 m × 100 m along ZERO-line (ZeroCalm 1 and 2). Furthermore, the modeled snow depth curves are validated by comparing the melt-out times for each grid cell and year to orthorec-15 tified images taken by an automatic camera system at 400 m a.s.l. at the slope of Zackenberg overlooking ZERO-line (Hinkler et al., 2002). For the study area, the modeled melt-out pattern matched the observations generally well, which suggests realistic maximum snow depths at least in case of realistic melt rates. 20 To run simulations of permafrost temperatures from 1958 to 2100, a continuous record of the driving data air temperature and snow depth was compiled from various sources. The method assumes that trends in air temperature and precipitation measured at one point or modeled by a medium-scale atmospheric scheme are representative for the trends along ZERO-Line.  for which an hourly record is available for the periods 1958-1975 and 1979-2011. 5 The gap was filled using linear interpolation between the averages from five years before and five years after for each hourly value.

Downscaling scheme from GCM to plot scale
-For present-day and future air temperatures, the near-surface air temperature from the HIRHAM5 5 km grid cell closest to the study area, which are available for three time slices, 1991-2010, 2031-2050 and 2081-2100. The gaps in between

Permafrost simulations along the ice edge to sea transect
To put the small-scale simulations along ZERO-line in a regional context, we perform first-order estimates for the permafrost development for sites at the ice edge (

Comparison to field data
To build confidence that the modeling is a satisfactory representation for the true ground thermal conditions, the model results are compared to available in-situ data sets. These comprise in particular thaw depth measurements at ZeroCalm 1 and 2 since 1996, as 5 well as temperatures from deep boreholes since 2012. Active layer thickness: The modeled and measured maximum thaw depths for seven years for which MicroMet/SnowModel was run and a realistic snow distribution can be assumed are shown in Fig. 3, with the areas selected for comparison equal to Elberling et al. (2013). Most importantly, CryoGrid 2 can capture the significant differences be-10 tween the three sediment classes Dryas, Cassiope and wetland caused by the different soil moisture contents. With a few exceptions, CryoGrid 2 can reproduce the measured thaw depth within the spatial variability in the validation areas (indicated by the error bars in Fig. 3), with the exception of the year 2006 which features stronger deviations from the measurements. The spatial variability within the target areas is significantly 15 smaller in the model runs than in nature, most likely since the sediment classification assumes constant soil properties within each class, while the soil composition can vary significantly within a class in reality.
Permafrost temperatures: In deeper layers, ground temperatures are influenced by the temperature forcing of an extended period prior to the measurement. Therefore, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | satisfactory agreement suggests that the statistical downscaling procedure (Sect. 3.4) employed to produce the forcing data for the model spin-up is adequate for the area.

Evolution of active layer thickness and ground temperatures
The modeled evolution of the temperature distribution at 1 m depth along ZERO-line is shown in Fig. 4. The modeled temperatures extend over a range of 2 to 5 • C from min-5 imum to maximum which is evidence of the significant spatial variability of the ground thermal regime in this landscape. According to the climate change scenario of the future projections (Sect. 3.2), ground temperatures will warm by about 4 • C until 2100, but permafrost will largely remain thermally sustainable along ZERO-line. However, the high-resolution simulations suggest a few localized "weak spots" where the yearly av-10 erage ground temperatures become positive in some years at the end of this century (Fig. 4). The future warming of air temperatures predicted by HIRHAM is not constant over the year, with the most pronounced warming of 0.4-0.6 • C per decade occurring in fall, winter and spring, while summer (June to August) temperatures increase by less 15 than 0.2 • C per decade. As a consequence, the annual maximum thaw depths increase only moderately until 2100, from 0.8-1.0 m to 1.1-1.4 m for Dryas, from 0.65-0.85 m to 0.8-1.1 m for Cassiope, and from 0.5-0.65 m to 0.6-0.8 m for the wetland class (Fig. 5).
The climate sensitivity of thaw depths is different between the classes, with a stronger increase for the classes with dryer soils than for the wetlands. It is remarkable that the 20 projected increase is only 0.1-0.2 m in the wetlands which can be related to the high ground ice content in the top permafrost layer. The biological activity in this high-Arctic ecosystem are strongly related to summer conditions. The simulations suggest a significant increase in average summer temperatures and thawing degree days (Fig. 6) within the effective root depth. The combination 25 of deeper active layer (Fig. 5) and warmer near surface (Fig. 6)  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | soil conditions allow plant to benefit from the additional nutrient and like to result in changes in the competition among plant species for light. This may lead to marked changes in vegetation over time, which however is beyond the scope of this study.

Ice edge to sea transect
The Zackenberg valley is located about halfway between the open Atlantic Ocean and 5 the edge of the Greenland ice sheet, which give rise to a strong climatic gradient from west to east over a distance of less than 100 km. To put the modeling results in this larger-scale perspective, permafrost temperatures at two additional sites along this west-east transect are estimated by using HIRHAM modeling results as input for Cryo-Grid 2 (Sect. 3.5). Figure 7 displays modeled ground temperatures for the ice edge, 10 Revet (located roughly half-way between Zackenberg and the ice edge), Zackenberg and Daneborg located close to the open ocean. Although driven by raw (i.e. not downscaled) coarse-scale modeling data, the Zackenberg temperatures are mostly within the temperature range provided by the high-resolution MicroMet/SnowModel-based simulations. While the temperatures at Revet and Daneborg are similar to Zackenberg, 15 the model suggests significantly colder ground temperatures close to the ice edge. Furthermore, the warming from present to 2100 is only about half at the ice edge, which suggests a different susceptibility to climate change at this site. The simulations suggest that permafrost further inland continue to be thermally sustainable until well into the next century, while coming close to the thaw threshold at least in localized spots 20 closer to the ocean. It must be emphasized that theses simulations of ground temperatures along this ice edge to sea transect must be considered first-order estimates due to simplified treatment of the snow cover (Sect. 3.5) and the lack of adequate representation of the complex topography in the 5 km forcing data. Furthermore, a thorough validation can 25 not be performed in the absence of in-situ data of ground temperatures at the ice edge, Revet and Daneborg. 3926 5 Discussion

Scaling strategies from GCM to plot scale
The presented simulations of ground temperatures and permafrost state variables are derived from a multi-step downscaling approach (Sect. 3.4), which bridges from the coarse spatial resolution of a GCM (hundreds of km) to the impact scale on the ground 5 (set to 10 m for this study). As such, the scheme is technically capable of bridging five to six orders of magnitude in space. The two main driving environmental variables for the thermal model CryoGrid 2 are surface temperature and snow depth.
The snow depth is assumed to be controlled by wind drift of snow at the target scale of 10 m which is modeled by the snow redistribution scheme MicroMet/SnowModel. MicroMet/SnowModel is a deterministic scheme, which is capable of predicting the snow depth for each model grid cell and thus e.g. reproduce the location of snow drifts and bare-blown spots. Such deterministic high-resolution modeling facilitates a better comparison and validation with ground observations, but is restricted to small model domains for computational reasons. However, SnowModel also includes the ability to 15 simulate snow distributions over large areas (e.g., the ice-free parts of Greenland, several 100 000 km 2 ) using, for example, subgrid snow distribution representations (e.g. Liston et al., 1999;Liston, 2004;Liston and Hiemstra, 2011). Gisnås et al. (2014) recently presented a statistical approach to account for the impact of the small-scale variability of snow depths on ground temperatures which is applicable on large spatial 20 domains.
The surface temperature is derived from air temperature for which the regional gradients are based on the RCM at a scale of 5 km. Within the target area along ZERO-line (a distance of 4 km), variations in air temperature are generally small. Further downscaling to 10 m is accomplished by using a high-resolution NDVI satellite image and 25 the NDVI vs. n factor relationship (Sect. 3.1) which is used to convert air to surface temperatures during the snow-free season. By this scheme, a high-resolution data set of surface temperatures is generated from comparatively low-resolution air tempera-3927 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ture data. More physically-based approaches make use of the surface energy balance (SEB) to compute surface temperatures from air temperature, wind speed, humidity and incoming radiation (e.g. Zhang et al., 2013;Fiddes et al., 2013). In addition to accounting for more complex topography through spatially distributed fields of incoming radiation, the surface energy balance and thus the surface temperature can directly be 5 connected to surface soil moisture and land cover/vegetation type which circumvents the use of n factors. On the other hand, SEB models require additional driving data sets, in particular incoming radiation, which must be compiled e.g. from large-scale atmospheric modeling (Fiddes and Gruber, 2014) and/or from sparse in-situ measurements (Zhang et al., 2012). Due to the potential for serious biases in such driving data 10 sets in remote locations (such as Zackenberg), it must remain unclear whether the capacity of SEB models in reproducing the true surface temperature is superior to the simple empirical concept employed in this study.

Model uncertainty
The presented model results must be considered a first-order approximation on the 15 future thermal state of the permafrost, which is subject to considerable uncertainty due to a variety of factors. Firstly, only one climate change scenario is considered, which does not account for the considerable spread in climate predictions. With permafrost approaching the thaw threshold at the end of this century for RCP 4.5 forcing, wide-spread permafrost degradation is e.g. likely for more aggressive climate change 20 scenarios. Secondly, the downscaling procedure from large-scale model data to high-resolution fields of temperature and snow depth is susceptible to uncertainties, since it assumes constant offsets between the two data sets based on the climate conditions of a sevenyear reference period, which may not be justified for a one-hundred-year period. This is 25 in particular critical since the temperature regime in the study area is characterized by strong inversion settings during a large part of the year (Meltofte et al., 2008). A modification of such inversions could lead to a different climate susceptibility of the study 3928 area compared to the large-scale trend, which cannot be captured during the reference period. Furthermore, the future snow distribution patterns are based on random years from the seven-year reference period, implying that the patterns are unchanged in a warmer future climate and that the reference period allows a representative sample of the interannual variability within the rest of the century. It is not unlikely that both 5 assumptions are violated at least to a certain degree. Furthermore, the CryoGrid 2 permafrost model assumes properties and relationships compiled and adjusted for the present state to be valid in the future. This concerns in particular the NDVI-based summer n factor relationship employed to derive surface from air temperatures (Sect. 3.1), as well the thermal properties of the snow and the ground stratigraphies. As an example, the snow density and thermal conductivity are likely to increase in a warmer climate, which would decrease the insulation of the winter snow cover and thus lead to colder temperatures as suggested by the model. Furthermore, constant soil water and ice contents are assumed, thus neglecting both seasonal and long-term changes in the water cycle. However, at least for the Cassiope class, our results for the future increase 15 in maximum thaw depth are in good agreement with the study of Hollesen et al. (2011) who used the coupled heat and water transfer model COUP (Jansson and Karlberg, 2004) to simulate the ground thermal and moisture regimes in this century. Finally, new processes not accounted for by the modeling scheme might become relevant in the course of climate warming, e.g. the occurrence of wintertime rain events, which can 20 lead to significant additional ground warming (Westermann et al., 2011).

From model results to permafrost landscape development
Most real-world applications of permafrost models assume non-interacting grid cells with spatially constant soil properties. Consequently, permafrost degradation in model studies (e.g. Westermann et al., 2013) is generally described as talik formation manifested in the temperature profile of a 1-D-grid cell. While such is indeed observed in instrumented boreholes, it can be accompanied by significant changes in the hydrological regime by thawing of hydrological barriers or the formation of new aquifers. operational permafrost models are not capable of predicting such developments, which is a significant limitation for sound predictions on the permafrost carbon feedback. For the study area, Elberling et al. (2013) demonstrated that the potential CO 2 emissions from carbon-rich wetland soils strongly depend on the future hydrological regime of the wetland, with a drying of the wetland leading to significantly faster carbon turnover. Fur-5 thermore, thawing of excess ground ice can entirely modify the landscape, e.g. through thermokarst or thaw slumps, which can be hotspots of greenhouse gas emissions and thus modify the carbon budget of an entire permafrost landscape.
From the perspective of model development, a simple increase of the spatial resolution seems a prerequisite to resolve such shortcomings in the next generation of per-10 mafrost models. At a 10 m resolution, this study could capture two important aspects which can be seen as part of the "thermal signature" of the permafrost landscape in Zackenberg: (a) the differences in maximum thaw depth between different ecosystem classes, and (b) the spatial variability of ground temperatures to a large extent caused by spatially variable snow depths. Compared to large-scale (as in Daanen et al., 2011) 15 or point-scale simulations (as in Hollesen et al., 2011), it provides a far more detailed (though still incomplete) assessment of the possible development of the Zackenberg permafrost landscape which can be better linked to studies on the future ecosystem carbon turnover (e.g. Elberling et al., 2013).

20
This study presents numerical simulations of present and future ground thermal conditions for a transect across the low-lying parts of the high-Arctic Zackenberg valley in NE Greenland. At the modeling scale of 10 m, the governing factors for the ground thermal regime are accounted for in a deterministic way. This involves snow depth (derived from a blowing snow model), ground properties (derived from an NDVI-based 25 ecosystem classification) and surface properties (derived from empirical correction factors for summer surface temperature based on NDVI). Past and future climate trends 3930 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | for the general area are derived from in-situ records and a Regional Climate Model. The following conclusions can be drawn from this study: -The model approach can capture the measured differences in maximum thaw depth between different ecosystem classes encountered in the area. The simulated ground temperatures are in agreement with available borehole measure-5 ments.
-The modeled average ground temperatures increase by about 4 K until 2100, but generally remain below 0 • C. However, a few model grid cells feature positive average 1 m temperatures in single years after 2060, which suggests that the thaw threshold could be reached at localized "weak" spots.

10
-The modeled maximum thaw depths increase moderately in all ecosystem classes, with the lowest value of 0.2 m for the wetland sites.
-The spatial variability of the average ground temperatures within a distance of a few kilometers is between 3 and 5 K and thus on the order of the projected increase of ground temperatures until the end of this century. 15 The study exemplifies that grid-based simulations at coarse scales with only one or few model realizations, as they are typical for the land-surface schemes of atmospheric models, cannot account for the spatial variability of the ground thermal regime in many permafrost areas. They are hence not capable of correctly predicting the onset of permafrost thaw which can trigger transformative landscape changes due to erosion and 20 hydrological processes.