Retrospective Analysis of Permafrost Landscape Evolution in Yakutia during the Holocene Warm Intervals

The observed global warming has significant impacts on permafrost. Permafrost changes modify landscapes and cause damage to infrastructure. The main purpose of this study was to estimate permafrost temperatures and active-layer thicknesses during the Holocene intervals with significantly warmer-than-present climates—the Atlantic (5500 years BP), Subboreal (3500 years BP) and Subatlantic (1000 years BP) optimums. Estimates were obtained using the ready-to-use models derived by G.M. Feldman, as well as mathematical modeling taking account of the paleogeography of the Holocene warm intervals. The data obtained were analyzed to reveal the regional patterns of warming impacts on different permafrost landscapes. The study results will be useful in predicting future permafrost changes in response to climate warming.


Introduction
The topic of this study is of relevance to assess the stability and optimal development of the natural environment under current climate warming. It is already apparent that the onset of climate warming is changing the evolution of permafrost landscapes as a result of thaw processes in disturbed areas [1][2][3][4]. Yakutia is among the permafrost regions most affected by the warming trend.
The main task of these studies is to assess the change in the permafrost temperature and active layer thickness (ALT) in significant periods of warming in the Holocene-in the Atlantic, Subboreal and Subatlantic optimum, respectively, 5500, 3500 and 1000 years ago. Knowledge of these data, when the main changes in permafrost landscapes occurred with the degradation of permafrost, can serve as a control characteristic of the loss of stability and become the basis for predicting the state of permafrost landscapes in the future. Fragmentary and sporadic data on the temperature of permafrost and ALT available in the literature are not sufficient to understand the processes occurring during the warming periods in the Holocene, therefore, the conditions for the development of permafrost landscapes at this time are important in assessing modern changes in permafrost landscapes and their future state. To achieve this goal, the following tasks were set to assess changes in the permafrost temperature and ALT during the Holocene warming, taking into account changes in air temperature and precipitation during these periods in the tundra, northern and middle taiga landscapes. The choice of plain landscapes-tundra, northern and middle taiga-is due to the fact that it is in them that the most unstable Yedoma landscapes in terms of climate warming and anthropogenic disturbances are developed.
We estimated the changes in ground temperature and ALT in the tundra, northern taiga and middle taiga zones of Yakutia during the Holocene warm intervals based on paleogeographic data reported in the literature [5][6][7][8][9][10][11]. Permafrost conditions are seldom documented by paleogeographic studies. Geocryologists study in detail the conditions and time of formation of syngenetic permafrost [12][13][14]. The oxygen isotope composition of ground ice is used to reconstruct winter temperature conditions and mean annual ground temperatures with removed snow and vegetation covers [15].
Understanding the Late Pleistocene-Holocene evolution of climate and its impact on permafrost landscapes can provide guidance for predicting future permafrost changes in response to anticipated climate warming. This is one of the major tasks in geocryology and physical geography. The retrospective analysis of permafrost and climatic conditions is only possible through predictive assessments and modeling of permafrost conditions. Such assessment studies have not been made before, except for a few permafrost temperature estimates for selected periods of the Late Pleistocene and Holocene [16].
Based on paleoenvironmental interpretations of diatom [17] and biogenic silica [18,19] records from Lake Baikal, Fotiev [16] related their contents and concluded that permafrost formation in eastern Siberia began 3.1 million years ago. According to Fotiev [16], the climate during the Sartan epoch was very cold, with mean annual air temperatures in the range of 8-14 • C lower than those at present, correlating with low contents of biogenic silica (SiO2bio) of 3-5%. The formation of thick ice wedges in the Sartan sediments suggests that permafrost temperatures were very low as well. By the end of the Boreal period, mean annual air temperatures were 10-12 • C higher than during the Sartan. Then, as Fotiev [16] assumed, permafrost temperatures in the Baikal area could have dropped to −5 • C during the longest and coldest early Subboreal period.
Diatom records show consistent warming in the Baikal region between 3800 and 2000 years BP and Fotiev [16] suggested that permafrost temperatures in Central Siberia varied from −3 • C in the lower Vilyui River area to −10 • C on the Taymyr Peninsula. The formation of alas basins during the Holocene warm intervals is indicative of changes in permafrost terrain occurring at that time. During the last 2000 years, air temperatures were 3-5 • C lower in the colder intervals and 1-2 • C higher in the warmer intervals compared to the present day. However, the relict ground temperature field in Central Siberia with minimum temperatures preserved at great depths and the existence of ancient ice wedges in Central Yakutia [16] are evidence that permafrost did not experience widespread thaw during the warmer intervals of the Late Pleistocene and Holocene. The map of alas distribution compiled by Bosikov [20] shows the maximum values of alas coverage of up to 19% of the area, also implying the partial rather than ubiquitous degradation of the Ice Complex.
The proposal of the research is to reconstruct permafrost conditions during periods of climate warming in the Holocene, which will make it possible to assess the conditions under which periglacial processes and the formation of modern landscapes took place. This will become the basis for the assessment of the level of changes in permafrost landscapes in the conditions of modern warming and forecasting their state in the future, which is absent in modern studies.

