A Hybrid Multi-Hazard Susceptibility Assessment Model for a Basin in Elazig Province, Türkiye

Preparation of accurate and up-to-date susceptibility maps at the regional scale is mandatory for disaster mitigation, site selection, and planning in areas prone to multiple natural hazards. In this study, we proposed a novel multi-hazard susceptibility assessment approach that combines expert-based and supervised machine learning methods for landslide, flood, and earthquake hazard assessments for a basin in Elazig Province, Türkiye. To produce the landslide susceptibility map, an ensemble machine learning algorithm, random forest, was chosen because of its known performance in similar studies. The modified analytical hierarchical process method was used to produce the flood susceptibility map by using factor scores that were defined specifically for the area in the study. The seismic hazard was assessed using ground motion parameters based on Arias intensity values. The univariate maps were synthesized with a Mamdani fuzzy inference system using membership functions designated by expert. The results show that the random forest provided an overall accuracy of 92.3% for landslide susceptibility mapping. Of the study area, 41.24% were found prone to multi-hazards (probability value > 50%), but the southern parts of the study area are more susceptible. The proposed model is suitable for multi-hazard susceptibility assessment at a regional scale although expert intervention may be required for optimizing the algorithms.


Introduction
Hazards are natural, technological, and human-induced events that potentially cause injury or loss of life, and damage to property, the socioeconomic well-being of society, the natural environment, and to historical and cultural resources (Bordbar et al. 2022). Earthquakes, floods, landslides, avalanches, wildfires, tsunamis, debris flows, and volcanoes are among the frequently observed natural hazards in the world. Most studies in the literature focus on a single hazard type (Wang et al. 2020a;Cetinkaya and Kocaman 2022;Feizizadeh et al. 2022;Karakas et al. 2022). However, multiple natural hazards often affect an area, interact with each other, and may occur as cascading events (Khatakho et al. 2021;Bordbar 2022). In such areas, disaster-related risks need to be predicted precisely to mitigate the adverse effects.
In recent years, a rise has been observed in the number of natural hazards associated with climate change and rapid urbanization (Dhar et al. 2017;Kocaman et al. 2020;Wang et al. 2020b) and it is expected to continue in the future. According to the reports published by the Disaster and Emergency Management Presidency of Türkiye (AFAD), approximately 344,560 natural hazards-including earthquakes, landslides, rockfalls, floods, and avalanchesoccurred in the country between 1950 and 2019 (AFAD 2020). Most natural hazard assessment studies in Türkiye have also focused on a single hazard type such as earthquakes (Erdik et al. 1999;Alpyurur and Lav 2022), landslides Akinci and Ozalp 2021), floods (Tiryaki and Karaca 2018), and wildfires (Arca et al. 2020). One natural hazard may trigger another, such as landslides induced by earthquakes or extreme rainfall.
Although susceptibility maps for a single natural hazard have been widely produced, there is no consensus on the method to be used when combining these maps. Therefore, this issue is still a subject that is open to development in the international natural hazards literature. For this reason, in this study, we proposed a novel multi-hazard susceptibility assessment (MHSA) approach for the regional scale and implemented it in an area prone to earthquakes, landslides, and floods, based on the Mamdani fuzzy inference system (FIS) (Mamdani and Assilian 1975). Yanar et al. (2020) presented the first example of employing Mamdani FIS for multi-hazard susceptibility assessment, but only two natural hazards were considered in their study. Thus, this study differs from previous studies as it is the first to use the Mamdani FIS to combine three natural hazards. The study site is located in Elazig Province, in the southeastern part of Türkiye, and is prone to multi-hazards due to its young and active tectonism, high seismicity, the weak shear strength of the lithological units, steep slopes, and climatic conditions. A 6.8 Mw earthquake in the province on 24 January 2020 caused over 30 fatalities, injured thousands, and triggered several landslides (Karakas, Nefeslioglu, et al. 2021). Extreme meteorological events also frequently cause floods in the region, yet no studies have appeared in the literature so far.
The methodology proposed here involves data-driven and expert-based methods for the production of univariate susceptibility maps for earthquakes, landslides, and floods individually, and combining them using an expertbased method. The expert inference required at various steps was provided by a coauthor, Gokceoglu, who has visited the region frequently. The methods preferred here were selected based on data availability, suitability for the target region, and their prediction performances. Thus, an ensemble machine learning (ML) method, random forest (Breiman 2001), was applied for the landslide susceptibility assessment (LSA) with optimized hyperparameters. Regarding the flood susceptibility assessment (FSA), an expert-based approach with a modified analytical hierarchical process (M-AHP) proposed by Nefeslioglu et al. (2013) was preferred and the score factors were determined based on the characteristics of the study area. For earthquake susceptibility, Arias intensity (AI) proposed by Arias (1970) was employed as the AI is a reliable parameter used to describe the earthquake shaking to trigger landslides. For the MHSA an expert-based approach, the Mamdani FIS, was used and the map was validated by the expert. Yanar et al. (2020) used the Mamdani FIS for the joint assessment of landslide and flood hazards in an urban expansion area in Ankara, Türkiye, but the present study includes earthquakes as well and proposes a more advanced strategy for defining the membership functions for each input. The proposed method is novel and produced an accurate multi-hazard susceptibility map (MHSM) for the region.
In the following, work related to the LSA, FSA, and MHSA is presented and discussed in Sect. 2. The methodological details, together with the data, are explained in Sect. 3, and the results of the LSA, FSA, and MHSA are presented in Sect. 4, followed by a brief discussion and conclusion in Sects. 5 and 6.

