*2.2. Landslide Surveying*

Landslide events were collected from di fferent sources. The CORONA KH-4B image (Mission ID: 1110-1089Fore) of 1970 (~2 m) (https://corona.cast.uark.edu/atlas) and aerial photos (1:20,000) from 1966 were used to survey old landslides (~230 samples) mostly in the protected forest. The new landslides (~430 samples) were obtained from field observations, the available database [83], and high-resolution images of Google Earth for 2016 (Figure 1). In addition to these samples, about 210 and 1650 polygons of old and new landslides, which were mapped using Sentinel images by Shirvani et al. [72], were used along with the landslide samples for the two study areas (Figure 1). The size of the landslide samples ranged between 0.095 and 239.6 ha in the protected forest and between 0.018 and 60.82 ha in the non-protected forest.

**Figure 1.** Location of the study areas in the Hyrcanian ecoregion in NE Iran: the Golestan National Park as a "protected forest" (**a**); and disturbed forests by mining, logging, and road building as a "non-protected forest" (**b**). Landslide events I were collected from different resources in the current research and the landslide events II were adopted from Shirvani et al. [72].

#### *2.3. Image Segmentation for Generating Objects*

A group of pixels that have similarity in their spectral and spatial properties is defined as an object [84]. The object-based paradigm has the ability to derive all possible features from the spectral, geometrical, contextual, and textural properties of either satellite images [84] or GIS (Geospatial Information Systems)-based data [85]. To obtain homogeneous objects, multi-resolution segmentation was implemented on the spectral bands of Landsat 8 of 2016 and SRTM (Shuttle Radar Topography Mission)-derived digital elevation model (30 m) [86]. After testing different scales by trial and error, the scale of 150 was selected with a higher weight for the near-infrared band (NIR) and the compactness and shape values of 0.8 and 0.2, respectively. The final segmentation was optimized using some digital elevation model (DEM) derivatives [72] such as slope, hillshade, terrain ruggedness index (TRI), and flow direction river (FDR). These image-segmented objects were assigned to calculate anthropogenic-induced deforestation and forest fragmentation from remote sensing data (Figure 2). Moreover, the summary statistics of conditioning and other triggering factors were calculated within each object segment.

#### *2.4. Conditioning and Triggering Factors*

We divided the driving forces of landslide susceptibility into conditioning and triggering factors based on the previous studies in the literature review [51,87] and the landslide characteristics of the study areas. Conditioning factors are cumulative events that show the potential of landslide occurrence, but do not necessarily trigger landslides [88], while triggering factors activate the landslides and increase the probability of occurrence by disturbing the balance between driving and resisting forces [89]. The selected conditioning variables included topographic, hydrologic, geology, soil, and vegetation layers. We generated topographic layers from the DEM; the hydrologic layers from the DEM and digital topographic maps (1:25,000); geological maps (1:100,000) from the Geological Survey and Mineral Explorations of Iran (GSI) and soil maps (1:100,000) from the Agriculture Department of Iran; and forest types were based on the thematic maps of forest managemen<sup>t</sup> plans (Table 1). Furthermore, triggering factors were classified into natural (i.e., rainfall, earthquake, and flood) and anthropogenic (i.e., forest fragmentation, forest loss, logging, and mining). Rainfall intensity and earthquake magnitude were created by kriging [90] and inverse distance weighted (IDW) interpolation methods [91] from the historical events (Table 2). Flood frequency was calculated based on the frequency of flood occurrence along the main rivers within a specific catchment during the last five decades (Table 2).

We created the forest layer of 1970 from CORONA images and aerial photos using object-based nearest neighbor classification (Figure 2). A multispectral image was created from the original band of the images and two created channels using the sharpening and embossing filters. The segmentation was implemented through the multispectral resolution algorithm on these images. The possible classes were defined and some objects from each class were selected as training samples. Varieties of ancillary spectral and textural features were created to improve the accuracy of the classification. After optimizing the dimensions of the features, the standard nearest neighbor algorithm was applied for classifying the images to the forest and non-forest classes [82]. Moreover, the forest layer of 2016 was mapped from Landsat 8 using the object rule-based classification (Figure 2). After image segmentation as described earlier in Section 2.3, some objects from each class were selected. Then, di fferent spectral, contextual, and textural features were derived from the main spectral bands of Landsat 8. The thresholds of non-forest classes from the forest class were determined by matching their values within the derived object features. The objects that had maximum di fference values less than 0.52 were classified as water, residential areas, and grasslands; objects that had green vegetation index (GVI) [92] values that were negative and standard deviations derived from the gray-level co-occurrence matrix (GLCM) less than 21 were classified as dry-farming; the objects that had enhanced vegetation index2 [93] values higher than 0.98 and entropy derived from the GLCM less than 7 were classified as irrigated-farming; and the remaining unclassified objects were classified as forest [82]. The accuracy of the classifications was validated using the provided ground truth samples and confusion matrix [94], as shown in Table A1.

We calculated a number of forest loss and forest fragmentation metrics by comparing these two forest layers such as the rate of forest loss [65], edge density (ED), mean patch size (MPS), mean shape index (MSI), mean patch edge (MPE), mean perimeter-area ratio (MPAR), and number of patches (NumP) within each object [95] (Table 2). These metrics were included as triggering of landslide susceptibility along with the other variables. The total volumes of logging and weights of mining materials were calculated within each object from 1970 to 2016 as indicators of logging and mining intensity in the two studied forests.

**Figure 2.** The process of discriminating forests from non-forests using object-based nearest neighbor and rule-based classifications by utilizing CORONA images/aerial photos of 1970 and Landsat 8 of 2016 over the protected and non-protected forests in NE, Iran (after Shirvani et al. [82]). Abbreviations: stdDev.: Standard deviation; Max. diff.: Maximum difference; GVI: Green Vegetation Index; EVI2: Enhanced vegetation index2; GLCM: Gray-level co-occurrence matrix.
