Modeling Shallow Landslide Runout Distance in Eocene Flysch Facies Using Empirical–Statistical Models (Western Black Sea Region of Türkiye)

Abstract

Uncertainties related to runout distances in shallow landslide analyses may not only affect lives but may also result in economic losses. Owing to the increase in shallow landslides, which are especially triggered by heavy rainfall, runout distances have been investigated to decipher whether applications of a functional runout distance are feasible. This paper aims to give insights into the modeling of the shallow landslide runout probability in Eocene flysch facies in the Western Black Sea region of Türkiye. There are two main stages in this study—which are dominated by empirical models, the detection of initiation points, and propagation—which help us to understand and visualize the possible runout distances in the study area. Shallow landslide initiation point determination using machine learning has a critical role in the ordered tasks in this study. Modified Holmgren and simplified friction-limited model (SFLM) parameters were applied to provide a good approximation of runout distances during the propagation stage using Flow-R software. The empirical model parameters suggested for debris flows and shallow landslides were investigated comparatively. The runout distance models had approximately the same performance depending on the debris flow and shallow landslide parameters. While the impacted total runout areas for the debris flow parameters were predicted to amount to approximately 146 km², the impacted total runout areas for the shallow landslide parameters were estimated to be about 101 km². Considering the inclusion of the RCP 4.5 and RCP 8.5 precipitation scenarios in the analyses, this also shows that the shallow landslide and debris flow runout distance impact areas will decrease. The investigation of runout distance analyses and the inclusion of the RCP scenarios in the runout analyses are highly intriguing for landslide researchers.

1. Introduction

Considering the increase in extreme rainfalls in the last decade due to climate change, it is possible to claim that shallow landslides are like any other type of natural hazard and that researchers should increasingly focus on them [1]. Runout distance determination in shallow landslide analysis is crucial because hazard evaluations involving runout distance information constitute one of the bases of sustainable development for societies living in mountainous areas [2]. Runout methods are yielding extraordinary insights into landslide runout analysis. For many years, researchers have studied landslides using existing runout methods, which are empirical–statistical methods, physical methods, and numerical models. The empirical–statistical model was applied in studies by Poltnig et al. [3], Kaafarani et al. [4], Mergili et al. [5], Zhou et al. [6], Apriani et al. [7], Di Napoli et al. [8], and Guthrie and Befus [9] to inspect the runout of landslides. A novel terrain-matching targeted machine learning model, which depends on statistical methods, was used to estimate landslide runouts in a study by Ju et al. [10]. Additionally, a data-driven method for the estimation of shallow landslide runout was used in a study by Giarola et al. [11]. Observation of flow-like landslide runouts was accomplished using the physical method, which was applied as a flume test, in studies by Baselt et al. [12] and Gao et al. [13]. The numerical model, which has frequently been used to scrutinize flow-like landslide runouts in recent studies, was also used in studies by Gao et al. [13], Yang et al. [14], Chae et al. [15], Calista et al. [16], Abraham et al. [17], Dash et al. [18], Vicari et al. [19], Mikoš and Bezak [20], Oh et al. [21], Zhou et al. [22], Alene et al. [23], and La Porta et al. [24]. Investigations into the probabilities of failures were based on merged models in some studies. For instance, the combination of the numerical and physical models yielded enhanced comprehension of the failure mechanism in a study by Zhu et al. [25]. The physical method and numerical method were also combined in order to probe the runout distances of debris flows in studies by Clark [26], Melo et al. [27], and Gan and Zhang [28]. In addition, numerical models aiming to calculate the runout of flow-like landslides were successfully verified with a physical model (flume test) in a study by Gao et al. [13]. Numerical and empirical models were combined in order to investigate landslide runouts in a study by Peruzzetto et al. [29]. By considering data-driven approaches, runout analyses can be successfully prepared by applying empirical–statistical, physical, and numerical methods. The mentioned runout methods are often evaluated as extremely advantageous, but, in fact, they also have disadvantages. For example, based on results obtained during an experiment, the physical model is more suitable for in situ conditions than the numerical model [30]. However, it is demanding for physical models to not only fully reconstruct 3D topography, due to scaling effects [31], but also reflect the dynamic real conditions of field observations [32]. The numerical model’s precision is contingent upon the parameters established within the model [25]. The numerical model, which is an expensive method, should be applied by experts who specifically possess experience in numerical analysis [33]. It is not suitable for rapid decision making because the preparation of all potential simulations is a time-consuming task [29]. In contrast, the empirical–statistical model enables quick assessments due to its lower computational times [34]. In fact, the empirical–statistical model can be used practically and realistically within a reasonable timeframe with a lower cost and minimal calculation demands. However, volume neglect and volume prediction based on past landslide events are challenging in complex geological conditions [34,35]. To decide on a suitable method for runout analysis, available data are important [36]. It may be claimed that the best runout model method is inconclusive, considering the advantages and disadvantages of both methods. Nevertheless, recent studies point out that although researchers have used various methods over the years to investigate methods of estimating runout distance that may contribute to landslide risk evaluations, empirical–statistical methods are still generally preferred because of computational costs [37]. The reach angle method was often considered as an empirical method by many researchers [3,4,38,39]. The multiple flow direction [40], Holmgren [41], modified Holmgren [42], cellular automata [43], and random walk [44,45] models provide opportunities to simulate possible runout paths. A combination of the empirical–statistical simplified friction-limited model (SFLM), which considers the reach angle, and flow direction algorithms, enables us to observe details such as feasible runout distances and paths using Flow-R software [42], as shown in recent studies [39,46,47,48,49,50,51,52,53,54]. Owing to the development of computer technology, recent studies on landslide runout distance prepared with software go well beyond the conventional techniques. DebrisFlow Predictor [9], DFLOWZ [55], and the Progressive Debris-Flow Routing and Inundation model [56] are alternative empirical model programs that currently give the ability to inspect not only debris flow propagations but also shallow landslide runout scenarios. Consequently, the main purpose of this study is to investigate the empirical models of the shallow landslide runout distance in Eocene flysch facies in the Western Black Sea region of Türkiye. Concentrated precipitations often contribute to the occurrence of shallow landslides [57]. Extreme rainfalls in the region containing Eocene flysch facies are on the rise, so residents have complained about high frequencies of shallow landslide occurrences. Considering shallow landslide occurrences and the spatial probability evaluation of the hazard, integration of landslide runout distance with a landslide susceptibility assessment seems to be indispensable [5,8,58,59]. Flow-R software, which was developed at the University of Lausanne [42], was used to empirically examine shallow landslide runout scenarios in this study. This software offers researchers the opportunity to try a variety of landslide runout simulations [39,46,47,48,49,50,51,52,53,54,60,61,62,63] to obtain runout distances by considering different parameter configurations.

