using a high spatial resolution dataset

Abstract. Climate projections for the 21st century indicate that there could be a pronounced warming and permafrost degradation in the Arctic and sub-Arctic regions. Climate warming is likely to cause permafrost thawing with subsequent effects on surface albedo, hydrology, soil organic matter storage and greenhouse gas emissions. To assess possible changes in the permafrost thermal state and active layer thickness, we implemented the GIPL2-MPI transient numerical model for the entire Alaska permafrost domain. The model input parameters are spatial datasets of mean monthly air temperature and precipitation, prescribed thermal properties of the multilayered soil column, and water content that are specific for each soil class and geographical location. As a climate forcing, we used the composite of five IPCC Global Circulation Models that has been downscaled to 2 by 2 km spatial resolution by Scenarios Network for Alaska Planning (SNAP) group. In this paper, we present the modeling results based on input of a five-model composite with A1B carbon emission scenario. The model has been calibrated according to the annual borehole temperature measurements for the State of Alaska. We also performed more detailed calibration for fifteen shallow borehole stations where high quality data are available on daily basis. To validate the model performance, we compared simulated active layer thicknesses with observed data from Circumpolar Active Layer Monitoring (CALM) stations. The calibrated model was used to address possible ground temperature changes for the 21st century. The model simulation results show widespread permafrost degradation in Alaska could begin between 2040–2099 within the vast area southward from the Brooks Range, except for the high altitude regions of the Alaska Range and Wrangell Mountains.


Introduction
According to State of the Climate in 2010 Report (Richter-Menge and Jeffries, 2011), the arctic cryosphere is undergoing substantial changes such as loss of sea ice and warmer ocean temperatures, melting of the Greenland Ice Sheet and glaciers as well as continuous increase of permafrost temperatures.Permafrost1 , one of the main components of the cryosphere in northern regions, influences hydrological processes, energy exchanges, natural hazards and carbon budgets.The World Meteorological Organization (WMO) identifies permafrost as one of the six cryospheric indicators of global climate change (Brown et al., 2008).Changes in the thermal state of permafrost in Alaska were reported recently by Romanovsky et al. (2010) and Smith et al. (2010), where they observed an increase in permafrost temperatures by 0.5-3 • C over the last 30 yr.Thawing of permafrost causes land surface changes, damaging forests, houses, and infrastructure.Changes in permafrost thermal state can have a significant impact on the state's economy due to the additional repair costs of public infrastructure (Larsen et al., 2008).
Mapping of permafrost distribution, especially its thermal state remains a challenging problem due to the sparsity of observed data.Despite the fact that geophysical surveys and boreholes are the most reliable sources of information E. E. Jafarov et al.: Numerical modeling of permafrost dynamics in Alaska about permafrost, they are extremely costly and are mostly available from relatively small areas such as oil fields and transportation corridors.
In the past, one of the most popular methods in permafrost mapping was the use of a frost number, or permafrost index methods (Nelson, 1986), which usually takes into account seasonal air temperature variations explicitly, and do not include other important factors (snow, soil moisture content, soil thermo-physical parameters) that affect the permafrost thermal state.Riseborough et al. (2008) divided permafrost models into three main categories: empirical, equilibrium and numerical.Empirical models relate permafrost occurrences to topoclimatic factors and use empirically derived landscape parameters representing the response of the active layer and permafrost to both climatic forcing and local factors, such as soil properties, moisture conditions and vegetation (Nelson et al., 1997;Shiklomanov and Nelson, 2002;Zhang et al., 2005;Shiklomanov et al., 2008).Equilibrium models employ transfer functions between the air and ground temperatures to define the active layer depth.The Kudryavtsev Model, N factor, TTOP and GIPL 1.0 models are classified as equilibrium models (Kudryavtsev et al., 1974;Romanovsky and Osterkamp, 1995;Romanovsky and Osterkamp, 1997;Shiklomanov and Nelson, 1999;Klene et al., 2003;Sazonova and Romanovsky, 2003;Wright et al., 2003).However, the applicability of equilibrium models are restricted to problems of limited complexity, and where the transient effects may be neglected or unimportant.Numerical solution methods are generally used to solve freezing and thawing problems over a short (engineering) time scale, where the transient effects of phase change are important (see Williams and Smith, 1989, p. 86).Therefore, transient numerical modeling with incorporated phase change effect is the most effective method to simulate and forecast the thermal regime of permafrost over a relatively short time interval.The transient numerical model (e.g., Gruber et al., 2004a,b;Nicolsky et al., 2007;Farbrot et al., 2007;Etzelmuller et al., 2008) has a higher computational cost, however, they are more accurate in determination of the temperature field and phase change boundary dynamics.
In this article we propose a method based on numerical modeling, which allows us to map the temporal dynamics and spatial distribution of permafrost with high spatial resolution.To simulate ground thermal regimes, we implement the GIPL2-MPI numerical transient model.The GIPL2 was developed by G. Tipenko and V. Romanovsky (Tipenko et al., 2004), and applied for first time to the entire Alaskian permafrost domain with 0.5 • spatial resolution by Marchenko et al. (2008).In order to quantify the socio-economical impact of permafrost degradation, the permafrost distribution maps with higher spatial resolution are required.In the current project, we map the permafrost thermal state with 2 by 2 km spatial resolution, which is 410 205 grid points, to cover the entire Alaskan region.Due to increase in the amount of spatial grid points, the computational load increases as well.To perform calculations efficiently, we developed the GIPL2-MPI version of the GIPL2 model, which runs in parallel on several processors.To run the parallel model, we used the Arctic Region Supercomputing Center facility at the University of Alaska Fairbanks.
The current paper is constructed in the following order.The mathematical description section gives more detailed information on the used method.In the methods section we outline all necessary input datasets, such as initial temperatures, snow, thermal properties of multi-layered organic and mineral soils, and geothermal heat flux.The model calibration and validation section provides a detailed description of the model calibration with observed monthly averaged temperatures available from shallow boreholes stations, and validation of the mean annual ground temperatures from deep boreholes, as well as validation of the active layer thickness values with CALM2 active layer observation stations in Alaska.The Optimization of Ground Thermal Parameters section illustrates the effect of the additional upper layer(s) of organic matter on the mean annual ground temperatures.In the results section, we illustrate the model outputs and provide the per decadal warming rates at different depths during the 21st century.The discussion section gives an overview on major physical factors that affect the permafrost thermal regime and outlines how the model can be improved.Finally, the conclusion section outlines major results and concludes the current work.

