A new approach to mapping permafrost and change incorporating uncertainties in ground conditions and climate projections

Spatially detailed information on permafrost distribution and change with climate is important for land use planning, infrastructure development, and environmental assessments. However, the required soil and surficial geology maps in the North are coarse, and projected climate scenarios vary widely. Considering these uncertainties, we propose a new approach to mapping permafrost distribution and change by integrating remote sensing data, field measurements, and a process-based model. Land cover types from satellite imagery are used to capture the general land conditions and to improve the resolution of existing permafrost maps. For each land cover type, field observations are used to estimate the probabilities of different ground conditions. A process-based model is used to quantify the evolution of permafrost for each ground condition under three representative climate scenarios (low, medium, and high warming). From the model results, the probability of permafrost occurrence and the most likely permafrost conditions are determined. We apply this approach at 20 m resolution to a large area in Northwest Territories, Canada. Mapped permafrost conditions are in agreement with field observations and other studies. The data requirements, model robustness, and computation time are reasonable, and this approach may serve as a practical means to mapping permafrost and changes at high resolution in other regions.


Introduction
About a quarter of the land in the northern hemisphere is underlain by permafrost (Zhang et al., 1999).The climate in northern high latitudes has warmed at about twice the rate of the global average during the 20th century, and continued warming is projected for the 21st century (ACIA, 2005).Climate warming causes ground temperature increasing, activelayer thickening, talik formation, and thawing of permafrost (Vaughan et al., 2013).These changes have significant impacts on infrastructure foundations, hydrology, ecosystems, and feedbacks to the climate system (ACIA, 2005).Mapping the distribution of permafrost and its possible changes with climate is important for land use planning, infrastructure development, ecological and environmental assessments, and modelling the climate system.
Different approaches have been used to map permafrost.Due to lack of field data in the early years, Wild (1882) delineated the southern boundary of permafrost in Russia using a simple model driven by air temperature.Later on, more field data were included in drawing the boundary of permafrost (see reviews by Nelson, 1989;Heginbottom, 2002;Shiklomanov, 2005).With more field experience and observations available, permafrost distribution was divided into continuous and discontinuous zones, and the latter was further divided into two to four subzones (Heginbottom, 2002).Permafrost zones were defined based on spatial continuity, i.e., the proportion of an area underlain by permafrost (Heginbottom et al., 1995;Brown et al., 1998).The most widely used permafrost maps are based on this concept, with boundaries delineated using long-term mean annual air temperature, field observations, and regional physiography Y. Zhang et al.: Permafrost mapping (Nikiforoff, 1928;Heginbottom et al., 1995;Brown et al., 1998).These maps are coarse in spatial resolution (usually without an explicit spatial resolution) with broad categories; they therefore have limited utility for practical land use planning and engineering applications or as a basis for assessing the impact of climate change on permafrost.
Permafrost distribution can be mapped at high spatial resolution by interpreting aerial photographs for some regions (Vitt et al., 1994;Reger and Solie, 2008).However, the association between land surface features and permafrost conditions has to be clear and interpretation needs good knowledge and experience.Several studies map permafrost distribution by classifying satellite images (e.g., Morrisey et al., 1986;Leverington and Duguay, 1997;Nguyen et al., 2009).Since permafrost cannot presently be directly imaged by optical satellite-based sensors, permafrost condition must be indirectly derived from permafrost-related and remote-sensing detectable geophysical or surface features.However, the relation between these features and permafrost conditions are complex and may change with space and time (Leverington and Duguay, 1997;Duguay et al., 2005).
Permafrost is a ground thermal condition rather than a substance.The impact of climate change on permafrost distribution was recognized a long time ago (Wild, 1882) and is now an important issue for northern land use planning, infrastructure development, and global climate projections (ACIA, 2005).Woo et al. (1992) estimated the shift in boundaries of continuous and discontinuous permafrost in Canada by assuming a 4-5 • C increase in surface temperature across the country.The Stefan solution, the frost index method (Nelson andOutcalt, 1983, 1987), Kudryavtsev's approach (Kudryavtsev et al., 1974), and the temperature at the top of permafrost (TTOP) model (Smith and Riseborough, 1996) are some of the major developments for mapping permafrost in recent years.These approaches integrate the effects of air temperature, snow and ground conditions in simplified ways and can be used to map permafrost spatially explicitly and to assess the effects of climate change (e.g., Anisimov and Nelson, 1997;Wright et al., 2003).For mountain regions, Haeberli (1973) developed a method to map permafrost using the basal temperature of snow (BTS) in the Alps and the method has been used in many highelevation areas.The method was further improved by including the effects of topography on air temperature, solar radiation, and snow conditions, and by incorporating field observations using permafrost probabilities for high resolution permafrost mapping (e.g., Lewkowicz and Ednie, 2004;Bonnaventure et al., 2012) and spatial sensitivity analysis to climate change (Bonnaventure and Lewkowicz, 2013).Recently, Gruber (2012) developed a global high-resolution (30 arcseconds, or < 1 km) permafrost map using relations between air temperature and probability of permafrost occurrence.Although these methods emphasized the impact of climate on permafrost, they assume that permafrost is in equilibrium with the atmospheric climate (therefore, they are categorized as equilibrium models by Riseborough et al., 2008).However, ground temperature observations show that current permafrost conditions are not in equilibrium (Osterkamp, 2005), and modelling studies show that the response of permafrost to climate warming during the 21st century will be transient and nonequilibrium (Zhang et al., 2008;Zhang, 2013).In addition, these methods use simple parameterization and statistic representation to consider the effects of ground conditions on permafrost (e.g., the seasonal n factor).They tend to underestimate the complexity and importance of ground conditions on permafrost dynamics, and the empirical parameters vary with time and space.
In the past two decades, process-based models have been used to quantify the impact of climate change on permafrost conditions and their spatial distributions.These models can integrate climate and ground variables based on energy and water transfer and can capture transient changes from seasonal to centennial, but they require detailed input data and lengthy computation time.Most spatial modelling studies, especially at national to hemispherical scales, have been conducted at half-degree latitude/longitude or coarser spatial resolutions (e.g., Lawrence and Slater, 2005;Marchenko et al., 2008;Zhang et al., 2008).At such resolutions, it is difficult to consider detailed spatial variations in vegetation and ground conditions, and the results are not suitable for land use planning and engineering applications.Recently, several studies have modelled and mapped permafrost at finer spatial resolutions, ranging from 2 km to 10 m (Duchesne et al., 2008;Jafarov et al., 2012;Zhang et al., 2012Zhang et al., , 2013)).In the North, however, the required maps for soil and ground conditions are coarse, with polygons in soil and surficial geology maps usually covering more than a hundred square kilometres.The lack of detailed soil and ground information is a major source of uncertainty for high resolution permafrost mapping (Zhang et al., 2012(Zhang et al., , 2013)).Another source of uncertainty is climate change projections.Projected permafrost conditions differ significantly under different climate scenarios (Anisimov and Reneva, 2006;Zhang et al., 2008Zhang et al., , 2012Zhang et al., , 2013;;Zhang, 2013), and consequently such results are difficult for decision-makers to use.
In this study we develop a new approach to mapping permafrost and change in an objective, replicable, and quantitative way while incorporating uncertainties in ground conditions and climate projections.This approach integrates remote sensing, field observations, and modelling and can map permafrost at high spatial resolution.The mapped permafrost probability is similar to the concept of the traditional permafrost zones.In this paper we describe the methodology of the approach and demonstrate its application for a large area in northern Canada.
Figure 1.The scheme of the new permafrost mapping approach.A land cover map classified from satellite images is used to represent the general land conditions and to improve the spatial resolution.For each land cover type, the probability of different ground conditions is estimated based on field observations.A process-based permafrost model (NEST) was used to calculate the evolution of permafrost for each ground type under three representative climate change scenarios generated by general circulation models (GCMs).From the model outputs, the probability of permafrost occurrence and the most likely permafrost conditions are determined.In the above scheme, X i,j is the model output based on ground type i and climate scenario j .p i is the probability of ground type i, and q j is the probability of climate scenario j .io is the most probable ground type.