Study Area
We examined lowland permafrost landscapes-tundra, northern and middle taiga-where the most characteristic Yedoma landscapes are unstable to climate warming and anthropogenic disturbances [3,21]. On the permafrost landscape map of the Republic of Sakha (Yakutia), these landscapes reflected as an inter-alas type of terrain [22]. The volumetric ice content of their deposits averages 0.6, and in the Arctic regions it reaches 0.8. The permafrost temperature in the tundra reaches up to −12 • C, and in the middle taiga, up to −2 • C. The ALT varied from 0.2 to 1.4 m in zonal types of vegetation.

The Data
For tundra landscapes, climatic data are given based on the results of long-term observations of meteorological stations Saskylakh, Tiksi and Chokurdakh, for the northern taiga-meteorological stations Olenek, Zhigansk, Verkhoyansk and Srednekolymsk, and for the middle taiga-meteorological stations Vilyuisk, Yakutsk, Churapcha and Ust-Maya ( Figure 1). Sections no. 29 and no. 30 of ice-rich permafrost with the physical and thermal properties considered to be representative of the Ice Complex sediments were taken as a lithologic substrate-loams with an ice content of 0.6, and the bulk density of dry loam is 980 kg·m 3 [23]. Climatic data-monthly air temperatures, thickness and density of snow cover-were taken from climatic reference books [24][25][26]. Measured ground temperatures were derived from the Permafrost Institute research reports, Turbina's [27] as well as our unpublished data. The thermal parameters of the active layer and permafrost were specified depending on the lithological section, soil moisture content, density, etc., in the generalizing work of Gavriliev [28] on the thermophysical properties of rocks and soils of the permafrost zone of Siberia. Soil moisture contents in the active layer were assumed to be 0.30 for the tundra, 0.25 for the northern taiga, and 0.13-0.20 for the middle taiga. The following soil thermal parameters were used [28]: from the surface to a depth of 3.0 m, the thermal conductivity, thawed is 0.96-1.51 W/(m • C), thermal conductivity, frozen-1.28-1.89 W/(m • C); heat capacity, thawed-768-872 W/m 3 • C, heat capacity, frozen-570-651 W/m 3 • C; freezing temperature is minus 0.2 • C. In the interval of depths of 3-15 m, thermal conductivity, thawed is 1.16 W/(m • C), thermal conductivity, frozen-2.15 W/(m • C); heat capacity, thawed-896 W/m 3 • C, heat capacity, frozen-605 W/m 3 • C; freezing temperature is minus 0.2 • C. The thermal conductivity and heat capacity of soils were set specifically for the composition and moisture content of soils, which are given in the intervals above.
Land 2020, 9, x FOR PEER REVIEW 3 of 11 meteorological stations Vilyuisk, Yakutsk, Churapcha and Ust-Maya ( Figure 1). Sections no. 29 and no. 30 of ice-rich permafrost with the physical and thermal properties considered to be representative of the Ice Complex sediments were taken as a lithologic substrate-loams with an ice content of 0.6, and the bulk density of dry loam is 980 kg·m 3 [23]. Climatic data-monthly air temperatures, thickness and density of snow cover-were taken from climatic reference books [24][25][26]. Measured ground temperatures were derived from the Permafrost Institute research reports, Turbina's [27] as well as our unpublished data. The thermal parameters of the active layer and permafrost were specified depending on the lithological section, soil moisture content, density, etc., in the generalizing work of Gavriliev [28] on the thermophysical properties of rocks and soils of the permafrost zone of Siberia. Soil moisture contents in the active layer were assumed to be 0.30 for the tundra, 0.25 for the northern taiga, and 0.13-0.20 for the middle taiga. The following soil thermal parameters were used [28]: from the surface to a depth of 3.0 m, the thermal conductivity, thawed is 0. The thermal conductivity and heat capacity of soils were set specifically for the composition and moisture content of soils, which are given in the intervals above.

