Modelling and mapping climate change impacts on permafrost at high spatial resolution for an Arctic region with complex terrain

Most spatial modelling of climate change impacts on permafrost has been conducted at half-degree latitude/longitude or coarser spatial resolution. At such coarse resolution, topographic effects on insolation cannot be considered accurately and the results are not suitable for landuse planning and ecological assessment. Here we mapped climate change impacts on permafrost from 1968 to 2100 at 10 m resolution using a process-based model for Ivvavik National Park, an Arctic region with complex terrain in northern Yukon, Canada. Soil and drainage conditions were defined based on ecosystem types, which were mapped using SPOT imagery. Leaf area indices were mapped using Landsat imagery and the ecosystem map. Climate distribution was estimated based on elevation and station observations, and the effects of topography on insolation were calculated based on slope, aspect and viewshed. To reduce computation time, we clustered climate distribution and topographic effects on insolation into discrete types. The modelled active-layer thickness and permafrost distribution were comparable with field observations and other studies. The map portrayed large variations in active-layer thickness, with ecosystem types being the most important controlling variable, followed by climate, including topographic effects on insolation. The results show deepening in active-layer thickness and progressive degradation of permafrost, although permafrost will persist in most of the park during the 21st century. This study also shows that ground conditions and climate scenarios are the major sources of uncertainty for high-resolution permafrost mapping.


Introduction
Climate warming at high latitudes was about twice the global average during the 20th century (ACIA, 2005).Observations have shown increases in near-surface ground temperature and active-layer thickness, and at some places, disappearance of permafrost (e.g., Vitt et al., 2000;Smith et al., 2010).Most climate models project that climate warming in northern high latitudes will continue at a rate higher than the global average during the 21st century (ACIA, 2005), which will further induce permafrost degradation.Permafrost thaw affects infrastructure, ecosystems, wildlife habitats, and has strong feedbacks on the climate system (ACIA, 2005).
Spatial modelling is an effective approach to understanding the distribution of permafrost and its changes with climate.Spatial permafrost models can be broadly classified as equilibrium models and transient models (Riseborough et al., 2008;Shiklomanov et al., 2007).Due to the slow response of ground thermal regimes to changes in atmospheric climate, permafrost conditions are and will continue to be in disequilibrium with the atmospheric climate in the 21st century (Osterkamp, 2005;Zhang et al., 2008).Therefore it is necessary to quantify changes of permafrost based on transient models.Most transient permafrost modelling and mapping studies have been conducted using half-degree latitude/longitude or coarser spatial resolution (e.g., Anisimov and Reneva, 2006;Marchenko et al., 2008;Zhang et al., 2006Zhang et al., , 2008)).These coarse resolution studies cannot accurately consider topographic effects on insolation since topographic conditions usually vary over shorter distances.In addition, coarse resolution results are difficult to validate by comparing with field observations and are not suitable for land-use planning and for ecological assessment.
Modelling and mapping climate change impacts on permafrost at high spatial resolution requires detailed input data, efficient computation schemes, and robust models.Recently, several studies modelled permafrost at higher spatial resolutions.Jafarov et al. ( 2012) mapped ground thermal conditions and changes with climate in Alaska at 2 km resolution using an implicit finite-difference numerical model.The effects of snow were considered explicitly but the topographic effect on solar radiation was not considered.Duchesne et al. (2008) mapped permafrost conditions and changes with climate in the Mackenzie basin at 30 m resolution based on a process-based heat conduction model.The model used seasonal n-factors to estimate the near-surface ground temperature from air temperature.Zhang et al. (2012Zhang et al. ( , 2013) ) mapped permafrost in the northwestern Hudson Bay Lowlands at 30 m resolution using a more detailed model that integrated the effects of other climate variables (e.g.precipitation, solar radiation, and vapour pressure) and changes in snow and soil moisture conditions.This region mainly consists of a peatland plain, where peat layer thickness can be estimated based on elevation.However, Arctic regions are usually not flat and many factors affect the distribution of organic layer thickness.More importantly, complex terrain has significant effects on soil, vegetation, and local climate, and thus has high spatial variations in permafrost distribution.Therefore, higher spatial resolution is needed to map permafrost in complex terrain, especially for land-use planning and for assessing the impacts of permafrost on hydrology, ecosystems and geohazards.The objectives of this study are to develop an approach to model and map permafrost at high spatial resolution for an Arctic region with complex terrain, to test the effects of spatial resolutions, and to identify the input data gaps for future studies.

The study area and field data sources
The study area is Ivvavik National Park (INP) located in the northern tip of Yukon Territory in Canada (Fig. 1).The park covers 10 170 km 2 and is composed of two major zones: the high-relief inland region which belongs to the British Mountains, and the coastal plain lying along the Beaufort Sea.The highest mountain reaches 1655 m.The park is an area of outstanding natural variety, containing arctic tundra, alpine tundra, treeline woodlands and arctic coastal habitats.Most of the park was not glaciated, and so contains ancient landforms that evolved to the present form over millions of years.The park is part of the calving ground and summering area of Porcupine Caribou Herd, numbering about 165 000 animals (Canadian Parks Service, 1993).Currently, permafrost is distributed continuously in the park (Heginbottom et al., 1995).
Three major field data sources were used for this study.One is the resources report for the park by Canadian Parks Service (1993), which provides general information about the climate, geology, geomorphology, soil, hydrology, vegetation and wildlife in the park.The second is a detailed ecosystem classification along the Firth River watershed, which covers about a third of the park (Mackenzie and MacHutchon, 1996).Based on this classification, we extended the ecosystem map to the entire park and defined soil and drainage conditions for each type.The third data source is our field observations conducted over four summers from 2008 to 2011.Field observations include location, topography, drainage, soil condition, plant composition and status, ecosystem type, and summer thaw depth.Summer thaw depth was measured at 162 sites across the park using steel probes and by digging soil profiles, of which 62 sites did not reach the frozen layer due to high stone content in the soil.Figure 1a shows the distribution of the 100 sites where summer thaw depths were recorded.

