Climate-induced landsliding within the larch dominant permafrost zone of central Siberia

Climate impact on landslide occurrence and spatial patterns were analyzed within the larch-dominant communities associated with continuous permafrost areas of central Siberia. We used high resolution satellite imagery (i.e. QuickBird, WorldView) to identify landslide scars over an area of 62 000 km2. Landslide occurrence was analyzed with respect to climate variables (air temperature, precipitation, drought index SPEI), and Gravity Recovery and Climate Experiment satellite derived equivalent of water thickness anomalies (EWTA). Landslides were found only on southward facing slopes, and the occurrence of landslides increased exponentially with increasing slope steepness. Lengths of landslides correlated positively with slope steepness. The observed upper elevation limit of landslides tended to coincide with the tree line. Observations revealed landslides occurrence was also found to be strongly correlated with August precipitation (r = 0.81) and drought index (r = 0.7), with June–July–August soil water anomalies (i.e., EWTA, r = 0.68–0.7), and number of thawing days (i.e., a number of days with tmax > 0 °C; r = 0.67). A significant increase in the variance of soil water anomalies was observed, indicating that occurrence of landslides may increase even with a stable mean precipitation level. The key-findings of this study are (1) landslides occurrence increased within the permafrost zone of central Siberia in the beginning of the 21st century; (2) the main cause of increased landslides occurrence are extremes in precipitation and soil water anomalies; and (3) landslides occurrence are strongly dependent on relief features such as southward facing steep slopes.