Modeling
The model derived by Feldman et al. [23], as well as the results of numerical modeling with the HEAT software system developed at Moscow State University [29,30] were used to obtain permafrost characteristics for the Holocene warm intervals. Preliminary work employing these methods was performed earlier in the application to the Ice Complex of Central Yakutia [31]. Transient heat transfer problems in the climate-substrate system were solved in the HEAT software using air temperature trends expressed as departures from the "regional normal", as well as other data from the reference weather stations in accordance with the generally accepted estimates for individual Holocene periods.
The predictive models of ground temperature in Yakutia [23] are not widely utilized now; however, we took them as a basis because mathematical modeling by G.M. Feldman covered the

Modeling
The model derived by Feldman et al. [23], as well as the results of numerical modeling with the HEAT software system developed at Moscow State University [29,30] were used to obtain permafrost characteristics for the Holocene warm intervals. Preliminary work employing these methods was performed earlier in the application to the Ice Complex of Central Yakutia [31]. Transient heat transfer problems in the climate-substrate system were solved in the HEAT software using air temperature trends expressed as departures from the "regional normal", as well as other data from the reference weather stations in accordance with the generally accepted estimates for individual Holocene periods.
The predictive models of ground temperature in Yakutia [23] are not widely utilized now; however, we took them as a basis because mathematical modeling by G.M. Feldman covered the entire range of Yakutian permafrost with ground temperatures ranging from −14 • C in permafrost to +3 • C in taliks. In this study, we exploited the ready-to-use forecast tables provided by G.M. Feldman, as well as the cartograms presented in the guidelines [23]. Estimates of the ground temperature and ALT changes were made for the Saskylakh, Tiksi and Chokurdakh weather stations located within the tundra zone; the Olenek, Zhigansk, Verkhoyansk and Srednekolymsk stations in the northern taiga; and the Vilyuisk, Yakutsk, Churapcha and Ust-Maya weather stations in the middle taiga. Sections no. 29 and no. 30 of ice-rich permafrost with the physical and thermal properties considered to be representative of the Ice Complex sediments were taken as a lithologic substrate [23]. Soil moisture contents in the active layer were assumed to be 0.30 for the tundra, 0.25 for the northern taiga, and 0.13-0.20 for the middle taiga. The following soil thermal parameters were used [28].