Mathematical model
The GIPL2-MPI numerical model solves the Stefan problem (Vasilios and Solomon, 1993) with phase change, which is the problem of thawing or freezing via conduction of heat.The enthalpy formulation used in the solution of Stefan problem is the most common method, which does not require explicit treatment of the moving freeze/thaw boundary (Caldwell and Chan, 2001).The core of the GIPL2-MPI numerical model is based on the 1-D quasi-linear heat conductive equation (Sergueev et al., 2003): where x ∈ (x u , x l ) is a spatial variable which changes with depth, x u and x l are the upper and lower boundaries of the vertical grid correspondingly, and τ ∈ (0, T ) is a temporal variable.t (x, τ ) is temperature, k(x, t) is thermal conductivity (Wm −1 K −1 ), H (x, t) is an enthalpy function.
where C(x, s) is volumetric heat capacity (MJm −3 K −1 ), and (x, t) is volumetric unfrozen water content (fraction of 1); L is the volumetric latent heat of freeze/thaw (MJm −3 ).The Eq. ( 1) is complemented with boundary and initial conditions.The upper part of the domain corresponds to the air layer which is at two meters height above the surface.The fictitious domain formulation (Marchuk et al., 1986) allows to incorporate a seasonal snow layer into the current air layer.The Dirichlet type boundary condition was used as an upper boundary condition: where t air is a monthly averaged air temperature.The geothermal gradient was set at the lower boundary: where g is geothermal gradient, a small constant number (K m −1 ).For initial temperature distribution we used an appropriate ground temperature profile based on the point location t (x, 0) = t 0 (x). (5) The formula for unfrozen water content (x, t) is based on empirical experiments and has the following form: Parameters a and b are dimensionless positive constants (Lovell, 1957), η(x) represents the volumetric soil moisture content.The right side of the Eq. ( 6) represents the liquid pore water fraction.The constant t * = (1/a) b is a freezing point depression, which from a physical point of view means that ice does not exist in the soil if t > t * .(x, t) changes with depth and depends on the soil type.The discritized form of the Eq. ( 1) can be found in Sergueev et al. (2003) and Marchenko et al. (2008).A detailed mathematical description of the model and numerical solution methods can be found in Nicolsky et al. (2007).

