The study was carried out first through the analysis and interpretation of remote sensing data, such as Google Earth satellite optical images (CNES/Airbus from 2011 to 2017), aerial photos of Geographical Mapping Organization of Iran (acquired in 1955 and 2002, scale 1:55,000 and 1:40,000), and vector topographic maps (National Cartographic Center of Iran, topographic map of Kuhdasht, scale: 1:25,000) and by a 10-m digital elevation model (DEM) obtained using ArcGIS 10^® software package with ANUDEM interpolation algorithm (Hutchinson et al., 2011). The remote study led to the first detection of the landforms used in the statistical analysis (block clusters, ridges, and gullies), and to a preliminary size measurement of the blocks outcropping in the landslide debris.

Field surveys were performed with the purpose of i) cross-checking the landforms previously mapped during the remote study; ii) defining the basal surface of the landslide deposit identifying bedrock outcrops, iii) identifying the location of springs, iv) sampling the landslide debris to perform permeability lab tests, and v) measuring the density of the shallow covers by volumetric techniques.

For validating the statistical model, permeability analysis of the fine portion of landslide debris sediments was performed.

To carry out the permeability test in the laboratory, in a first step, a total number of 63 debris matrix samples were taken to classify them according to the Unified Soil Classification System (USCS) standard (ASTM, 2017a). The outcrop of each sampling point was cleaned to a depth of 5 cm and then 1 kg samples were taken.

With this purpose, the grain size analysis was performed, and Atterberg limits were defined. The grain size analysis was carried out following ASTM (2007) recommendations at the Geotechnical Laboratory of Sapienza University. The ASTM 200 sieve passing was characterized with X-ray sedigraphy. The X-ray sedigraph consists of an analyzer unit (SediGraph III by Micromeritics) and a multisampler (MasterTech 52) interfaced by management software. Atterberg limits were measured following ASTM (2017b) recommendations at the Geotechnical Laboratory of Sapienza University.

In order to proceed with the test, it was necessary to measure the in situ density of the debris sediments in order to obtain a reconstituted laboratory sample. Therefore, the rubber balloon method (ASTM, 2015) for measuring the in situ density was used. This test was performed for 28 locations on the landslide debris surface with the values ranging between 1.2 and 1.7 kg/m³ for the finer part of the landslide debris (ASTM 40 sieve passing).

Then, the ML group samples were aggregated into five subgroups based on their relative proximity for performing the permeability tests. To determine the permeability of ASTM 40 sieve passing, various falling head tests were performed using a variable load permeameter, with the laminar flow directed from the bottom to top (ASTM, 2006). It was necessary to reconstitute the samples according to different densities corresponding to different depths for recognizing the variation of permeability with respect to depth. Therefore, assuming a linear relation between stress and density (where the stress at the surface is equal to zero corresponding to the density value of 1.2 kg/m³, increasing up to a maximum of 600 kPa, linearly correlated to the depth and corresponding to the density value of 1.7 kg/m³ and maximum permeable depth).

Moreover, to reconstruct the paleo-topography buried below the debris inferring insights on the emplacement dynamics of the debris, the geological map of the neighboring areas and 13 cross sections (7 along and 6 across the overall run-out direction) were performed. In particular, for the geological sections we referred to the stratigraphic column proposed by James and Wynd (1965) and Alavi (2004) and to the detailed geological mapping of Kabir-Kuh carried out by the National Iranian Oil Company (NIOC, Setudehnia and Perry, 1967; Takin et al., 1970; Macleod 1970). In addition, field surveys were carried out to define the geometry of the boundary between the bedrock and the landslide debris. For this purpose, bedrock outcrops were identified on which the cross sections were then built. Furthermore, to accurately estimate the landslide debris volume, stratigraphic logs were sampled using AutoCAD from the sections along and across the run-out direction, each 350 and 500 m, respectively, obtaining a total of 524 stratigraphic logs.