Table 1. Summary of studies that utilized Flow-R.

Authors	Research	Landslide Initiations	Runout Method	Flow Direction Algorithm
Horton et al. [42], Jiang et al. [51], Xu et al. [54], Park et al. [61], Sharma et al. [63]	Debris flow runout susceptibility map	Initiations detected by Flow-R	SFLM	Modified Holmgren
Pastorello et al. [46], McCoy [47], Paudel [48], Putra et al. [53]	Debris flow runout susceptibility map	User-defined initiations considering flow accumulation and slope [46], landslide susceptibility map [48], and Sentinel images [53]	SFLM	Modified Holmgren
Ali et al. [37], Bera et al. [62]	Debris flow runout susceptibility map	User-defined initiations based on remote sensing, slope angle distribution [37], and GPS and multi-temporal satellite images [61]	SFLM	Holmgren
Giano et al. [2]	Debris flow runout susceptibility map	Initiations detected by Flow-R	Perla	Modified Holmgren
Polat and Erik [60]	Debris flow runout susceptibility map	User-defined initiations considering landslide susceptibility map	Perla	Modified Holmgren
Charbel and El Hage Hassan [50]	Mudflow runout susceptibility map	Initiations detected by Flow-R	SFLM	Holmgren
Do et al. [49]	Landslide runout susceptibility map	Initiations detected by Flow-R	SFLM	Modified Holmgren
Liu et al. [52]	Landslide runout susceptibility map	User-defined initiations based on previous studies and D-InSAR technology	SFLM	Modified Holmgren

provides a summary of previous Flow-R studies.

If the runout distance can be estimated, the resilience to shallow landslides will be stronger. Accordingly, this shallow landslide investigation focuses on the development of more accurate empirical models to obtain better estimations of shallow landslide runout distances. Using machine learning to detect landslide initiations is the critical stage of this study; however, accurate model parameter configuration also constitutes one of the most prominent stages in the algorithm provided by Flow-R. Modified Holmgren and SFLM parameters are used to estimate runout distances in the propagation stage in Flow-R. Considering the difficulty in distinguishing and the connectedness of shallow landslides and debris flow, the decision of whether a movement is a shallow landslide or debris flow depends on a researcher’s judgment [64]. It is also impossible to estimate which types of landslides will occur (whether they are debris flows or shallow landslides) in Eocene flysch facies. Therefore, runout results are comparatively presented with respect to two different parameter configurations, which are shallow landslide and debris flow. Additionally, the question of whether the empirical–statistical method is suitable for our large study area is also discussed deeply in the later sections of this paper. Moreover, it is expected that the changes in precipitation trends caused by climate change will lead to changes in landslide occurrence in future years [65,66,67,68,69]. Thus, this study also aims to show that it is possible to analyze the near and distant future runout distances of shallow landslides considering the climate RCP scenarios’ precipitation differences. The empirical model parameters suggested for debris flow and shallow landslides are also investigated comparatively by including RCP 4.5 and 8.5 monthly precipitation values.

Finally, the key research points of this paper can be listed as follows: (i) select the most suitable threshold values for detecting shallow landslide initiations for scenarios with and without RCP values, (ii) determine the appropriate propagation parameters, and (iii) estimate the probable maximum shallow landslide runout distances for scenarios with and without RCP values.

2. Materials

2.1. Study Area