Methods
We implemented the Scenarios Network for Alaska Planning (SNAP) dataset as a baseline input for the GIPL2-MPI numerical model.The dataset is composed of five GCMs, which according to SNAP, performs the best for Alaska (Walsh et al., 2008).It includes monthly averaged temperatures and precipitation data for the years 1980-2099, using A1B carbon emission scenario.The output from the selected five models were downscaled to 2 by 2 km resolution by SNAP using the knowledge-based system PRISM 3 .PRISM uses 3 Parameter-elevation Regressions on Independent Slopes Model climate mapping system (http://www.prism.oregonstate.edu).
a digital elevation model that contains information describing Alaska's topography (slopes, aspects, elevations) and observed precipitation measurements to determine variations in precipitation as functions of elevations.
To calculate snow depth and its thermal conductivity we developed the following method outlined below.Snow density is calculated according to (Verseghy, 1991) where τ f = 0.24 corresponds to an e-folding time of about 4 days, ρ s is the snow density in unit of kg m −3 , and the value of ρ s ranges from minimum snow density ρ smin to maximum snow density ρ smax , with a time step τ of one month.We chose ρ smin and ρ smax from the corresponding snow class according to Sturm et al. (1995).Snow depth h s is calculated by extracting snow water equivalent (SWE) from the downscaled five GCM composite precipitation dataset: Snow thermal conductivity k s was calculated according to Sturm et al. (1997): Snow water equivalent was extracted from the precipitation dataset by comparing monthly mean temperatures (MMT) with the water freezing point temperature.If the MMT is less than 0 • C then we accumulate the SWE on monthly basis (i.e. if snow stays on the ground after the first snow fall we add existing SWE to the SWE for the current month).When the SWE is obtained, we calculate density, depth and thermal conductivity of the snow by employing Eqs. ( 7)-( 9).We analyzed the ground temperature profiles (ground temperature distribution profiles are available online at Geophysical Institute Permafrost Laboratory4 and Advanced Cooperative Arctic Data and Information Service5 websites) in more than 25 relatively deep boreholes from 29 to 89 m in depth (Osterkamp andRomanovsky, 1999, Osterkamp, 2003) along the Trans-Alaskan transect.This analysis revealed a ground temperature zonality in Alaska with generally lower permafrost temperatures in the north and higher ground temperatures in the south.Based on this zonality, we extrapolated available initial ground temperature profiles to the wider areas and classified them into 18 ground temperature zones (Fig. 1).The 18 ground temperature zones represent the 18 classes of temperature distribution with depth and were used as initial conditions for simulation.
The thermo-physical properties (volumetric soil ice/water content, unfrozen water curve parameters, soil heat capacity and thermal conductivity, thickness of soil layers, etc.) for 18 ground temperature zones might be different and depend on many factors including surficial geology.The number of soil type classes we used in these simulations was 26 and each class had its own number of soil and bedrock layers with different thermal properties (e.g.peat, silt, bedrock, gravel etc.).The multilayered soil columns were assigned for each soil class according to the Modified Surficial Geology Map of Alaska (Karlstrom, 1964).The thermo-physical properties were assigned to each ground mineral layer according to surficial geological (soil type) class.The model was calibrated against the ground temperature measurements from the shallow boreholes, which were specific for each soil class and geographical location (the method used and its limitations were described in more detail by Nicolsky et al., 2007).Organic layer in the model was introduced as a separate layer(s) which could be added at the top of mineral soil column.For upper organic soil layers, we used the data obtained from the numerous field observations and Ecosystem Map of Alaska from the National Atlas of the United States of America6 .To further optimize the number and the thermal properties of the organic layers we developed an algorithm described in the Optimization of ground thermal parameters section.
The active layer thickness (ALT) in the model is calculated according to the unfrozen water function Eq. ( 6) and unfrozen water saturation coefficient (UWSC).The UWSC is equal to 95 % and serves as a threshold between frozen and thawed state (i.e.freeze/thaw moving front corresponds to state when the degree of soil layer saturation passes UWSC).The ALT correspond to the deepest upper freeze/thaw front over the hydrological year cycle.
Each grid point on the map uses a one-dimensional multilayer soil profile down to the depth of 700 m.The verti-cal grid has fine resolution between nearby points at the near surface ground layer (0.01 m) and becomes coarser towards the bottom boundary (100 m).The geothermal heat flux was assigned as a lower boundary condition.The values for the geothermal heat flux were generated using Pollack's geothermal heat model (Pollack et al., 1993).

Model calibration and validation
For the initial model calibration, we used measured data from more than fifteen shallow boreholes (1-1.2 m in depth) across Alaska.These high quality ground temperature measurements (precision generally at 0.01 • C) are available from the mid-1990s to 2010.Soil water content and snow depth measurements were also available at most of these boreholes.Figures 2 through 4 illustrate the results of the model calibration for the three shallow borehole stations.The West Dock site (Fig. 2) is located on the outer Arctic Coastal Plain within the Prudhoe Bay oil field.The polygonized "uplands" and drained thaw-lake basins constitute the primary relief at this site.Landcover units include wet nonacidic graminoidmoss tundra.The site was described by Osterkamp (1987) and Romanovsky and Osterkamp (1995).The SagMat site (Fig. 3) is located on a north facing slope of about 2 degrees.The vegetation cover is a moist acidic tundra (Walker et al., 2008).The Galbraith Lake site (Fig. 4) is located in the previously glaciated mountain valley.Landcover units include graminoid-moss tundra and graminoid, prostratedwarf-shrub and moss tundra (wet and moist nonacidic).Further site descriptions can be found in Ping et al. (2003).All these sites were instrumented by at least ten thermistors arranged vertically at depths from 0 m to 1 m.The detailed description on thermistor set up and installation can be found in Nicolsky et al. (2007).For all the compared stations, correlations between GCMs and the observations are higher than 90 % for monthly mean air temperatures (MMAT).Despite the high correlation between downscaled and observed MMATs, the freezing and thawing periods are not always well satisfied with the downscaled GCMs composite.As it can be seen from Fig. 2, the variances between simulated and measured ground temperatures during the 1999 freezing period increased with depth.The simulated monthly averaged ground temperatures decrease sharply when the actual ground temperatures freezing period is longer.The same pattern can be observed during the thawing periods 2004, 2005 in SagMat (Fig. 3).Furthermore, there were winter periods where the actual monthly averaged ground temperatures were colder than the simulated (see SagMat thawing period 2002, 2004, 2005, Fig. , Fig. 3, Galbraith Lake winter periods 2004, 2005, Fig. 4), which could be due to snow conductivities and snow depths biases over those winter periods.The simulated ground temperatures were, in general, smoother than measured temperatures and did not represent the seasonal effects, which could be an another from 29 m to 89 m in depth (GTN-P7 ) were available for the model calibration in terms of permafrost temperature profiles.For most of the borehole stations, the deviation between observed and measured mean annual ground temperatures (MAGT) was less than 1 • C (Fig. 5).There were four stations where the deviation was greater than 1 • C, three of them located in the tussock area.The MAGTs for the tussock like areas are usually sensitive to the snow depth.If snow depth does not exceed the height of tussock, then cold air could penetrate deep into the ground; otherwise, when snow depth exceeds the height of tussock it isolates the ground from the cold air.The permafrost observation station with the measured −3.81 • C and simulated −6.32 • C MAGTs at a depth of 20 m corresponds to the site in the continuous permafrost zone with MAATs around −10 • C. The fact that measured MAGT was almost 3 degrees warmer might be due to the site's location, which was around small lakes.The convective heat transfer due to ground water movements or heat from the open water reservoirs might be essential factor in producing warmer MAGT.
In addition to deep borehole stations, we validated the MAGTs from the permafrost observation stations from US Schools project8 (Fig. 6).The measurements from these stations were taken at relatively shallow depths, ranging from 1 to 6 m.Most of the stations were located in close proximities to public schools, rivers or lakes.There were three stations where differences between measured and simulated temperatures were greater than 3 • C. One of them was located in the city, another one close to the lake, and the other one was close to one of the branches of the Yukon river.The first permafrost station experienced the influence of anthropogenic warming and the other two could experience the laterial heat exchange due to subsurface ground water movements.Finally, we validated the simulated active layer thicknesses with observed ALTs from 43 CALM 9 observation stations in Alaska (Fig. 7).For the ALT comparison test we compared averaged active layer thicknesses over the available time pe-  riods with the corresponding simulated ALTs.During ALTs simulation the soil moisture content was specified for each of the simulated stations and held constant throughout the entire simulation period.The major restriction of this approach reflects the limitation of available data on soil moisture content and its dynamics over time for each of the compared stations.There were two main uncertainties while comparing simulated and measured ALTs.First coming from the active layer measuring method, where there are several methods on measuring ALT and all of them have their own limitations (Nelson and Hinkel, 2003).Second comes from the simulated AL depth, which was driven by monthly averaged climate data and by the amount of prescribed soil moisture content.
During model validation, the values for several observation stations were adjusted by assigning additional organic layers.
To evaluate the measure of the overall model performance and the model bias, we calculated the mean absolute error (MAE), root mean square error (RMSE) and mean bias error (MBE) according to the following series of equations (Willmott and Matsuura, 2005): where e i is a difference between simulated and observed MAGTs, and ALTs and n is the number of stations.
The MAE shows an overall error for all compared stations, when the RMSE emphasizes an error variation within the individual stations and the MBE shows that the model underestimates or overestimates the observed data.The MBE in Table 1 shows that our simulations were mainly underestimated.

Optimization of ground thermal parameters
To ensure that the initial temperature conditions did not influence the results, we ran a spin up with assigned initial ground temperature profiles (Fig. 1) corresponding to the middle of August 1980, until the soil temperature profile reached equilibrium with the upper and lower boundary conditions.The equilibrium temperature profile was determined when the maximum difference of soil temperatures at all levels between two successive annual cycles was less than 0.01 • C and used as the initial condition.
After 30 yr of simulations, the MAGTs at 1 m depth in the western and southwestern regions of Alaska were slightly warmer than the measured ground temperature data from those areas (Fig. 8).These parts of Alaska correspond to subarctic oceanic and continental sub-arctic climate.During spring, when the Bering Sea is ice free, the moderating influence of the open water helps to melt the snow early for some areas adjacent to the sea, when winter temperatures are more continental due to presence of sea ice (Serreze and Barry, 2005).The mean annual average air temperatures range from −1 • C to 2 • C from west to southwest of Alaska.These regions correspond to ecosystem protected permafrost zones that have formed under colder climate conditions and currently persist only in undisturbed late-successional ecosystems (Shur and Jorgenson, 2007).The MAGTs for the areas with a sufficient amount of organic cover and soil moisture content usually experience gradual increase even when MAATs are slightly higher than 0 • C (Jorgenson et al., 2010).Therefore, an additional layer of the organic matter might provide permafrost with neccessary resilience.However, high variability of precipitation and thermal properties of mineral soils in western and southwestern parts of the region makes the choice of the appropriate additional organic layer (AOL) non-trivial.
To address this issue, we developed an algorithm, which assigns optimal additional organic layer(s) for every grid point, based on deviation coefficient between the modeled equilibrium temperature profile and the assigned initial ground temperatures.During model calibration we tested the effects of varying climatic and ecological conditions on ground temperatures, and developed nine classes of addi- tional organic layers.The additional organic layer(s) varied from thinner to thicker, had different amounts of soil moisture and slightly different thermal properties.
Our algorithm is based on the following principle.If the initial temperature profile represents actual ground temperature distribution for the year 1980, then equilibrated ground temperatures should not deviate significantly from its initial temperature distribution profile.Otherwise, we assign additional organic layer(s) in successive manner.The layer(s) corresponding to the smallest deviation between equilibrated and initial ground temperature profiles was assigned as the final additional organic layer(s) at the corresponding location.
The obtained additional organic layer mask (Fig. 9) excludes lakes, rivers and mountain areas, and shows places where differences between initial temperature profiles and equilibrated temperatures vary significantly.The places required an additional organic layer(s) located in areas where we do not have many observation stations, mostly in the discontinuous permafrost zone.The amount of AOL in the northwestern part is thin and not so extensive in comparison to the southwestern territories.The thickness of the AOL is becoming more diverse in south and southwest parts of the region.The AOL mask indicates the places where the initial temperature profiles should be adjusted.As mentioned in the methods section, the initial temperature profiles were assigned according to the measured data, which were available for a limited number of places and does not cover the entire region.
The difference map between MAGTs at 1 m depth simulated with and without additional AOL (Fig. 10) showed colder annual ground temperatures in the western and southwestern parts of Alaska, which might more closely represent current ground temperature distribution.This method does not necessary guarantee that assigned AOL is going to cool the 2010 MAGT after spin up and 30 yr of simulation.Nevertheless, we were able, overall, to cool the MAGT for year 2010 for south and southwest areas (Fig. 10).This approach emphasized the importance of improved maps of soil organic layers and improved initial temperature profiles, and the necessity of further development of the permafrost temperatures observation system, especially in western and southwestern parts of Alaska.

Results
According to the decadal average of MAGTs histogram over twelve decades from 1980 to 2099, the overall area covered by MAGT less than 0 • C is going to decrease at different rates.The model projects the areal increase of the averaged decadal MAGT warmer than 0 • C at 2, 5 and 20 m depth with rates 3.7 %, 3.5 % and 2.4 % per decade, respectively (Fig. 11).
The spatial snapshot of MAGTs at 2 m depth below the surface (Fig. 12) shows pronounced warming almost everywhere including current continuous permafrost areas (e.g.Seward Peninsula and the south part of Brooks Range).As it was mentioned in Optimization of Ground Thermal Parameters section, the amount of ground surface organic material may retard the permafrost thaw in the discontinuous permafrost zone.This effect can be carried onto continues permafrost zone via changes in vegetation.Therefore, discontinues and sporadic permafrost areas with small or no organic layer and little soil moisture content will be more vulnerable to the rapid permafrost thaw.High altitude areas such as the Chugach Mountains, the Wrangell Mountains, the Alaska Range and the Brooks Range maintain relatively stable MAGTs during first half of the 21st century due to cold annual air temperatures.
The MAAT dynamics for 120 yr from downscaled five GCMs composite showed larger positive temperature trend during the 21st century for the northern region.This high positive trend can be observed almost everywhere north of the Brooks Range.The simulated mean annual snow depth showed that the amount of annual snow fall decreases in the south and increases in the north, when the number of snow-free days increase in the whole region due to warmer MAATs.The increase of snow-free days and at the same time, the increase of the thickness of snow greatly affects the mean annual surface temperatures in the northern part of Alaska.Eventually, this effect propagates further to the ground and transfers to the MAGTs.Figure 13 illustrates the higher trend of the MAGT for the two northern sites in comparision with central and south-central locations.

Discussion
The GIPL2-MPI transient numerical model with proper input parameters is a valuable tool for mapping the thermal state of permafrost, and its future dynamics with high spatial resolution.However, it is important to understand the limitations of the current model and the downscaled GCMs composite dataset.
The composite of five downscaled GCMs simulate well the seasonal cycle variations of near-surface temperature with a correlation between models and observations of 90 % or higher (Figs.2-4).However, the precipitation bias still remains high, correlation between GCMs and observations is 50 to 60 % (Bader et al., 2008).
The stations where GIPL2-MPI model showed high discrepancies with observed MAGTs (Fig. 6) were established recently (2005)(2006)(2007)(2008)(2009) and do not have long term data for more comprehensive analysis.The significant number of those stations are located along the rivers and in populated areas.The strong differences between measured and simulated MAGT might be caused by changes in surficial geology (due to flooding), ground water movement (convective heat transfer) or anthropogenic disturbances.
Besides the factors described above affecting the permafrost ground temperature simulations, forest fires are also an important factor.As a result of forest fires, the surface albedo decreases and the soil thermal conductivity of the surface soil layer increases (Hinzman et al., 1991).The areas where wildfires removed the upper organic layer are vulnerable to the active layer increase and warming of permafrost.If the burned area is ice-rich then the increase of active layer might melt the buried ice, which would cause the soil to collapse and form thermokarst depressions (Yoshikawa et al., 2002).In order to address the dynamics of the mean annual ground temperatures of the burned area, it is important to include the development of the vegetation and soil organic layer after the fire event.At the current stage the model does not include the effect of forest fires.
The MAGTs in continuous permafrost zone are mostly climate-driven, as apposed to the discontinuous permafrost zone, where the effect of the ecosystem on the permafrost thermal state is more pronounced (Shur and Jorgenson, 2007).With climate warming, present day continuous permafrost will turn into discontinuous, and vegetation change can develop enough upper organic layer providing additional resilience to thaw.Therefore, an introduction of the dynamic vegetation layer might decrease current modeling bias.The impact of humans and wildfires have to be taken into consideration where it is necessary.Further work is needed to improve parametrization of soil properties for each type of surface and soil conditions.Development of the spatial soil moisture map for Alaska will improve understanding of the soil moisture distribution and its dynamics as well as the results of permafrost modeling.Methods for calculating snow depth and snow thermal properties require further improvement.The 2 km distance between adjacent points is sufficient to neglect lateral heat transfer.However, for modeling watershed areas with grid resolution substantially finer than 1 km, the convective heat transfer by ground water movement, most likely, needs to be taken into account.Therefore, coupling the current model with a hydrological model could be an important step towards better simulation results for watershed and wetland areas.

Conclusions
The increase in mean annual air temperatures and the amount of precipitation for the northern part of the region predicted by the five GCM composite with A1B emission scenario could accelerate permafrost warming in the North (Fig. 13).The central part of the region, depending on the corresponding ecosystem, would experience permafrost degradation at different severity levels.The effect of the upper organic layer, soil water saturation and soil thermal properties play a significant role, providing necessary resilience for permafrost to sustain, even when MAAT are close or slightly above 0 • C. To provide an estimate of the permafrost degradation severity level for a specific geographic location, the effects of varying climatic and ecological conditions require more detailed investigation.
According to the modeling results, average areal decadal permafrost degradation at 20 m depth can be pursued with 2.4 % rate per decade, which concludes that for the next 90 yr Alaska is going to lose about 22 % of its currently frozen ground.Further analysis and development of the model is required to improve the MAGTs and ALTs simulations.In order to cover more accurately the entire region and to decrease uncertainty in the model predictions, more permafrost temperature observation stations are necessary.The increase in the grid spatial resolution requires high resolution maps of surficial geology, precipitation, vegetation, surface organic layer and soil moisture content.

Fig. 1 .
Fig. 1.Permafrost observation station locations and 18 ground temperature zones.Each ground temperature zone corresponds to its own initial ground temperature distribution profile.

Fig. 8 .
Fig. 8. Simualted mean annual ground temperatures at 1 m depth for the year 2010.

Fig. 10 .
Fig. 10.Projected differences between MAGTs simulated with and without additional organic layer(s) at 1 m depth for year 2010.

Fig. 11 .Fig. 12 .Fig. 13 .
Fig. 11.The amount of area over entire State of Alaska occupied by colder and warmer than 0 • C MAGTs averaged over ten years time interval from 1980 to 2099 at different ground depths.

Table 1 .
The model error statistics obtained by comparing MAGTs and ALTs from 3 different datasets: 60 deep borehole stations compared for 2007-2009 years; 77 US School project stations from relatively shallow depths compared for 2009; 43 averaged active layer thicknesses from CALM stations compared over entire available time periods.