Based on the remote analysis and on the field cross-check, we mapped the main landforms on the landslide debris (Figure 2). It was possible to point out the evidence of gravity-induced landforms such as blocks, ridges, extensional features, denudational scarps as well as fluvial and runoff landforms such as gullies, badlands, and a suite of three orders of strath terraces including a flood plain and springs. As regards the primary and secondary landforms, 1,129 ridges, 1,162 blocks, and 1,162 gullies were mapped (Table 2). The mean lengths are 256, 205, and 15 m, respectively.

Figure 3 shows some examples from the field of primary and secondary zones deduced on the base of landform associations described in the Definition of Primary and Secondary Landforms.

In Table 2, the density value ranges for each landform are reported. The ridge density range is between 0 and 8 km/km², the gully density range is between 0 and 18.29 km/km², and the block density range is between 0 and 254.64 km/km². The aggregate density range is between 0 and 320.72 km/km², 0–550.18 km/km², and 0 to 4,937 km/km², respectively. Figure 4 shows the map of distributions of the landform density (Figures 4A–C) and the map of probability of primary and secondary zones based on the criteria described in the Definition of Primary and Secondary Landforms (Figure 4D). As can be seen in Figure 4A, the highest density of ridges is concentrated in the middle zone of the debris, while at the perimeter, the lowest density occurs. Figure 4B shows that the density of gullies is higher close to the Seymareh River gorge and in the northern part of the debris, while it tends to lower values in the rest of the debris. In Figure 4C, the highest density of blocks occurs along the perimeter of the debris, especially in the northern part, where they outcrop in a wider area, including the middle part. The lowest density is recorded in the area near the Seymareh River gorge. Based on the density distribution of each landform type, we reclassified the values into four classes according to the quantile criterion. In detail, the density classes for ridges, gullies, and blocks are reported in Table 2.

Combining the aggregate density zonation (Figures 4A–C) according to the aforementioned conditions, it was possible to compute the probability model of primary and secondary zones (Figure 4D). As reported in Figure 4D, the probability of primary zones is located along the perimeter and in the northern part where the highest density values of blocks and ridges occur. The probability of secondary zones is recorded along the Seymareh River gorge, in correspondence with the highest density of gullies.

The results of USCS classification are shown in Figure 5, while the in situ density results are reported in Figure 6. Based on the obtained results shown in Figures 5, 6, groups of matrix types can be distinguished into high- and low-compressibility silt (MH and ML), clay of low plasticity (CL), sandy silty clay (CL-ML), clayey sand (SC), and silty sand (SM). Fine-grained soils (ML, MH, CL, and CL-ML) are mainly located in the Asmari Formation debris area along the middle and distal portions of the landslide debris, while the coarse-grained (SC and SM) soils are mainly distributed within the debris in the proximity of the detachment area of the Seymareh landslide which is the Pabdeh-Gurpi Formation debris area. However, the ML group constitutes 68% of identified soil types with 42 samples of 62, and for this reason, we chose to focus the analysis on them. In detail, the 42 samples of the ML group were aggregated into five subgroups based on their relative proximity (Figure 7). Moreover, five samples were discarded because of the high amount of matrix likely due to their location close to the scar area.

Regarding the in situ density analysis, field measurements were performed in 28 sites (Figure 6). The density values range from 1.2 to 1.62 kg/m³ and on average 1.43 kg/m³. However, a density of 1.2 kg/m³ was considered as surface density with (σ = 0 KPa) for the permeability analysis because we assumed it was more representative than measurements on other samples where the presence of rocks with higher density biased the total density result.