The study area was the Eocene flysch facies located in the Western Black Sea region of Türkiye, covering an area of about 877 km² (Figure 1). Bartin and Zonguldak are the main cities located on the Eocene flysch facies. The research area was determined by considering the natural geological boundary consisting of Eocene-aged clastics and carbonates [70]. The Eocene flysch facies contain weak and weathered sandstone, claystone, and silt and have low topographic values. Their altitude varies between 4 and 1571 meters (m). The altitude in the northeast is low, while in the southwestern part it is relatively high. Rainfall is almost equally drastic at all meteorological stations in the study area. According to the Turkish State Meteorological Service [71], extreme rainfall of 252.8 mm was observed on 11 August 2021 in only three hours at the Bartin Ulus Ceyupler station, which is very close to the study area. Therefore, it is not surprising to posit that shallow landslides abound in this formation owing to heavy rainfall. In the recent past, this region was frequently exposed to heavy rains. For instance, Can et al. [72] stated that due to heavy rainfall in May 1998, shallow landslides were observed in Eocene flysch facies in this region. In addition, the Turkish State Meteorological Service’s historical meteorological precipitation data indicate that 93.2 mm of daily rainfall occurred in Bartin on 21 May 1998 [73]. During this meteorological event, 91.8 mm of daily rainfall was also observed in Zonguldak. As mentioned, shallow landslides are frequently observed in Eocene flysch units due to heavy rainfalls [74] (Figure 2). Figure 2a clarifies the shallow landslides triggered after the May 1998 precipitations, while Figure 2b–d show examples of shallow landslide occurrences after heavy rainfall in November 2023. In total, 60 and 111 mm of daily precipitation were observed in the cities of Bartin and Zonguldak on 19 November 2023, respectively [75]. Figure 2 shows the types of landslides detected in the study area according to the Varnes (1978) classification frame [76]. The study area is predominantly covered by complex cultivation patterns and broad-leaved forests [77]. Mixed forests and land principally occupied by agriculture with significant areas of natural vegetation are also observable in the study area [77]. Shallow landslides are frequently observed in both mixed forests and broad-leaved forests in this study area, according to Corine Land Cover [77].

Although shallow runout distances have not yet been investigated in Eocene flysch facies, some landslide analyses have already been conducted in the formation. Duman et al. [78] prepared a landslide inventory of Northwestern Anatolia involving Eocene flysch facies in the region. Akgun et al. [74] evaluated the general characteristics of the landslide sizes in the formation by considering magnitude and frequency relations. In their study, they emphasized that while rollover effects were noticed at 0.048 km², the fractal dimensions of the distributions were noticed as 1.41 [74] in the flysch facies.

Studying large areas makes parameter selection difficult. For this reason, the study area was divided based on the basin boundaries. ArcGIS software enabled the division of the study area with respect to the basin boundaries by considering the watersheds in the region. The Eocene flysch facies were divided into three sub-basins in which runout analyses could be performed (Figure 3). To perform appropriate parameter selection by decreasing the size of the investigation area, smaller sub-basins could also be considered. However, studying larger areas is important for researchers investigating shallow landslide runout distances. The areal distributions of the sub-basins are as follows: Egerci basin in the south is 269 km², Beycuma basin in the middle is 311 km², and Ihsanoglu basin in the north is 297 km².

2.2. Data Source and Pre-Processing

It should primarily be noted that the preparation of a susceptibility map was especially critical to determine the initiation points that were used during the runout analyses. Therefore, the following data were obtained to prepare a shallow landslide susceptibility map and to calculate runout distances in this study: First, Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) data with a spatial resolution of ~30 m were obtained. They were converted to a 25 m resolution using the nearest resampling technique in ArcGIS. The coordinate system was set as WGS 1984 UTM zone 36N. Second, a DEM with a spatial resolution of 25 m was also implemented to obtain topographic features such as the slope gradient and slope aspect using the System for Automated Geoscientific Analyses (SAGA) GIS. In addition, the slope length factor, valley depth, topographic wetness index (TWI), plan curvature, profile curvature, convergence index, closed depressions, channel network base levels, and channel network distance were obtained in the SAGA GIS environment to identify the geomorphometric features of the study area. Third, the geological information was evaluated based on [38]. Subsequently, Corine Land Cover [77] was used to assess the land cover in the study area with a 100 m spatial resolution. Fifth, NDVI data were obtained from Sentinel-2 images with a 10 m spatial resolution.

NDVI = (NIR − Red)/(NIR + Red)

(1)

Here, NDVI is the Normalized Difference Vegetation Index, NIR is the Near-Infrared Band (band 8), and Red is the Red Band (band 4) of a Sentinel-2 image. Finally, the spatial distributions of recent shallow landslides were mapped using Sentinel-2 images and Google Earth. Short field surveys were also carried out to perform ground checks of the mapped shallow landslides. Furthermore, available landslide inventories published in national databases and the recent literature [72,78,79,80] were also evaluated, particularly to identify shallow landslides during mapping using satellite images. When an available landslide inventory was examined considering historical satellite images from Google Earth, 262 shallow landslides were mapped. The landslides were geometrically classified in terms of failure depths, reach angles, and observed runout distances. The volume and area relations proposed by Hovius et al. [81] based on Jaboyedoff et al. [82] and the semi-ellipsoid approach were used to determine the landslide depths. The obtained values proved that the failures mapped in this study can be classified as shallow landslides since their depths are less than 5 m [83]. The descriptive statistics of the failure depths, reach angles, and observed runout distances are given in

Table 2. The properties of the landslide inventories prepared for the sub-basins.

Sub-Basin	N	Statistics	Area (m2)	Failure Depth (m)	Travel Angle (°)	Observed Runout Distance (m)
Egerci	111	Min.	21	0.2	1	7
		Max.	4116	3.2	49	122
		Mean	536	1.0	23	53
		Median	344	0.9	24	47
		Std. Deviation	572	0.5	12	29
Beycuma	15	Min.	23	0.2	4	7
		Max.	1843	2.1	28	83
		Mean	376	0.9	13	36
		Median	258	0.8	13	31
		Std. Deviation	438	0.5	8	20
Ihsanoglu	136	Min.	7	0.1	1	5
		Max.	1710	2.1	32	122
		Mean	140	0.5	11	25
		Median	65	0.4	10	17
		Std. Deviation	226	0.3	7	21

