2.1.1 Landslide catalog

Landslide events used in this study come from two sources: the Global Landslide Catalog (GLC) (Kirschbaum et al., 2010, 2015) and the Global Fatal Landslide Database (GFLD) (Froude and Petley, 2018). The GLC has been compiled by NASA since 2007 and contains all types of mass movements triggered mostly by rainfall. The sources of the GLC are mainly media reports, disaster databases, and scientific reports. The GFLD only includes landslide events that caused fatalities and is obtained from media reports. It currently covers the period from 2004 to 2017. These two landslide inventories were chosen because, to the best of our knowledge, they are the only ones with the exact landslide dates available for the study region.

Figure 1Landslide events from 2004 to 2018 extracted from the GLC (white points) and the GFLD (black points). Background contour is topography from digital elevation model (DEM) data from Shuttle Radar Topographic Mission (SRTM).

We selected landslide events triggered by atmospheric factors in Kyrgyzstan and Tajikistan from 2007 to 2018 from the GLC and 2004 to 2017 from the GFLD. Then we merged these two data sets and deleted duplicate events that occurred on the same day and came from the same source link, resulting in 96 landslide events for Kyrgyzstan and Tajikistan from 2004 to 2018 (Fig. 1).

2.1.2 Atmospheric data

Rainfall and snowmelt data are extracted from the High Asia Refined analysis version 2 (HAR v2). The HAR v2 is a newly developed regional atmospheric data set. It was generated by dynamical downscaling of the ERA5 reanalysis data using the Weather Research and Forecasting (WRF) model. It provides atmospheric data with high resolution and accuracy over High Mountain Asia (Hamm et al., 2020; Wang et al., 2021). Detailed modeling strategies of the HAR v2 are described in Wang et al. (2021). The HAR v2 has a grid spacing of 10 km and is available in hourly, daily, monthly, and yearly aggregations. Daily products were used in this study to determine the climatic trigger of each landslide event (Sect. 2.2.1) and to define thresholds for landslide triggering (Sect. 2.2.2). Rainfall was calculated as the difference between total precipitation and snowfall. Snowmelt is not a standard output of the WRF and was calculated using the surface energy balance (SEB). The SEB in the HAR v2 is resolved by the Noah land surface model (LSM) (Tewari et al., 2004):

where R_n, H_s, H_l, and H_g are net radiation, sensible heat flux, latent heat flux, and ground heat flux (in W m⁻²), respectively. These four variables are directly available in the HAR v2. H_m is the heat flux for melting and refreezing (in W m⁻²). H_m>0 indicates melting process, while H_m<0 refers to refreezing process. When H_m>0, snowmelt h_m (kg m⁻² s⁻¹) is calculated as

where λ_m is the latent heat of fusion. When the calculated h_m is greater than snow water equivalent, then h_m is set to be equal to snow water equivalent.

2.2.1 Determine the atmospheric trigger of landslide events

The atmospheric trigger of a landslide event is determined by the co-occurrence of the landslide event with rainfall and snowmelt event. If a landslide event only occurred within or 1 d after a rainfall (snowmelt) event, then this landslide event is defined as rainfall (snowmelt) triggered. If there are both a rainfall event and a snowmelt event on the day or 1 d before the landslide occurrence day, then the atmospheric trigger of this landslide event is mixed.

To define a rainfall (snowmelt) event, the daily time series of rainfall (snowmelt) were extracted from the grid cells where landslides occurred. For each time series, an independent rainfall (snowmelt) event is defined as a series of consecutive days in which more than 0.2 mm d⁻¹ of rainfall (snowmelt) is simulated. The value of 0.2 mm d⁻¹ is chosen because it is the traditional precision of daily precipitation measurement (Jarraud, 2008) and can be applied to separate dry and wet conditions (Rodwell et al., 2010).

2.2.2 Threshold model for atmospheric triggers

The threshold model developed in this study contains three steps: (1) define LTEs and NLTEs; (2) define the thresholds for rainfall, snowmelt, and the sum of rainfall and snowmelt (hereafter referred to as rainfall + snowmelt) based on maximizing the predictive performance using 2×2 contingency tables; and (3) validate and assess the uncertainties of the defined thresholds. The methods for the first two steps were adopted from Leonarduzzi et al. (2017). Only the landslide events for which the atmospheric triggers could be determined were used for threshold modeling.

