Multisensory satellite observations of the expansion of the Batagaika crater and succession of vegetation in its interior from 1991 to 2018

On-site monitoring in large areas located in inaccessible regions can be difficult and costly. Thus remote sensing is an essential tool for mapping and monitoring changes in such regions. Therefore, this paper describes long-term multisensory satellite observations of the expansion of the Batagaika crater in Northern Siberia and natural succession of vegetation in its interior from 1991 to 2018. Landsat 5 TM, Landsat 7 + ETM, Landsat 8 OLI/TIRS imageries were mainly used as a data source for analyses, although Sentinel-2A imagery and DEM image from ASTER satellite were also employed for calculating a vegetation index and expansion in the crater area. The observations were conducted in years 1991–2018 and were made in a summer season. The results reveal that the crater area increased by almost three times during these 27 years and that the fastest expansion took place between 2010 and 2014 with 22.7% increment. The analysis of elevation of the crater revealed that in 2018 its maximum depth was ca 70 m and that depth was decreasing towards its north-east tail. Additionally, the satellite imagery of land surface temperature which is a driving force of crater expansion was visualized for chosen hot days within the time frame 2010–2018. The study of temporal and spatial changes in NDVI spatial distributions inside the crater revealed also a high rate of the succession of vegetation, which may reduce melting of permafrost inside the Batagaika crater and its further expansion.