with respect to the sub-basins. Figure 2b,d, which were taken after heavy rainfall in November 2023, provide good evidence that shallow landslides have continued to occur due to excessive rainfall in this region, and their observed runout distances are compatible with Table 2.

3. Methods

The method of this study was divided into four critical interconnected steps: (i) prepare a shallow landslide susceptibility map using the machine learning logistic regression (LR) method; (ii) determine shallow landslide initiations by considering the shallow landslide susceptibility map and the RCP scenarios’ precipitation values; (iii) assess the shallow landslide runout distance methodology by taking into consideration detected shallow landslide initiations; and (iv) select suitable model parameters for runout distance assessment.

3.1. Preparation of Shallow Landslide Susceptibility Map

Considering the constant lithological information of the Eocene flysch facies, fifteen features, including the altitude, slope gradient, slope aspect, landslides, land use, NDVI, slope length factor, valley depth, topographic wetness index (TWI), plan curvature, slope profile curvature, convergence index, closed depressions, channel network base levels, and channel network distance, were used as the landslide conditioning factors for landslide susceptibility mapping (

Table 3. Descriptive statistics of the topographic parameters with respect to the sub-basins.

Sub-Basin	Statistics	A	SG (°)	SA (°)	PLC	PRC	CI	CD	TWI	SLF	CNBL (m)	CND (m)	VD (m)
Egerci	Min.	41	0	0	−0.01	−0.01	−92	0	2	0	41	−49	−3
	Max.	1574	60	6.28	0.01	0.012	96	20	21	38	1382	449	549
	Mean	444	15	3.25	0	0	0	0	7	5	381	64	236
	Median	425	14	3.04	0	0	0	0	6	4	369	49	229
	Std. Deviation	259	8	1.91	0.002	0.002	8	1	2	3	233	60	133
Beycuma	Min.	13	0	0	−0.01	−0.01	−100	0	3	0	13	−27	0
	Max.	658	52	6.28	0.008	0.009	100	19	20	24	444	304	277
	Mean	170	11	3.27	0	0	0	0	7	3	129	41	140
	Median	158	9	3.33	0	0	0	0	7	2	118	32	146
	Std. Deviation	91	6	1.91	0.001	0.001	9	1	2	2	69	38	61
Ihsanoglu	Min.	4	0	0	−0.01	−0.02	−98	0	2	0	6	−39	0
	Max.	633	70	6.28	0.017	0.021	92	32	20	30	284	468	339
	Mean	121	11	3.35	0	0	0	0	7	3	79	42	117
	Median	104	10	3.43	0	0	0	0	6	2	69	32	117
	Std. Deviation	71	7	1.83	0.001	0.001	10	1	2	2	43	42	51

Abbreviations: altitude: A, slope gradient: SG, slope aspect: SA, plan curvature: PLC, profile curvature: PRC, convergence index: CI, closed depression: CD, topographic wetness index: TWI, slope length factor: SLF, channel network base level: CNBL, channel network distance: CND, valley depth: VD.

).

An open-source machine learning library in Python provided the opportunity to produce a landslide susceptibility map at this stage. The LR method is valuable, given its high prediction accuracy, in machine learning research, such as the preparation of landslide susceptibility maps [84,85,86,87,88,89,90]. For instance, Polykretis and Chalkias [91] prepared a landslide susceptibility map for shallow rotational and translational slides and the flow types of landslides using this method. The authors stated that the LR method is the better option. Nhu et al. [86] stated that the LR method was the second most successful and trustworthy model in their comparison of the production of shallow landslide susceptibility maps. Nwazelibe et al. [90] indicate that LR is the better method to produce landslide susceptibility maps. Landslides, which are shallow translational debris slides, slumps, and occasional debris topples, frequently occur in the Eocene-aged Ameki-Bende formation [90]. Taking these examples from the literature into account, in this study, LR was used to determine shallow landslide susceptibility.

P = 1/(1 + e^−y) = e/(1 + e^y)

(2)

Here, P is the landslide probability occurrence, which changes from 0 to 1, while e is the exponential constant [92]. In Equation (3), y is a dependent variable. When b₀ is a constant intercept from the model, b_n is the triggered landslide factor of a_n, which represents an independent variable.

y = b₀ + b₁a₁ + b₂a₂ +⋯… + b_na_n

(3)

The ratio of the training/testing samples was selected as 80:20. The data were divided into training (80%) and testing (20%) sets [91,93,94,95,96,97]. The “train_test_split” and “cross_val_ score” functions of the sklearn model selection were used in the open-source Python libraries for training and testing the data.

Shallow landslide susceptibility maps obtained using the LR method for the sub-basins are given in Figure 4. There are three classes indicating low (0–0.4), moderate (0.4–0.7), and high (0.7–1) probabilities of shallow landslide occurrences.

3.2. Determination of Shallow Landslide Initiations

Possible shallow landslide initiation determination is a significant concern for landslide researchers. Even though it is difficult to make precise predictions about initiations, a susceptibility map can provide a basis to estimate possible shallow landslide initiations. Similar considerations were also given by Xu et al. [54]. Xu et al. [54] indicated that more susceptible areas on susceptibility maps have higher initiation abilities. Can et al. [72] marked that due to heavy rainfall in May 1998, shallow landslides were observed in very highly susceptible zones in the Eocene flysch facies in this region. In this study, highly susceptible areas, which are called shallow landslide initiations, were implemented for runout distance estimation using Flow-R software. In fact, the determination of the susceptibility threshold of each basin was key to identifying the shallow landslide initiations in this study because they were used as pre-defined initiations in the runout distance assessment stage. The QGIS environment was used to visualize the shallow landslide susceptibility maps and apply the necessary threshold value of susceptibility to determine the initiations (Figure 5).