The first step is to define LTEs and NLTEs for rainfall, snowmelt, and rainfall + snowmelt. Here, we take rainfall as an example to describe the procedure. First, the method used in Sect. 2.2.1 is applied to define rainfall events for each time series extracted from grid cells where landslides occurred. Next, if a landslide event occurred during or 1 d after a rainfall event, then this rainfall event is classified as a landslide-triggering event. Given the uncertainty in timestamps of landslide events, the day after is also considered as a temporal relaxation. Otherwise, if a rainfall event is not associated with any landslide events, it is classified as a non-landslide-triggering event. For each rainfall event, we calculated three event properties: mean intensity I_mean, maximum intensity I_max, and the accumulated amount of rainfall for the entire event Q. For triggering events, we also calculated these three properties by only considering the period up to the day of the landslide occurrence (hereafter referred to as UTL, meaning up-to-landslide). Note that not all the landslide events co-occurred with a rainfall event. For these events, we set I_mean, I_max, and Q to zero. The same procedure for defining LTEs and NLTEs was conducted for snowmelt and rainfall + snowmelt as well.

The second step is to define thresholds of rainfall, snowmelt, and rainfall + snowmelt for entire events and UTL events using I_mean, I_max, and Q. No single threshold can perfectly separate LTEs from NLTEs since their distributions overlap. We applied 2×2 contingency tables to select the threshold that yields the best predictive performance. Using a certain threshold as a binary classifier, LTEs and NLTEs were categorized into true positive (TP), true negative (TN), false positive (FP), and false negative (FN). The Peirce skill score (PSS) (Hanssen and Kuipers, 1965) was applied as the measure of the predictive performance because it is trail-independent, which means it is unbiased even when the numbers of LTEs and NLTEs are not equally presented (Woodcock, 1976). The PSS is also known as the Hanssen–Kuiper skill score and the true skill statistic. It is calculated as the difference between hit rate (HR) and false alarm rate (FAR):

We chose the threshold that maximizes the PSS. We also computed the Euclidean distance (d) to the optimal point (HR = 1, FAR = 0), which is another commonly used skill score in this application (e.g., Gariano et al., 2015; Piciullo et al., 2017; Postance et al., 2018; Zhuo et al., 2019). Additionally, the receiver operating characteristic (ROC) curve was used to determine the general predictive power of a certain predictor by calculating the area under the ROC curve (AUC) (Fawcett, 2006).

The last step is to validate the threshold model and assess uncertainty. For the calibration of thresholds, all landslide event samples were utilized, and corresponding statistic measures were calculated; i.e., the threshold model was trained and tested on the same data set. To test the model's predictive ability on an unseen data set, we performed k-fold cross-validation. Landslide events were randomly split into k folds with k=8. Then for each unique fold, the fold was taken as the testing set, and the remaining k−1 folds were taken as the training set. Mean values of thresholds, the corresponding statistic measures, and their uncertainties represented by standard deviations were reported.

2.2.3 Mean annual exceedance

Mean annual exceedance (${\stackrel{\mathrm{‾}}{N}}_{\mathrm{th}}$) is calculated for each HAR v2 grid cell. It is defined as the number of events that exceed a certain threshold over a certain period (N_th) divided by the total number of years (N_a):

The unit of ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{th}}$ is the number of events per year. Mean annual exceedance transforms weather-scale triggering conditions to climate-scale disposition. It depicts where landslides are likely to occur from the climatic aspect.

4.1.1 Uncertainty of landslide inventories