A general description of the approach
Figure 1 shows the scheme of the new permafrost mapping approach.A land cover map from satellite images is used to represent the general land conditions and to improve the spatial resolution of the final products.For each land cover type, the probability of different ground conditions is estimated based on field observations.The evolution of permafrost was simulated using a process-based model for each ground type under three representative climate change scenarios (low, medium, and high warming).From these model outputs, the probability of permafrost occurrence and the most likely permafrost conditions are determined for each land cover type.The following is the rationale and feasibility of the approach.
Climate, especially air temperature, is the dominant factor controlling the spatial distribution of permafrost at large scales (from a hundred kilometres to continental).However, in a small area within which the climate is somewhat similar, permafrost can be present or absent; if present, the permafrost conditions can be significantly different from place to place (Shur and Jorgenson, 2007).Such local variation depends primarily on the distribution of vegetation and ground conditions including soil composition, snow, topography, and drainage (e.g., Brown, 1973;Shur and Jorgenson, 2007;Morse et al., 2012).High resolution land cover maps developed from satellite imagery can capture some general features of soil and hydrological conditions, therefore they can explain some of the major differences in permafrost conditions (Nelson et al., 1997;Nguyen et al., 2009;Jorgenson et al., 2010) especiallyopticalimages, only contain information about vegetation and near-surface conditions.Therefore ground conditions within a land cover type developed by satellite images can vary substantially.For example, hollows and hummocks are common at local scales and there are also topographically controlled variations at large scales.Such differences can cause significant variations in permafrost conditions (Duguay et al., 2005;Zhang et al., 2013).Field observations show that organic layer thickness (OLT) on the top of the mineral soil, including mosses and lichens, is a dominant factor affecting active-layer thickness and ground temperature (e.g., Harris, 1987;Kasischke and Johnstone, 2005;Johnson et al., 2013).Model sensitivity tests also indicate that OLT is the most sensitive factor affecting active-layer thickness and ground thermal conditions, especially when OLT is less than 50 cm (Yi et al., 2007;Zhang et al., 2012;Riseborough et al., 2013).Therefore we used OLT to divide soil conditions into several ground types within a land cover type.As current remote sensing technologies have difficulties measuring OLT, we estimated the probability distribution of OLT based on field observations for each land cover type.
Other input parameters for ground conditions, including mineral soil texture, organic matter content in mineral soils, gravel fractions, drainage parameters, and snow-drifting parameters, can be estimated based on land cover features, field observations, and soil and surficial geology maps.These parameters may be unique for each ground type or the same for all ground types within a land cover type.Land cover maps can be developed from satellite remote sensing images.Leaf area indices (LAI) and wildfire history can also be estimated using optical satellite images.Spatially detailed monthly climate data, necessary for high resolution spatial modelling, have been developed by interpolating station observations (Wang et al., 2006;McKenney et al., 2011).The temporal patterns of air temperature and precipitation can be estimated using observations at climate stations (Zhang et al., 2012(Zhang et al., , 2013) ) or from reanalysis of climate observations.Water vapour pressure can be estimated based on daily minimum air temperature; daily total insolation without topographic effects can be estimated using latitude, the day of the year, and the diurnal temperature range (Zhang et al., 2012).Topographical effects on solar radiation can be calculated using a digital elevation model (Zhang et al., 2013).To consider the uncertainty of future climate projections, three scenarios (low, medium, and high warming) can be selected to represent the medium and general range of the probable scenarios based on climate model projections.If we assume that the low and high warming scenarios represent the lower and upper quartiles of the probable scenarios and the medium warming scenario represents the two middle quartiles, the probabilities of the low, medium, and high warming scenarios would be 0.25, 0.5, and 0.25, respectively.
To map permafrost and its evolution, permafrost dynamics for each grid cell are modelled under the possible ground types and the three representative climate scenarios.The probability of permafrost occurrence in a grid cell can be determined based on the probabilities of the ground types and climate scenarios.In turn, the distributions of permafrost probability and the most likely permafrost conditions can be mapped.The probability of permafrost within a grid cell can be interpreted as a probability due to the uncertainties in ground conditions and climate projections, or can be interpreted as permafrost proportion or extent due to spatial heterogeneity within the grid cell, especially when the grid cell is large.