When shallow landslide initiations are identified according to a susceptibility threshold, it is a good idea to be on the safe side. The threshold value was selected according to the frequency histograms of the shallow landslide susceptibility values of the sub-basins (Figure 5). The threshold value for the shallow landslide initiations in the sub-basins was determined to be 0.70. Therefore, shallow landslide susceptibility grid cells with values equal to or higher than 0.70 and recent landslides mapped during the inventory stage were chosen as the shallow landslide initiations for the sub-basins (Figure 6). According to Corine Land Cover [77], the initiations in the Beycuma and Ihsanoglu basins are frequently abundant in broad-leaved forests, while those in the Egerci basin are abundant in broad-leaved forests and mixed forests.

In this part of the study, Community Climate System Model (CCSM4) projections based on RCP 4.5 and RCP 8.5 monthly precipitation scenarios for both cities were also analyzed using an IDW heterogeneous environment method in ArcGIS based on data obtained from NCAR [98]. Geological hazards, such as shallow landslide initiation points, may exhibit a spatially heterogeneous distribution. IDW interpolation demonstrates greater success in heterogeneous environments [99]. The IDW method is considered the best method due to its reasonable processing timeframe and adequate precision for heterogeneous data [100]. The IDW method was also used to analyze the spatial distributions of the RCP scenarios’ rainfalls in studies by Mohamed Yusof et al. [101] and Nasidi et al. [102]. Although 30 years is the most widely used reference period, 20 years has recently been used in some works [103]. Therefore, the time periods were classified as 2025–2044, 2044–2063, 2063–2082, and 2082–2100. In addition, 20-year monthly median values were used to analyze the general precipitation trend in order to apply the IDW method. Selected precipitation threshold values were utilized in order to detect critical shallow landslide initiations by considering the IDW results of the three sub-basins. Historical meteorological data for two stations, Bartin and Zonguldak, and RCP 4.5 and 8.5 scenario precipitation values for both cities were evaluated by separately using cross-correlation and choosing the reference period of 1995–2022. This step aimed to determine how extreme precipitations are represented in the obtained RCP data. Therefore, extreme precipitation values that resulted in shallow landslides in the study area were considered using this cross-correlation equation in order to reach predicted precipitation values. Considering these values, finding the median precipitation values was a critical stage. Both scenarios gave 81 mm monthly critical precipitation threshold values. Although the probability of this precipitation amount falling in one month or one day will be uncertain, it is possible to expect this amount of rainfall in one day by considering previous daily rainfall data and recent extreme daily rainfall data in this area. Hence, it was revealed that for the two scenarios the critical precipitation threshold of 81 mm would be applied in the analyses. Nevertheless, in this condition, shallow landslide initiations should supply both conditions together. These conditions are that the shallow landslide susceptibility is equal to or more than 0.70 and the RCP precipitation scenarios are equal to or more than 81 mm. In addition, while including the RCP precipitation scenarios’ effects in the runout distance analyses, the impacts of the other triggering landslide parameters were neglected, and it is admitted that they might stay the same except for the RCP precipitations. Shallow landslide initiations, which were detected by considering the RCP precipitation scenarios with respect to divergent time periods, are shown in Figure 5. Despite the reduction in shallow landslide initiations, the Egerci sub-basin and the Beycuma sub-basin may be exposed to shallow landslides, while the Ihsanoglu sub-basin may not be affected by shallow landslides in the presented time range for both scenarios (Figure 7 and Figure 8). Therefore, the Ihsanoglu sub-basin is not presented in Figure 7 and Figure 8. Although there is no big difference between the results of the other three scenarios apart from the RCP 4.5: 2025–2044 scenarios in the Egerci sub-basin, differences can be observed in the northwestern part of the Egerci sub-basin in terms of shallow landslide initiations. The differences between in the presented RCP 4.5 and 8.5 scenarios are considerably larger in the period of 2025–2044 for both the Egerci and Beycuma sub-basins. In addition, the RCP 8.5 scenario only affects shallow landslide initiations in the period of 2025–2044 for both sub-basins. In the Beycuma sub-basin, shallow landslide initiations can be observed for the RCP 4.5 scenario in the periods of 2025–2044 and 2082–2100, while in the Egerci sub-basin shallow landslide initiations can be examined in the periods of 2025–2044, 2063–2082, and 2082–2100. All results show that shallow landslide initiations decrease if the RCP precipitation scenarios are considered; therefore, this study also has the same opinion about landslide susceptibility decrement depending on RCP scenarios as the previous studies by Ageenko et al. [65] and Park and Lee [104].

3.3. Shallow Landslide Runout Distance Assessment Methodology

The empirical–statistical methods were applied in this study to analyze the runout distances of shallow landslides. Models are generally prepared using the empirical method for large-scale studies [46]. It is important to choose a suitable method and algorithm to determine how to implement the software [42]. The flow direction algorithm and the simplified friction-limited model (SFLM) are two main stages that need to be applied to obtain probability of runout from the software during the stages of spreading and energy calculation. Essentially, the combination of the spreading algorithm and energy calculation gives the possible runout paths and distances. First, it is indispensable to choose the flow direction algorithm to be used in a study. Therefore, the question of which flow direction algorithm method is preferable for use in landslide runout studies should be pondered. Holmgren [41] added the “x” parameter to the multiple flow direction approach of this algorithm, allowing studies to be carried out in more detail. In Equation (1), while i and j represent flow directions, ρ_i^fd is the sensitivity ratio in the i direction, tanβ_i is the gradient between the center cell and the cell in the i direction, tanβ_j is the slope gradient between the center cell and the cell in the j direction, and x is the variable parameter that controls the variable deviation. As can be seen in Equation (4), x can take values between one and infinity.