Validation
To verify the data on ground temperature and ALT changes, geocryological forecasts were made for Saskylakh for the tundra, Srednekolymsk for the northern taiga and Yakutsk for the middle taiga for 100 years with mean annual air temperature (MAAT) trends calculated for the considered periods of warming. MAAT trends were calculated using the data suggested by paleogeographers [5,6,8].
The snow cover data included in the modeling was adapted to their proposed precipitation data. For Srednekolymsk and Yakutsk, our measurements at the sites were taken as the basic characteristics of ground temperature verification. For Saskylakh, due to the lack of data, we used data on the ground temperature in the Tiksi tundra area. The calculations were carried out using the HEAT software system developed at the Department of Geocryology, Moscow State University [29,30]. The program considers monthly air temperatures, changes in the density and depth of snow cover, and trends in MAAT. When carrying out the calculations, we used the procedures for setting predictive modeling problems, which were successfully applied to the permafrost of western Siberia [32][33][34].

Tundra (Saskylakh, Tiksi and Chokurdakh Weather Stations)
The present (before 1980) mean annual air temperatures in the tundra zone range from −14 to −14.5 • C, and annual precipitation is from 205 to 340 mm. Permafrost temperatures vary from −6 to −8 • C in the southern shrub tundra and from −10 to −12 • C in the Arctic herb tundra [22]. ALT values range from 0.4 to 0.6 m and from 0.2 to 0.4 m, respectively.
During the Atlantic Optimum (ca. 5500 y BP), spruce-birch-larch forests expanded into the shrub-moss tundra [35]. An increase in temperature by 3 • C and an increase in precipitation by 70 mm caused an almost critical state of permafrost. Permafrost temperatures during this period were 5-6.5 • C higher and the active layer was 0.5-0.7 m deeper than present (Table 1). Modeling using the Saskylakh station records and geotechnical section no. 30 according to Feldman et al. [23] has shown that the permafrost temperature during the Atlantic Optimum was 4.8 • C higher than today.
In the Subboreal period (ca. 3500 y BP), increasing MAAT in the Tundra zone were 2 • C warmer and annual precipitation was 100 mm higher than today, which caused the significant warming of permafrost. Permafrost temperatures might have increased by 3.5-5 • C, and active-layer thicknesses by 0.3-0.5 m. The results of modeling for Saskylakh suggest a 3.3 • C increase in permafrost temperature.
During the Subatlantic period (about 1000 years ago), the increasing mean annual air temperatures in the Tundra zone were 1.3 • C warmer and annual precipitation was 50 mm higher than at present, thus the permafrost temperatures may have been increased by 2-3 • C (according to modeling in Saskylakh by 2.3 • C), while ALT was up to 0.1 m deeper.

Northern Taiga (Olenek, Zhigansk, Verkhoyansk and Srednekolymsk Weather Stations)
The present (before 1980) mean annual air temperatures in the northern taiga zone range from −12 to −15 • C. Annual precipitation is from 175 to 310 mm. Permafrost temperatures range from −4 to −6 • C, while the average ALT values are from 0.6 to 0.8 m.
During the Atlantic Optimum (ca. 5500 y BP), with the predominance of birch-alder associations with the participation of conifers and the undergrowth of shrub birches [36], an increase in the mean annual air temperature by 2 • C and an increase in precipitation by 50 mm increased permafrost temperatures by 2.5-4 • C and ALT by 0.35-0.40 m. Modeling for Srednekolymsk suggests a 3.2 • C higher permafrost temperature compared to the present value.
In the Subboreal period (ca. 3500 y BP), the mean annual air temperature in the northern taiga zone, increased by 1.5 • C compared to the present, and the annual precipitation, which was 75 mm greater, were associated with permafrost temperatures 2.5-4.5 • C higher and an ALT 0.4-0.5 m deeper than present. According to the modeling, the deviation of permafrost temperature from the present value is +2.7 • C.
During the Subatlantic period (about 1000 years ago), due to an increase in the mean annual air temperatures in the northern taiga by 1.5 • C and annual precipitation 50 mm higher than that at present, the permafrost temperatures may have been increased by 1-3.5 • C (according to the modeling by 2.2 • C), while the ALT could have been 0.35 m deeper.