Introduction
Permafrost is defined as a frozen ground that remains at or below 0 °C for two or more years and it underlies about a fifth of the land surface of the Earth (Schirrmeister et al. 2011). Permafrost is vulnerable to degradation and thaws due to climate warming (Hinzman et al. 2005, Babkina et al. 2019).
Significant amounts of permafrost deposits are found at Batagay in the Yana Uplands of northern Yakutia, which is commonly known as Batagay mega-thaw slump or the Batagaika crater (Ivanova et al. 2003, Günther et al. 2015. The gigantic crater was formed in the 1960s of the twentieth century as a result of deforestation. The transpiration of trees was cooling the surface but later, the ground was no longer shaded by the trees in the summer. Consequently, the sunlight slowly warmed the ground. As it was mentioned above, the Batagay region was abundant with the thawing permafrost, which led to the quick expansion of the crater (Nace 2017). Due to the high thaw rate of permafrost, the mega-thaw slump reached a width of up to 80 m in 2014 (Günther et al. 2015;Murton et al. 2017). It exposes a profile of yedoma deposits, reaching a thickness ranging from 7 to 22 m (Kunitsky et al 2013) and underlying ice-rich periglacial alluvial sand deposits of around 60 m thickness (Slagoda 2004). Similar craters are commonly created due to global warming over vast areas of difficult access in Northern Siberia (Osipov et al. 2019;Osawa et al. 2010); hence satellite observations of such phenomena allowing for monitoring large areas with appropriate time resolution are of particular importance.
The aim of the study was systematic satellite observations of the expansion of the Batagaika crater over the three last decades (precisely, for the timeframe from 1991 to 2018), including observations of the succession of vegetation in its interior. In this work, optical remote sensing data acquired from Landsat and Sentinel satellites were used to observe the expansion of the crater between 1991 and 2018. The detailed analysis of the depth of the crater was also an important part of the study. For this purpose, the profile graph method using digital elevation model (DEM) was employed as input data. Measurement of land surface temperature gives an insight into the ground temperature changes in the Batagaika crater between 2010 and 2018. Landsat satellite imagery was used to observe the land surface temperature. The determination of spatial and temporal distributions of the Normalized Difference in Vegetation Index (NDVI) inside the crater made it possible to understand the succession of vegetation at different stages of the development of the crater. The NDVI distributions were calculated mainly from Landsat imagery. However, since the imagery of Sentinel-2A has become recently accessible, it was used to check some analyses, for instance for more precise delineating NDVI distributions at the studied area. It is worth to stress that it was the first time when such systematic analysis of vegetation cover over the Batagaika crater was made. The beginning of the paper is the master thesis written by one of the authors (Vishakh 2018).

Study area
Located roughly 10 km southeast of Batagay, in the centre of the Verkhoyansk district, Sakha Republic (Yakutia) on the left bank of the Batagay River, the Batagaika Crater has a characteristic shape resembling a tadpole (Fig. 1). The center geographical coordinates of Batagaika crater are 67° 34′ 54″ N, 134° 46′ 44″ E. Batagay is characterized by a continental subarctic climate, according to the climate classification of Köppen (2011). According to Lydolph (1985) the temperatures in the region can vary considerably during the year: the mean air temperature ranges from + 15.5 °C in July to − 44.7 °C in January and the temperature range from difference between the absolute minimum (− 67.8 °C), recorded in 1892 to the absolute maximum (+ 37.3 °C), equals 105.1 °C. The mean annual precipitation is only 181 mm, with the lowest rate during the winter (13% of the annual precipitation) and the highest rate during the summer months (51% of the annual precipitation). In such climatic conditions, soils were formed which are characteristic for zones of the cryoarid steppe (Elovskaya et al. 1979). The Batagay region is covered by an open forest dominated by larch (Larix cajanderi) with ericaceous shrubs and a lichenmoss ground layer (Murton et al. 2017). In openings and above the tree line in the uplands, shrubs of birch (Betula middendorffii), alder (Duschekia fruticosa) and P. pumila occur. Chosenia arbutifolia and Populus suaveolens form stands on river floodplains (Murton et al. 2017).

Landsat processing
The spectral information from Landsat images was used to analyse the expansion of the Batagaika crater and vegetation distribution inside it. The visible and near-infrared bands were used to extract vegetation information, while the thermal band was used to extract the land surface temperature. Minimum cloud images of Landsat 5 TM (Thematic  Table 1. Since the satellite data from 1992 to 1998 were not available in the USGS website, the expansion of the crater could not be observed during the whole period of time. This was due to the unavailability of cloud free images and because of its high spatial resolution in the thermal band. Moreover, the Global Digital Elevation Map, version 2 (GDEM-2) was used to determine elevation changes around the crater. GDEM-2 is generated from imagery taken by Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) instrument onboard Terra satellite. According to literature, the GDEM-2 absolute vertical accuracy assessment expressed in root mean square error (RMSE) in most cases, even in mountainous regions, was in the range 6-9 m. (Tachikawa et al. 2011, Gómez et al. 2012, Satgé et al. 2015, Morais et al. 2017).

Estimation of area changes and expansion rate from 1991 to 2018
ArcGIS 10.3 software was used to calculate the area and the perimeter of the Batagaika crater. To do so, six datasets of the Landsat and Sentinel 2 were collected, namely from 1991, 1999, 2005, 2010, 2014 and 2018. Hence, in this paper, the shape of crater over three decades has been discussed. In the head part of crater, where expansion was the fastest, seven (out of ten) transverse interpolate lines through the crater were created to measure its width at each observation time. This part of the analysis was carried out using 3D Analyst Tool in ArcGIS 10.3 software.

Estimation of the elevation of the crater
To estimate the elevation of the Batagaika crater the Profile Graph method was used together with the ASTER DEM image taken on 17th October 2011. Figure 2 shows the DEM image of the Batagaika crater. The exact shape of the Batagaika crater using a Landsat image which was taken on 12th August 2011 was traced and the shape on the ASTER DEM image was marked using the interpolate line tool. Then, ten interpolate lines through the crater using the 3D Analyst Tool were created. Figure 2b shows the DEM image of the crater with ten interpolate lines. The first longitudinal line (number 1) has a length equal to the length of the crater, namely, c.a. 2.1 km. The other transverse nine lines indicate the width of the crater. Using the Profile Graph option in ArcGIS software, the elevation along all interpolate lines was determined, which helped to obtain precise and easy to interpret information about geometry of the Batagaika crater.
The NDVI values range approximately between − 1 and 1. The higher the value of this index the denser and healthier the vegetation is, whereas the negative values are commonly attributed to the presence of water (Carlson and Ripley 1997).
For the calculation of land surface temperature (LST) Landsat imagery from 2010 to 2018 was used. Landsat 7 images for the timeframe 2010-2013 were collected and Landsat 8 images were used for the timeframe 2014-2018. Only minimum clouded images which give the clear vision of the study site were chosen. Simultaneously, to clearly show the changes in LST in the Batagaika crater, specific hot days in the hottest month of each year were chosen. All LST images were created using ArcGIS 10.3 software, using corrections dependent on the emissivity, NDVI and land cover type, as it was described in detail by Jeevalakshmi et al. (2017), which assured that error in LST should be not greater than 3 °C. To retrieve the LST with better accuracy, accurate in situ emissivity measurements should be made. However, for the purpose of understanding of the expansion of the Batagaika crater the accuracy of LST which was used, was sufficient.

Expansion of the crater from 1991 to 2018
In Fig. 3, which shows one of the most essential results of the remote sensing study, the expansion of the Batagaika crater over the last 27 years (from 1991 to 2018) is delineated in detail. This figure has been created using the editor tool in ArcGIS 10.3 software. It is clearly visible that the crater is growing faster towards south-west, where it occupies a substantial area with a round shape, while its northwestern extremity remains relatively narrow. The area of slump obtained on 13-08-1991 and on 29-06-2018 was 0.19 km 2 and 0.78 km 2 , respectively, i.e. the area of the crater increased three times during the studied period since 1991. All assessed values of the area and the perimeter of the crater for the timeframe between 1991 and 2018 are given in Table 3. The increase in area of the Batagaika crater in the period from 1991 to 2018 is confirmed by satellite images shown in Fig. 4. The image in 2018 was taken by Sentinel-2A while the other images were taken by Landsat satellite. It may be clearly seen in Fig. 4 that Sentinel-2A satellite provides satellite images of better quality than Landsat one. We have also calculated the increase rate of the crater during the studied period Table 4 indicates the increase rate from 1991 to 2018. It turned out that it underwent significant temporal changes. For example, the increase rate of the Batagaika crater between 1991 and 1999 was ca 0.016 km 2 /year. During the next years, the increase rate of the crater was higher, namely ca 0.039 km 2 /year during 2010-2014, and 0.022 km 2 /year during 2014-2018. The average expansion rate between 1991 and 2018 was ca 0.026 km 2 /year.
Using 3D analyst tool in ArcGIS, the image of the head part of the Batagaika crater with seven interpolate lines was created. It is shown in Fig. 3. These lines are used to figure out the expansion of crater along each line within the period 1991-2018. Table 5 indicates the changes in the width of seven interpolate lines in the head part of the Batagaika crater within five different subperiods covering the timeframe 1991-2018. The obtained results show that the head part of the crater expanded the fastest during subperiod 1991-1999. In Fig. 5 it is also apparently a characteristic parabolic shape of the head part of the Batagaika crater. The maximum width of the Batagaika crater in 2018 was 979 m.

Profiles of the elevation of the Batagaika crater
In Fig. 5 all calculated profiles of the elevation of the Batagaika crater above sea level (a.s.l.) are shown. In Fig. 5a the span of the abscissa of the elevation graph is equal to the length of the crater (2.1 km), while in Figs. 5b-j these spans are equal to the widths of the crater corresponding to transversal lines 2-10 in Fig. 2. The ordinates of the diagrams in Fig. 5 indicate the elevation a.s.l. in meters. From Fig. 5b-j it can be easily concluded that crater is clearly the deepest and widest in its main southwestern part (the head of the "tadpole") and quickly narrows and shallows approaching the north-east end (the tail of the "tadpole"). It is also easy to notice that the terrain is lowered towards the north-east, where water from the crater flows down to neighbouring rivers. Thanks to this the crater is not filled completely with water.

The distributions of LST inside the Batagaika crater on the hottest days of the year
According to Table 4, the highest crater expansion occurred for the timeframe 2010-2014, but it is visible that during the next four years (from 2014 to 2018) there was not any significant expansion. Air temperature during this period was high in Batagay region (World Weather Online), so that land surface was also intensively heated. This is also in accordance with measurements of thaw layer depths in Northern Siberia (Global Cryosphere Watch), which show that in the period 2010-2014 thaw layer depths were higher than after 2014.
Increasing ground temperature is the main driving factor for the thawing of permafrost, so it was decided to get deeper insight into the LST spatial distribution inside the crater between 2010 and 2018. Figure 6 shows the examples of such LST distributions on low clouded, and hottest days of the year during the timeframe 2010-2018. Figure 6 clearly shows the high-temperature gradients inside the crater in those days. All of the images show that the warmest ground areas are located mainly in the centre and in the northern part of the crater, whereas the coldest ground areas, where permafrost is still melting are located mainly in its southern part, as well as on the edge of the "tadpole's" head. Right there the expansion of the Batagaika crater was the fastest in most of the studied periods. Thus, the LST observations confirm the fact that the expansion of the crater is driven by the permafrost thawing due to ground heating during warm months of the year (comp. Figs. 3 and  6).

Variation in the NDVI value over three decades
The satellite observations of the spatial distribution of vegetation inside the Batagaika crater allow to understand the overgrowing process related to the crater expansion within the last three decades. Figure 7 shows the changes in NDVI spatial distributions inside the Batagaika crater in the period from 1991 to 2018, whereas in Table 6 the main descriptive statistics of NDVI value inside the Batagaika crater are given.
Both from Table 6 as from Fig. 7 it may be concluded that noticeable and constant increase in the mean value of NDVI started in 2010. Before this year the mean values of NDVI inside the crater were negligibly small and fluctuated around zero, which is characteristic of bare soil, barren areas of rock, sand, or snow (Carlson and Ripley 1997). The changes in NDVI after 2010 may be interpreted as the result of natural succession of vegetation inside the crater. Comparing Fig. 7 with Fig. 5c-e it can been seen that the vegetation succession was the most intense in the northern part of the crater and its centre which has maximum depth where most of the permafrost has already melted.
The mean value of NDVI inside the crater in 2018, was equal ca 0.24, significantly higher compared with those in all previous years, which could imply that relatively dense vegetation can be already found in some parts of the Batagaika crater with no permafrost. The obtained mean values of NDVI in 2014of NDVI in , 2010of NDVI in , 2005 were ca 0.09, 0, and ca − 0.16, respectively. It is interesting that the expansion rate of crater  between 2014 and 2018 was lower than from 2010 to 2014 when the highest expansion rate was observed. Taking into consideration many possible causes of this slowdown of expansion rate in recent years such as climatic, geological, etc., it cannot be excluded, that one of them is a succession of vegetation, the presence of which reduced the heating of the soil (Osawa et al. 2010). Although this issue requires further research, the results already suggest the possibility of reducing the extent of permafrost thawing by appropriate artificial afforestation of newly created craters, or by activities aimed at maintaining or strengthening the forest cover on the areas prone to permafrost thawing.

Conclusions
Optical satellite observations have great potential for mapping and monitoring of permafrost disturbances over vast and hard-to-reach Siberian areas resulting in potentially threatening methane emissions. As an example, using multisensory satellite imageries, the expansion of the Batagaika crater over 27 years, namely between 1991 and 2018 was studied in detail. A systematic and comprehensive observation of the crater was of great importance because it may give deeper insight into processes forming such craters on the Northern Hemisphere. It was shown that during this period the area of the crater increased by almost three times from 0.19 km 2 , on 13-08-1991, to 0.78 km 2 on 29-06-2018. The study showed that major expansion took place between 2010 and 2014, but then the expansion unexpectedly slowed down. Using the profile graph method, the geometry of the crater was precisely analysed. In particular, it was assessed that the maximum length, width and depth of the crater in October 2011, were 2104 m, 979 m and 70 m, respectively. From the elevation analysis of the crater it was also concluded that its elevation was decreasing towards the tail of the crater. The expansion in two different parts of the crater, namely, in its head and tail is clearly different, probably due to different geological conditions. It was concluded that the parabolic shape of the head of the crater can be the result of thawing permafrost on the slope of the terrain. This has been confirmed due to the systematic observations of the shape of the head part of the crater as well as by the precise measurements of the crater width increment along with seven different directions since 1999. It was found that the highest expansion occurred during the hottest days between 2010 and 2014 and that high values of LST inside the crater are one of the main driving factors for this expansion. The observed LST in the northern part of the crater was much higher than in the southern one. This can be attributed to the fact that the southern part of the crater is abundant in permafrost and it is still thawing. At the same time, the northern part with high LST is deficient in the thawing permafrost (thermokarst), which means that huge amount of permafrost had already thawed inside the crater in this part. It is fully consistent with calculated profiles of the elevation of the Batagaika crater as well as with satellite vegetation observations. Detailed study of the spatial distributions of NDVI inside the Batagaika crater over three decades and its statistical results allowed us to conclude that the increase in NDVI which was observed after 2010 may be interpreted as the result of natural growing succession of vegetation inside the crater. In particular, it was observed in the Northern part of the crater in the period from 2010 to 2018 after the thawing of permafrost. It is important to note that the study also suggests that succession of vegetation inside the crater interior may delay further process of permafrost melting. It might be therefore possible that appropriate man-made afforestation could reduce the rate of the process of craters formation end expansion of Siberian permafrost.
We believe that the study may be helpful for scientists to better understand the formation and expansion of craters in permafrost, and possibly to build more advanced models of this expansion. The results obtained may also be supportive in planning future satellite observations of such craters.