ρ_i^fd = (tanβ_i)^x/(∑⁸ _(j=1) (tanβ_j)^x) ∀ {(tanβ > 0@ × ε[1;+∞])}

(4)

The modified Holmgren algorithm proposed by Horton et al. [42] is one of the most favorable flow direction algorithms because it includes the “dh” parameter, which varies the gradient values by varying the height of the central cell. It smooths DEM roughness. Since low-resolution data were used, the modified Holmgren algorithm [42] was used in this study. Putra et al. [53] stated that they also chose this algorithm because it is less sensitive to the DEM resolution. The SFLM is the second key part of the propagation stage. In this stage, the energy should be sufficient when it travels to another cell. Energy loss is caused by friction [42,48].

E_kinⁱ = E_kin⁰ + ∆E_potⁱ − E_fⁱ

(5)

In Equation (5), E_kinⁱ, which is the kinetic energy (cell in direction i), is calculated using ∆E_potⁱ and E_fⁱ. While ∆E_potⁱ is the variability in the potential energy (cell in direction i), E_fⁱ is energy loss due to friction (cell in direction i).

E_fⁱ = gΔx tanφ

(6)

The E_fⁱ calculation multipliers are g, Δx, and tan φ (Equation (6)). Δx represents the horizontal displacement increment, while tan φ is the energy line gradient. g is gravitational acceleration and is also used to calculate energy loss. The reach or travel angle is significant for the SFLM during the energy calculation. Past landslide runout travel angles are the key to understanding the region’s landslide travel angle trends. Since the study area was divided into three sub-basins, different travel angles were determined for each sub-basin separately. The properties of the landslide inventory prepared in this study, involving the depth of the estimated failure surface, runout distances, and travel angles, were used to identify the parameter sets that achieved ideal model performance across the models of the sub-basins. The determination of the median travel angle for the sub-basins was important for the assessment of the probabilities of runout distances in larger areas.

The selection of an appropriate velocity and energy limit was also important for this study to remain within the acceptable energy levels and to avoid an endless energy loop during the analysis (Equation (7)).

V_i = min{√(V₀² + 2g∆h − 2g∆x tanφ), V_max}

(7)

Here, V₀ and V_max indicate the initial speed and the maximum speed limit, respectively.

3.4. Selection of Model Parameters

Appropriate parameter selection considering the region is important during model preparation. In this study, Flow-R software, which is a distributed empirical model, was used to model possible runout distances for the sub-basins. The software enabled comparison of the runout distance results not only for debris flows but also for shallow landslides by considering proposed parameters. Flow-R software promoted the utilization of a range of parameters to enhance the accuracy of the research. Many simulation scenarios can be considered to investigate the effects of the selected parameters because the parameters can change more rapidly in actual conditions. The parameter x can become valuable for understanding the effect of the flow direction. The value of x is equal to 4 for debris flows, while it can be set as 22 to 26 for shallow landslides [105]. The value of x has also been set as 4 for debris flows in recent studies [37,47,51,61]. In this study, comparative analyses were performed depending on whether failures occurred as debris flows or shallow landslides. Therefore, considering the suggestions already given in the literature, the value of x was set as 4 or 25 for the debris flows and shallow landslides in the analyses, respectively. The value of dh is also a significant parameter that affects results by changing the reduced rate of roughness. The value of 1 was suggested for dh in recent studies [2,49,53]. Horton et al. [105] also proposed that dh can be equal to 1 for both debris flows and shallow landslides. Since slope morphology can change in various ways, travel angles can also change to various degrees. For instance, slope gradients are higher in the western part of the Eocene flysch facies. Travel angles were also evaluated separately for the three sub-basins by considering back analyses. Accordingly, median travel angle degrees were determined for the sub-basins separately as follows: 24° for the Egerci sub-basin, 13° for the Beycuma sub-basin, and 10° for the Ihsanoglu sub-basin. Horton et al. [105] suggested that travel angles of 8°–13° can be used for debris flows, while they may be 16° to 20° for shallow landslides. Different researchers considered different travel angles in the literature, e.g., 5° [38,49], 9° [51], 10° [47], 11° [42], 12° [37], and 13° [61]. Back analyses also show that if the other parameters are kept the same and only the travel angle is decreased, runout distances will increase. The last parameter that should be considered during analyses is velocity. The estimation of velocity is difficult. Back-calculated debris flow velocities were observed in the range of 5 to 15 m/s in a study by Prochaska et al. [106]. Similarly, different researchers considered different velocity values in the literature, e.g., 5 m/s [49], 10 m/s [37], 15 m/s [42,47,54], 16 m/s [107], 20 m/s [46], and 25 m/s [51]. By considering both back analyses and the local information obtained by residents having memorials related to the shallow landslides that occurred in May 1998 in this region, it was concluded that the velocity parameter could be set as 15 m/s for shallow landslides in the Eocene flysch facies. The same velocity value was also used for debris flow during the comparative analyses. This parameter selection was compatible with the suggestions that were given by Horton et al. [42]. Apart from the velocity parameter, energy limitation should also be considered to obtain logical runout distance results. Consequently, the DEM and pre-defined landslide initiations were imported into Flow-R. Files should be converted to ASCII format. It is also necessary to be aware that files should have the same spatial resolution and coordinate system. The model parameters implemented in Flow-R for shallow landslides and debris flow are summarized in

