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

Introduction Conclusions References


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 Figures 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., 2008a).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., , 2008a)).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 management and for ecological monitoring and 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 conducted high spatial resolution permafrost modelling.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)  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 conditions and soil moisture (Zhang et al., 2012).This region mainly consists of a peatland plain, where peat layer thickness can be estimated based on elevation.However, northern regions are usually not flat and organic layer thickness can not be estimated readily from elevation.Complex terrain has significant effects on climate and ecosystems, and thus has high spatial variations in permafrost conditions.Therefore, spatially detailed information about permafrost conditions in complex terrain is very important 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 a region with complex terrain, and to identify the major controlling factors and input data gaps for future studies.
2 Methods and data

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.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 in the north and northeast of the park (Fig. 1).The park straddles the transitional zone between the arctic tundra and subarctic boreal forest.
Currently, permafrost is distributed continuously in the park (Heginbottom et al., 1995).Only a small northeastern portion of the park was covered by ice during the last glacial period.
Three major field data sources were used for this study.One is the resources report for the park by Canada Parks Service (1993), which provides general information about the climate, geology, geomorphology, soil, hydrology, vegetation and wildlife in the park.Figures

Back Close
Full 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 conditions, vegetation composition and status, ecosystem types, 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 measured summer thaw depths were used to validate the model results.

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-vegetationatmosphere 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 liquid water, are determined based on energy conservation.Detailed description of the model and validations can be found in Zhang et al. (2003Zhang et al. ( , 2005Zhang et al. ( , 2006Zhang et al. ( , 2008b)).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 (Appendix A).Figures

Back Close
Full

2.3
The input data and processing

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 interpolate monthly climate based on station observations (e.g., Wang et al., 2006;McKenney et al., 2011).Since the temporal patterns are usually affected by the changes in climate stations with time used for the interpolation, we estimated the spatial distribution of long-term monthly averages for air temperature and precipitation based on Wang et al. (2006) (Fig. 2) and estimated the temporal patterns using observations at climate stations.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 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 = 15 706).However, the monthly total precipitation is correlated between the stations (R = 0.50, n = 516).Therefore we used station observed daily climate data to represent the temporal patterns of the climate for a region surrounding the climate station.
The coastal region of the INP neighbouring the Arctic Ocean has a marine climate while the inland region has a continental climate (Canada 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 at nearby stations.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 0m,g (M) is the 30-yr (1971-2000) mean air temperature in the month M for the grid cell, interpolated from climate stations based on Wang et al. (2006) using a spatial resolution of 30 m by 30 m in INP.∆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 monthly air temperature difference between the grid cell and the climate station in month M year Y .Monthly total 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 (CCCma (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  simulated by the same model, and then we used these relative changes to estimate the future monthly climate for all the grid cells Figures 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, respectively (Fig. 3).The monthly air temperature and precipitation were converted to daily data based on the above described method (Eqs. 1 and 2) using the daily observations at the two climate stations from 1968 to 2011.
Vapour pressure was estimated based on minimum air temperature (Zhang et al., 2012).Daily total insolation without topographic effects (an input to consider the 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).

Classifying the average climate and topographic effects on insolation
Since the average spatial patterns 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 variations.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 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. 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.Average annual insolation was calculated for each grid cell (Appendix A).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 equation.Although we did not use annual mean air temperature and the annual total degree-days when air temperature is less than 0 • C (TDD T<0 ) in clustering, their errors were small as well (Table 1).
To reduce the errors in estimated insolation, 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

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 by Mackenzie and MacHutchon (1996) using air photos and detailed field observations at 367 sites.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 peat thickness, the texture of mineral soils, organic mater content in the mineral soils, coarse fractions, 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 coarse materials (pebbles and even boulders).We estimated the coarse fractions based on the field observations and the biophysical and terrain features of the ecosystem types described by Mackenzie and MacHutchon (1996).The thickness of the subsoil was estimated assuming that the Introduction

Conclusions References
Tables Figures

Back Close
Full coarse fraction increases 2 % for every 10 cm until the coarse 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 moisture regime 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 water table is above the land surface (5 % of the water depth above the land surface each day).
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 Full where L 10m,i is the estimate 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).
The equation was developed based on Landsat imagery and field observations.B 0 (E j ) is the average foliage biomass (g cm −2 ) for an ecosystem type E j 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 ).LAI was discretized into 18 classes (increasing exponentially from 0 to 4.0) in modelling.

Procedure for modelling and mapping permafrost
The erally became thinner from the south to the north.Active-layer was thinner in most of the valleys and in north and northeast coastal regions.Thicker active-layer was evident in upper mountains and along drainage channels in some valleys.The results also show significant difference among ecosystem types (Fig. 6c).Alpine slopes and rocklichen areas have thick active-layers because these ecosystems 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) have thicker active-layer as well due to their wet conditions and high gravel content.Pediment slopes (types 14-17) and graminoid wetland (type 21) have thin active-layer due to the thick peat layer and wet conditions.
We calculated the average active-layer thickness for each climate type and LAI class (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 is thinner due to the peaty and wet soil conditions and the marine climate in this region.Active-layer thickness also increased with annual insolation (Fig. 6b), but the correspondence is scattered due to the effects of Figures 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 types, 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 types and then by LAI classes.This result indicates that the ecosystem types were 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, the ecosystem types were still the most important factor controlling the spatial distribution of active-layer thickness in the park.The importance of the LAI was close to that of the climate because the variation of LAI was associated with the ecosystem types and climate.Ecosystem types were the dominant factor mainly because the corresponding ground conditions are the most important factors controlling active-layer thickness in this region (Wang et al., 2010).

Changes in active-layer thickness and permafrost distribution
The model results suggest that active-layer became thicker from the 1970s to the 2000s (i.e., [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009], especially in higher elevations in the southern part of the park (Fig. 5a, b).Average active-layer thickness in INP increased by 13 % (0.1 m) from the 1970s to the 2000s.Figure 7 shows the temporal changes in average active-layer thickness and permafrost conditions in INP.The average active-layer thickness increased steadily under the CCCma scenario and increased significantly under the ECHAM scenario.From the 2000s to the end of the 21st century, the modelled average active-layer thickness increased by 34 % and 99 % (0.3 m and 0.9 m) under the scenario of CCCma Introduction

Conclusions References
Tables Figures

Back Close
Full and ECHAM, respectively.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 with relatively small spatial differences.Under ECHAM scenario, however, the relative changes in active-layer thickness was larger and with a wider range (mainly from 70 % to 150 %).Large inter-annual variations existed due to climate fluctuations, especially under the ECHAM scenario.
The model results showed that permafrost currently underlies almost all 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 account for only 0.0007 % and 0.004 %, respectively).From the mid 21st century, more areas were predicted to develop taliks and permafrost disappears in some southern areas (Fig. 7).By the end of the 21st century, taliks develop in 0.1 % and 13.7 % of the land area in the park under CCCma and ECHAM, respectively.Permafrost disappears in 0.9 % and 5.1 % of the area in the park under these two scenarios (Fig. 7).Talik development occurs mainly in southern areas in the park (Fig. 5).Taliks could develop in most of the ecosystem types, except in pediment slopes (Types 14-17) and graminoid wetland (type 15) (Fig. 8).Permafrost could disappear in forests (types 5-8) and in some southern valleys, mainly in floodplains (types 24-26) and inactive alluvial terraces (types 19 and 20) (Fig. 8).

Comparison with observations and results from other studies
Figure 9 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. 9, 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 ground conditions or mistaking rocks for frozen ground.Greater than one meter summer thaw depths were observed in some alpine slopes although the top of the permafrost were not reached in most of the sites.Model Introduction

Conclusions References
Tables Figures

Back Close
Full tests showed that active-layer thickness is thick in alpine slopes due to lack to organic layer, high gravel content in soil and sparse vegetation.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 (ground material and drainage 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 active-layer 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 to 300 m, which is similar with the observations in this region (Smith and Burgess, 2002).

Discussion and conclusions
This study mapped permafrost distribution and changes with climate in a region with complex terrain at 10 m spatial resolution using a process-based model.The high resolution results show both the large-scale transitions across the park and landscapescale variations due to topography, ecosystems and vegetation conditions.Topography induced differences in climate and insolation are very important for permafrost distribution and changes with climate warming.The mapped active-layer thickness and permafrost distribution were generally comparable to field observations and other studies.
This study demonstrates that it is practical to model and map permafrost distribution and changes with climate at high spatial resolution for regions with complex terrain.

TCD Figures Back Close
Full Such high resolution results are more useful for land management and ecosystem assessment than previously developed coarse resolution permafrost maps (e.g., Zhang et al., 2006Zhang et al., , 2008a)).For example, in our previous half-degree latitude/longitude modelling study (Zhang et al., 2006), INP was covered by about ten grid cells and the modelled active-layer thickness was similar among the grid cells.
Several features of this study are noteworthy.First, this study used a higher spatial resolution than previously published modelling and mapping work, such as the 30 m resolution modelling work (Duchesne et al., 2008;Zhang et al., 2012), the 2 km resolution work (Jafarov et al., 2012), the half-degree latitude/longitude or even coarser resolution modelling work (Zhang et al., 2006(Zhang et al., , 2008a)).A high spatial resolution can more effectively use satellite remote sensing data, DEM, and field data, and can provide more detailed spatial results.Secondly, although several studies have mapped permafrost in mountainous regions at high spatial resolution and considered topographic effects on insolation, they assumed that permafrost conditions are in equilibrium with the 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., 2008a).The modelled talik development in this study is an example of the disequilibrium response of permafrost to climate change.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 observations and the clustering approach significantly reduced the computation time for high resolution modelling.Such methods could be useful for other process-based high resolution spatial modelling, such as for ecological and biogeochemical processes.
The results of this study show significant differences in projected permafrost degradation under the two climate scenarios.Under the warmer scenario ECHAM, the projected changes in active-layer thickness not only showed a larger relative increase but also with a wider range, and the areas with taliks or even disappearance of permafrost Figures are much larger than that of the cooler scenario CCCma.This result shows that the uncertainty in climate scenarios has significant impacts on permafrost projections.Another information gap is about ground conditions.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 only represent permafrost conditions under typical ecosystem conditions.Improving the information about soil and drainage conditions in the north is important not only for permafrost mapping but also for other applications, such as ecosystem assessment and infrastructure development under a changing climate.
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 water bodies.
Although the effect of LAI on the spatial distribution of permafrost seems small in this area, temporal changes in vegetation conditions and LAI could affect snow drifting and energy exchanges between the land surface and the atmosphere, therefore affect permafrost conditions.In this study, however, we did not consider changes in LAI caused by disturbances or plant growth with climate change.
Appendix A

Calculating topographic effects on insolation
Topography affects insolation due to the slope and aspect of the local land surface and the shading effects of the surrounding mountains and hills.The shading effect was calculated based on the viewshed approach (Fu and Rich, 2000).A viewshed is the angular distribution of sky obstructed by the surrounding mountains and hills.The maximum obstruction zenith angle in an azimuth direction can be determined based on the elevation of the cells in that direction.Introduction

Conclusions References
Tables Figures

Back Close
Full the slope and aspect of the surface of the grid cell.R dir0 is the direct solar radiation accepted on a horizontal surface without obstruction from the surroundings (W m −2 ).
The irradiance received in a grid cell at any time in a day is the sum of the direct and diffuse solar radiation received in that grid cell (R diff plus R dir ).
It is difficult and usually unnecessary for the users to provide hourly climate data for long-term permafrost modelling.Therefore we estimated the diurnal variations of both direct and diffuse solar radiation using cosine functions based on Wang et al. ( 2002) where θ n is the zenith angle of the sun at noon (in radians), F 1,θn is the ratio between direct solar radiation at noon and the daily mean direct solar radiation on a horizontal surface, and F 2,θn is the ratio between diffuse solar radiation at noon and the daily mean diffuse solar radiation on a horizontal surface.R day0 is daily average insolation received on a horizontal surface (W m −2 ), and F diff0,day is the fraction of daily diffuse solar radiation in the daily total insolation received on a horizontal surface.R day0 can be calculated by where R day0 is the daily total insolation received on a horizontal surface (MJ m −2 day −1 ), which is the input for the model.D is day length (s) determined based on latitude and the day of the year.Wang et al. (2002) developed regression relationships to determine F 1,θn and F 2,θn based on θ n for a location.Our numerical simulation showed that the relationships Introduction

Conclusions References
Tables Figures

Back Close
Full changes with latitude, especially when the day length reaches 24 h.Therefore we determined F 1,θn and F 1,θn numerically for each day of a year based on the original equations from Wang et al. (2002).
where t is time (s).The daily fraction of diffuse radiation in daily irradiance received on a horizontal surface (F diff0,day ) can be estimated using a logistic function based on Boland et al. (2008).We slightly modified the parameters based on the daily insolation data observed at Inuvik climate station  Full  3. Ground conditions defined for each ecosystem types.They are estimated based on field observations and the soil features by Mackenzie and MacHutchon (1996).
Discussion Paper | Discussion Paper | Discussion Paper | 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.We mapped permafrost in the northwestern Hudson Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 longest continuous observations within and nearby INP.Observations from these two stations are available from 1958-2010 and 1968-2010, respectively.Data gaps are filled based on observations Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Slope, aspect and viewshed are the topographic attributes affecting insolation.Slope and aspect were calculated using 10 m resolution digital elevation model (DEM), which were re-sampled from 30 m DEM data.The DEM data are from Topographic Data Centre of Natural Resources Canada.The viewshed of a grid cell is the angular distribution Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 5a, b shows the modelled spatial distribution of active-layer thickness in the1970s and 2000s (i.e. from 2000 to 2009).Across the park, active-layer thickness gen- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

a
Drainage classes 0 to 6 represent very poorly, poorly, imperfectly, moderately well, well, rapidly, and very rapidly, respectively, defined based on National Soil Database (http://sis.agr.gc.ca/cansis/nsdb/detailed/name/drainage.html).b Soil moisture regimes 0 to 8 represent very xeric, xeric, sub-xeric, sub-mesic, mesic, sub-hygric, hygric, sub-hydric, and hydric, respectively, defined based on Walmsley et al. (1980).c Soil texture symbols are S for sand, L for loam, Si for silt, and C for clay.d SOM is the soil organic mater content in the top mineral soil layer.Discussion Paper | Discussion Paper | Discussion Paper |

a
In percentage of the average active-layer thickness of the type or class.Discussion Paper | Discussion Paper | Discussion Paper |

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

Table 4 .
The average absolute intra-class deviation and the average absolute inter-class deviation of the active-layer thickness in the 1970s and their ratios for the input data of ecosystem types, climate types, and LAI classes.