Then, for each of the five subgroups, four permeability tests were performed with increasing density values that were related to the depth correlated stresses (σ). The maximum value of density reconstructed in the permeameter was 1.7 kg/m³ for σ equal to 600 kPa, which corresponds to a depth of 35.29 m. From this depth, the landslide debris was assumed to be impermeable. Figure 7 represents the plots of the permeability analysis results of Table 3 as a function of the σ, γ, ρ, and H. Figure 7A shows the relationship between γ and σ used for obtaining the relative values of H (Figure 7B). Figures 7C–E report K as a function of γ, σ, and H for each permeability subgroup. In detail, the obtained K values range from 1.25 × 10⁻⁷ to 3.27 × 10⁻⁴ cm/s, from 3.37 × 10⁻⁶ to 4.52 × 10⁻⁴ cm/s, from 4.75 × 10⁻⁶ to 1.56 × 10⁻³ cm/s, from 1.25 × 10⁻⁵ to 2.3 × 10⁻³ cm/s, and from 2.40 × 10⁻⁵ to 6 × 10⁻³ cm/s, respectively, for subgroups 1, 2, 3, 4, and 5. The graph shows a similar behavior in relation to γ, σ, and H among subgroups 1, 2, and between subgroups 4 and 5. Indeed, although K decreases with γ, σ, and H for all the samples, subgroups 1, 2, and 3 intercept the K axis with a lower value in comparison with subgroups 4 and 5.

To reconstruct the paleo-topography beneath the landslide debris, 13 geological cross sections (7 along and 6 across the run-out direction) were built. Among them, two geological cross sections, whose traces are drawn in the map of Figure 1, are reported in Figure 8. The other geological sections and the details related to the landslide debris volume computation are reported in the Supplementary Material.

As can be observed in section A of Figure 8, oriented SW-NE, the topography beneath the landslide basal surface does not follow the anticline–syncline structure. Indeed, the topographic depocenter, as the relative lowest elevation point along the section, is not located in correspondence of the synclinal axis. Instead, it lies right next to the current course of the river, 2 km toward SW with respect to the same axis. Moreover, at the core of the syncline, there is the Agha Jari Formation with thickness greater than 500 m. From section B of Figure 8, oriented NW-SE, the Agha Jari Formation outcrops in the most NW part of the profile, reaching a thickness of almost 500 m.

The mechanisms of paleo-landslides are generally difficult to be studied and interpreted since there are not enough data available from the period before the landslide event (e.g., Shoaei 2014; Pánek and Klimeš, 2016) and on damming landslide debris, original and reshaped surface landforms can outcrop extensively (Dewitte and Demoulin, 2005; Magnarini et al., 2019; Pánek et al., 2020; Magnarini et al., 2021). Original primary landforms are linked to the landslide emplacement and kinematics, while secondary reshaped landforms are associated with the run-off and fluvial erosion (Crozier, 2010). Different types of surface features can be recognized on the landslide debris, such as ridges, gullies, denudation scarps, and blocks (Parise, 2003; Rouhi et al., 2019; Pánek et al., 2020). The spatial associations of ridges, gullies, and blocks on the Seymareh landslide debris offer the opportunity to study the relationship between these features and the processes that formed them. Through the statistical model proposed in this research, it was possible to distinguish primary (original) and secondary (reshaped) zones on the landslide debris (Figure 9). By analyzing the primary landforms (ridges and blocks) on the original zones of the landslide debris (Figure 9), it was possible to understand the landslide emplacement kinematics, as they provide insights for the zonation of debris in sectors characterized by different stress configurations. So far, little attention has been given to the Seymareh landslide debris as a possible indicator of its emplacement kinematics after the slope collapse (Roberts and Evans, 2013; Shoaei, 2014). To this aim, the ridge curvature and azimuthal direction were measured. As shown in Figure 9A, kinematic arrows, deduced from the ridge orientation analysis only in the primary zones, show two main orientations, NE-SW and NW-SE, which agree with the proposed emplacement kinematics of the Seymareh landslide debris by Shoaei (2014). The block distribution also revealed that a back and forward motion involved the landslide debris during its emplacement, as the huge blocks are distributed not only along the front of the landslide mass but also in its middle zone, where the Chenareh-Kuh ridge is expected to be buried. This evidence allows us to confirm what has already been stated in the previous studies (Harrison and Falcon, 1937; Watson and Wright, 1969; Roberts, 2008; Roberts and Evans, 2013; Shoaei, 2014) that a part of the Chenareh-Kuh and Halush-Kuh anticlines in front of Kabir-Kuh are buried under the landslide debris. It means that during landslide debris emplacement a part of the debris stopped against the topographic obstacles sweeping back into the Seymareh paleo-valley, and another part with higher energy overpassed these two anticlines.