Introduction
Landslides are a widespread phenomenon within Eurasian and North American permafrost areas Studies have shown that in recent years warming in permafrost areas has resulted in an increase of landslide incidents, with landslides expected to be more frequent with continued warming temperature and an increase in precipitation (Montrasio and Valentino 2008, Blunden andArndt 2011, Shan et al 2015).
Substantial reduction in the range of the geographical limits of permafrost has been observed since 1975 in Russia (IPCC 2013). During the last four decades, an increase of permafrost temperatures of 0.3°C-2.0°C has also been observed in Siberia (Romanovsky et al 2010). Temperature increases of 2°C or greater may impact local industrial infrastructure, including the gas and oil industries (Anisimov and Reneva 2011).
The region of interest for understanding landslides on permafrost is large, remote and well suited for satellite scenes analysis and GIS techniques (Huscroft et al 2003, Chau et al 2004, Lyle et al 2004, Booth et al 2009. Remote sensing of landslides in forested areas is based on change-detection using vegetation indices (e.g., EVI, NDVI), or detection of the denudation of landslide beds (Chau et al 2004, Booth et al 2009. Starting from 2002, Gravity Recovery and Climate Experiment (GRACE) satellite mission has provided estimates of the Earth's gravitational field anomalies resulting, in particular, from change of water mass. GRACE data have been analyzed for water mass changes in the Arctic and Antarctic The goal of the work reported herein was to (i) estimate spatial pattern of landslides and their dependence on terrain elevation and slope azimuth and steepness, (ii) discover if there is a dependence of landslide occurrence on climate variables (air temperature, precipitation and drought index) and gravimetric soil moisture as estimated from GRACE data.

Study area
The study area is located within the northern part of central Siberia and includes the watershed of the Kochechum River including the Tembenchi and Embenchime tributaries. The area lies to the north and west of the small town of Tura and encompasses about 62 000 km 2 (figures 1 and 2). This site is with a portion of the Siberian Traps, a vast basaltic plateau dissected by many rivers. This is a hilly area with elevations ranging from 100 to 1100 m a.s.l. and permafrost thickness about 200-400 m. Active layer thickness is about 0.5-1.5 m within sediments and about 5 m within bedrocks (Ershov 1989). Sediments are composed of sandy and clay loams and contain about 10%-30% of fragmented debris. Ground ice content within sediments on the slopes reached 10%-15%. Kurums (rock fields) were found at the upper parts of the slopes. Solifluction processes are widespread within the area and observed mainly within the middle and lower part of slopes. Regularly, solifluction rates are low (about 1 mm yr −1 ) (Ershov 1989). Trees within the solifluction zone deviate from the vertical direction and form so called 'drunken forest'.
The mean permafrost temperature was about −2°C to −4°C on the middle and lower part of valleys, reaching −5°C to −7°C at the highest elevations locations (Ershov 1989). Forests are formed by larch (Larix gmelinii Rupr.) with rare birch (Betula pendula Roth) admixture. The mean crown closure for larch stands was about 0.2. Mean height, diameter at breast height and age were 8.5 m, 12.5 cm and 250 years, respectively, based on sample plot measurements acquired in 2007 and 2012. Ground cover was composed of small shrubs (Betula nana, Salix sp, Ribes sp, Rosa sp., Juniperus sp, Vaccinium sp), lichen and moss. Soils are cryogenic brown soils (Ershov 1998).

Climate
Climate within the study area is strongly continental with long cold winters and short warm summers The study area covers about 62 000 km 2 near Tura, Russia in northern Siberia as shown by the box. Insets: on-ground photo (above) and high-resolution satellite scene of a typical landslide.
(table 1, figure 3). Recorded maximum July temperatures have reached +39°C. Snow melting is observed regularly during the month of May. Stable snow cover is formed beginning in early October. Mean snow depth is about 40-50 cm.
We used daily, monthly, summer and annual air temperature, and also sum of positive temperatures in June-July-August (JJA), May-June-July-August (MJJA) and May-June-July-August-September (MJJAS), mean summer temperatures, and the number of thawing days N t>0 (i.e., number of days with positive temperatures: t max >0). Precipitation was analyzed as monthly and summer (i.e., JJA) values. In addition, a correlation of landslide occurrence with drought index SPEI was analyzed. The Standardized Precipitation-Evapotranspiration Index (the SPEI) can measure drought severity according to its intensity and duration (Vicente-Serrano et al 2010). The SPEI uses the monthly difference (D i ) between precipitation and potential evapotranspiration (PET): where T is the monthly mean temperature in°C; I is a heat index, which is calculated as the sum of 12 monthly index values, m is a coefficient depending on I, and K is a correction coefficient computed as a function of the latitude and month which takes into account number of sun hours in a day. SPEI data were obtained from (http://sac.csic.es/spei/database. html) and averaged for a cell size 0.5°×0.5°( ∼33×56 km 2 at the study location).
Data (monthly and daily air temperature and precipitation) from the weather station at the nearby town of Tura (coordinates: 64.27 N. 100.23 E.) were used in the analysis (http://aisori.meteo.ru/ClimateR).
Within the study area positive annual temperature trends were observed since the 1980s ( figure 3(a)). An increase of August precipitation was observed since the 1990s, whereas annual precipitation decreased (figure 3(b)). A positive August drought index SPEI trend (i.e., drought decrease) has been observed since the 1990s (figure 3(c)).

Satellite data
WorldView-1, -2 and QuickBird-2 high-resolution (pixel size 0.5-0.6 m) scenes (Neigh et al 2013), Landsat-5, -7 panchromatic band (pixel size 15 m), and gravimetric measurements (GRACE; http://www.csr. utexas.edu/grace) were used in this study. A time series of WorldView, QuickBird and Landsat scenes were compiled for landslide detection. High-resolution data (N=110 summer-acquired scenes) covered the period 2004-2012. Each scene was corrected (i.e. geometry and radiometry corrections) and covered an area of 320 km 2 (with total analyzed area about 17 500 km 2 ). These scenes (both black/white and spectral format) were used for landslide detection by manual photointerpretation. Landslides were identified based on texture and spectral characteristics and contextual information. The high quality of the  WorldView-1, -2 and QuickBird-2 scenes used allowed very accurate detection of landslides (e.g., insert on figure 1, and figures A1-A5 in appendix). A digitizing tablet was used to measure length width and area of landslides. There were no misclassification with separating landslides from burned areas and other disturbances. On the other hand, precise dating of landslides was complicated by lack of high-resolution data for the period 2005-2008. That problem was solved by using Landsat scenes for the analysis, since for each year we have at least three Landsat scenes. Summer acquired, good quality Landsat scenes (N=50) covered the period 1989-2012. For the period 1989-1999 only two scenes were available with the other 48 covering the period 2000-2012. Landsat scenes were obtained from USGS GloVis (http:// glovis.usgs.gov). GRACE data were used for soil water content estimation. We used annual and summer minimum and maximum gravimetric values, and equivalent of water thickness anomalies (EWTA, measured in cm). EWTA accuracy was approximately 10-30 mm month GRACE spatial resolution was 1°×1°degree (∼112×44 km 2 at latitude 66°). The data were processed using ERDAS Imagine software (http:// geospatial.intergraph.com) and ESRI ArcGIS software (http://www.esri.com).

GIS analysis
The distribution of landslides with respect to relief features (elevation, azimuth, slope steepness) was analyzed based on the ASTER global digital elevation model (DEM; http://earthexplorer.usgs.gov/). The ASTER DEM horizontal and vertical accuracy were ±20 m and ±30 m, respectively (https://www. jspacesystems.or.jp/ersdac/GDEM/E/index.html). The elevation range was quantized to 50 m strata. Aspect and slope steepness data were calculated using DEM and ArcGIS tools. The aspect data were quantized to sixteen directions (i.e., by 22.5°referenced clockwise from north). Because the distribution of relief features within the analyzed area was uneven, it could lead to bias. To avoid this, the data were normalized by the following procedure. The analyzed area with a given azimuth, slope steepness and elevation (shown by boxes on figure 2) was referenced to the total study area (rectangle on figure 2) with similar parameters: given on-ground class within the ith category of the topographic feature c, and A c(i)I is the area of the ith category of topography feature c over the entire analyzed territory. Statistical analysis of the data was carried out with Microsoft Excel and Statsoft Statistica (http://statsoft. ru) software. We used linear regression and Pearson correlation (r) analysis and Akaike information criterion (Akaike 1974) to determine significant relationships between landslides occurrence and climate variables and soil water anomalies (EWTA).

Landslides statistics
Analysis of satellite imagery found that a total 145 landslides occurred during the time period from 2000 to 2012 within our study area. All of the observed landslides were used in the spatial analysis of landslides occurrence. Out of the total, 31 landslides were excluded from the temporal analysis because they could not be dated with one-year precision.  4(b)). The landslide statistics including slope azimuth (aspect) and slope steepness are presented in table 2.

Landslides and relief features
The spatial distribution of landslides is strongly uneven with respect to azimuth as seen in figure 6(a).
The majority of landslides occurred at slopes with southern and south-western exposures (figure 5(a)). With respect to elevation, most landslides had an upper limit mainly within the range of 200-300 m a.s.l. Occurrence of landslides declined exponentially as elevation increased (figure 5(c)). The maximal elevation of landslide initiation was found near the upper tree line, at about 650 m a.s.l. Landslides also exponentially increased with an increase in slope steepness (figure 5(b)). In figure 5(d) slope steepness is presented as sin(α), where α is the maximum slope angle along the extent of the landslide scar. There is a weak (but significant) positive correlation between landslide length and sin(α) (figure 5(d)).

Landslides and soil water content anomalies
Temporal dynamics of EWTA coincided with landslide occurrence (figure 6(a)). A significant (r 2 =0.48) correlation between landslides occurrence and June-August EWTA was observed (table 3, figure 6(b)). Note that a significant temporal increase of EWTA dispersion (i.e., soil water anomalies variation) extremes is observed (figure 6(a)). Drought index SPEI and soil water anomalies EWTA are correlated (figure 6(c)). Table 3 lists the correlation coefficients of climate variables and landslide occurrence. Landslide occurrence was found to be significantly correlated with annual mean temperatures (r=0.51), but not annual precipitation. A higher correlation was observed with the N t>0 , duration of thawing period (i.e., number of

Landslides and relief features
The spatial pattern of landslides is uneven with respect to azimuth as all landslides occurred on slopes with the highest insolation and warmest temperatures, i.e. south and south-west facing slopes. Not a single landslide was found on the shadowed northern exposures. That effect is likely related to permafrost active layer depth. The latter varies from several cm to 1.0 m depending on exposure (Kharuk et al 2008). Landslide occurrence is also strongly dependent on slope steepness. The number of landslides increased exponentially with increases in slope steepness ( figure 5(c)). Landslide lengths varied within a wide range-from short (50 m) to very long (>400 m) with a mean value about 170 m. The majority of landslides began at elevations between 200 and 250 m, with number decreasing exponentially with elevation increase ( figure 4(b)). The maximum elevation of landslide headscarps approximates that of the upper tree line (which within the study area is about 650 m a.s.l.). The presence of trees may promote landslide activation, because (1) the weight of trees provides a downslope driving force, and (2) tree roots help bind together the active layer. It is known that larch roots exist partly within the frozen soil horizon even during summer (Abaimov et al 2002). This occurs because (1) in anomalously warm years roots penetrate to deeper soil horizons and then are frozen in cold years, and (2) the active layer decreases from the moment of tree establishment. The latter is caused due to moss and lichen ground cover that acts as a thermal insulator (Kharuk et al 2008). Warming causes the active layer to increase, which releases the roots from the frozen soil. These, together with increased soil water content leads to the soil layer sliding over the permafrost while precipitation increases. The estimated landslide hazard area is about 30% of total area.

Landslides, climate variables and soil water anomalies
Landslides are significantly and positively correlated with July-August drought index and June-August soil water anomalies; thus, the probability of landslides increases as soil water anomalies increase (figure 6(b)). Globally, landslides were reported most frequently from July to September (Kirschbaum et al 2015), which coincided with our data with the exception of September. Meanwhile, there is a trend of increased variance of soil water anomalies ( figure 6(a)). The latter indicates that the probability of a landslide occurring will increase even given a stable mean precipitation level. Landslides occurrence are significantly correlated with August precipitation (table 3, figure 6(a)). Notably, August precipitation increased during the last decade whereas the annual precipitation decreased ( figure 3(b)). Moreover, precipitation itself increased active layer thickness due to high heat capacity of water (about four times in comparison with air). Along with precipitation, active layer thickness is significant for landslide occurrence. The active layer captures rainfall, and when pore water pressure is sufficient to reduce normal friction to a critical level, landslides can occur. The deeper active layer provides a higher rainfall trapping, increasing active layer weight over permafrost. The increasing probability of landslides triggering obeys Newton's Second Law (Iverson 2000). It is known that due to shallow active layer and no underlying permafrost permeability the majority of rainfall goes directly to the rivers. Along with rainfall, water seepage from thawing permafrost also increases the landslide probability. Significantly, landslides occurrence correlated positively not only with precipitation, but also with SPEI drought index ( figure 7(b)). Thus, with a decrease in drought conditions landslides occurrence also increased.
No correlation was found between landslides occurrence and summer air temperatures. Meanwhile, landslides occurrence was significantly correlated (r 2 =0.67) with the number of days with t max >0°C (N t>0 ) during the May-August period. Thus, the annual period of warming is a significant determinant of landslides occurrence. The main variability of N t>0 was observed in May (r 2 =0.55) table 3, i.e. N t>0 increase or decrease occurred during May depending on the year. Similarly, landslides occurrence was correlated with the sum of positive temperatures during May (r 2 =0.52). This coincides with the general trend of climate warming, i.e. earlier snow melting.
A weak (and significant) correlation was observed with annual temperatures (table 3). That correlation may be a consequence of 'permafrost warming'. When permafrost temperature is increasing there is a decrease in shear and normal stresses of frozen ground due to less ice-bonding (there is particularly strong decrease from −3°C to 0°C (Streletskiy et al 2012). Although there are no data on the permafrost temperature increase within the study area, Romanovsky et al (2010) showed an increase of permafrost temperatures of 0.3°C-2°C in Siberia during the last four decades.
Thus, the main cause of observed increase in landslide occurrence is an increase of precipitation and soil water anomalies. Climate scenarios forecast an increase of air temperature in the Arctic from 7°C to 11°C by the end of the 21st century (Sillmann et al 2013, Vaks et al 2013. This warming may lead to increases in the permafrost active layer thickness and ultimately more landslides. The data obtained supports the hypothesis that landslides occurrence will be more frequent with warming and an increase in precipitation (Montrasio and Valentino 2008).
Along with heavy rainfall and permafrost thawing, forest fires can trigger landslides too. There is evidence of warming-induced higher fire frequency within the Siberian permafrost zone (Kharuk et al 2011). According to predictions, a future increase of forest fires within the boreal zone is expected (e.g., Flannigan et al 1998), which in combination with permafrost thawing should increase the occurrence of landslides.

Post-landslides vegetation growth
Landslides create patches of disturbed soil that are the initiation of succession of forest species, including potential establishment of new species into the larch habitat. This problem has been addressed in only a few papers (e.g, Abaimov et al 2002). In particular, landslides scars present opportunities for establishment of less cold-tolerant species into larchdominated forests. There is evidence of Siberian pine (Pinus sibirica) and fir (Abies sibirica) migration into larch-dominated communities (Kharuk et al 2005). In general post-landslide vegetation growth in permafrost is poorly understood and needs more investigation.

Conclusion
Based on the analysis of high-resolution satellite images and climate data for the period from 2000 to 2012, landslides have increased within the study area, an area of continuous permafrost in central Siberia. This phenomenon correlates with August precipitation, drought decrease, and soil water anomalies. The main cause of the observed landslide occurrence increase is an increase of precipitation and soil water anomaly extremes. Landslide occurrence is strongly dependent on relief features and were found in the study to be located on southward facing slopes only on steeper slopes. The area studied represents the vast larch forests of the central Siberian Plateau. We will expand our analysis to other parts of the Arctic forest and tundra to more fully understand the impacts of landslide dynamics on Arctic ecosystems and carbon balance.