Table 4. The model parameters implemented in Flow-R in this study.

Sub-Basin	Failure	Flow Direction Algorithm	x	dh (m)	Travel Angle (°)	Velocity (m/s)
Egerci	Debris flow	Modified Holmgren	4	1	24	15
	Shallow landslides	Modified Holmgren	25	1	24	15
Beycuma	Debris flow	Modified Holmgren	4	1	13	15
	Shallow landslides	Modified Holmgren	25	1	13	15
Ihsanoglu	Debris flow	Modified Holmgren	4	1	10	15
	Shallow landslides	Modified Holmgren	25	1	10	15

with respect to the different sub-basins.

4. Results

4.1. Shallow Landslide Runout Distance Assessment

Effective determination of the runout distances of shallow landslides is significant to understand and assess the importance of the negative effects of shallow landslides. As runout research has been developing rapidly, it is difficult to agree on which methods will become the most practical in the future. This study has pointed out that the empirical–statistical method can determine the runout distances of shallow landslides. Runout distance models have approximately the same performance depending on the debris flow and shallow landslide parameters. Runout distances in the Eocene flysch facies can be estimated at a coefficient of determination level of 0.64 with debris flow parameters, while this value was obtained as 0.62 for the models developed using the shallow landslide parameters (Figure 9). Therefore, it is possible to bring the hazard regions up to date by integrating the resultant runout zones into landslide susceptibility maps. The results of the runout distance evaluations for the sub-basins are represented in Figure 10, Figure 11 and Figure 12 with respect to the model parameters considered for debris flows and shallow landslides.

The uncertainty about whether observed landslides are shallow landslides or debris flows continued in the statistical–empirical evaluations. Therefore, models were prepared for both parameter options to determine where shallow landslides or debris flows may reach because it was difficult to make decisions about the exact detection of whether landslides are shallow landslides or debris flows on a regional scale. It was possible to compute that the runout distances did not differ significantly when the simulations with both parameters were examined. It can also be stated that it was hard to distinguish differences between the debris flow analysis results and the shallow landslide analysis results. Changes in distances generally happened gradually in various areas for the debris flow analysis parameters. Moreover, the runout distances reached farther in the debris flow analyses. However, highly probable runout distances were not as common in the debris flow analyses as in the shallow landslide analyses. The highly probable range (0.5–1) indicates that a grid cell is involved in a runout zone with a probability value higher than 50%, while the low-probability range (0–0.5) means that a grid cell is included in a runout zone with a probability value lower than 50%. An assessment of the runout areas indicated that while the models with debris flow parameters impacted a total area of 146 km², the models with shallow landslide parameters impacted a total area of 101 km².

4.2. Shallow Landslide Runout Distance Assessment Considering RCP Scenarios

Including the RCP precipitation scenarios in the runout distance analyses in the study area adequately fulfilled the need for necessary precautions related to shallow landslides in the future. Simulations of possible runouts were carried out for shallow landslides and debris flow to investigate whether their parameters affect the runout distances in the Eocene flysch facies for detected landslide initiations, which were determined considering the RCP 4.5 and RCP 8.5 precipitation scenarios. Possible future runout areas that will be affected by debris flows or shallow landslides will decrease because of a reduction in shallow landslide initiations (Figure 13, Figure 14 and Figure 15). It is also not surprising that due to a lack of shallow landslide initiations, shallow landslides will not occur in the Egerci, Beycuma, and Ihsanoglu sub-basins during some periods because shallow landslide initiations were factored into the estimation of the runout distances in this study. Although shallow landslide initiations will decrease in the Egerci sub-basin in the period of 2025–2044 in the RCP 4.5 scenario (Figure 14), it is possible to claim that there are no big differences between these runouts and the runouts shown earlier, which did not include the RCP precipitation scenarios. However, in the same period (2025–2044), there will be discernable declines in shallow landslide initiations and their runout distances in the Beycuma sub-basin. It is notable that dramatic declines in shallow landslide initiations will also be clearly visible in the Egerci and Beycuma sub-basins in the period of 2025–2044 for the RCP 8.5 scenario. In addition, Figure 13, Figure 14 and Figure 15 give the opportunity to compare the effects of the RCP 4.5 and 8.5 precipitation scenarios for the Egerci sub-basin and the Beycuma sub-basin in the period of 2025–2044 in terms of the divergence of their runout distances. Compared with the initial analyses, there is a trend of decreases in shallow landslide initiations and their runout distances in the Egerci and Beycuma sub-basins. However, it is noteworthy that fluctuating changes in shallow landslide initiations and runout distances will be observed in both sub-basins. For instance, they will show an increasing trend after the period of 2044–2082 in the Beycuma sub-basin, despite the fact that they will not reach the levels of the initial analyses (Figure 13).

It should be underlined that all obtained runout model results depend on the RCP scenarios. If the RCP scenarios are changed or do not turn out as projected, results such as landslide initiations and their runout distances might change.

5. Discussion