The model
The Northern Ecosystem Soil Temperature model (NEST) was used to compute and map permafrost conditions in the park.NEST is a one-dimensional transient model that considers the effects of climate, vegetation, snow, and soil conditions on ground thermal dynamics based on energy and mass transfer through the soil-vegetation-atmosphere system (Zhang et al., 2003).Ground temperature is calculated by solving the one-dimensional heat conduction equation.The dynamics of snow depth, snow density and their effects on ground temperature are considered.Soil water dynamics are simulated considering water input (rainfall and snowmelt), output (evaporation and transpiration), and distribution among soil layers.Soil thawing and freezing, and associated changes in fractions of ice and water, are determined based on energy conservation.Detailed descriptions of the model and validations can be found in Zhang et al. (2003Zhang et al. ( , 2005Zhang et al. ( , 2006)).Lateral flows and snow drifting are parameterized in a simplified way (Zhang et al., 2002(Zhang et al., , 2012)).We improved the model to address the topographic effects on insolation so that the NEST model could be used for areas with complex terrain.The algorithms are presented in Appendix A, and the source code is presented as Supplementary Material.

Climate data
There are seven climate stations within or nearby INP.Three stations near the northern shoreline have about 50 yr of observations while the other four stations, including all the three inland stations, have less than 15 yr of observations.Most of the data have some gaps.Several studies developed methods to spatially interpolate monthly climate based on station observations (e.g., Wang et al., 2006;McKenney et al., 2011).Since temporal patterns are usually affected by the changes in climate stations with time used for the interpolation, we estimated the temporal patterns using observations at representative climate stations and estimated the average spatial distributions of air temperature and precipitation based on the monthly spatial data from Wang et al. (2006).Observations show close correlations of air temperature among the stations in the park.For example, the correlation coefficient of the daily air temperature from 1968 to 2010 at Komakuk Beach Station (69.62 • N, 140.20 • W) and Old Crow Station (67.34 • N, 139.50 • W) is 0.92 (N = 15706) although these two stations are about 250 km apart across the park from the coast to inland (Fig. 1).The correlation of daily precipitation between these two stations is poor (R = 0.16, N = 15706).However, the monthly total precipitation is correlated between the stations (R = 0.50, N = 516).Therefore we used station observed climate data to represent both the long-term patterns of the monthly climate and the daily fluctuations within a month for a region surrounding the climate station.
Monthly air temperature in a year for a grid cell was estimated by where T m, g (Y, M) is the monthly mean air temperature for grid cell g in month M year Y .T 0 m, g (M) is the 30 yr  mean air temperature in the month M for the grid cell g, interpolated from climate stations based on Wang et al. (2006) using a spatial resolution of 30 m by 30 m (Fig. 2).T m, s (Y, M) is the deviation of the monthly air temperature in year Y month M from the 30 yr average for that month estimated based on the station observations.The daily air temperature for a grid cell was estimated based on the daily observations at a climate station where T d, g (Y, M, D) and T d, s (Y, M, D) are daily air temperatures (daily maximum or minimum) on day D month M year Y for a grid cell and a climate station, respectively.T m, gs (Y, M) is the difference of monthly air temperatures between the grid cell and the climate station in month M year Y .Precipitation was estimated in a similar way but using ratios instead of differences so that daily and monthly precipitation can be zero.
For the climate conditions in the future (2012-2100), we selected two scenarios generated by two climate models (CC-Cma (CGCM3.1)and MPI ECHAM5).The two scenarios will be referred to as CCCma and ECHAM, respectively.We selected these two scenarios because the two GCMs had a better behaviour in modelling 20th century climate in North America (Walsh et al., 2008), and the two scenarios represented the general range of air temperature change under the intermediate emission scenarios (A1B) for this region.The monthly climate scenario data were downloaded from World Data Center for Climate (http://mud.dkrz.de/wdc-for-climate/, accessed in April 2011).We first calculated the monthly relative changes between the projected climate in the future and the 30 yr averages during 1971-2000 simulated by the same model, and then we used these relative changes to estimate the future monthly climate for all the grid cells in the park since the park is mainly covered by one climate model grid cell.From the 30 yr average  to the 2090s, annual mean air temperature was projected to increase 4.2 • C and 8.8 • C under the scenarios CCCma and ECHAM, respectively, and annual total precipitation was projected to increase 26 % and 37 % under the two scenarios,  2 for the names of the ecosystem types).respectively (Fig. 3).The monthly air temperature and precipitation were converted to daily data using the above described method (Eqs. 1 and 2) based on the historical daily observations at climate stations.
Vapour pressure was estimated based on minimum air temperature (Zhang et al., 2012).Daily total insolation without topographic effects (an input to consider topographic effects on insolation as described in Appendix A) was estimated based on latitude, the day of the year, diurnal temperature range and vapour pressure (Zhang et al., 2012).The parameters were estimated based on observations at Inuvik climate station (68.32 • N, 133.52 • W).
The coastal region of the park neighbours the Arctic Ocean and has a marine climate while the inland region has a continental climate (Canadian Parks Service, 1993).We delineated the southern boundary of the coastal region based on the 300 m elevation contour.We used the daily climate observations at Komakuk Beach Station and Old Crow Station to represent the temporal patterns of the climate in the coastal and the inland regions, respectively; as these two stations have the longest continuous observations within or nearby the park.Observations from these two stations are available from 1958-2011 and 1968-2011, respectively.Data gaps are filled based on observations at nearby stations.

Classifying the average climate and topographic effects on insolation
Since the spatial distributions of the climate were estimated based on long-term monthly means, which do not change with time, we discretized the average climate conditions of the grid cells into different clusters to reduce computation time.Permafrost conditions are mainly dependent on seasonal and annual climate conditions and are not very sensitive to day-to-day climate fluctuations.Active-layer thickness, for example, is mainly determined by the annual total degree days when daily air temperature is above 0 • C (TDD T >0 ) according to the Stefan's equation (Lunardini, 1981).Therefore we used TDD T >0 and annual total precipitation to cluster the average climate.The 30 m resolution TDD T >0 and annual total precipitation were re-sampled to 10 m resolution using the nearest neighbour method.Slope, aspect and viewshed are the topographic attributes affecting insolation.Slope and aspect were calculated using 10 m resolution digital elevation model (DEM), which was re-sampled from 30 m DEM data using bilinear interpolation method.We did this re-sample to match the resolution of the ecosystem map.The DEM data are from the Topographic Data Centre of Natural Resources Canada.The viewshed of a grid cell is the angular distribution of sky visibility versus obstruction, similar to the view taken by upward-looking hemispherical (fisheye) photographs from the centre of the grid cell (Fu et al., 2000).Since most valleys in the park are less than 3 km wide between the peaks, we calculated viewshed for each 30 m by 30 m grid cell using a window of 201 by 201 grid cells of the DEM (or 3 km from the grid cell to the sides of the window).The viewshed angles were calculated for each 22.5 • azimuth direction, with a total of 16 azimuth directions for each grid cell.The viewshed angles were then interpolated to 10 m spatial resolution.Using slope, aspect and the viewshed angles, average annual insolation was calculated for each grid cell for clustering analysis.
We clustered the average climate and topographic features based on TDD T >0 , annual precipitation, average annual insolation, slope and aspect using PCI Isodata unsupervised classification method (PCI Geomatics Enterprises Inc. Canada).Since the software allows a maximum of 255 clusters, we first clustered the inland region into nine large clusters and then further clustered each of them into 254 clusters.The coastal region was directly clustered into 254 clusters and the clustering errors were smaller than the errors for the inland region.Figure 4 shows an example of the climate, insolation and topographic conditions of the clusters in the coastal region.Table 1 lists the average clustering errors.The effects of the clustering errors on the active-layer thickness were generally less than 0.6 % based on Stefan's equation.
To reduce the errors in estimated insolation for each cluster, the topographic attributes of a cluster were assigned using that of a grid cell within the cluster where the difference of the annual total insolation between the grid cell and the cluster average was the smallest.The latitude of that grid cell was used for calculating the insolation of that cluster as well.Similarly, monthly air temperature and precipitation of a cluster were assigned using that of a grid cell in the cluster where the difference of TDD T >0 between the grid cell and the cluster average was the smallest.

Ecosystem types and distribution
Figure 2d shows the ecosystem types in INP mapped by Fraser et al. (2012) using 10 m-resolution SPOT imagery acquired during 13 July to 22 August 2006.The ecosystem types are described in Table 2.The map was developed using a new image-based predictive ecosystem mapping method, which integrates remote sensing-based vegetation map with predictive terrain attributes from a DEM (Fraser et al., 2012).A decision tree classifier was trained based on the ecosystem map of the Firth River basin, developed using air photos and detailed field observations at 367 sites by Mackenzie and MacHutchon (1996).The overall accuracy of the classification is 85 % (Fraser et al., 2012).

Ground conditions and hydrological parameters
The ground conditions and hydrological parameters needed for the model include top organic layer thickness, the texture of mineral soils, organic matter content in the mineral www.the-cryosphere.net/7/1121/2013/The Cryosphere, 7, 1121-1137, 2013  .TDD T >0 and TDD T <0 are annual total degree days for daily air temperature above and below 0 • C, respectively.soils, fraction of gravel (from pebbles to boulders), lateral inflow and outflow parameters, and the snow-drifting parameter.Since no detailed maps are available for these input data, we estimated their spatial distributions based on the ecosystem map.For each ecosystem type, the texture of the mineral soil was determined based on the typical conditions reported by Mackenzie and MacHutchon (1996) and our field observations (Table 3).Peat thickness was estimated based on our field measurements and the land features of the types, especially the drainage and moisture regime described by Mackenzie and MacHutchon (1996).Soil organic matter content in the top mineral soil was estimated according to the typical field measurements.The changes in soil organic content with depth was estimated from the field data and the general patterns described by Hossain et al. (2007).Most areas, especially in deeper ground, contain some gravel.We estimated the fraction of gravel based on the field observations and the biophysical and terrain features of the ecosystem types described by Mackenzie and MacHutchon (1996).Gravel reduces the volume of soil to hold and exchange water.The thickness of the subsoil was estimated assuming that the fraction of gravel increases 2 % for every 10 cm until the gravel fraction reaches 100 %.The ground conditions for the ecosystem types were listed in Table 3. Lateral flow parameters were estimated based on the drainage class and moisture regime of the ecosystem types reported by Mackenzie and MacHutchon (1996) (Table 3).We assumed that the lowest water table depth for lateral outflow (ranged from 0 m to 0.4 m) is proportional to the soil moisture class (ranging from 0 to 8 representing moisture regimes from hydric to very xeric), and the rate of lateral outflow is exponentially related to the drainage class number (ranging from 0 to 6 representing drainage types from very poor to very rapid).We also assumed that there is no surface lateral inflow and water will flow away gradually when the water table is above the land surface (water table decreases 5 % each day when it is above the land surface).
The snow-drifting parameter was estimated considering the effects of land exposure, plant forms and density.We defined the land exposure of a grid cell based on the typical terrain features of the ecosystem types described by Mackenzie and MacHutchon (1996).We divided vegetation into six plant forms (coniferous trees, deciduous trees, medium to tall shrubs (≥ 0.5 m), low shrubs (< 0.5 m), sedges or grasses with low shrubs, and sedges or grasses).The effects of plant density were estimated based on leaf area indices (LAI).
Since geothermal observations were sparse in Arctic regions, all the grid cells in the park were assumed to have the same geothermal heat flux (0.08 W m −2 ) and the same thermal conductivity for the bedrock (2.28 W m −1 • C −1 ), estimated based on site observations (Pollack et al., 1993;Jessop et al., 1984;Majorowicz et al., 2004).

Leaf area indices
Leaf area index was downscaled to 10 m spatial resolution using the 30 m resolution foliage biomass from Landsat imagery where L 10 m, i is the estimated LAI for a 10 m pixel i (cm 2 leaf cm −2 land), E i and E j are the ecosystem types for the 10 m pixels i and j , respectively.j ranges from 1 to 9 for the nine 10 m pixels in the 30 m Landsat pixel.B 30 m is the foliage biomass of the 30 m pixel (g cm −2 ).B 30 m is calculated using a regression equation of simple ratio vegetation index (band-4/band-3) of Landsat imagery (acquired on 30 August 2001) (Chen et al., 2013).B 0 (E j ) is the average foliage biomass (g cm −2 ) for an ecosystem type E j , which was estimated as the average of all the 30 m pixels in the park where E j is the dominant ecosystem type.S(E i ) is the specific leaf area (cm 2 g −1 ) of the ecosystem type (122 cm 2 g −1 for sedges; 154 cm 2 g −1 for deciduous shrub; 74 cm 2 g −1 for evergreen shrubs; and 50 cm 2 g −1 for black spruce) (Kattge et al., 2011).Some steep areas were blocked by shadows in the Landsat imagery.We estimated LAI in these areas as the average LAI of the ecosystem type B 0 (E i )S(E i ). Figure 2c shows the distribution of LAI in the park.LAI was discretized into 18 classes (increasing exponentially from 0 to 4.0) in modelling.

Computation procedure
The model was run for all the different combinations of the input-data classes (ecosystems types, climate-topography clusters, and LAI classes) for the coastal and inland regions, respectively.One combination of the input-data classes usually represents the input-data conditions of more than one grid cell.Except water bodies and aufice, there are 49 441 and 258 633 different combinations of input-data classes for the coastal and inland regions, respectively.The total number of these input-data combinations was only 0.3 % of the total number of the land grid cells in the park at 10 m resolution (There are 107 million land grid cells at 10 m resolution in the park excluding water and aufice).Thus, the computation  time was reduced significantly for this high-resolution spatial modelling.The daily climate data observed at the Old Crow and Komakuk Beach stations were from 1968 and 1958, respectively.For spatial consistency, we modelled the permafrost in the park from 1968.To initialize the model, we selected the daily climate data in a year (1972 and 1964 for the Old Crow and Komakuk Beach stations, respectively) which had a climate close to the climate extrapolated for the year 1967.We ran the model iteratively using these years' climate until the modelled ground temperature was stable, then ran the model from 1968-2011 based on the observed daily climate.
After that, we ran the model for the two climate scenarios to 2100.Across the park, active-layer thickness generally became thinner from the south to the north.The active layer was thinner in most of the valleys and in coastal regions as well.A thicker active layer was evident in the upper mountains and along drainage channels in some valleys.The results also show significant differences among ecosystem types (Fig. 6c).Alpine slopes and rock-lichen areas had thick active layers because they have no organic layer on the surface and usually contain high fractions of gravels in the soil.Active floodplains (types 23-26) and drainage channels (types 11-13) had thicker active layer as well due to their wet conditions and high gravel content.Pediment slopes (types 14-17) and graminoid wetland (type 21) had thin active layer due to the thick peat layer and wet conditions.

Spatial distribution of active-layer thickness
We calculated the average active-layer thickness for each climate-topography cluster (Fig. 6).The results show that the average active-layer thickness increased steadily with summer air temperature.Figure 6a also shows a slightly different trend for the coastal region, where the active layer was thinner due to the peaty and wet soil conditions in this region.
Active-layer thickness also increased slightly with annual insolation (Fig. 6b), but the correspondence is scattered due to the effects of other factors, especially the ecosystem types and air temperature.Active-layer thickness decreased with LAI due to the shading effects of plants (Fig. 6d).
To compare the relative importance of the input variables on the modelled spatial distribution of active-layer thickness in the park, we calculated the average absolute intra-class deviation and inter-class deviations based on the input variables of ecosystem types, climate-topography clusters, and LAI classes (Table 4).The average inter-class deviation is a measure for the differences among classes, while the average intra-class deviation is a measure for the variability within each class.When the averages were not weighted by the land areas of the classes, the inter-class deviations and the ratios between the inter-class and intra-class deviations were the largest for the ecosystem types, followed by  climate-topography clusters, and then by LAI classes.This result indicates that ecosystem type was the most important factor differentiating the modelled active-layer thickness in the park whereas LAI was the least important factor.When the averages were weighted by the land areas of the classes, ecosystem type was still the most important factor control-ling the spatial distribution of active-layer thickness in the park.The importance of the LAI was close to that of the climate-topography because the variation of LAI was associated with the ecosystem type and climate.Ecosystem type was the dominant factor mainly because the corresponding ground conditions are the most important factor controlling active-layer thickness in this region (Wang et al., 2010).

Changes in active-layer thickness and permafrost distribution
The model results show that the average active-layer thickness in INP increased by 0.1 m (13 %) from the 1970s to the 2000s (i.e., 2000-2009).The increase was more significant in higher elevations in the southern part of the park.From the 2000s to the end of the 21st century, the modelled average active-layer thickness increased by 0.3 m and 0.9 m (34 % and 99 %) under the scenarios of CCCma and ECHAM, respectively (Fig. 7a).The difference between these two climate scenarios was significant.Under the CCCma scenario, the relative changes in active-layer thickness were mainly between 20 % and 40 % across the park.Under the ECHAM scenario, however, the relative changes in active-layer thickness was larger and with a wider range (mainly from 70 % to 150 %) (Fig. 5d).Large inter-annual variations existed due to climate fluctuations, especially under the ECHAM scenario (Fig. 7a).With progressive deepening in summer thaw depth, the deep layer thawed in summer may not be frozen back in the following winter, thus forming a year-round unfrozen layer, called talik, between the permafrost and the seasonal thawing/freezing layer.Currently permafrost underlies almost all the land areas in the park and very small areas contain taliks (in the 2000s, the land area with taliks and the area without permafrost in the park accounted for only 0.0007 % and 0.004 %, respectively).From the mid-21st century, more areas were predicted to develop taliks and permafrost disap-peared in some southern areas.By the end of the 21st century, taliks developed in 0.1 % and 13.7 % of the land area in the park under CCCma and ECHAM, respectively, and permafrost disappeared in 2.1 % and 5.2 % of the area in the park under these two scenarios.Talik development occurred mainly in the southern areas of the park (Fig. 5c and d).Taliks could develop in most of the ecosystem types, except in pediment slopes and graminoid wetland.Permafrost could disappear in forests and in some southern valleys, mainly in floodplains and inactive alluvial terraces.

Permafrost degradation stages
Zhang (2013) divided permafrost degradation process induced by climate warming into five stages: gradualthawing stage, increased-thawing stage, frequent-talik stage, isothermal-permafrost stage, and permafrost-free stage.We determined the stages for each grid cell in the park based on the model results.From 1967 to 2011, 8.7 % of the permafrost area in the park transferred from gradual-thawing stage to the increased-thawing stage, with very small areas (< 0.01 %) in frequent-talik and permafrost-free stages (Fig. 7b).By the end of the 21st century, the model predicted that more areas would be in higher degradation stages, although significant differences existed between the two climate scenarios.Under scenario CCCma, the area with gradual-thawing permafrost would be reduced to 84.8 % of the area in INP by the end of the 21st century, while the area with increased-thawing stage and permafrost-free area would expand to 12.8 % and 2.1 % of the area in INP, respectively.More significant degradation was predicted under the ECHAM scenario.By the end of the 21st century, only 19.9 % of the area in INP would contain gradual-thawing permafrost.Permafrost in increased-thawing stage would expand to 53.0 % of the area in INP.Permafrost in other degradation stages would also expand under this scenario (Fig. 7).Spatially, permafrost in increased-thawing stage and permafrost-free areas are mainly located in some valleys and south-facing slopes under scenario CCCma.Under scenario ECHAM, most of the northern area in the park would be in increased-thawing stage while most of the southern area would be in frequent-talik and isothermal-permafrost stages, and some areas would be permafrost free.

Comparison with observations and results from other studies
Figure 8 shows a comparison between the modelled and measured summer thaw depth at 100 sites.The modelled summer thaw depth was comparable to the observations for most of the sites, although large differences exist for some sites, especially the seven sites in the circle in Fig. 8, where the model significantly over-estimated the observations.These seven sites are in alpine slopes containing widespread rocks on the land surface and in the soil.The observed thin summer thaw depths at these sites are probably due to local variations in soil conditions or mistaking rocks for frozen ground.Summer thaw depths greater than one meter were observed in some alpine slopes although most observations did not reach the frozen layer (therefore they are not included in Fig. 8).
Model tests also suggested that the active layer was thick in alpine slopes due to lack of organic layer, high gravel content and sparse vegetation (Wang et al., 2010).The correlation coefficient between the modelled and observed summer thaw depth is low (R = 0.56, N = 93), even excluding the seven sites.The low correlation is mainly due to variations of ground conditions within an ecosystem type.Therefore, more spatially detailed maps for ground conditions are needed for high-resolution mapping of permafrost.
Based on field observations, Tarnocai (1986) created a polygon map of active-layer thickness for the INP area.The mapped median active-layer thickness is 0.4 m (0.3-0.5 m) in most northern coastal regions and 0.75 m (0.7-1.0 m) in the northern deltas.The median active layer is 0.7 m (0.5-0.9 m) along the large valleys and 1.35 m (0.5-2 m) in inland mountain regions.This spatial pattern and the range of the activelayer thickness are similar with our modelled results although his map was represented by only 15 polygons within the park.The modelled permafrost in this region is continuous, which is in agreement with the permafrost map of Canada (Heginbottom et al., 1995).The modelled permafrost thickness ranged from 150 m to 300 m, which is similar with the observations in this region (Smith and Burgess, 2002).Burn and Zhang (2009).The modelled active-layer thickness in 2012 is the average modelled based on the two climate scenarios for that year.Burn and Zhang (2009) reported changes in active-layer thickness and ground temperature in Herschel Island just outside of the coastal region.By calibrating a model using detailed field observations, they estimated that annual mean temperature at 1 m and 20 m depths had increased by 2.6 • C and 1.9 • C, respectively, from the beginning of the 20th century (1899)(1900)(1901)(1902)(1903)(1904)(1905) to 2006, and the perturbation of ground temperature had reached about 120 m depth although the change was less than 0.1 • C. We checked the modelled ground temperature profiles for similar land cover conditions (sedge tussock) in the coastal region.From 1967 (assuming in equilibrium) to 2006, the modelled ground temperature increased by 3.4 • C and 2.5 • C at depths of 1 m and 20 m, respectively.These values are larger than that of Burn and Zhang (2009) although we considered a shorter time.This is probably because the climate data selected for initializing the model have a lower precipitation.And the increasing trends in precipitation and snow depths in the last 45 yr could not represent the trends in the entire 20th century since higher precipitation was recorded during 1958-1967 (Fig. 3c) (Fig. 3c).Our model results show that ground temperature at 120 m has increased by 0.04 • C, which is similar to the result of Burn and Zhang (2009).Burn and Zhang (2009) also measured summer thaw depth in Herschel Island from 2003 to 2007.Its magnitude and inter-annual fluctuations are comparable with our modelled active-layer thickness except in 2007, in which the observation shows a significant increase (The squares in Fig. 9).Burn and Kokelj (2009) reported a long-term record of summer thaw depth since 1983 measured at Illisarvik in Richards Island (69.48 • N, 134.59 • W) (The red dot in Fig. 1b), probably the longest continuous measurements of summer thaw depth available.The modelled active-layer thickness is closely correlated with these measurements (R = 0.71, N = 30) although the absolute values are somewhat different due to the differences in climate and ground conditions between the observation sites and the coastal region in INP (Fig. 9).The modelled deepening trend of active-layer thickness was 2.9 cm per decade from 1983 to 2012, which was similar to that of the observations (deepening 2.3 cm per decade).

Discussion
This study mapped permafrost distribution and changes with climate in an Arctic region with complex terrain at 10 m spatial resolution using a process-based model.Compared to previous studies, several features of this work are noteworthy.First, this study used a higher spatial resolution than previously published modelling and mapping work, such as the 30-m resolution modelling studies (Duchesne et al., 2008;Zhang et al., 2012), the 2 km-resolution study (Jafarov et al., 2012), or the half-degree latitude/longitude or even coarser resolution modelling studies (Zhang et al., 2006(Zhang et al., , 2008)).In this high-resolution study, satellite remote sensing data were used more effectively; topographic effects on solar radiation can be calculated more accurately; the plot size of field observations (large uniform plots are usually hard to find in regions with complex terrain) is comparable to the grid size, therefore they can be used directly for model calibration and validation; and the results reveal more spatial details, and are therefore more suitable for land-use planning and ecological monitoring and assessment.
To test the effects of spatial resolution on the modelled permafrost conditions, we selected an area of 20 km by 20 km (comparable to the grid size of half-degree latitude/longitude, which is about 20 km by 55 km in this region) and run the model at 100 m, 1 km, 5 km and 20 km resolutions based on the dominant land cover type in each grid cell.Coarser spatial resolution not only revealed less spatial details (Fig. 10), the modelled result also had a bias towards the permafrost conditions of the dominant land cover type.In this testing area, the dominant land cover type is rock-lichen, which has a deeper active layer than most of other types.With the increase in grid size, the modelled average active-layer thickness in the area increased by 22 %, and the spatial variations among the grid cells decreased from 0.4 m to 0 m from 10 m resolution to 20 km resolution (Fig. 11).The temporal patterns of the modelled average active-layer thickness for the area were similar among different spatial resolutions (R ranged from 0.97 to 1.00.N = 134).
Secondly, although several studies have mapped permafrost in mountainous regions at high spatial resolution and considered topographic effects on insolation (e.g., Bonnaventure et al., 2012), they assumed that permafrost conditions are in equilibrium with climate (see reviews by Riseborough et al., 2008).Ground temperature observations and modelling studies have shown that the current and future permafrost conditions are not in equilibrium with the atmospheric climate (Osterkamp, 2005;Zhang et al., 2008).The modelled talik development and different degradation stages in this study are examples of the disequilibrium response of permafrost to climate change.
Thirdly, although permafrost will persist in most of the park in the 21st century, the model results show changes in permafrost degradation stages with climate warming.Compared to fluctuating active-layer thickness, the degradation stages clearly categorized the degradation processes of permafrost with climate change (Fig. 7).Our sensitivity tests also show that the effects of snow conditions on active-layer thickness were dependent on degradation stages.When permafrost is in the gradual-thawing stage, changes in snow drifting condition had almost no effects on active-layer thickness.When permafrost is in increased-thawing stage or taliks developed, reducing snow drifting (i.e., with deeper snow accumulated) would increase summer thaw depth.
And finally, this study estimated daily climate variations based on station observations and used a clustering method to consider climate distribution and topographic effects on insolation.The estimated climate is directly linked to station observations and the clustering approach significantly reduced computation time for high-resolution spatial modelling.Such methods could be useful for other process-based high-resolution spatial modelling, such as for ecological and biogeochemical processes.
This study also indicates several sources of uncertainties or information gaps for permafrost mapping.One major source of uncertainty is the information about ground conditions.Although satellite data now can provide high spatial resolution maps for land cover and vegetation conditions, soil maps are still very coarse, especially for Arctic regions.We estimated soil and drainage conditions based on ecosystem types.Soil and drainage conditions, however, could differ widely within an ecosystem type.Therefore the results can only represent permafrost conditions under typical ground conditions of the ecosystem types.Improving the information about soil and drainage conditions is important not only for permafrost mapping but also for other applications, such as ecosystem assessment and infrastructure development under a changing climate.Another major uncertainty is the future climate.The results of this study show significant differences in projected permafrost conditions under the two climate scenarios.Under the warmer scenario ECHAM, the projected changes in active-layer thickness not only showed a larger increase but also with a wider range of fluctuation, and the areas with taliks or even disappearance of permafrost are much larger than that of the cooler scenario CCCma.The sparse and relatively short climate records in Arctic regions also affect the accuracy of the model results.The climate of the first year selected to initialize the model could affect the simulated permafrost thickness and the long-term trends of ground thermal conditions.
Although we considered snow drifting and lateral water flows, the NEST model is still one-dimensional, assuming each grid cell to be uniform and the lateral heat flux can be ignored.Therefore the results do not represent areas with strong lateral heat fluxes, such as sharp mountain peaks or edges, and the areas very close to waterbodies.Climate change and fire disturbances have significant impacts on vegetation conditions and LAI, which will affect permafrost through their effects on snow drifting and energy exchanges between the land surface and the atmosphere.In this study, however, we did not consider changes in LAI caused by disturbances or plant growth with climate change.Our sensitivity tests suggest that increasing LAI without changes in the snow drifting factor will reduce active-layer thickness due to shading effects.Reducing snow drifting or even capturing some snow from surroundings will cause faster thawing from the bottom of permafrost.Its effects on active-layer thickness depend on the stages of permafrost degradation.Therefore changes in vegetation and associated snow drifting could affect permafrost conditions in different ways.

Conclusions
This study mapped climate change impacts on permafrost at high spatial resolution for an Arctic region with complex terrain using a process-based model.Ecosystem types and leaf area indices were mapped using satellite data.Soil and ground conditions were estimated based on ecosystem types and field observations.Solar radiation was calculated explicitly based on topographic attributes.Spatial distributions of monthly average climate were interpolated based on station observations and elevation.The temporal variations of climate were estimated based on the patterns of station observations.The computation time was reduced significantly by clustering climate distribution and topographic effects on insolation into discrete types.The high-resolution modelling results show both the large-scale transitions across the park and landscape-scale variations due to topography, ecosystems and vegetation conditions.Such high-resolution results are more useful for land-use planning and ecosystem assessment than previously developed coarse resolution permafrost maps, which not only show less spatial details and variations, but also have a bias towards the permafrost conditions of dominant land cover types.The mapped active-layer thickness and permafrost distribution were generally comparable to field observations and other studies.The results show deepening in active-layer thickness and progressive degradation of permafrost, although permafrost will persist for most The Cryosphere, 7, 1121-1137, 2013 of the region during the 21st century.This study also shows that ground conditions and climate scenarios are the major sources of uncertainty for high-resolution permafrost mapping.

Fig. 1 .
Fig. 1.(a) Relief map of Ivvavik National Park and summer thaw depth observation sites (red dots), and (b) location of the park.Green squares are the locations of climate stations used in the study.The red dot in (b) is the location of Illisarvik in Richards Island.

Fig. 2 .
Fig. 2. Distribution of (a) mean annual air temperature, (b) mean annual total precipitation, (c) leaf area indices, and (d) ecosystem types in Ivvavik National Park.Shaded relief was added for easy interpretation.The air temperature and precipitation were averages from 1971 to 2000 derived based on Wang et al. (2006).The leaf area indices were estimated based on Landsat imagery.The ecosystem map was from Fraser et al. (2012) (see Table2for the names of the ecosystem types).

Fig. 3 .
Fig. 3. Changes of (a, b) annual mean air temperature and (c, d) annual total precipitation at the Komakuk Beach and Old Crow climate stations.The observations at these two stations were used to represent the temporal patterns of climate in the coastal and the inland regions, respectively.

Fig. 4 .
Fig. 4. The conditions of the clusters in the coastal region.(a) annual total insolation of the clusters, (b) slopes and aspects of the clusters corresponding to the annual total insolation of the clusters, (c) annual total degree days when air temperature is above 0 • C and annual total precipitation of the clusters, and (d) viewsheds of some selected clusters (the legend shows their corresponding annual total insolation (MJ m −2 yr −1 )).Aspect is 0 • for a slope facing north and it increases anti-clockwise.In (b), aspect was shown as 360 minus the aspect if it is larger than 180 • .

Figure
Figure 5a and b show the modelled spatial distribution of active-layer thickness in the 2000s (i.e. from 2000 to 2009).Across the park, active-layer thickness generally became thinner from the south to the north.The active layer was thinner in most of the valleys and in coastal regions as well.

Fig. 5 .
Fig. 5. Modelled distribution of (a) average active-layer thickness in the 2000s (i.e., 2000-2009), (b) an enlarged area for the black square in (a), and relative changes of active-layer thickness from the 2000s to the 2090s under the climate scenarios of (c) CCCma and (d) ECHAM.Shaded relief was added for easy interpretation.

Fig. 6 .
Fig. 6.Correspondence of the modelled average active-layer thickness in the 1970s to (a) the average summer air temperature (June-September), (b) annual solar radiation, (c) ecosystem types, and (d) leaf area indices.The range bars in (c) and (d) are the average absolute deviations within the classes.

Fig. 7 .
Fig. 7. Changes in active-layer thickness and degradation stages.(a) Mean active-layer thickness.(b) and (c) are the percentages of land areas in different degradation stages in the park under the climate change scenarios of CCCma and ECHAM, respectively.The degradation stages are based on Zhang (2013).

Fig. 8 .
Fig. 8.Comparison between the modelled and measured summer thaw depth.The model significantly over-estimated the seven sites in the circle, which are rock-lichen and alpine slopes, probably due to spatial heterogeneity or mistaking rocks for frozen ground.The solid line corresponds to the linear regression equation shown in the panel.The dashed line is a 1-to-1 line for reference.

Fig. 9 .
Fig. 9. Comparison between the modelled active-layer thickness in the coastal region and measured summer thaw depth at Illisarvik in Richards Island (1983-2012.The circles) and in Herschel Island (2003-2007.The filled squares).The measurements at Illisarvik in Richards Island during 1983-2008 are from Burn and Kokelj (2009), and the unpublished data during 2009-2012 are kindly provided by Christopher Burn.The measurements in Herschel Island are fromBurn and Zhang (2009).The modelled active-layer thickness in 2012 is the average modelled based on the two climate scenarios for that year.

Fig. 10 .
Fig. 10.The modelled spatial distributions of active-layer thickness (averages during the 2000s) at different spatial resolutions for a 20 km by 20 km area.The spatial resolutions are indicated in (a-g).The location of the area in the park (the red square) is shown in (h).To show the difference between (a) and (b), a small area (the black squares in (a) and (b)) was enlarged and shown in (d) and (e).

Fig. 11 .
Fig. 11.The effects of spatial resolution on the modelled activelayer thickness for a 20 km by 20 km area shown in Fig. 10h.The average is calculated for all the grid cells in the area during 1967-2100.The spatial standard deviation was calculated based on all-year averages (1967-2100) of the grid cells.The future years (2012-2100) was modelled based on scenario ECHAM.The patterns are very similar using the scenario CCCma.

Table 2 .
Mackenzie and MacHutchon (1996)ystem types in Ivvavik National Park.More detailed descriptions of the terrain features, vegetation compositions, and soil conditions can be found inMackenzie and MacHutchon (1996).

Table 3 .
Ground conditions defined for each ecosystem type.Type Drainage a Moisture b Sites c Org. layer (cm) d Mackenzie and MacHutchon (1996)sent very poorly, poorly, imperfectly, moderately well, well, rapidly, and very rapidly, respectively, defined based on Agriculture Canada Expert Committee on Soil Survey(1987).The data were fromMackenzie and MacHutchon (1996).bSoilmoistureregimes0 to 8 represent veryxeric, xeric, sub-xeric, sub-mesic, mesic, sub-hygric, hygric, sub-hydric, and hydric, respectively,  defined based on Walmsley et al. (1980).The data were fromMackenzie and MacHutchon (1996).cThenumber of sites observed during 2008-2011 for estimating the ground conditions.dThetoporganiclayerthickness.eObservedrange of top organic layer thickness during 2008-2011.fThetoporganiclayerthicknessdefinedforthemodel input.gSoiltexturesymbolsareSforsand,Lfor loam, Si for silt, and C for clay.hObservedsoiltexture was fromMackenzie and MacHutchon (1996)and from our fieldwork during2008-2011.iThesoiltexturedefinedfor the model input.jThesoilorganicmattercontent in the top mineral soil.It was defined based on observations during2008-2011.kThedepth(cm)wheregravels first appears and the content of the gravel (%) at this depth.Gravel content gradually increases with depth after this layer.The data were estimated fromMackenzie and MacHutchon (1996)and our fieldwork during2008-2011.

Table 4 .
The average absolute intra-class deviation and the average absolute inter-class deviation of the active-layer thickness in the 1970s calculated based on the ecosystem types, climate-topography clusters, and LAI classes.
*In percentage of the average active-layer thickness of the type or class.