Related Work
Based on the literature analysis, it was apparent that most susceptibility assessment studies considered a single hazard type. Recent reviews with bibliometric analyses of landslide susceptibility research were provided by Liu et al. (2022), who analyzed a period between 1999 and 2021, and by Lima et al. (2022), who analyzed literature on data-driven methods utilized for LSA between 1985 and 2020. According to Liu et al. (2022), the trends in the LSA methods have evolved from expert-based, statistical, or probabilistic analysis to supervised ML in recent years. The number of LSA studies, especially with data-driven methods, is much greater thanks to the availability of comprehensive inventories, the increasing availability of geospatial datasets, and novel ML methods. Among the supervised ML methods, logistic regression (LR), artificial neural networks (ANN), support vector machines (SVM), and RF appear most frequently in the LSA literature (Lima et al. 2022).
The LSA studies mainly differ from each other based on the data sources utilized for the assessment, the geographic locations, the size of the assessed area (local or regional), the spatial resolution of the data, and the predictors (that is, the conditioning factors). The data used for the landslide investigations can be sourced from photogrammetric and remote sensing techniques-such as using passive (satellite and aerial optical images) or active sensors, for example, light detection and ranging (LiDAR) or synthetic-aperture radar (SAR) (Colesanti and Wasowski 2006;Jaboyedoff et al. 2010;Nava et al. 2022)-or can be obtained from existing geodatabases.
The accuracy assessment approaches are either qualitative (expert inspection) or quantitative based on various metrics, such as precision, recall, receiver operating characteristic (ROC) curve, the area under the ROC curve (AUC) value, Kappa index, F-1 score, overall accuracy 1 3 (OA), and so on, or both approaches at the same time. In LSA studies, highly accurate results can be obtained depending on the data quality and the suitability of the learning method. Merghadi et al. (2020) compared different ML algorithms for the LSA, and found that the RF method was more robust than the other tested methods. Can et al. (2021) obtained an AUC value of 0.96 in an area over 2700 km 2 using the decision-tree (DT) based extreme gradient boosting (XGBoost) method. Karakas et al. (2022) obtained an AUC value of 0.93 using the RF method with very high-resolution aerial photogrammetric datasets. But although very accurate results can be achieved, the LSA methods and the predictors have not yet been standardized, mainly due to geological and geographic differences, and spatially and temporally incomplete inventories.
Considering the conditioning factors for the FSA and LSA, topographic (for example, altitude, gradient slope, slope length factor, topographic wetness index (TWI), stream power index (SPI), plan curvature), hydrologic (for example, distance to rivers), land use and land cover (LULC), and geological (for example, soil drainage) parameters have often been utilized for both types of studies. In this study, before selecting the factors, the most frequently utilized parameters for the LSA and the FSA in the last 5 years were evaluated after filtering with the number of citations (> 15). The topographic derivatives (altitude, slope, aspect, curvature), lithology, distances to hydrologic elements (for example, rivers), infrastructure (for example, roads), and geological structures (for example, faults), and LULC were found the most frequently utilized factors for the LSA. For the FSA, topography, lithology, hydrology, and LULC were commonly used variables.
The peak ground acceleration (PGA) has been frequently used for the probabilistic assessment of seismic hazards in site selection and design of engineering structures. Erdik et al. (1999) assessed the probabilistic seismic hazard in Türkiye and its surrounding areas, and produced a PGA map according to specific return periods. Alpyurur and Lav (2022) evaluated the seismic hazard for the southwest of Türkiye and produced a new seismic susceptibility database.
Only few studies on MHSA appear in the literature. The majority of these are stepwise studies and combine univariate susceptibility maps with weighted overlay analysis or the AHP. Mukhopadhyay et al. (2016) used the multiple-criteria decision analysis (MCDA) for the MHSA considering coastal erosion, storm surge, sea-level rise, coastal floods, tsunamis, and earthquakes in India, and the weights were defined by experts involved in the study. Khatakho et al. (2021) applied the AHP for evaluating flood, landslide, earthquake, and urban fire hazards in Nepal and aggregated them again by using the AHP. Aksha et al. (2020) also utilized the AHP for MHSM production in Nepal. Rusk et al. (2022) used an ML method, the maximum entropy (Maxent), for producing the FSM, LSM, and a fire susceptibility map in the Hindu Kush Himalaya. In addition, recent studies by produced an MHSM using a new ensemble model named stepwise weight assessment ratio analysis (SWARA) utilizing adaptive neuro-fuzzy inference system (ANFIS) and grey wolf optimizer (GWO) for the LSA, FSA, and earthquake assessment in Iran by using a PGA map. Yanar et al. (2020) produced an MHSM for a part of Ankara, Türkiye by integrating an FSM and an LSM with the Mamdani FIS. In the present study, we integrate the FSM, LSM, and AI map for the Elazig region using a Mamdani FIS with membership functions defined using the data characteristics. Therefore, our study is the first example of combining the landslide, flood, and earthquake hazards for MHSA using a FIS.