Estimating the probabilities of ground types within each land cover type
The probability distribution of ground types in all the areas of a land cover type was estimated using a modified logistic function based on field sampling observations of OLT at different locations.We did the modification because OLT cannot be negative and is usually thicker than a certain value for peatlands and some other northern land types. where where x is OLT (cm) and F (x) is the cumulative distribution probability of OLT for a land cover type.x 0 is the minimum OLT (cm) of the land cover type, µ and s (in centimetres) are the average and the standard deviation of the logistic distribution.Equation (1) was derived by modifying logistic distribution function considering x ≥ x 0 and F (x 0 ) = 0.
where a is a parameter.It equals 1 + 1/e x 0 based on the fact that F (x) equals 1 when x is infinite.The probability density can be derived as where f (x) is the probability density function.The probability of a ground type i can be determined from Eq. ( 5) where p(i) is the probability of a ground type i for OLT ranging from a i to b i .e ai and e bi are calculated from Eq. ( 2) when x equals a i and b i , respectively.i ranges from 1 to N , where N is the number of ground types.The thinnest and thickest OLT ground types can be determined from Eq. ( 6) using a predetermined probability level (e.g., 0.1).Then, the OLT ranges of the other ground types can be defined.The three parameters in Eq. ( 1) (x 0 , µ and s) can be determined by comparing the fitted cumulative probability with the observed relative cumulative frequency.Figure 2 shows an example of distribution of ground types in a land cover type.
For the model input, the OLT of a ground type was represented by 0.5(a i + b i ), except for the thinnest and thickest OLT ground types, which were determined as the OLT at which the cumulative probability F (x) equals a predetermined probability level (0.05 and 0.95 for the thinnest and thickest OLT ground types, respectively).Sensitivity tests showed that the modelled average permafrost conditions converged quickly with the increase in the number of ground types.Permafrost conditions were similar when the difference of the OLT was less than 5 cm, especially when OLT was thicker than 30 cm.Permafrost conditions were not sensitive to the changes in OLT when OLT was thicker than 50 cm (Zhang et al., 2012;Riseborough et al., 2013).Based on this response we can select the number of ground types and the thickness of each ground type.

The model
The permafrost condition and its changes with climate for each ground type were quantified using the Northern Ecosystem Soil Temperature model (NEST).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 water 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 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. ( , 2006)).Lateral flows and snow drifting are parameterized in a simplified way (Zhang et al., 2002;2012).The model can also consider the effects of topography on solar radiation (Zhang et al., 2013).

Study area and field observations
We applied this approach to a 94 km by 94 km area in Northwest Territories, Canada, centred at 62.77 • N, 114.27 • W (Fig. 3a).The area is located in the Great Slave geological province, a region of high mineral resource development in northern Canada where seasonal and all-weather roads and other infrastructure are important.With the effects of the last glacial ice sheet and the glacial lake McConnell, most of the surficial geology consists of fine-grained glaciolacustrine sediments and wave-washed bedrock (Wolfe et al., 2011).The study area includes the Great Slave lowlands and uplands with vegetation ranging from dense forest to arctic tundra (Ecosystem Classification Group, 2007).The Great Slave lowlands are poorly drained low relief terrain characterized by numerous water bodies separated by fens, white birch (Betula papyrifera) and black spruce (Picea mariana) forests, whereas the uplands are areas of higher relief and bedrock outcrops.Permafrost in the region is discontinuous and highly variable, usually with abrupt transitions from frozen to unfrozen ground over short distances (Wolfe et al., 2011).Therefore it is important to provide spatially detailed information about permafrost conditions and possible changes with climate warming.We selected this area also because it covers a large spatial environmental and vegetation gradient and many field observations are available.
In this area, the 30-year  mean annual air temperature ranged from −6.9 to −4.9 • C (Fig. 3f), and 30-year mean annual precipitation ranged from 288 to 309 mm.We conducted fieldwork in this area in 2005 and 2011.Data were collected on land cover types, vegetation conditions, surface OLT, and soil profile conditions (depths of the horizons, texture, and organic matter content), drainage, and summer thaw depths.We also collected observations made by other investigators in this area (Gaanderse, 2011;Brown, 1973;Karunaratne et al., 2008;Roujanski et al., 2012).In total we compiled 124 observation sites (Fig. 3b).Most sites were located near roads in the south providing easy access.Wolfe et al. (2014) also identified 1777 lithalsas in the Great Slave lowlands and uplands using 1 : 60 000 scale aerial photographs acquired between 1978 and 1980, of which 955 sites were within our study area (Fig. 3b).Lithalsas are permafrost mounds formed by ice segregation within mineral soil.These field data and lithalsa observations were used to validate the modelled permafrost distribution.