Uncertainties of the GLC and the GFLD are comprehensively discussed in Kirschbaum et al. (2010), Kirschbaum et al. (2015), and Froude and Petley (2018). The first major problem of these two data sets is that they underestimate the total number of landslides. This is because these two data sets' primary sources are media reports, which are biased towards events with human casualties (Carrara et al., 2003). The second issue is that the spatial distribution of landslides is biased towards populated areas. In our study area, landslide events also tend to cluster in areas with high population density, e.g., the eastern rim of the Fergana Basin and the Tajik Basin. Landslide number over remote areas is much more likely to be underreported. In addition, there is large uncertainty in landslide location because most media reports do not contain the exact location where landslides were initiated but rather just the name of the village, road, or city affected by landslides. An example in our case is the landslide event in the Issyk-Kul Basin (Fig. 1), the location of which is in a flat area, and the location accuracy provided by the GLC is “exact”. This landslide event's initial zone must be different from the reported location and somewhere nearby with slopes. We also failed to determine the climatic trigger of this landslide event using the HAR v2. Last but not least, landslide timing was also reported with a certain degree of uncertainty. Although it is more typical that a landslide was reported after its actual occurrence (positive errors), negative errors are also possible depending on the interpretation of historical landslide information by an analyst (Peres et al., 2018). Our results show that using the entire weather event period leads to a better predictive performance than just using the UTL period (Table 1). This could be an indication of negative errors in the landslide timing.

Despite these known limitations, the GLC and the GFLD still provide the lower boundary of landslide number and are proven to be valuable in global and regional landslide studies. For example, the GLC has been successfully applied to detect the initiation of rainfall-induced landslides globally (Jia et al., 2020), to investigate the spatiotemporal distribution of potential landslide-triggering factors (Stanley et al., 2020), to explore the synoptic-scale precursors of landslides (Hunt and Dimri, 2021), and to evaluate the Global Landslide Hazard Assessment Model (Kirschbaum and Stanley, 2018). Although the landslide number is known to be incomplete, our results show that they can still present the seasonal distribution of landslide occurrence reasonably well (Fig. 4). This was also concluded by Kirschbaum et al. (2015), who stated that the reason for the unbiased seasonal distribution of landslide occurrence is that the compilation method depends on media alerts, which are consistent throughout the year. Additionally, even though location uncertainty exists, we could determine atmospheric triggers of 91 % of landslide events (87 out of 96). The reason could be that landslide-triggering rainfall and snowmelt events generally have a large spatial extend (Leonarduzzi et al., 2017).

4.1.2 Uncertainty of atmospheric data

Extracting weather data that can represent the exact weather conditions at landslide sites is always a challenge in studies investigating rainfall thresholds for landslide triggering. Rain gauges are the main source of rainfall information (Segoni et al., 2018), and it is very seldom that landslide initial locations are gauged. Due to the highly heterogeneous spatial distribution of precipitation, especially over complex terrain, there exists great uncertainty when rainfall is not directly measured from landslide initial points. Additionally, Marra et al. (2016) found that the initial points of shallow landslides and debris flows generally correspond to the local peak of rainfall. Rain depth decreases with distance, causing an underestimation when rainfall is measured away from the landslide initial point. Traditionally, the nearest gauge is used to represent the weather condition at the landslide site, which sometimes can be kilometers away. Nikolopoulos et al. (2015) examined other more complicated interpolation methods, such as inverse distance weighting and ordinary kriging, and concluded that these methods did not bring any particular added value to the simplest nearest neighbor method.

Using gridded data can avoid this allocation problem (Leonarduzzi et al., 2017). But uncertainties still exist since gridded data only represent the grid-mean value but not the “true” weather conditions at landslide sites. Nevertheless, it is still essential that the gridded data used in our study can accurately represent the grid-mean value. The WRF model configurations of the HAR v2, such as the forcing strategy and physical parameterization schemes, were carefully chosen to ensure its quality (Wang et al., 2021). Several studies (Pritchard et al., 2019; Li et al., 2020) indicate the high accuracy and quality of the old version of the High Asia Refined Analysis (HAR) (Maussion et al., 2014). Wang et al. (2021) compared the performance of the two versions of the HAR against in situ observations from 57 GSOD stations over the High Mountain Asia in terms of daily precipitation and air temperature at 2 m. It was concluded that compared to the old version, HAR v2 generally produces slightly higher precipitation amounts with a mean bias of 0.36 mm d⁻¹. Furthermore, Hamm et al. (2020) compared the HAR v2 with other gridded precipitation data sets at different spatial resolutions, including reanalysis data and satellite-based precipitation retrieval, over a rugged terrain of the central Himalayas and the southwestern Tibetan Plateau. It was concluded that the HAR v2 is the only product that can resolve orographic precipitation, which is a fundamental process over complex terrain. Simulation of air temperature at 2 m in the HAR v2 is better than the old version due to the snow depth correction approach (Wang et al., 2021). Snowmelt in the HAR v2 is resolved by the Noah LSM, which only considers a single layer of snowpack (Koren et al., 1999). Several studies found uncertainty of the Noah LSM in reproducing the snow-related process, e.g., the overestimation of snow albedo (e.g., Chen et al., 2014; Minder et al., 2016; Tomasi et al., 2017). Nevertheless, the snow-related process is the major weakness of LSMs and needs further improvement in the future (Chen et al., 2014).

