Snow cover and permafrost evolution in Siberia as simulated by the MGO regional climate model in the 20th and 21st centuries

An approach to ground thermodynamics description using a coupled system of atmospheric regional climate and ground heat transfer models is improved by accounting for the time varying snow density. The simulations are compared to available observational analyses, and the sensitivity of the ground thermal regime to variable snow density is analysed. Projected changes of the ground thermal regime in the 21st century relative to the late 20th century are shown and compared to earlier estimates.


Introduction
Within the last few decades climate-permafrost interactions have received thorough scientific attention due to the possibility of permafrost degradation under climate warming. The future of permafrost remains a major point of discussion , Arzhanov et al 2007, Demchenko et al 2002, Euskirchen et al 2006, Nadyozhina et al 2008. The monitoring data of subsurface warming and active layer depth increasingly exhibit significant changes in many regions (Romanovsky et al 2007, Sazonova et al 2004, Pavlov 2003, Pavlov and Malkova 2005, Pavlov 2008). Permafrost and seasonal snow cover evolution are considered as cryospheric indicators of global climate change. It is known that permafrost properties are largely constrained by surface air temperature and snow cover. The influence of snow cover and surface temperature variations on permafrost properties has been extensively studied using observational and modelling data (Stieglitz et al 2003, Zhang 2005, Kitaev 2006, Sherstyukov 2008, Sosnovsky 2006. However, the studies were mostly focused on local analyses of snow cover and ground temperature. Recent studies (Dye and Tucker 2003, Konstantinov et al 2006, Groisman et al 2006 show that global warming may substantially affect the spatial distributions of snow cover and its duration; the changes in snow exhibit large interregional differences. Here we evaluate the spatial distributions of permafrost coverage and thermal properties as well as permafrost evolution in northern Eurasia using the Main Geophysical Observatory Regional Climate Model (MGO RCM) (Shkolnik et al 2000(Shkolnik et al , 2006 output. The RCM incorporates a description of all fundamental physical processes and climate system feedback. It is driven by the MGO general circulation model (Shneerov et al 2001).
Approaches to ground warming simulation are discussed in , , Romanovsky et al (2007) and Malevsky-Malevich et al (2001), (2005. There exist a number of modelling efforts to describe spatial variations of permafrost response to climate change. The climate forcing for such models can be provided by either a coarse resolution GCM with a typical grid cell size of the order of 10 2 km or a higher resolution RCM (grid step of about 10-50 km). Regional interactions between snow cover and permafrost reveal a great deal of spatial inhomogeneities and thus are best investigated using RCM, due to its more detailed representation of such regional features as orography, land surface and deep soil properties, coastlines, etc. The MGO permafrost model (PM) (Malevsky-Malevich et al 2001, 2005, 2007 was developed to be driven in off-line mode by GCM/RCM output. The model is used to investigate the snow cover influence on thawing and freezing processes, and ground temperature. The permafrost model employs a multi-level one-dimensional heat transfer scheme; vegetation and snow cover over land are treated as additional model levels. In previous studies , Nadyozhina et al 2008 we used the monthly mean data of surface temperature and water equivalent of snow at a constant density (200 kg m −3 ) produced by GCM/RCM.
In this study the evolution of ground temperature distribution using daily RCM output data is analysed, and the snow thermal properties and time varying density of snow are parameterized according to Mocko and Sud (2001). Our study aims at improving permafrost modelling with regard to snow physics. The daily output allows us to more correctly assess snow duration using model produced data instead of applying interpolation procedures to roughly estimate the duration using monthly mean data.
There exist various approaches to analysis of the observational snow data (Groisman et al 2006, Kotlyakov et al 2002, Kopanev 1982. In section 1 a comparison of the simulated snow cover with that obtained using different observation analyses is made, including those derived from satellite observation. The study is mostly focused on snow duration. In section 2 different approaches to snow density evolution are considered and compared against each other. In section 3 the analysis of interactions between various impacts on the permafrost evolution for climatic regions in the permafrost zone is performed. The model projection of the near surface permafrost degradation was developed for several slices in the 21st century using MGO RCM output data (section 4). The changes in snow cover and permafrost are evaluated as arithmetic differences between the future scenario and the late 20th century simulation. The modelling climate is assumed trend-free for each time slice. The RCM was run at 50 km resolution over the territory of Russia under the IPCC A2 scenario of greenhouse gas and aerosols increase (IPCC 2007). The permafrost model was used in an off-line mode, so potential feedback between permafrost and climate system are not accounted for. It is assumed that the snow properties are homogeneous within each RCM grid cell and vegetation distribution will remain unchanged due to climate warming. One has to use multiple IPCC scenarios to assess a range of uncertainties in future changes in the cryosphere due to different equally probable emissions. In spite of using a single model and particular IPCC scenario (A2), the key conclusions of our study on permafrost sensitivity to snow cover parameterization are expected to hold for other forcing scenarios and climate models.

Snow cover distribution: simulations and observational data
The complete and detailed description of current snow cover temporal evolution over Russia is available from Groisman et al (2006). In that paper the long-term in situ observational data are presented in summary fashion as the mean climatic surface characteristics. As a first step we used the RCM simulated monthly mean characteristics from the baseline period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). The simulated snow cover and that derived from observations (Groisman et al 2006) for  appear to be in reasonable agreement. In Climate of Russia (2001) the snow climatology is provided for a different period spanning 1961-1995. When comparing data averaged over different time slices it is assumed that modelling discrepancies are larger than those due to different representation of snow inter-annual variability in respective 10, 30 and 35 yr climatologies. The observation analysis over the vast territories of Russia (e.g. Siberia) with a very sparse observational network is prone to a great deal of uncertainty in snow cover estimation. It is also assumed that the uncertainty associated with snow density variation has the same order as the uncertainty in calculations, measurements and periods of data averaging.
In figure 1 the annual mean number of days with snow on the ground is shown for 1961-1990 in comparison to simulated snow cover duration. The snow duration has been assessed as an annual number of days with snow depth greater than zero. A reasonable agreement can be found between RCM simulated snow cover duration and that derived from the data provided by Groisman (2009). This is, however, not evident for the snow area boundaries that appear to be shifted with respect to observation analysis. The shift is most pronounced in the southern part of the Siberian permafrost area, where the model simulates about 15-20% longer snow durations. The discrepancies between the model and analysis can be explained by larger modelling errors and observation uncertainties over the complex regional orography. Another comparison has been conducted with a different observational analysis (Kopanev 1982). Despite the differences between observational estimates and their averaging periods, the results of both comparisons are generally similar. One can see the same model bias in the peripherals of the southern permafrost boundary. In the mountainous regions one can find some considerable differences between simulations and observational estimates. The model simulates little increase of the snow period duration in the southern regions of Siberia.
In Dye (2001) and Dye and Tucker (2003) the snow-free period is evaluated using satellite data. The segment of the map from Dye and Tucker (2003) is shown in figure 2(a). Figure 2(b) displays the simulated number of snow-free days within the same area. The satellite data are averaged over 1972-2000. The model satisfactorily portrays the picture obtained by satellite, however, the modelling distribution is characterized by substantially larger spatial variability. One can find the model simulated maximum snow extent slightly overestimates, according to the adopted classification, the extent derived from satellite measurements.

The role of seasonal snow density in cryosphere modelling
There are a number of papers concerning snow density (ρ s ) parameterization with respect to snow cover duration (Essery et al 1999, Osokin et al 2001, Pavlov 1980  2006). The duration of this period is estimated as the percentage of time in a year between the first and last occurrence of the snow amount exceeding some threshold. All the parameterizations typically include region dependent constants and appear not to be universal. There exist many empirically based parameterizations, one of which for ρ s was used in the SSiB (Simplified Simple Biosphere) Model (Mocko and Sud 2001). It has been sufficient to fit our goal in relating the changes of permafrost characteristics to snow cover distribution.
The water equivalent of snow depth, calculated by the RCM, was used to drive the permafrost model with different snow density properties. The relationship between snow density and snow cover duration is described using an exponential function, as has been done in other studies. The snow density is assumed to increase, starting at some minimum value ρ s min at the day of first snow fall up to maximum value ρ s max near the end of snow season. The change of ρ s is given by the following expression: Snow heat conductivity (λ s ) depends upon snow density. Firstly, the dependency ρ s = f (t) has been used to reveal the its effect at a regional scale. Afterwards, the dependency λ s = f (ρ s ) has been added and the combined effect has been estimated. The snow heat conductivity used is according to Mocko and Sud (2001) where λ 0 = 2.22 W m −1 K −1 , ρ a = 1000 kg m −1 . Figure 3 illustrates the spatial pattern of climatic snow depth h s distribution averaged over the 1991-2000 reference period. Figure 3 presents the snow depth distribution for February as simulated using three different approaches to setting snow density, with ρ s = 200 and 400 kg m −3 , and ρ s calculated during cold season according to expression (1). All three numerical experiments reveal similar patterns of snow distribution, with well pronounced areas of maximum snow accumulation. In terms of accumulated snow mass, the density calculation approach using expression (1) leads to moderate snow cover, as compared to snow cover simulated using constant values of ρ s that represent extremely high (ρ s = 200 kg m −3 ) and low (ρ s = 400 kg m −3 ) snow cover. Let us investigate the influence of snow cover on permafrost conditions at different soil levels.
Among permafrost indicators there are the depth of seasonal thawing and temperatures at different levels. Figure 4 displays the contemporary influence of snow density evolution and time varying heat conductivity on the spatial distribution of the ground temperature at 1.6 m depth (T 1.6 ) for Autumn (September) and Spring (May). It should be noted that September and May correspond to the annual maximum and minimum ground temperatures. The 1.6 m level has been chosen because it corresponds to a standard measurement level in the soil and thus is widely used in the investigations of snow-ground interactions (Kitaev 2006, Romanovsky et al 2007, Vasilyev 2009). It can be seen that the differences δT = [T 1.6 (ρ s = const) − T 1.6 (ρ s = const)] are positive. The effect of snow density is mostly pronounced in T 1.6 in March. A combined effect of both dependencies (ρ s = f (t) and λ s = f (ρ s )) on the temperature calculation revealed significant changes in the spatial structures of δT in September and March. The correspondence between the model simulated and observed boundaries of snow cover in Autumn and Spring depends weakly on the approach used for snow density parameterization. This feature can also be found in the spatial distribution of seasonal thawing in figure 3(c). Figure 5 also displays the dependence of the active layer depth in August on cold season snow density evolution.
The effect of variable snow density on maximum seasonal thawing is mostly noticeable within the areas of discontinuous permafrost. In some regions of Central and Western Siberia ignoring the influence of variable ρ s and λ s on seasonal cycle of temperature leads to larger biases in thawing depth by 0.5-1 m. It is important to notice that an increase in thawing of about 1-1.5 m is equivalent to seasonal freezing and, therefore, a considerable northward shift of the permafrost boundary.

RCM and PM simulated regionally averaged climatologies
The RCM output at 50 km resolution is a step forward to explicitly describe the effects of meso-scale forcing (topography, vegetation, coastal lines). The RCM allows us to provide realistic spatial details important to distinguish interregional snow cover evolution and regional ground thermodynamics.
Several climatic regions are included into the current RCM domain, e.g. Arctic and Siberia. In our study the selection of climatic regions is based on classification (Alisov 1957, Groisman et al 2005 with a focus on terrestrial permafrost areas. Three climatic regions were considered: one northward adjacent to Arctic Circle, and the two others in Western and Central Siberia (the regions 1, 5, 8 according to Groisman et al (2006)). The area averaged warm season air temperature sums, maximum depth of snow, snow duration and ground temperature at 1.6 m depth are presented in tables 1 and 2 for different climatic regions. The mean climatic values correspond to the 10 yr averaging period. The ground temperatures, as well as snow depth, were produced using constant and time variable values of ρ s and λ s .
The estimates of snow day number agree reasonably with that obtained in Groisman et al (2006) for all regions considered (table 1), although simulated values are slightly overestimated. Both model and observation analysis agree that the snow duration is most extensive in the far North of Eurasia. In Western Siberia modelled snow cover appears to be largest as compared against snow cover over the rest of the region. It can be seen from the estimates in table 2 that the combined effect of the time varying ρ s and λ s substantially affects ground temperature normals. Relative changes in the temperature normals are greater in September than in March. In tables 1 and 2, for each region, inter-annual variations of the area averaged ground temperature, summer temperature sums, and snow cover properties are shown. Owing to the relatively short simulations and vast areas of averaging it is difficult to separate the particular contributions from snow and air temperature impacts from the total impact on thermal state of the ground. One can only notice that the relatively high ground temperatures in region 5 are explained by the shortest snow duration, with the largest accumulated snow mass accompanied by the greatest temperature sums in summer.
The permafrost area is relatively small in this region. As noticed above, the influence of time varying snow density causes a shift in the permafrost boundary (most pronounced in Western Siberia).
The comparison between thawing depths averaged over different regions is uncertain due to the large spatial variability and variety of permafrost areas. However, the influence of different parameterizations of snow cover on thawing depths exhibits considerable interregional differences. Figure 5 displays the distribution of thawing depth change depending on constant and time varying snow density. Also shown are climatic region boundaries, as estimated by Groisman et al (2006). The effect of variable density is noticeable northwards to the Arctic circle and in the regions of discontinuous permafrost.

Permafrost evolution during the 21st century
The evolution of permafrost has been evaluated using RCM output for the mid-and late 21st century (2041)(2042)(2043)(2044)(2045)(2046)(2047)(2048)(2049)(2050)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098)(2099)(2100) according to the IPCC A2 scenario. In earlier studies (Malevsky-Malevich et al 2001, 2005 the changes in seasonal thawing depths were assessed using coarser resolution AOGCM output. The study (Nadyozhina et al 2008) presented RCM estimates of future permafrost evolution under an assumption of constant snow density. Figure 6 displays simulated seasonal thawing and freezing depths over Russia by the end of 21st century using time variable and constant snow density during the cold season. It can be seen that the approach to setting snow density significantly affects the projection.  Given a time varying ρ s that includes both intra-and interseasonal variability along with heat conductivity variations results in a shrinkage of the intermediate area between the seasonal thawing and freezing areas (see figure 6). This leads to a significant widening of both the seasonal thawing and freezing areas. Additionally, the seasonal freezing area tends to shift northwards. The case of constant snow density is characterized by an excessive southward shift of the seasonal freezing area. Thus, time variable snow density leads to a less dramatic scenario in shifting the permafrost boundaries but has a more pronounced effect in shifting the seasonal freezing boundaries northwards.

Summary
In this study the evolution of snow cover and permafrost conditions is simulated using a sophisticated model. The model employs climatic forcing at daily resolution (RCM simulated near surface temperature) and incorporates the calculation of a time varying snow density approach. RCM snow simulation has been compared against different observation analyses and satellite data. The thermal state of the ground strongly depends on snow cover duration and its depth. Most of observation analyses estimate spatial distributions of snow duration at reasonable levels of confidence, while great uncertainties exist for evaluating snow depth. In this study the snow duration data is used for comparison with RCM simulated snow cover. It has been found the RCM simulated snow duration is in reasonable agreement with the observation analyses. The spatial variability of simulated snow cover reveals slightly better agreement with satellite data analyses. The sensitivity of modelled permafrost to snow cover has been investigated. With considerable uncertainties, the thawing depth and ground temperature depend substantially on the approach adopted to parameterize snow cover density, other than the seasonal cycle of snow boundaries that does not reveal significant changes. However, the soil temperature is only weakly constrained by snow cover variations at the surface, as temperature changes Figure 6. Seasonal thawing/freezing depths (cm) as projected for the end of 21st century (2091-2100) relative to the late 20th century for monthly mean climate forcing and ρ s = const (a) and daily mean climate forcing and ρ s = f (t) (b). Isolines denote the respective thawing and freezing depths (cm). appear to be within the range of temperature inter-annual variability.
The area averaged soil temperature at different levels and active layers, as well as snow cover in different climatic regions, have been analysed. It has been shown that the interregional differences in the above characteristics are well captured by the RCM-PM system in regions with continuous permafrost.
The scenario of changes by the end of the 21st century in permafrost and thawing/freezing depths largely depends on snow density parameterization. It has been found that the permafrost boundaries migrate northwards more slowly than seasonal freezing boundaries, significantly affecting the local wetland ecosystems.