Materials and Methods
In this section, we present the study area, the geospatial datasets employed in the study, and the methods utilized for the production of the LSM, FSM, AI, and the MHSM.

Study Area Characteristics
The study area is located in the southeastern part of Türkiye ( Fig. 1a) in Elazig Province and covers part of the Keban Dam Lake basin. It covers approximately 1010 km 2 and the altitudes obtained from the EU-DEM v1.1, published by the Copernicus Land Monitoring Service (CLMS 2021), range from 822 to 2146 m (Fig. 1b). A worldwide land cover map published by the European Space Agency (ESA-WorldCover 2020) with 75% overall accuracy was used to represent the environmental conditions as LULC classes, which consist of tree cover (1) (1.23%), shrubland (2) (0.07%), grassland (3) (39.20%), cropland (4) (27.94%), built-up area (5) (5.50%), bare/sparse vegetation (6) (15.32%), and permanent water bodies (7) (10.72%) (Fig. 1c). The average annual temperature of Elazig Province is 13.1 °C, and the average annual precipitation is 415.1 mm (MGM 2022). The region is in the East Anatolian Fault Zone (EAFZ), tectonically active, with high seismicity Karakas, Nefeslioglu, et al. 2021), and has historically been subject to frequent devastating earthquakes (AFAD 2020; METU 2020). The geological features consist of lithological units with weak shear strength characteristics (Fig. 1d). Avci and Sunkar (2018) investigated the relationship between lithological units and the distance to faults with landslides in Elazig and the neighboring province in the north, and observed that most of the landslides in the region were triggered by earthquakes.
The lithological units in the area (Fig. 1d) were obtained by digitizing the 1:100,000 scale geological map published by Keskin (2011) and Herece (2016) and are listed in Table 1 in the order of their age. The most common units are terraces (unconsolidated gravel, sand, silt, clay), the Kırgeçit formation (conglomerate, sandstone, siltstone, claystone, limestone), and basalt. The landslide inventory was provided by the General Directorate of Mineral Research and Exploration (GDMRE) (MTA 2022), Türkiye. The inventory is incomplete and covers an area of 1.42 km 2 , with the smallest and the largest landslide areas of 0.01 km 2 and 0.34 km 2 , respectively. The landslides in the region are mostly present in lithological units that are not firmly attached under the thick basalts (Herece 2016), and in the Kırgeçit formation.