4.1.3 Impact of spatial resolution of atmospheric data

Previous studies have shown that the spatial resolutions of gridded rainfall data have impacts on identifying landslide-triggering thresholds (Marra et al., 2017; Nikolopoulos et al., 2017). To investigate the influence of spatial resolution of rainfall + snowmelt data on the event properties of landslide-triggering weather events and the triggering thresholds, we resampled the rainfall + snowmelt data from HAR v2 to lower resolutions (20, 30, and 40 km). Then, we repeated the procedure described in Sect. 2.2.2 to determine the event properties of LTE UTL events and their associated thresholds. The results are presented in Fig. 7. There are nine “not detected” events when using the original HAR v2 10 km data (Fig. 3), which means the rainfall + snowmelt amounts at these landslide grid points are near zero (≤0.2 mm d⁻¹) on the day and 1 d before landslide occurrence. By lowering the spatial resolution, more events can be detected. This implies the uncertainty in the reported landslide location since resampling of rainfall + snowmelt encompasses rainfall + snowmelt information from nearby grid points. In general, I_mean and I_max decrease with the increase in grid size, which is in line with the findings of Hamm et al. (2020) that higher-resolved products generally capture more extreme events than coarser products. I_mean and I_max thresholds defined by coarser products are also generally lower. The impact of grid size on Q is the opposite: larger grid size leads to higher Q and threshold value. This is closely associated with the increase in event duration with the increase in grid spacing, resulting from the fact that the resampling process can blend several localized events temporally together. However, lowering the spatial resolution does not lead to worse predictive performance. This, on the one hand, implies again that lower resolution can partly compensate for the uncertainty in landslide locations. On the other hand, it indicates that although landslide initiation itself is a highly localized phenomenon, the weather processes that ensure sufficient water input into the system and that trigger landslides can be clearly identified at the mesoscale (Prenner et al., 2018).

Figure 7Boxplots demonstrating the impact of spatial resolution of atmospheric data on I_mean, I_max, Q, and duration of LTE UTL events, as well as the associated landslide-triggering thresholds (blue stars). The yellow line denotes the median, and the green triangle indicates the mean. Outliers are not shown for a better intercomparison, and n denotes the number of landslide events detected by rainfall + snowmelt.

Download

Based on the above analysis, it can be expected that a convection-permitting-scale (<10 km) downscaling simulation would provide a more realistic representation of weather events that initialized landslides. Compared to such a high-resolution simulation, the HAR v2 10 km data would underestimate the intensity and overestimate the duration of landslide-triggering rainfall + snowmelt events. Moreover, the 10 km resolution of the HAR v2 is not able to explicitly resolve convection processes. Convection-permitting-scale simulations show improvement over simulations applying cumulus parameterization schemes in several aspects, such as more accurate reproduction of the timing of precipitation peaks (Ou et al., 2020; Zhou et al., 2021). However, a finer resolution has a lower tolerance for uncertainty in the landslide location. The potential of a kilometer-scale simulation cannot be realized if the landslide location uncertainty is larger than the grid size. Thus, for our study region, future studies should not focus only on acquiring high-resolution and high-quality atmospheric data but also on developing landslide inventories with higher location accuracy.