The recognition of the ridges on the surface of the Seymareh landslide debris allowed the identification of three distinct directions of the rock avalanche debris propagation during its emplacement (Figure 9B). A possible correlation among the different movement directions and zones of the landslide debris related to a specific regime of mass emplacement introduced in the literature can also be identified (Dufresne et al., 2010; Dufresne, 2012; Longchamp et al., 2016). The first zone (ZA) is in the middle area of the debris and reflects the extension regime during mass emplacement. The ridges in this zone have longer distances from each other and are perpendicular to the run-out direction. The second zone (ZB) of the landslide mass testifies the compression affecting this area during debris emplacement. The ridges in this zone are close to each other with a tendency to become parallel to the run-out direction. The third zone (ZC) is located in the distal margins of the landslide debris which reflects shearing and has ridges both parallel and perpendicular to the run-out direction. The results show that the largest ridges are in the medial (and eastern) zone, changing into smaller ridges in the distal part (Dufresne et al., 2010).

The morphology of the paleo-valley of the Seymareh landslide, reconstructed based on the field evidence collected for this research, allows locating the topographic depocenter of the valley close to the actual position of the Seymareh River course. Moreover, the distribution of the kinematic indicators extracted from ridges direction and the block location is consistent with the position of the depocenter of the Seymareh valley before the landslide emplacement. This paleo-morphology indicates that the paleo-Seymareh River was flowing to the NE of the present river gorge which agrees with Roostai et al. (2018). In detail, the confluence between Seymareh and Kashkan occurred at the toe of the collapsed slope, as suggested also by Shoaei (2014). Furthermore, the estimated volume of the landslide debris (43.97 km³, details are reported in the Supplementary Material) calculated by the reconstruction of paleo-valley morphology is in good accordance with the estimated volume of the debris by Roberts and Evans (2013), which shows the high accuracy of reconstructed paleo-morphology of the Seymareh landslide debris area in the current research.

As it can be seen in Figure 10, the landslide debris is divided into three zones with different permeability values which fit very well with the primary and secondary zonation. It also confirms the deduced kinematics of the debris emplacement. The NE zone where the debris overpassed Chenareh-Kuh and swept back to the paleo-valley has a higher permeability value with respect to the zone where the debris emplaced in the NW direction. Moreover, we can interpret the dynamics of the lake overflow. Through the statistical approach proposed in this research, it was possible to recognize the zones with reshaped surface features. As an example, the results show that the statistical approach correctly recognized the secondary zones on the landslide debris corresponding to the current Seymareh River location. By comparing the direction of the gullies in the secondary zones with the results of the permeability test, it was possible to validate the statistical model and provide a better interpretation to the movement process which terminated the landslide dam cutting. Indeed, the distribution of gullies in the SW part of the debris mainly has a direction between 0 and 90°. In particular, that part corresponds to the area of the landslide debris with the Pabdeh–Gurpi Formation; while gullies on the left side of the river are between 90–180° which are perpendicular to the other side of the river. This area is the boundary between Pabdeh–Gurpi and Asmari Formations debris with several springs distributed along the contact between the landslide debris and its basement. Based on the results of permeability laboratory tests, a non-negligible change in the debris permeability can be observed within the landslide debris and related to the upper part of the debris in the Asmari Formation area. There is another change of the permeability related to the intersection of the Asmari Formation and Pabdeh Gurpi Formation. Moreover, it can be deduced that the surface erosion of the debris by gullies in combination with the change of the permeability and springs distribution along with the contact between the landslide debris basement from basement led to the cutting of the Seymareh landslide dam.