Geospatial Datasets and Conditioning Factors
The types, sources, and scales/spatial resolutions of the data are summarized in Table 2. The topographic parameters were derived from the EU-DEM v1.1 tile E60N20 with a grid size of 25 m and a height accuracy of 7 m (CLMS  Table 1 1 3 2021). The distances to rivers and drainage channels were computed with data from the topographic vector geodatabase (TopoVT) of the General Directorate of Mapping (GDM), Türkiye, which was also partly updated by manual delineations from a recent orthoimage. The lithology and the distance to faults were used as geological parameters.
The topographic and geological features were computed using the ArcGIS software from ESRI and the open-source SAGA tool (Conrad et al. 2015). Slope is among the most important parameters for the FSA and the LSA. Steep slopes are more susceptible to landslides. Lower slopes and flat areas are more prone to flooding since the surface inundation level may increase in these areas (Sozer et al. 2018;Khatakho et al. 2021). The aspect is the line direction of the steepest descent and denoted in degrees clockwise from the north. The sunlight exposure times, and the freezing and thawing events that affect the decomposition and erosion of the material of the slopes facing different directions, can be related to the aspect (Dai and Lee 2002). The plan and profile curvatures are the second derivatives of the DEM and denote the magnitude of the change in slope and aspect. With respect to the FSA, while the plan curvature explains the flow acceleration and erosion/deposition rate, the profile curvature denotes the flow velocity variation (Kalantar et al. 2018). The TWI describes the hydrological conditions of the topography and is used to explain the water-saturated areas in a catchment. The SPI represents the erosional strength of running water and is related to the discharge to a catchment. The drainage density represents the proportion of the lengths of all streams and rivers to the total catchment area. Distances to permanent rivers and dry drainage channels were used to analyze the effect of rivers in the event of flooding by assuming that closer proximity increases susceptibility. Lithology is essential in landslide, flood, and earthquake hazard assessments. Rock type, permeability, and surface runoff are affected by lithological properties (Kalantar et al. 2018;Khatakho et al. 2021). Since the study area is in the East Anatolian Fault Zone and is prone to landslides triggered by earthquakes, the distance to fault parameter was included. The faults were digitized from the web portal of the GDMRE (MTA 2022). The different types of LULC in a region exhibit various influences in natural hazard events and were also used as conditioning factors.

Methodology
In this study, we produced univariate susceptibility maps for landslides, floods, and earthquakes, and integrated all three maps into a single multi-hazard map as explained below. We applied a data-driven method (the RF) for the LSA. For the flood and multi-hazard susceptibility assessments, expertbased methods were employed and the required inputs (expert opinion) were provided by a coauthor of this article, Gokceoglu.