The runout distances of shallow landslides are important outputs that may facilitate improved resilience against shallow landslides. According to the results of empirical evaluations, runout model outputs indicate that shallow landslides can damage transportation networks and residential areas in this region. For example, Nie et al. [108] assessed the landslide risk along national highways based on information quantity and random forest coupling methods. Flow-R, a distributed empirical model, was implemented with great success not only in this study but also in many studies in the literature. Giano et al. [2], Polat and Erik [60], Charbel and El Hage Hassan [50], and Jiang et al. [51] discussed whether Flow-R and observed landslide runout results are compatible. The researchers concluded that Flow-R results are well matched with the observed results. Despite not considering rheological parameters, Flow-R model results are coherent [51,53]. Do et al. [49] stated that the Flow-R model may give exaggerated results. However, the results can be used to evaluate all possibilities. Putra et al. [53] also put forward that although Flow-R does not provide complete hazard mapping, it does enable the prediction of possible runout debris flow paths. In terms of analyzing shallow landslide runout distances in the three sub-basins, the selection of model parameters and data quality had significant roles. Despite employing a limited dataset, Flow-R was proven to attain accurate results, especially for the Ihsanoglu sub-basin, with an R-squared level of about 79%. Applying a machine learning algorithm was advantageous in this study during the preparation of a susceptibility map to obtain quite reliable shallow landslide initiations. Flow-R software can also be implemented to distinguish landslide initiations by applying a grid-based overlay approach [2,50,51,54]. However, the pre-determination of landslide initiations via machine learning is thought to have a positive effect, increasing the success of the analyses with respect to the accuracy of the runout distance probability. Shallow landslide initiations will decline in the future RCP scenarios for the sub-basins, which means the landslide runout distances will decrease [65,104]. If a more detailed examination is conducted on shallow landslide initiations, they are abundant in broad-leaved forests in the three sub-basins in the initial analyses. The literature review also showed that shallow landslide initiations are often observable in coniferous forests and broad-leaved forests [109,110]. However, climate change scenarios indicate that the susceptibility of broad-leaved forests will decrease [111,112]. If shallow landslide initiation detection analyses, including the RCP precipitation scenarios, are considered, shallow landslide initiations show a decreasing trend in broad-leaved forests, especially in the Beycuma sub-basin. The period between 2025 and 2044 will be the most affected period for both sub-basins in the RCP 4.5 scenario in terms of both shallow landslide initiations and their runout distances. These analysis results are critical for estimating the hazard rates of shallow landslides in the study area in terms of their contributions to disaster management strategies.

Limitations of this Study

Improving the spatial resolution of the DEM is important to obtain better results [50]. Horton et al. [42] stated that DEMs with a 10 m spatial resolution provide higher accuracies. However, a spatial resolution of 25 m can also be implemented successfully in analyses. The models in this study were also sensitive to DEM effects. Nevertheless, considering the size of the study area, the computational cost of higher resolutions could not be managed. Supercomputers may be a solution to this situation. Moreover, despite many attempts at changing the DEM resolution to analyze runout distances, propagations were not triggered for some initiations in datasets examined considering both 10 m and 25 m resolutions. The main reason for this peculiarity may be that the DEM stopped or was deflected by natural environment features and building features [53]. Furthermore, the selection of the velocity parameter also may have affected the results. Increasing runout distances should not be surprising in view of an increasing velocity. For instance, Horton et al. [105] suggest that the velocity limitation is 8 m/s for mudflows. Therefore, it has become clear that runout distances are affected by the maximum velocity limitation and should be considered cautiously during model parameter determinations. In this study, a back analysis was carried out just for this purpose. In addition, the landslide initiation process is primarily influenced by geological conditions [113]. Conducting further research that integrates lithological and geological investigations along with hydrogeological conditions would highlight the significance of these data in the initiation process of landslides [114]. Therefore, an additional focus on hydrogeological conditions may be required to improve the accuracy of probable shallow landslide initiations.

6. Conclusions

Runout distances have been ruled out in many studies. Considering the successful results, runout will be considered more in future studies because runout distance probabilities in susceptible areas enable us to make more appropriate decisions in hazardous regions. The goal of this study was to determine the shallow landslide runout distance probabilities in the Eocene flysch facies in the Western Black Sea region of Türkiye. This research offers evidence that the empirical–statistical method can successfully evaluate the runout distances in the study area. Despite the different empirical model parameters suggested for shallow landslides and debris flow, the models give nearly identical runout distance results in the Eocene flysch facies. Runout distances can be predicted at R-squared levels of 0.64 and 0.62 for the debris flow and shallow landslide models, respectively, in this study area. The total impacted runout areas were estimated to be approximately 146 km² and 101 km² for the debris flow and shallow landslide parameters, respectively. This study also reveals that machine learning algorithms can be remarkably fast and accurate in providing shallow landslide initiations in the sub-basins. The successful inclusion of the climate change scenarios in the runout distance analyses was ensured. The RCP 4.5 and RCP 8.5 scenarios also give the possible future runout distances in the Eocene flysch facies. The simulated runout results, including the RCP 4.5 and RCP 8.5 scenarios, prove that areas affected by shallow landslides or debris flows will decrease in the future due to the decrement in shallow landslide initiations. The analysis results also show that the time period of 2025–2044, in accordance with the RCP 4.5 scenario, persists as a challenging interval for the Egerci and Beycuma sub-basins due to shallow landslide occurrences and their runout distances. In conclusion, this study reveals that if runout analyses can be carried out successfully, they contribute to the preparation of effective disaster management policies to increase the resilience to landslides for societies living in mountainous areas. Using more advanced machine learning and statistical techniques to predict failure initiations for shallow landslides is recommended to improve the accuracy of the runout distance results in future work. It is foreseen that relations between the RCP precipitation scenarios and the runout distances of shallow landslides will increasingly seize attention in enhanced landslide studies in the future.