Different morphological parameters control the emplacement mechanism, kinematics, and behavior of landslide dams (Ermini and Casagli, 2003; Tacconi Stefanelli et al., 2016; Rouhi et al., 2019). The parameters describing the geomorphology of the valley such as valley width, slope angle, and the height of the release area exert control on the emplacement mechanism and kinematics of the landslide deposit (Ermini and Casagli, 2003; Tacconi Stefanelli et al., 2016; Rouhi et al., 2019), and indirectly control the likelihood of landslide dam formation (Nicoletti and Sorriso-Valvo, 1991).

Landslide dams can be characterized by morphological parameters such as the size and shape of the depositional feature (Tacconi Stefanelli et al., 2015). In the case of cut landslide dams, the size, shape, and distribution of features on their surface help to forecast the cutting process. Some authors, by extracting the geomorphic parameters of landslide dams, conducted inventories for some countries in the world. The interpretation of landslide dam failure is a complex topic because of the numerous variables involving the hydrodynamic interference between landslides, rivers, and dam creation at the same time.

Numerous attempts have been made to investigate the complex behavior of the landslide dams. Some authors using landslide dams’ inventories, proposed some geomorphic indexes (Swanson et al., 1988; Ermini and Casagli, 2003; Coico et al., 2013; Dal Sasso et al., 2014) to predict the formation and evolution of landslide dams. These indexes result from the composition of two or more morphological attributes that characterize the landslide (e.g., landslide volume or length) and the involved river valley (e.g., valley width). They proved that these indexes could be a useful tool for making accurate predictions concerning the stability of landslide dams. Nash (2003) took into consideration the block size and proposed Modified Dimensionless Blockage Index (MDBI) and showed that the average block size (D50) of the dam material exerts significant influence on dam longevity and evolution. Recently, Tacconi Stefanelli et al. (2016) to overcome some limitations of other indexes, introduced two new ones: The Morphological Obstruction Index (MOI) and the Hydromorphological Dam Stability Index (HDSI) and analyzed about 300 landslide dam events extended all over Italy. The newly proposed indexes (Tacconi Stefanelli et al., 2016) show an improvement in the forecasting effectiveness and have the advantage of being based on morphometric input parameters that can be easily and quickly assessed in a distributed way even over large areas. These indexes consider mainly geomorphic variables characterizing both the landslide (e.g., landslide volume or length) and the river channel (e.g., valley width) to understand the key parameters controlling the stability and evolution of landslide dams, and to use them for providing a better modified geomorphic index. Generally, the stability and evolution of a landslide dam should be assessed based on numerical simulations or an experimental model (Chen et al., 2011). In most cases, numerous accurate parameters and materials are not easily accessible; therefore, for an immediate evaluation of landslide dams’ stability, geomorphic valley parameters and morphometric dam parameters should be analyzed (Jakobsen, 2016). A comprehensive formula or approach for assessing the behavior of landslide dams should consider all parameters related to the landslide, landslide dam, and river.

Furthermore, the analysis of morpho-structural features related to the debris of rockslides and rock avalanches provides new insights into their emplacement kinematics and failure mechanism (Dufresne and Davies, 2009; De Blasio, 2014; Esposito et al., 2014; Longchamp et al., 2016; Wang et al., 2018; Rouhi et al., 2019). Although the scientific community has focused heavily on landslide dams in various parts of the world, their evolution and behavior have not been completely understood. In this regard, the new statistical approach presented in this study represents a new tool for the recognition and distinction of original and reshaped regions, allowing to investigating the emplacement kinematics of landslide debris.