Production of Univariate Susceptibility Maps
The RF algorithm (Breiman 2001), which is frequently used for regression and classification problems, was utilized for the LSA. Consisting of multiple DTs with different subsamples, the RF is relatively immune to outliers and noise in the data. For hyperparameter optimization, the Rand-omizedSearchCV method (Bergstra and Bengio 2012) was employed. The Gini index, which measures the variance in the classification, was preferred here. The number of estimators and the maximum depth allowed in a tree were 418 and 27, respectively. The bootstrap was used as the sampling strategy for training each tree. A balance of class weights was selected considering the imbalance problem in the classification and the parameter is usually introduced to give focus on the minority class.
The predictors used for the LSA were altitude, slope, aspect, plan and profile curvatures, lithology, drainage density, distance to fault, SPI, and TWI. The features were selected based on previous LSA studies performed for the area (Can et al. 2021;Karakas, Kocaman, et al. 2022;Karakas, Nefeslioglu, et al. 2021;Karakas et al. 2022), and the expert opinion. The landslide samples used for the training were selected from the pixels in the inventory. The landslide/ non-landslide classes include 2264 and 3396 pixels, respectively, yielding to a total of 56,600 pixels for all features. An 80/20 ratio was used to split the training and test samples. The AUC value, precision, recall, and F-1 score metrics were used for validation with test samples. In addition, the SHapley Additive exPlanations (SHAP) methodology (Lundberg and Lee 2017) was used to investigate the predictive performances of the input features.
The M-AHP was applied for the FSA due to the lack of learning data. The M-AHP eliminates the expert subjectivity in pairwise factor comparisons in AHP (Nefeslioglu et al. 2013). A normalized factor score difference matrix is prepared with the maximum weights given and a factor comparison matrix is generated by considering the modified importance value chart to ensure consistency. The factor percentage importance distributions at the decision points are computed and linear distances between the normalized factor score and the decision points are measured. A new weight vector is computed for each grid cell. Seven conditioning factors were considered in the FSA (altitude, slope, lithology, TWI, LULC, distance to permanent rivers, and distance to dry drainage channels) based on literature analysis and expert opinion. The class weights were defined by the expert. Both the LSA and the FSA methods were implemented in a Python programming environment and the workflows are presented in Fig. 2.
As stated by Rathje et al. (1998), it is often useful in earthquake engineering practice to characterize the frequency content of an earthquake ground motion with a single parameter. However, it is still a difficult problem to illustrate earthquake hazard with a single map when preparing the regional multi-hazard susceptibility maps. Due to the significant role of local site conditions on earthquake shaking, the effect of the current seismic code provisions is described through appropriate elastic design spectra based on different site categories (Kotha et al. 2018).
The main parameter proposed for soil categorization is the Vs30 (Borcherdt and Glassmoyer 1992). Despite its common use, there is no universal agreement that this is a valid proxy for seismic amplification (Castellaro et al. 2008). Chousianitis et al. (2018) emphasized that not all of the main properties of a ground motion can be captured through a single parameter and proposed several different engineering parameters to reflect one or more ground motion characteristics concurrently. However, there is still an inadequacy of standard methods for site classification to evaluate the expected amount of peak horizontal acceleration amplification (Del Gaudio et al. 2019).
Another parameter used in ground motion prediction equation is the horizontal component of cumulative absolute velocity. Campbell and Bozorgnia (2010) used cumulative absolute velocity employed as an index to indicate the possible onset of structural damage to nuclear power plant facilities and liquefaction of saturated soils. Foulser-Piggot and Goda (2015) stated that Arias intensity (AI) and cumulative absolute velocity are ground motion measures that have been found to be well suited to application in a number of problems in earthquake engineering, and these parameters reflect multiple characteristics of the ground motion (for example, amplitude and duration), despite being scalar measures. Due to the availability of reliable data, Arias intensity is employed in the present study to classify the study area relatively in terms of earthquake effects, that is, the Arias intensity map classifies the area relatively in terms of earthquake impact. This classification does not express the absolute effect of a possible earthquake, but rather demonstrates the relative influence in the study area. Therefore, the Arias intensity used to show the earthquake effect in the production of the multi-hazard susceptibility map is expressed by the susceptibility map. The Arias intensity proposed by Arturo Arias (1970) expresses the energy discharge considering the earthquake shaking duration and time-dependent variation of the frequency content (Kayen and Mitchell 1997). The intensity of the shaking is determined by measuring the acceleration of transient seismic waves and the summation of the horizontal and vertical components of the acceleration record (Eq. 1). In Eq. 1, g is the acceleration of gravity, t is time, and td is the total recording length. Among the ground motion parameters that measure ground shaking, AI is well related to the landslide distribution and density (Keefer 1984;Wang et al. 2007;Lee et al. 2008;Karpouza et al. 2021), and a reliable In this study, the AI was used to produce the earthquake susceptibility map with the inverse distance weighting (IDW) interpolation based on the records of 25 earthquake stations in the study area and its surroundings provided by the AFAD earthquake database.