Middle Taiga (Vilyuisk, Yakutsk, Churapcha and Ust-Maya Weather Stations)
The present day (before 1980) mean annual air temperatures range from −9 to −11.5 • C and annual precipitation is in the range of 250-300 mm. Permafrost temperatures vary from −2 to −4 • C, while ALT values are between 1.2 and 1.4 m in forests and between 1.8 and 2.2 m in meadows.
During the Atlantic Optimum (ca. 5500 y BP), with the predominance of larch and birch forests with pine, with an increase in the mean annual air temperature in the middle taiga by 2 • C and an increase in the annual precipitation by 50 mm, the permafrost temperature was 1.5-2.5 • C higher, and the ALT 0.1 m deeper than today. The modeling performed for the Umaibyt site near Yakutsk, Central Yakutia, using the data presented by Turbina [27], indicates that the permafrost temperature was then 1.8 • C warmer compared to the present day.
During the Subboreal period (ca. 3500 y BP), under the development of larch forests with the participation of pine and birch, as well as dwarf birch trees, with an increase in the mean annual air temperature in the middle taiga higher than the modern one by 1.5 • C and an increase in annual precipitation by 50 mm, the permafrost temperatures increased by 1-2 • C (1.2 • C according to the simulation results) and an increase in the ALT by 0.1 m.
The Subatlantic period (ca. 1000 y BP), with an increase in the mean annual air temperature in the middle taiga by 1 • C and an increase in annual precipitation by 20-40 mm from modern values, the permafrost temperature increased by 0.5-1.5 • C (0.5 • C according to the simulation results), and an increase in the ALT-up to 0.10 m.