Land cover types and the probabilities of ground types
The land cover types (Fig. 3c) were classified using SPOT high-resolution visible and infrared images.The area is covered by a SPOT 4 image acquired on 7 August 2008 and a SPOT 5 image on 4 August 2011, both at 20 m spatial resolution.Radiometric normalization was performed using robust Theil-Sen regression on the overlap between the two scenes by adjusting the 2008 radiometry to match the 2011 image.The enhancement classification method (Beaubien et al., 1999) was used to improve image contrast  (Cihlar et al., 1998) was used to manually merge and label the clusters to 16 land cover types considering spectral similarity and spatial proximity to depict land cover features visible in high resolution images in Google Earth TM and our field data.Validation of a national land cover product (Olthof et al., 2014) that includes this study area suggests an overall accuracy of about 70 %.We modelled permafrost conditions for 13 land cover types, excluding water, roads, and recently disturbed areas (Table 1).We estimated the probability distribution of the ground types in each land cover type using field observations of OLT.Most of the land cover types had five to eight ground types (Table 2).Rock outcrop was assumed as one ground type without organic layer.The parameters were determined by fitting the observed relative cumulative frequency with the cumulative distribution probability of the modified logistic function (Table 2).Our field observations in the study area did not include land cover types of erect shrubs, shrub-herb mixture, and herbs (types S, SH, and H, respectively).To fill this data gap we used field observations of these land cover types that were conducted in 2013 near Tuktoyaktuk (69.44 • N, 133.03 • W) and Fort McPherson (67.44 • N, 134.88 • W).Due to the differences in geological setting, these observations may overestimate the OLT in our study area, especially when bedrock is near the surface.

Leaf area indices and fire disturbances
Leaf area indices in midsummer (Fig. 3d) were mapped using Landsat-5 imagery (acquired on 21 July 2002) based on Abuelgasim and Leblanc (2011), who developed an equation of LAI using field observations and the reduced simple ratio index of Landsat channels 3, 4, and 5.As they normalized their Landsat scenes based on SPOT VGT composite in the summer of 2005, we first normalized our Landsat images to the same SPOT VGT composite and then calculated LAI using their LAI equation.
A digital database of historically burned areas and year of burning (Fig. 3e) was obtained from the Government of Northwest Territories Centre for Geomatics  et al., 2011).For simplicity, we assumed that fires consumed all the foliage, and the remaining standing woody parts had little effect on radiation and surface energy exchanges.The albedo of the land surface during snow-free periods was reduced by 50 % immediately after the fire (Yoshikawa et al., 2003;Mack et al., 2011).We assumed that fire could consume a maximum of 13 cm of mosses and top organic matter excepting wetland and fen (type W), for which we used half of this depth due to their wet conditions (Johnestone et al., 2010;Turetsky et al., 2011).
Although trees can take several decades to re-establish, the regeneration and re-establishment of sedges and shrubs occur rapidly (Bond-Lamberty et al., 2002).Based on the LAI map for this area (representing the conditions in 2002), the average LAI in the area burned after 1994 (mostly in 1998) was about 60 % of the average LAI in the unburned area, and the average LAI in the area burned before 1973 was about 120 % of the unburned area.Based on these results and other studies (Bond-Lamberty et al., 2002), we assumed that land cover types dominated by sedges and shrubs would regenerate in the year following a fire and LAI would reach pre-burn levels in 5 and 15 years, respectively.Treedominated land cover types would regenerate in 5 years after a fire and LAI would reach pre-burn level in 50 years (Bond-Lamberty et al., 2002).We also assumed that the surface organic matter and mosses consumed by fire would recover in 50 years (Mack et al., 2011).The albedo of the land surface would also return to the pre-burn level.The increases in LAI, albedo, and surface organic layer were treated as linear patterns with time after regeneration.After recovery, the land cover types will be the same as pre-burn conditions, and summer LAI will not change with time (LAI varies seasonally).We did not consider new fires after 2013.2011) developed a spatial data set of 30year monthly averages  for air temperature and precipitation in North America by interpolating station observations.The data set has a spatial resolution of 5 arcmin (about 10 km).We interpolated these monthly data to 20 m resolution using bi-linear interpolation (Fig. 3f) and calculated annual total degree-days when daily air temperature was above 0 • C. Based on annual total degree-days and annual total precipitation, the climate was clustered into 20 classes.The average relative errors were 0.25 and 0.45 % for the annual total degree-days and annual total precipitation, respectively.Temporal climate changes were estimated using daily observations at the Yellowknife airport climate station as a template.A detailed description of the method can be found at Zhang et al. (2013).From the 1940s to the recent past decade (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), annual mean air temperature increased by 2.1 • C. The increase mainly occurred after the mid-1960s.Annual total precipitation increased by 99 mm from the 1940s to the past decade.