Multi-Hazard Susceptibility Assessment with the Mamdani-Fuzzy Inference System
Mamdani FIS (Mamdani and Assilian 1975) aims at creating a control system using linguistic rule sets based on expert opinion. It is an intuitive and easy-to-use method for solving complex and nonlinear problems. Here, the LSM, the FSM, and the AI maps were aggregated for producing the MHSM using Matlab Fuzzy Logic Toolbox. The Mamdani FIS method consists of four main steps-fuzzification, rule evaluation, aggregation, and defuzzification. (1) In the fuzzification step, the membership classes of crisp inputs are evaluated by the expert. Here, three membership classes-low, moderate, and high-were defined for each input (that is, LSM, FSM, and AI map) and the output was classified in five linguistic classes (very low, low, moderate, high, and very high). The membership boundaries were defined considering the severity of each hazard in the area based on the statistics (for example, loss of lives and damages to properties) in the AFAD report (AFAD 2020).The most harmful hazard type has been earthquakes, followed by landslides and floods. The histograms of the univariate susceptibility maps were also analyzed for probability distribution of each susceptibility map for the purpose of class boundary definition. Triangular functions define each membership (Fig. 3), with the x and y axes indicating the probability levels and membership degrees, respectively. Figure 3 shows a susceptibility level expressed with a probability between 0 and 1 (normalized) that falls within two membership classes (for example, 0.4 in the landslide membership function belongs to the moderate and high susceptibility levels at different degrees).
Linguistic if-then rules were defined by the expert (Gokceoglu, coauthor of this article) to be evaluated in the second step (Table 3). The 27 rules reflect the opinion of the expert, Fig. 3 The membership functions defined in the Mamdani fuzzy inference system (FIS). MHS multi-hazard susceptibility who has long experience in the region. In the third step, the aggregation is carried out and a fuzzy output of the model is obtained using the maximum operator (Osna et al. 2014). In the final step (defuzzification), the fuzzy output was transformed into crisp values using the centroid technique as the fuzzy output includes the results of multiple rules for every susceptibility level, which belongs to two membership classes each time. The fuzzy output is in fact an area, and the centroid technique computes the center of this area to obtain the crisp value. The crisp outputs were geocoded for the MHSM production.
The LSM, FSM, AI map, and the MHSM are presented in Figs. 5a−d. The AUC value obtained from the RF is 0.98, which indicates high classification accuracy confirmed by the OA (92.3%) and F-1 score (0.91) as well. The feature importance results indicate that slope, lithology, and altitude are the most predictive features, and slope and altitude values have a positive correlation with the susceptibility levels, whereas smaller distance values to faults yield an increase in susceptibility level. The similarity between the high landslide susceptibility levels (Fig. 5a) and higher slopes (see Fig. 4a) is visually apparent.
The FSM (Fig. 5b) indicates that the areas around the Keban Dam Lake and in the center of Elazig City are highly susceptible. In addition, Elazig Airport was classified as highly susceptible to flooding. The AI map (Fig. 5c) indicates higher intensity values in the southern parts of the