Discussion
Ice wedges occur from the depths of 0.5-0.9 m in the tundra, 1.7 m in the northern taiga, and 2.1 m in the middle taiga [37]. The estimated results of an increase in the ALT up to 0.5-0.7 m in the Atlantic and up to 0.4-0.5 m in the Subboreal periods obtained by us show that under natural zonal conditions only in the tundra and northern taiga, the natural risk of the degradation of ice wedges is possible. At levels of change in permafrost conditions in the Subatlantic period, no disturbances occur in natural landscapes. Changes in permafrost conditions in the Subboreal and Atlantic thermochrones of the Holocene have significantly affected the development of permafrost landscapes in eastern Siberia, which is confirmed by the data of Mathias Ulrich [38,39].
In the middle taiga, under natural zonal conditions, an increase in ALT cannot ensure the development of negative cryogenic processes, since the seasonal thawing does not reach the ice wedge. Here, this requires an additional external factor, which includes the violation of the vegetation cover.
There is, as of yet, no adequate knowledge of active-layer dynamics and permafrost warming in boreal forests in relation to the biomass growth of root and aboveground ecosystems [40] or to successional processes [41]. These factors may be responsible for the relative stability of permafrost temperature and active-layer thickness in the boreal forests of Central Yakutia reported by Varlamov et al. [42]. At present, this process should be further taken into account in modeling and forecasting the development of permafrost landscapes.
Research on the Late Pleistocene-Holocene evolution of periglacial landscapes in Central Yakutia is important for understanding the future trajectories of permafrost change triggered by current climate warming. Recent Russian-German studies in Central Yakutia [38,39] have shown that the modern appearance of alas landscapes is a complex series of expanding thermokarst features which have formed a single thermokarst basin. Detailed studies of alases on the Tyungyulyu (Khara-Bulgunnyakh alas) and Abalakh terraces (Yukechi alas) of the Lena River suggest thermokarst initiation on the Yedoma between 11,000 and 9000 y BP. The dating of the Yedoma deposits from the study sites shows ages from 17,000 to 41,000 y BP.
Intensive thermokarst activity occurred during the Holocene optimum from 7000 to 5000 y BP, when permafrost temperature on the Yedoma of the middle taiga of Central Yakutia was 2 . . . 2.1 • C higher than modern values. At this time, separate small thermokarst basins began to unite into one alas. Then, about 3.5 thousand years ago, in the Subboreal period, the permafrost was 1.4 -1.5 • C higher than modern values; according to paleogeographic data, lateral expansion was characteristic due to the next activation of thermokarst and the final formations of the configuration of Central Yakutian alases.
The three warm intervals of the Holocene discussed here are considered to be the most likely future warming scenarios [5,9]. The strong modification of the landscape structure in permafrost regions during these periods was largely due to thermokarst. Thermokarst and thermoerosional basins on the Yedoma in Yakutia are evidence of climate warming which caused permafrost degradation and surface ponding [43]. It is estimated that the area of the Yedoma affected by thermokarst and thermal erosion during the periods of warmings periods in the Holocene was up to 10-19% in the middle taiga in Central Yakutia [20] and up to 60-75% in the tundra in coastal lowland [44][45][46]. In Central Yakutia, forest fires are thought to have played a significant role in thermokarst activity, based on charcoal records [47,48], which is further evidenced by ALT simulation results that show slight increases in undisturbed taiga.
The risk of natural hazards in permafrost regions is rising, as emphasized by many researchers [49][50][51][52][53][54]. The thickening of the active layer can cause rapid thermokarst and landslides. Increases in permafrost temperature can reduce the bearing capacity of foundations, increasing the risks for buildings and infrastructure.
Balobaev et al. [55] made a prediction of permafrost temperature fields in Yakutsk for the years 1960-2200 based on the combination of harmonics with periods of 300, 110, 75 and 14 years, with no account for anthropogenic air pollution. Although the authors admit that their results are not very accurate, they predict that the stability of the cryogenic system can reach a critical threshold, beyond which it will start to degrade. The authors assume an air temperature of −8 • C in 2017, lowering to −11 • C by 2054, rising to −8.5 to −9 • C by 2098, and decreasing again to −11.5 • C by 2200. Such temperature changes occurred during the warmer and colder intervals in the Holocene. Ground temperatures at the depth of 20 m are predicted to be −1.2 • C in 2025, −2 • C in 2050, −1.5 • C in 2105, −1.5 • C in 2150, and −2.3 • C in 2190. The prediction was made based on the ground temperature of −2.6 • C in 1960 and −1.7 • C in 2000. If we use the 1960 ground temperature as a reference, the maximum deviation within 200 years will be +1.4 • C, which is roughly equal with our modeling to the Subboreal warming (3500 y BP) in the middle taiga of Central Yakutia.
There have been numerous modeling studies across the world to predict spatial permafrost changes due to climate warming. For example, a model developed at the University of Alaska [56] predicts that by the end of the 21st century, permafrost in eastern Siberia will have crossed the threshold of "its stability". In Russia, the most severe permafrost degradation is projected for northwest Siberia and the European north. In Yakutia, permafrost will be thawing by the mid-21st century in its southwestern and southern parts; by the end of the 21st century, widespread permafrost thaw is predicted along the latitudinal course of the Vilyui River to Yakutsk on the Middle Lena and to Ust-Maya on the Middle Aldan River. By using the Late Pleistocene and Holocene permafrost degradation patterns as a guide for future conditions, however, particular attention should be given to the fact that permafrost landscapes are under anthropogenic pressure in the present warming period.

Conclusions
The retrospective analysis of the climatic and landscape conditions in Yakutia based on published paleogeographic studies has shown the main range of past changes. Predictive models indicate the greatest changes in permafrost landscapes in the Arctic zone. Warming during the Atlantic interval (5500 y BP) caused ground temperatures to increase by 5-6.5 • C in the Tundra, by 2.5-4 • C in the northern taiga and by 1.5-2.5 • C in the middle taiga compared to the present. These differences in response are attributed to the climatic characteristics of the landscape zones. The magnitude of warming during the Subboreal and Subatlantic intervals was less than during the Atlantic period.
Increases in ground temperature observed presently in Yakutia vary from 2-3 • C above the long-term mean in the Arctic [57] to 0.5 • C in Central Yakutia [3] during the 1980s, which allows us to state that the level of modern warming remains high.