McKenney et al. (
The climate scenario data were downloaded from the World Data Center for Climate (http://mud.dkrz.de/wdc-for-climate/, accessed April 2011).First, we selected 18 climate projections of six climate models (CCCma, ECHAM, HadCM, GFDL, MIROC, and NRCAR) under three greenhouse gas emissions scenarios (B1, A1B, and A2); we then selected three climate change scenarios to represent low, medium, and high warming scenarios based on their temperature projections.They were generated by CC-Cma (B1), CCCma (A1B), and MIROC (A1B) (Fig. 4).We 4 F ig is the rate of ground inflow (day −1 ). 5 F os and 6 F og are the rates of surface and ground outflows, respectively (day −1 ).Detailed description of these parameters can be found in Zhang et al. (2002Zhang et al. ( , 2012)).-is assuming no lateral flows.LAI is for leaf area indices in peak growing season.
did not select CCCma (A2) because the projected air temperature is lower than several other projections and it is similar and sometimes lower than that of CCCma (A1B) before the 2050s.Under the selected three scenarios, annual mean air temperature is projected to increase 0.4, 1.7 and 2.7 • C, respectively, from the recent past decade to the 2050s, and to increase 1.0, 3.4 and 4.9 • C, respectively, from the recent past decade to the 2090s.The projected changes in precipitation are not very significant (Fig. 4b).

Other input parameters
In addition to OLT, the model also requires inputs about mineral soil conditions and hydrological parameters.Field surveys and coring indicate differences in sediment types above and below an elevation of about 205 m a.s.l.(above sea level), likely related to sedimentation within the glacial lake Mc-Connell (Stevens et al., 2012).Therefore, we first divided the study area into two regions using the 205 m a.s.l.isoline (Fig. 3b) and then defined the mineral soil conditions for each region.This delineation also closely corresponds to the separation of the Great Slave lowlands and uplands (Ecosystem Classification Group, 2007).At elevations below 205 m, mineral soils are mainly composed of clay (Wolfe et al., 2011).Above 205 m., mineral soils were assigned as sandy loam for the low vegetation or barren land (type L) and low-to medium-density coniferous forest (type C3), and as silty clay for other land cover types based on the general patterns of field observations (Stevens et al., 2012;Wolfe et al., 2011).The organic matter content in mineral soils was estimated based on the general patterns observed by Hossain et al. (2007).Except for low vegetation and bedrock (Types L and R), the depth of all surficial deposits (including peat and mineral soils) was assumed to be 7 m for bogs and wetlands (types C5 and W) and 5 m for other land types based on most of the borehole observations (Wolfe et al., 2011).The soil thickness for the low vegetation type was assumed to be 2.5 m and we assumed no soil on bedrock outcrops.The thermal conductivity of bedrock was 0.026 W m −1 • C −1 (Brown, 1973).The depth of the ground profile simulated was 120 m, and heat flux at the bottom of the ground profile was assumed as 0.007 W m −2 (Majorowicz and Grasby, 2010).
The model also requires several parameters to estimate lateral water flows (Zhang et al., 2002(Zhang et al., , 2012)).We defined these parameters based on the general drainage conditions and the range of water table for different land cover types (Table 3).The model also has a parameter to consider the effects of snow-drifting due to wind (Zhang et al., 2012).A positive value denotes the fraction of snowfall drifted away from the site and a negative value indicates drifting snow captured from the surroundings.We estimated this parameter based on the vegetation types and their density (represented by the LAI in summer) (Table 3).Since this region is relatively flat, we did not consider the effects of topography on solar radiation.

Spatial modelling
To reduce the computation time, LAI was divided into 18 levels based on the effects of LAI on solar radiation (i.e., exp(−0.5LAI)).The latitude in the area was divided into five grads (0.2 • difference between two adjacent grads).Together with the land cover types (13 types), elevation levels (2 types), climate (20 classes), and wildfire years (13 types), the land pixels in the study area had 31 724 unique www.the-cryosphere.net/8/2177/2014/The Cryosphere, 8, 2177-2194, 2014 combinations of the input types, which was 0.2 % of the land pixels in this area (we ran the model only for these unique input types rather than running the model pixel by pixel).For each land cover type, the model was run for all the ground types under the three climate scenarios.Model computations took about three weeks using three personal computers.
3 Result and analysis

Permafrost distributions
Figure 5 shows the distributions of permafrost probability modelled for the 1950s (1950-1959), the 2000s (2000-2009), and the 2050s (2050-2059).These results show significant reductions in permafrost probability due to climate warming from the 1950s to the 2050s.Permafrost is predicted to disappear completely in most of the area by the end of the 21st century.The spatial distribution of permafrost in this area was mainly related to land cover types and associated ground conditions.For example, the non-permafrost area in the 2000s mainly occurred in rock outcrops, low vegetation, and deciduous forest (land cover types R, L, and D, respectively), while sparse coniferous forest (types C4 and C5) and herb dominated areas and wetlands (types H and W) had high permafrost probability.Leaf area indices and fire disturbances also affected permafrost conditions, but their effects were mixed with that of land cover types.Permafrost probability in the northeast corner of the study area was higher, probably due to the relatively cooler climate in this area.However, the overall effects of the climate gradient on spatial distribution of permafrost were not so obvious because of the strong effects of ground conditions and relatively small climate gradient in this area (Fig. 3f).For example, although the southern part is slightly warmer than the northern part, LAI is higher and soils contain more clay in the south, which compensated for the effects of warmer temperature on permafrost development.
Figure 6a shows the modelled temporal change of average permafrost probability (or permafrost extent) in this area from 1942 to 2100.On average, permafrost extent was reduced from 72.0 % in the 1950s to 52.0 % in the 2000s.The model predicted that permafrost extent would be reduced to 12.4 % in the 2050s and to 2.5 % in the 2090s.The difference in permafrost extent between the low and high warming scenarios was relatively small, usually less than 5 % (in permafrost extent).
Different land cover types show different decreasing patterns in permafrost extent (Fig. 6b).The most rapid reduction in permafrost extent was in low vegetation or barren land (type L) due to their exposed land conditions.The reduction was rapid in deciduous forest area (type D) due to its shallow OLT.Spruce lichen bogs and wetlands (types C5 and W) were the last land cover types experiencing permafrost reduction.With time, the reduction of permafrost for a land cover type generally shows a slow-rapid-slow decline pattern.The initial slow reduction is due to the lower warming rate of climate (especially before the mid-1980s) and the time required to warm the ground to a critical level.The later slow reduction of permafrost is related to ground types with deep organic layers and high LAI, which can protect permafrost from thawing (Shur and Jorgenson, 2007).Model results show that some forest areas could still maintain permafrost by the end of the 21st century mainly due to their deep OLT and high LAI, assuming they were not disturbed by wildfires during the 21st century.

The most likely permafrost conditions
Permafrost conditions under varying ground types differ significantly due to the effects of OLT.Permafrost does not typically exist where OLT is very shallow.Therefore it is difficult to calculate the average permafrost conditions.A meaningful way is to present the most likely permafrost conditions, which is determined from the model results of the most likely ground types under the medium warming climate scenario.The most likely permafrost extent averaged in this area (dashed curve in Fig. 6a) was slightly different from the average of all the ground types under the three climate scenarios (solid curve in Fig. 6a).The average and the most likely spatial patterns were similar in the 1950s and 2000s (the most likely permafrost distribution only has two types, with and without permafrost, shown in Fig. 7 as coloured and white, respectively).However, the average and the most likely permafrost distributions were very different in the 2050s (comparing Figs.5c and 7c) as permafrost was predicted to disappear in most areas under the most probable ground and climate conditions.
Figure 7 shows the distribution of the most likely activelayer thicknesses in the 1950s, 2000s, and 2050s.The active layer deepened from the 1950s to the 2000s, and permafrost in most areas was predicted to disappear in the 2050s.The red patches (deeper active layer) in the northeastern corner in Fig. 7b are due to the effects of wildfires that consumed vegetation and some of the top organic matters.For the pixels with permafrost in all the years before 2012, the average active-layer thickness increased from 0.7 m in the 1940s to 1.4 m in the 2000s.

Result verification
According to the Canadian permafrost map (Heginbottom et al., 1995), the study area is within the extensive discontinuous permafrost zone.Based on the 124 observation sites in this area, permafrost was detected at 87 sites.Frozen ground was not detected at the other 37 sites.The permafrost sites account for 70.2 % of the observation sites, which confirm that this area is within the extensive discontinuous permafrost zone (50-90 %).The modelled average permafrost extent in this area was 64.2 % during 1942-2012, well within the range of the extensive discontinuous zone.Based on the ground types most closely matching local conditions, the model results show that permafrost occurred at 71 of 87 observed permafrost sites.For the 37 sites where no frozen ground was detected in the field, the model results show that 21 sites did not have permafrost, six sites had permafrost but their thaw depths were deeper than the measurement depths, and the model results were incorrect for the other ten sites.Overall, the model correctly simulated permafrost occurrence for 82 % of the observed permafrost sites and correctly simulated 73 % of the sites where frozen ground was not detected in the field, either because thaw depths were deep or permafrost did not exist.Coincidently, the modelled and observed total numbers of sites with permafrost (87 sites) were the same.Figure 8a shows a comparison between the modelled and observed percentage of sites with permafrost at each land cover type.The model results were close to the observations for most of the land cover types except for the low vegetation sites (L).Permafrost was observed at three of the four low vegetation sites although the model did not indicate permafrost at these sites.Field data show that two of the three sites are in disturbed areas with some peat layers added in the soil profile, and another site is near the road with some embankment material on the surface.These modified conditions may have actually promoted permafrost development.Wolfe et al. (2014) identified 955 lithalsas in our study area (Fig. 3b).By definition, permafrost exists at these locations.We calculated the average permafrost probability during 1978-1980, at which time the aerial photographs were acquired to identify the lithalsas.Eighteen locations were classified as water since they are very close to water bodies.The model results show that all the 937 non-water locations had a ≥ 14 % probability of permafrost occurrence.The probability of permafrost occurrence was ≥ 50 % at 838 locations, and was ≥ 90 % at 501 locations.The average probability was 76.9 % for all the locations.The lithalsa mounds are usually well drained due to their local topography, and contain segregation ice.The model did not consider these factors and thus they might contribute to the underestimation of the probability of permafrost occurrence for these locations.
For the observed sites, the modelled mean summer thaw depth for each land cover type was comparable with the observations, although the variation range was usually large among the sites (Fig. 8b).The modelled average thaw depth of all the sites with permafrost was 0.78 m, which was very close to the average of the observations (0.81 m). Figure 9 shows a scatter graph comparison for the 71 sites with permafrost.The magnitudes of the modelled summer thaw depths were similar to the observations for most of the sites, although significant differences existed for some sites due to their specific soil, vegetation, and hydrological conditions.The correlation coefficient is 0.365 (n = 71).The land cover types S, SH, and H were not shown because there were no field observations for in the study area for these three types.
The modelled dynamics of snow depth, especially the timing of the accumulation and disappearance of snow and the maximum snow depth, compared well with the observations at the Yellowknife climate station.The correlation coefficient is 0.892 for the 19 794 observations from 1955 to 2012.The modelled permafrost thickness for recent decades ranged from several metres to 130 m, which was comparable to the observations reported by Smith and Burgess (2002) and Roujanski et al. (2012).We also tested the model offline using observations in this area measured by Brown (1973), Karunaratne et al. (2008), andRoujanski et al. (2012).Based on their local ground and vegetation conditions, the model could capture their ground temperature regimes and permafrost conditions.

Discussion
This study proposes a new approach for mapping permafrost and change with climate considering uncertainties in ground conditions and climate projections.We apply this approach at 20 m resolution to a large area in northern Canada.The data requirement and the cost of computation are affordable, and the mapped permafrost conditions are in reasonable agreement with field observations and other studies.Compared to previous mapping methods, several features of this approach are noteworthy.
First, compared to the traditional zonal permafrost mapping methods (Nikiforoff, 1928;Heginbottom et al., 1995;Brown et al., 1998), this new approach is more objective, replicable, and quantitative in integrating most of the related factors, and can be used to map permafrost at higher spatial resolution and to assess the impact of climate change.Unlike the equilibrium models (as reviewed by Riseborough et al., 2008), our process-based modelling approach can quantify the transient responses of permafrost to changes in climate and ecosystems.
Second, this approach integrates satellite remote sensing data, field observations, and a process-based model.Satellite data can greatly improve the spatial resolution of existing permafrost maps and can provide various parameters, including land cover, LAI, wildfire, and topographic features.As current soil and surficial geological maps are coarse, the heterogeneous ground conditions are statistically represented based on field observations.Third, the permafrost maps produced by this approach are more useful for land use planners and decision-makers due to their higher spatial resolution and effective incorporation of incomplete and uncertain information about ground conditions and future climate projections.Previous studies usually presented different permafrost maps under assumed different ground conditions and climate change scenarios (e.g., Zhang, 2013;Zhang et al., 2008Zhang et al., , 2013;;Daanen et al., 2011).Such results are difficult to use due to their wide differences.The probability concept is similar to the traditional permafrost maps, but this new approach can develop permafrost maps with a much higher degree of precision and spatial resolution than the traditional ones.Figure 10 shows the importance of spatial resolution for permafrost zones.At 20 km resolution, our results show the same permafrost zone as in the traditional permafrost map (extensive discontinuous).At finer resolutions, however, other permafrost zones appeared.In addition, this approach also provides most likely permafrost conditions, such as active-layer thickness and permafrost thickness.
Finally, this approach may serve as a practical way to map permafrost evolution at high spatial resolution for other regions.Satellite remote sensing data are routinely available across the globe.To consider the heterogeneity of ground conditions, we estimated their statistical distributions rather than mapping ground conditions explicitly, which is difficult at present.Organic layer thickness and other general ground conditions can be surveyed in the field and the data can be accumulated gradually as they do not change significantly over time without disturbances.In this study we estimated the probability of ground types based on OLT.With more field data available, other ground features can be included to better define ground types.The spatial input data can be organized and used to run the model on multiple computers simultaneously, thus greatly reducing computation time for high resolution spatial modelling.In addition, this approach allows for modifications in spatial developments in model and remote sensing technologies.
Permafrost is a ground thermal condition influenced by many factors of geological, climatic, hydrological, and ecological processes at different temporal and spatial scales (Jorgenson et al., 2010;Shur and Jorgenson, 2007;Bonnaventure and Lamoureux, 2013).Although our new mapping approach quantitatively integrated most of the factors based on energy and water dynamics, it had several limitations.First, the NEST model is one-dimensional and assumes each grid cell to be uniform without lateral heat exchange.Therefore, the results could not represent areas with strong lateral heat fluxes, such as areas between very different land cover types, especially close to water bodies.The model did not consider permafrost thaw due to lateral heat flow, which could be important for patchy permafrost.Second, we considered fire disturbances that occurred in the past, but we did not consider future fires.Such a treatment might overestimate future permafrost extent, because high LAI and deep organic layer can effectively keep permafrost in disequilibrium from the warming of the atmospheric climate (Camill and Clark, 1998;Shur and Jorgenson, 2007).Ground subsidence and associated hydrological changes could have significant impacts on permafrost degradation as well.Finally, this study did not consider changes in land cover types and LAI caused by climate warming.Observations show that northern high latitudes are becoming greener and shrubbier with climate warming (e.g., Tape et al., 2006).These changes can be very important for the evolution of permafrost.

Conclusions
Spatially detailed information about permafrost evolution with climate warming is important for land use planning, infrastructure development, and environmental and ecological assessments.Although satellite remote sensing can provide high resolution maps of vegetation and land cover types, data regarding soil and ground conditions are much coarser, which greatly limits our capacity to map permafrost at higher resolutions.The wide range of projected climate scenarios also brings along uncertainties when mapping permafrost change.Our proposed new approach incorporates these uncertainties for permafrost mapping by integrating remote sensing, field observations, and a process-based model.Land cover types from satellite remote sensing are used to capture the general land features and to improve the spatial resolution of existing permafrost maps.Other vegetation features can also be derived from remote sensing imagery.The probabilities of different ground conditions within land cover types are estimated from field observations.A process-based model is used to quantify the dynamics of permafrost for each ground type under three representative climate scenarios of low, medium, and high warming.From the model results, the probability of permafrost occurrence and the most likely permafrost conditions are determined.We apply this approach at 20 m resolution to a large area in Canada.The mapped permafrost conditions are in agreement with field observations and other studies.This demonstrates that the data requirements, model robustness, and computation time are reasonable, and that this approach may serve as a practical means to map permafrost evolution at a high resolution in other regions.

Figure 2 .
Figure 2.An example of determining the probability distribution of organic layer thickness for a land cover type (C2, coniferous forest -medium density).(a) Fitting the observed relative cumulative frequency (the dash curve) with a modified logistic function (the bold curve).The inset shows the scatter graph comparison and linear regression for the fitting.(b) The probability density (the bold curve), and the distribution of the ground types and their probabilities (the rectangle boxes).

Figure 3 .
Figure 3. (a)The location of the study area in Canadian permafrost zones(Heginbottom et al., 1995), (b) the observation sites and the two regions divided based on elevation 205 m a.s.l., and the distributions of (c) land cover types, (d) leaf area indices, (e) wildfires, and (f)30-year  (1971-2000)  mean annual air temperature in the study area.

Figure 4 .
Figure 4.The historical climate trends and the selected three scenarios for the 21st century for the (a) annual mean air temperature and (b) annual total precipitation.The historical climate is observed at the Yellowknife airport climate station.The low, medium, and high warming scenarios were generated by CCCma (B1), CC-Cma (A1B), and MIROC (A1B), respectively.More details about the data and the climate models can be found at the World Data Center for Climate (http://mud.dkrz.de/wdc-for-climate/).

Figure 6 .
Figure 6.The modelled temporal changes of permafrost probability (a) for the entire study area and (b) for five selected land cover types in this area.

Figure 7 .
Figure 7. Modelled most likely active-layer thickness in the 1950s, 2000s and 2050s (a, b, and c, respectively).(d) shows an enlarged area in (b).This figure also shows the most likely distribution of permafrost (white for non-permafrost areas and the other colours for areas with permafrost).

Figure 8 .
Figure 8. Comparisons of observed and modelled (a) percentage of the sites with permafrost and (b) the average and standard deviation of summer thaw depths for each land cover type.The numbers in brackets below the x axis are the numbers of observation sites.The average or the standard deviation of summer thaw depths were calculated based on the sites with permafrost (observed or modelled).The land cover types S, SH, and H were not shown because there were no field observations for in the study area for these three types.

Figure 9 .
Figure 9.Comparison between the modelled and measured summer thaw depths.There are 71 observation sites where both observation and the model show the existence of permafrost.The correlation coefficient is 0.365 (n = 71).

Figure 10 .
Figure 10.The effects of spatial resolution (shown in the panels) on permafrost zonal distributions.These maps were produced by dividing permafrost probability in the 2000s (Fig. 5b) into five zones and then re-sampled to different spatial resolutions.

Table 1 .
Land cover types classified using SPOT.stretching the dynamic range representing land features to the full 8 bit range while compressing the dynamic range of water, clouds, and shadows.A large number of spectral clusters were generated from the enhanced multispectral data using fuzzy K means clustering and a pseudo-colour table was applied representing the enhanced colours.The progressive generalization classification approach by

Table 2 .
Parameter values of the probability distribution equation for organic layer thickness and the number of ground types for each land cover type estimated using field observations.
* R 2 is the square of the correlation coefficient between the fitted and observed relative cumulative frequency.

www.the-cryosphere.net/8/2177/2014/ The Cryosphere, 8, 2177-2194, 2014 (
http://www.geomatics.gov.nt.ca/).In the study area, about 26.2 % of the land area experienced wildfire since 1967.The major fires were in 1973 and 1998, which accounted for 71.6 and 19.3 % of the total burned area, respectively.Fires consume vegetation and soil organic matter.The amount of consumption and the post-fire regeneration process depend on fire intensity, vegetation condition, soil moisture, and local topography (Turetsky

Table 3 .
Hydrological parameters defined for different land cover types.WT ig is the lowest water table (cm below the surface) for beginning lateral ground inflow. 2 WT os and 3 WT og are the highest water table (cm below the surface) for beginning lateral surface and ground outflows, respectively. 1