If (LS is high) and (FS is high) and (AI is high) then (MHSL is very_high) 2 If (LS is high) and (FS is high) and (AI is moderate) then (MHSL is very_high) 3 If (LS is high) and (FSis high) and (AI is low) then (MHSL is high) 4 If (LS is high) and (FS is moderate) and (AI is high) then (MHSL is very_high) 5 If (LS is high) and (FS is low) and (AI is high) then (MHSL is high) 6
If (LS is high) and (FS is moderate) and (AI is moderate) then (MHSL is high) 7 If (LS is high) and (FS is low) and (AI is low) then (MHSL is moderate) 8 If (LS is high) and (FS is low) and (AI is moderate) then (MHSL is moderate) 9 If (LS is high) and (FS is moderate) and (AI is low) then (MHSL is moderate) 10 If (LS is moderate) and (FS is high) and (AI is high) then (MHSL is very_high) 11 If (LS is moderate) and (FS is high) and (AI is moderate) then (MHSL is high) 12 If (LS is moderate) and (FS is high) and (AI is low) then (MHSL is moderate) 13 If (LS is moderate) and (FS is moderate) and (AI is high) then (MHSL is high) 14 If (LS is moderate) and (FS is low) and (AI is high) then (MHSL is moderate) 15 If (LS is moderate) and (FS is moderate) and (AI is moderate) then (MHSL is moderate) 16 If (LS is moderate) and (FS is low) and (AI is low) then (MHSL is low) 17 If (LS is moderate) and (FS is low) and (AI is moderate) then (MHSL is moderate) 18 If (LS is moderate) and (FS is moderate) and (AI is low) then (MHSL is moderate) 19 If (LS is low) and (FS is high) and (AI is high) then (MHSL is high) 20 If (LS is low) and (FS is high) and (AI is moderate) then (MHSL is moderate) 21 If

(LS is low) and (FS is high) and (AI is low) then (MHSL is moderate) 22 If (LS is low) and (FS is moderate) and (AI is high) then (MHSL is moderate) 23 If (LS is low) and (FS is low) and (AI is high) then (MHSL is moderate) 24 If (LS is low) and (FS is moderate) and (AI is moderate) then (MHSL is low) 25 If (LS is low) and (FS is low) and (AI is low) then (MHSL is very_low) 26 If (LS is low) and (FS is low) and (AI is moderate) then (MHSL is low) 27
If (LS is low) and (FS is moderate) and (AI is low) then (MHSL is low)  Keefer and Wilson (1989). The gradient MHSM obtained from the Mamdani FIS (Fig. 5d) shows that the southern parts of the study area have very high susceptibility levels mainly due to the higher AI and landslide susceptibility values. Figure 5 shows that the landslide susceptibility is very high in most parts of the study area, although the flood susceptibility is lower in areas with higher slope values.

