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

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² (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⁻¹) (Ershov 1989). Trees within the solifluction zone deviate from the vertical direction and form so called 'drunken forest'.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 within the study area is strongly continental with long cold winters and short warm summers (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.

Zoom In Zoom Out Reset image size

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): D_i = P_i − PET_i.

PET (mm) is obtained by:

$$\begin{eqnarray*}&&{\rm{PET}}=16\times K\times {(10\times T\times {I}^{-1})}^{m},\end{eqnarray*}$$

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² 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)).

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² (with total analyzed area about 17 500 km²). 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⁻¹ (Riegger et al 2012, Long et al 2014). GRACE-derived EWTA values are caused by both water anomalies in the soil and in rivers and ponds. Because we used a watershed approach in this study, there was no incoming water flow from outside of the watershed. Because of that, EWTA is proportional to the soil water content anomalies. GRACE data are available since launch and activation in 2003 (http://grace.jpl.nasa.gov; Swenson and Wahr 2006, Landerer and Swenson 2012). GRACE spatial resolution was 1° × 1° degree (∼112 × 44 km² at latitude 66°). The data were processed using ERDAS Imagine software (http://geospatial.intergraph.com) and ESRI ArcGIS software (http://www.esri.com).

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:

$$\begin{eqnarray}&&{K}_{c(i)}=({A}_{c(i)f}/{A}_{c(i)I})*100/\displaystyle {\sum }_{i=a}^{b}({A}_{c(i)f}/{A}_{c(i)I}),\end{eqnarray} \tag{ 1 }$$

where K_c(i) is the coefficient of normalization, c(i) is ith category of landscape feature c, A_c(i)f is the area of the 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).

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.

Figure 4(a) indicates that the number of landslides increased since 2006 with 80% occurring after 2006. The majority (60%) of landslides had lengths within the range of 75–225 m (figure 4(b)). The landslide statistics including slope azimuth (aspect) and slope steepness are presented in table 2.

Zoom In Zoom Out Reset image size

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)).

Zoom In Zoom Out Reset image size

Temporal dynamics of EWTA coincided with landslide occurrence (figure 6(a)). A significant (r² = 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)).

Zoom In Zoom Out Reset image size

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 days with t_max > 0 °C) (r = 0.67). The highest correlation was observed with August precipitation (r = 0.81). When averaged over the MJJA time period the correlation with precipitation was lower (r = 0.58). Correlations between SPEI drought index and landside occurrence were significant for July and August (table 3). Figures 7(a) and (b) shows observed positive linear relationships for August precipitation (r² = 0.66) and SPEI (r² = 0.48), respectively.

Zoom In Zoom Out Reset image size

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 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² = 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² = 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² = 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.

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 larch-dominated 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.