Discussion
In this study, a new approach to combine landslide, flood, and earthquake hazards with different machine learning and expert-based methods was proposed for multi-hazard assessments at a regional scale. The RF method provided high performance (OA = 92.3%) for the LSA as expected.
Although the number of landslides in the inventory is low, and they cover a small part of the study area, the accuracy was in line with the results of Karakas et al. 2022), thus indicating the adequacy of the landslide inventory for the study purposes. Yet, further investigations and understanding of the relationships between the landslides and earthquakes in the region are needed similar to the study of Karakas and others (Karakas, Nefeslioglu, et al. 2021a, b), in which the landslides triggered by the 2020 Mw. 6.8 Elazig earthquake were reported.
The FSM obtained from the M-AHP proved the usability of the method based on expert validation. This finding is in line with the study of Sozer et al. (2018), which applied the M-AHP in Ankara, Türkiye for the FSM at the regional scale. Hammami et al. (2019), Swain et al. (2020), andKhatokho et al. (2021) also used the AHP for the FSA and confirmed the suitability of this approach. But since the M-AHP method eliminates expert subjectivity and creates a weight matrix for The conditioning factors used for the LSA and the FSA were selected after a review of the studies carried out between 2017 and 2021 and considering the study area characteristics. Since high quality maps could be obtained for both the FSA and the LSA, the selected factors were found suitable. Although the factors may involve multicollinearity, this issue was not considered here as the DT-based methods are relatively immune to this issue (Piramithu 2008;Can et al. 2021;Karakas et al. 2022). Among the factors used for the LSA, the most predictive features were slope, lithology, and altitude. In addition, the EU-DEM was found suitable for the study purposes.
The AI, which was utilized to produce the earthquake susceptibility map, has been frequently preferred in seismic hazard assessment studies in recent years (Costanzo 2018;Skilodimou et al. 2019;Karpouza et al. 2021;Gupta and Satyam 2022). As the parameter used to assess the earthquake shaking degree for triggering landslides, it is widely used in the literature and was preferred here.
The use of the Mamdani FIS for the MHSA is a major contribution of this study. The production of a multi-hazard susceptibility map (MHSM) is still an extremely complex problem due to regional complexities and the lack of ground truth data that prevents the use of data-driven machine learning methods. In the literature, the univariate susceptibility maps have usually been integrated with raster arithmetic operations or basic spatial analysis methods (Pourghasemi et al. 2019;Pourghasemi et al. 2020;Pouyan et al. 2021;Bordbar et al. 2022;Ullah et al. 2022), such as the cumulation of the individual hazard levels per pixel. The use of Mamdani FIS by creating rules for each hazard has increased the sophistication level and usability of the final maps. Yet, the method cannot be carried out simply in any desktop geographic information system (GIS) software but its application requires specific software and a certain level of expertise. In addition, the membership functions and the rules must be defined by an expert who is familiar with the area as there are no widely accepted procedures or standards. Further applications and case studies would help to increase the knowledge, which may eventually lead to standardization of the method. Figure 5 shows a large portion of the study area is susceptible to landslides, floods, and earthquakes. In addition, a MHSM can be more useful than univariate susceptibility maps for site selection efforts for highways, railways, settlements, and so on.

Conclusion
In this study, a novel approach for the combined assessment of landslide, flood, and earthquake hazard susceptibility was proposed and evaluated in a region in the Elazig Province of Türkiye, which is tectonically active and vulnerable to multi-hazards that affect each other. The RF was powerful for producing an accurate LSM (OA = 92.3%) even when the inventory is incomplete. The FSA was carried out with the M-AHP as no flood inventory was available. In addition, the AI values were used considering the earthquake susceptibility of the region. The univariate susceptibility maps were combined with the Mamdani-FIS to obtain the MHSM, which is the first example in the literature. The FSM and the MHSM were validated visually both in 2D and 3D views by the expert. The results exhibited high performance of the proposed methodology for the MHSA and can be recommended for similar efforts. The produced maps can be used by local authorities for identifying areas prone to multi-hazards, planning sustainable land use, site selection of engineering structures, and efficient disaster management. Yet, regional MHSA is still a highly difficult problem, and it is almost impossible to solve it with data-driven approaches. For this reason, an expert-based fuzzy inference system (FIS) was preferred. Although successful and promising results were obtained in the present study, there is still a long and difficult way to go in producing multi-hazard susceptibility maps. Multi-hazard susceptibility assessment studies need to be carried out in different regions and with different methods to increase the accuracy and reliability of such assessments. The application of FIS can be complex as it requires specific tools and expert involvement for defining the membership functions and the rules. The future work of the study will include the evaluation of the proposed method in the neigboring basins in terms of prediction performance and usability. In addition, risk assessments by considering the vulnerable elements are also planned.