2019 | 72/1 | 81–92 | 5 Figs. | 4 Tabs. | www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society 1. INTRODUCTION Plenty of landslide susceptibility models have been developed by using GIS technology within the past three decades. These mod- els have been prepared for various purposes and stakeholders, always intended to serve as an assessment of where landslides are to be expected in future. Landslide susceptibility models are cat- egorized according to various criteria, and always denoted as: quantitative or qualitative, direct or indirect, heuristic, probabi- listic or deterministic, for differentiated or undifferentiated types of landslide mechanism. Landslide susceptibility models and maps are being prepared on global (e.g. LIN et al., 2017), continental (e.g. GÜNTHER et al., 2013a), national (TRIGILA et al., 2013), (KOMAC & RIBICIC, 2006), (BĂLTEANU et al., 2010), (GLADE at al., 2012) and other scales, but mostly on the regional scale (e.g. KOMAC, 2005). Regional landslide mapping is largely accepted among scientists, bearing in mind that the spatial distribution of land- slide conditioning parameters changes over larger areas. In this Preliminary regional landslide susceptibility assessment using limited data Igor Peshevski1, Milorad Jovanovski1, Biljana Abolmasov2, Jovan Papic1, Uroš Đurić3, Miloš Marjanović2, Ubydul Haque4 and Natasha Nedelkovska5 1 University Ss. Cyril and Methodius, Faculty of Civil Engineering, Blvd. Partizanski odredi No.24, 1000 Skopje, Republic of Macedonia; (pesevski@gf.ukim.edu.mk) 2 University of Belgrade, Faculty of Mining and Geology, Đušina 7, 11000 Belgrade, Serbia 3 University of Belgrade, Faculty of Civil Engineering, Bulevar Kralja Aleksandra 73/I, 11000 Belgrade, Serbia 4 University of North Texas Health Science Center, Department of Biostatistics and Epidemiology, Fort Worth, TX 76107, United States of America 5 Geohydroconsulting Ltd, Manapo 7, 1000 Skopje, Republic of Macedonia doi: 10.4154/gc.2019.03 Abstract In this paper a heuristic approach for preliminary regional landslide susceptibility assessment using limited amount of data is presented. It is called arbitrary polynomial method and takes in- to account 5 landslide conditioning parameters: lithology, slope inclination, average annual rain- fall, land use and maximum expected seismic intensity. According to the method, in the first stage, a gradation is performed for each of the carefully selected conditioning parameters by assigning so called rating points to the grid cells on which the region is divided. Values of the rating points vary from 0 to 3 and depend on the parameter’s character and importance for land- slide development within the region of interest. A so called Total Landslide Susceptibility Rating (TLSR) model is obtained by summing the individual rating points of each parameter and divid- ing the region into five susceptibility zones according to Jenks natural breaks classification. Ver- ification of the TLSR model is then performed by overlaying the landslide inventory map of the selected region over the prepared susceptibility map. The sensitivity of the model can be addi- tionally tested by multiplying the conditioning parameter’s rating points by sensitivity coefficients. In this way, additional landslide susceptibility models are obtained, named Weighted Total Land- slide Susceptibility Rating (WTLSR) models. As a practical example of the method, two TLSR models are presented here for the Polog region in Republic of Macedonia, for return periods of maximum expected seismic intensity for 100 and 500 years. With over 74% of mapped landslides falling in zones of high and very high susceptibility, the results are considered satisfactory for regional scale landslide modelling and are comparable with more advanced quantitative meth- ods. Additional WTLSR models were prepared, and their correlation identified the best model. The presented approach is considered to be very convenient for conducting preliminary region- al landslide susceptibility assessments with the ability to fine-tune the results. Due to its simplic- ity, it can be applied to additional landslide conditioning parameters other than the one present- ed in the paper, depending on the region of interest and available data sources. It is especially practical for use in developing countries, where various organizational, technical and economic constraints prevent application of more advanced data driven methods. Limitations and restric- tions of the approach are also discussed. context, prominent landslide researchers (SOETERS & VAN WESTEN 1996; MIHALIĆ, 1998; HERVÁS, 2007; FELL et al., 2008; POSCH-TRÖZMÜLLER & BÄK, 2013; COROMI- NAS et al., 2014) have presented ideas of the minimal require- ments for performing landslide inventory, susceptibility, hazard and risk zoning in relation to the scale of mapping. One of the most applied methods in regional landslide mod- elling is the qualitative experience (knowledge) based – heuristic approach, which is considered a rather controversial topic in land- slide assessment (BARREDO et al., 2000; ARDIZZONE et al., 2002; ERCANOGLU et al., 2008). It is based on the use of an index or parameter maps (SOETERS & VAN WESTEN, 1996), and it is accepted as being applicable only for the preliminary levels of research (GÜNTHER et al., 2013b) if not combined with more advanced approaches (MARJANOVIĆ, 2014). Alterna- tively, according to DRAGICEVIĆ et al., (2012), heuristic meth- ods are empirical and subject to various levels of uncertainty, but Article history: Manuscript received February 13, 2018 Revised manuscript accepted October 26, 2018 Available online February 15, 2019 Keywords: regional landslide susceptibility, arbitrary-polynomial method, conditioning parameters, sensitivity, Polog region G eo lo gi a C ro at ic a Geologia Croatica 72/182 have proven to be a reliable and cost-effective method that allows detailed and comparable assessments of landslide susceptibility. There are many different possibilities to heuristically com- bine various conditioning factors and produce a susceptibility map at any scale, but despite this, it is rarely the case that a par- ticular model can be applied that gives the same results for dif- ferent study area. Also, using a large number of conditioning fac- tors in a particular analysis is not a guarantee for successful susceptibility modelling. Factor selection will always be subjec- tive, so the assessment of regional factors by the expert conduct- ing the analysis is of crucial importance. Statistical methods, which are generally considered as most appropriate for landslide susceptibility mapping at regional scales because they are objective, reproducible and easily updatable (NARANJO et al., 1994; HE & BEIGHLEY, 2008; VAN WESTEN et al., 2006), also have some limitations when com- pared to heuristic approaches. In this regard, THIERY et al. (2007) reveal major limitations that concern all statistical meth- ods. These are as follows: (1) their significant sensitivity to the quality and accuracy of the input thematic data and particularly to the landslide inventory used to train the model, (2) the absence of expert opinion which may result in a satisfactory statistical output in terms of degree of fit, but may not be realistic in terms of physical meaning, (3) the number of landslide events to incor- porate in the statistical model should be appropriate to the size of the study area, resulting in increased data requirements in large study areas, and (4) the use of oversimplified factors assumed to control or trigger landslides, by only considering information that is relatively easily mapped or derivable from a digital elevation model (DEM). Furthermore, GIS based statistical landslide sus- ceptibility assessment is often more focused on the tool than on the input data and frequently involves an extreme simplification of landslide controlling factors (VAN WESTEN et al., 2006). Therefore, in cases where a limited amount of data is avail- able, even with a simple heuristic approach such as that presented here, and for example similar ones, such as that presented by ZHU et al., (2014), it is possible to obtain satisfactory results, particularly at the regional level. Landslide susceptibility models sensitivity is also the fre- quent subject of scientific discussion. Despite the large number of recent advances and developments in landslide susceptibility modelling, there is still a lack of studies focusing on specific as- pects of the model sensitivity (CATANI et al., 2013). It is sug- gested by ROSSI et al., (2010) that optimal susceptibility predic- tions might be obtained through the combination of suitable basic landslide susceptibility models generated by different methods rather than by the application of a single prediction. In this con- text, most recent research on landslide susceptibility model sen- sitivity usually comprises a comparison of different models of the same study area. Examples of such correlations can be found in numerous publications (see CHANG & KIM, 2004; OH et al., 2009; MAGLIULO et al., 2009; ROZOS et al., 2010; MARJA- NOVIĆ, 2010; CONSTANTIN et al., 2011; AKGUN, 2012; DEVKOTA et al., 2013; LEOPOLD et al., 2013; KAVZOGLU et al., 2014; DEMIR et al., 2015; RAMESH & ANBAZHAGAN, 2015; SHAHABI et al., 2015; WANG et al., 2016), to name a few. In some studies, for example in LEE et al., (2012), relatively com- plex algorithms are being engaged to check the sensitivity of a landslide susceptibility model. Extraction of a particular land- slide conditioning parameter and repetition of the validation pro- cess is one approach to test model sensitivity (LEE & TALIB, 2005; GUZZETTI et al., 2006). Other approaches for testing the sensitivity of a particular landslide susceptibility model (and not by comparison with other models) or evaluating its quality, reli- ability and prediction skill, are presented in (GUZZETTI et al., 2006; FRATTINI et al., 2010; ROY & OBERKAMPF, 2011; PETSCHKO et al., 2014; FEIZIZADEH & BLASCHKE, 2014). Regardless of the approach used in the model preparation, it can be stated that every modelling result is highly dependent on the quality of the input data, the assumptions made to set up the model design and the selection of the appropriate performance metrics (PETSCHKO et al. 2014). Furthermore, the evaluation of reliability consists of changing some of the model integral parts and obtaining numerous variations of the outcome susceptibility maps, up to 350 in GUZZETTI et al., (2006). CATANI et al., (2013) suggest that careful sensitivity analysis which takes into consideration both conventional and innovative tools should al- ways be performed before producing final susceptibility maps at all levels and scales. When there are various constraints for some regions of inter- est such as: different scales of landslide conditioning parameters maps; limited amounts of geotechnical data of historical land- slides; low numbers of registered landslides in databases; no ac- cess to multi temporal aircraft imagery; no resources to perform field surveys, etc., the heuristic approach is possibly the most suit- able modelling option for the researchers. In this context, in order to present the possibility of producing respective preliminary landslide susceptibility models of a region by using a limited amount of datasets due to the previously mentioned constraints, a specific heuristic approach called the arbitrary-polynomial method is presented here. The approach takes into account the findings of past research (e.g. by ABOLMASOV & STOJKOV 1994, ABOLMASOV & STOJKOV 1995; JOVANOVSKI et al., 2013a; PESHEVSKI et al., 2013; PESHEVSKI, 2015a). The Po- log region in the Republic of Macedonia was used as a case study area. The produced models are further tested for sensitivity, with the aim of presenting that fine tuning of heuristic landslide sus- ceptibility maps obtained from poor datasets is feasible and is a practical tool in landslide susceptibility zonation practice. 2. THE ARBITRARY – POLYNOMIAL METHOD The approach used in this paper is named the arbitrary polyno- mial method, and engages five landslide conditioning parameters: lithology, slope inclination, average annual rainfall, maximum expected seismic intensity and land use. These parameters are selected as the most representative for landslide susceptibility as- sessment, and their selection is based on analysis of limited his- torical data of landslides that had a destructive effect on the in- frastructure in the Republic of Macedonia (PESHEVSKI et al., 2013; JOVANOVSKI et al., 2013b). It can be argued that rainfall and seismicity conditions are usually considered as triggering factors for landslide occurrence and in most cases related to landslide hazard rather than suscep- tibility assessment. In our case, because of their importance for the analyzed region, they were considered as landslide condition- ing parameters. For the same reasons, numerous landslide sus- ceptibility studies exist where earthquakes and/or /rainfall have been considered as landslide conditioning parameters. These studies have been performed on both local and regional scales (e.g. WANG et al., 2015; ERENER & DÜZGÜN, 2010; ILAN- LOO, 2011; CHALKIAS et al., 2014; HE & BEIGHLEY, 2008; SHIBIAO et al. 2013; WANG et al. 2015; PAMELA et al. 2018; SHAHABI & HASHIM, 2015, UMAR et al., 2014). There are also some case studies where both of these conditioning param- G eologia C roatica Peshevski et al.: Preliminary regional landslide susceptibility assessment using limited data. 83 eters are taken into account when assessing global landslide sus- ceptibility (e.g. LIN et al. 2017). According to the method, a heuristic gradation is performed by assigning rating points to the landslide conditioning param- eters. The value of the ratings is defined by expert assessment on the basis of particular conditioning factor importance for land- slide development in a given region. In order to investigate the sensitivity to changes of these rat- ings, analysis is further extended by adding so called sensitivity coefficients for each of the conditioning parameters. The general framework of the method is presented in Fig. 1 and details are given in the following subsections. The procedure is quite different to the statistical correlation of landslides and conditioning parameters which is the most com- mon approach in landslide susceptibility mapping practice. It is an alternative which brings expert judgment to a central role, and not translating the problem to a more simple or complex statisti- cal exercise. 2.1. Rating maps for landslide conditioning parameters 2.1.1. Lithology (LT-R) Analysis of reports and thematic engineering-geological maps of historic landslides in the Republic of Macedonia showed that bed- rock and superficial geology are among the most important fac- tors for development of landslide processes (PESHEVSKI et al., 2013). This is usually the case with such diverse lithologies as those present in the Western Balkans. In our case study, when preparing a rating map for lithology, all known data for the con- ditions and properties of the rock masses of known landslides in Macedonia were considered. After analysis of all the available thematic geological, engineering geological, geotechnical maps and cross-sections, and geotechnical reports, lithological units were grouped into six engineering-geological groups with simi- lar physical-mechanical behavior. The rating values were as- signed to all lithological units, wherein higher values indicate a higher potential for landslide occurrence. 2.1.2. Slope Inclination (SI-R) The second significant conditional parameter is slope inclination, which is an essential element for every landslide susceptibility analysis. It is closely related to all gravitational processes, and rates of material movement are largely dependent on it. However, in the sense of measuring its potential to develop landslides, rat- ing slope inclination is not a trivial problem, especially on re- gional scale landslide susceptibility and for models that should encompass different types of landslide mechanisms. In general, a more reliable connection between slope inclination and sliding processes can be found for rock fall and shallow landslides, while for deep seated landslides and especially those in the presence of groundwater, it is a very unpredictable property. Each suscepti- bility analysis needs to consider slope inclination, with gradation done in the most logical manner, and in relation to the region of interest. Slope inclination as a parameter can be appropriately obtained from existing Digital Elevation Models (DEM) of the region of interest. A polynomial interpolation is a convenient way to define ratings for the slope inclination parameter in order to encompass all types of landslides. Minimal rating values are as- signed even at the lowest value of 3° inclination (translational and rotational landslides in Pliocene sediments on relatively flat ter- rain), and for all terrain over 25° inclination (mostly rock falls and rockslides), the highest rating of 3 is assigned. This type of inter- polation was used for the analyzed case study region. In more detailed statistical models, the slope inclination is treated sepa- rately according to the slope movement types (landslides, rock- falls, debris flows, etc). Since the presented approach is intended to give a general presentation of a region’s susceptibility, the in- terpolation technique is considered as very practical. In cases where detailed databases on the slope movement type exist, dif- ferent interpolation curves should be applied for each of the slid- ing mechanisms analysed. 2.1.3. Rainfall (AP-R) Rainfall is governed by many factors in a region, especially ele- vation, mountain range orientation, climate zone, etc. In this con- text, severe rainstorms or periods of prolonged rainfall have been reported to trigger many shallow, or deep-seated slides and debris flows worldwide, with devastating consequences. R. Macedonia is no exception. In particular, the Gradot, Velebrdo landslides (JOVANOVSKI et al., 2013b), and more recent ones, (Probistip, Makedonska Kamenica, Delcevo and Tetovo in 2015, 2016, 2017, 2018 respectively) are clear examples of such scenarios. Apart from the economic loss, some of these landslides have been life- threatening and human losses have been reported (PESHEVSKI et al., 2017). In the absence of more precise data from particular rain gauge stations or correlations of rainfall quantities and land- slide movement rates, rating values for this conditional parameter can be assigned arbitrarily according to, for example, average an- nual rainfall records. 2.1.4. Seismic intensity (EI-R) Earthquakes (considered as a conditioning factor) are also one of the main precursors of landslides. Landslides in Veles, Skopje- Sopishte, Berovo, and Pehcevo in R. Macedonia have been re- ported to occur after an earthquake between 1900–2015 (PE- SHEVSKI, 2015a). The earthquake importance according to the approach presented herein is taken into consideration through available maps of the maximum expected seismic intensity of the country. The Medvedev–Sponheuer–Karnik (MSK) scale is used as a reference in the presented case study. If other types of seis- mic data are available for a particular region, it can be rated in a similar manner. Figure 1. Flow-chart of the arbitrary polynomial method. G eo lo gi a C ro at ic a Geologia Croatica 72/184 2.1.5. Land Use (LU-R) Land use ratings are assigned by comparison of the available landslide inventory of the region of interest. In the particular case presented here, CORINE Land Cover of the European Environ- mental Agency, ver. 2006 was used as reference. A landslide in- ventory map is overlain? with the land use map to determine in which land cover units landslides mainly occur. Based on this, appropriate ratings can then be applied. 2.2. Calculation of Total and Weighted-Total Landslide Susceptibility Ratings Each specific combination of conditioning parameters is con- nected with a different susceptibility to sliding. The sum of the individual ratings for each conditional parameter in a grid cell (for regional scale is suggested grid cell size of 100x100m) gives its Total Landslide Susceptibility Rating (TLSR). The maximum possible theoretical value of the TLSR is 10.0 and the minimum possible value is 0.3. The following relationship is adopted: TLSR100-500 = LT-R + SI-R + AP-R + EI-R + LU-R (1) Where: TLSR – total landslide susceptibility rating (for return period of maximum expected seismic (earthquake) intensity for 100 and/or 500 years) LT-R – value of rating for lithological type SI-R – value of rating for slope inclination AP-R – value of rating for annual precipitation EI-R – value of rating for maximum expected seismic (earth- quake) intensity (for return period of 100 and/or 500 years) LU-R – value of rating for land use After calculation of the TLSR value for each grid cell, the terrain is divided into 5 classes of landslide susceptibility using Jenks natural breaks classification (very low, low, medium, high, very high). By overlaying the generated TLSR model map with the available inventory map of landslides of a particular region, the validation process can then be conducted. In order to investigate the sensitivity of the basic landslide susceptibility models (TLSR models) obtained, each conditional parameter rating can be additionally weighted by so called sen- sitivity coefficients. In this way, a Weighted Total Landslide Sus- ceptibility Rating (WTLSR) is obtained. A sensitivity coefficient of 0.3 is assigned to the conditional parameter with the highest preference, and 0.175 for the remaining four (these all sum up to 1). This is repeated four times, wherein each time a different conditional parameter receives the 0.3 coefficient. Other propor- tions/values of the sensitivity coefficients can also be tested (e.g 0.35 vs. 0.1625, 0.4 vs. 0.15). In our case study, other sensitivity coefficients gave less realistic outcomes. It should be noted that depending on the analyzed region and conditioning parameters used, the proportions and values of the ratings and sensitivity co- efficients need to be carefully selected. This selection should be done by an expert team and then checked by running several combinations of ratings and sensitivity coefficients. It can be ar- gued that this approach translates the susceptibility analysis into a trial and error exercise, making it highly subjective. However, with careful analysis of all the available data, weightings can be reasonably well selected, which will ultimately lead to produc- tion of a satisfactory and representative preliminary landslide susceptibility map of a particular region. Following this procedure, 8 additional variants of the basic TLSR models were prepared. The maximum possible value of the WTLSR that can be obtained with the proposed sensitivity coefficients is 2.125, whereas the minimum possible value of the WTLSR is 0.0525. Equations for calculation of the WTLSR used to test the TLSR models sensitivity are as follows: WLT-TLSR100 = LT-R*0.3 + SI-R*0.175 + AP-R*0.175 + EI-R*0.175 + LU-R*0.175 (2) WLT-TLSR500 = LT-R*0.3 + SI-R*0.175 + AP-R*0.175 + EI-R*0.175 + LU-R*0.175 (3) WSI-TLSR100 = LT-R*0.175 + SI-R*0.3 + AP-R*0.175 + EI-R*0.175 + LU-R*0.175 (4) WSI-TLSR500 = LT-R*0.175 + SI-R*0.3 + AP-R*0.175 + EI-R*0.175 + LU-R*0.175 (5) WAP-TLSR100 = LT-R*0.175 + SI-R*0.175 + AP-R*0.3 + EI-R*0.175 + LU-R*0.175 (6) WAP-TLSR500 = LT-R*0.175 + SI-R*0.175 + AP-R*0.3 + EI-R*0.175 + LU-R*0.175 (7) WLU-TLSR100 = LT-R*0.175 + SI-R*0.175 + AP-R*0.175 + EI-R*0.175 + LU-R*0.3 (8) WLU-TLSR500 = LT-R*0.175 + SI-R*0.175 + AP-R*0.175 + EI-R*0.175 + LU-R*0.3 (9) The validation process of these 8 additional WTLSR models is performed in the same manner as for the TLSR models using the landslide inventory maps. In the final stage, by detailed performance analysis of all maps generated by the models, an optimal one can be chosen as being the most representative of the particular region of in- terest. 3. LANDSLIDE SUSCEPTIBILITY MODELLING FOR THE POLOG REGION IN THE REPUBLIC OF MACEDONIA 3.1. Study area The Polog region is located in the northwest part of the Republic of Macedonia (Fig. 2). This region covers ~2420 km2 including the densely populated towns of Tetovo and Gostivar (parts of which were developed on rugged hilly terrain). It includes many villages on the steep Mt. Shar Planina, and important infrastruc- ture including railways, a highway, well developed network of regional and local roads (mostly in the mountains), ski centers, and a very important hydro energy system consisting of 130 km of water distribution channels accompanied by 167 km of service roads. In a geological context, the study area belongs to a larger re- gional tectonic unit called the Western Macedonian Zone (WMZ). In this unit rock masses from the Palaeozoic, Mesozoic, Pliocene and Quaternary periods are represented. Igneous rock masses include granodiorites, granites, diorites, rhyolites, ser- pentinites, gabbros, diabases etc. The Palaeozoic is represented by a thick complex of metamorphic rocks, rarely igneous rocks. The rocks from the Devonian age are the commonly occurring ones in the area, and here belong to the phyllitic schists, meta- conglomerates, metasandstones, quartzites, quartz-chlorite schists, carbonate schists and marbles. It is important to note that most landslides in the study area have been reported to occur at the contact of the weak schist type rocks and the soil debris which covers them. G eologia C roatica Peshevski et al.: Preliminary regional landslide susceptibility assessment using limited data. 85 The study area belongs to seismic zones where the maximum expected seismic intensities of 7, 8 and 9 are likely, (according to the MSK scale) for return periods of 100 and 500 years. The strongest recorded earthquake in the study area in the past 100 years was registered at its southwesternmost border on 30.09.1967 near the town of Debar (magnitude 6.5, intensity 9), Another strong earthquake was recorded in the northern part of the study area on 12.03.1960 in the town of Tetovo (magnitude 5.6, inten- sity 8). It is worth mentioning that many other earthquakes rang- ing in magnitude from 4 – 5 have occurred in the past. Despite its importance, the problem of earthquake related landslides in the region has not been well investigated. The climate in the study area is warm continental with many sunny days throughout the year. The summer is warm – hot and the winter is cold and snowy, while the spring and fall are char- acterized by rainfall. The medium air temperature in the region is around 11 °C, while the annual range is –20 to +38 °C. In the winter months, the snow cover reaches a thickness of 50 cm, and in the higher altitudes up to a metre. In relation to rainfall, the study area is characterized by the highest average annual precip- itation in the whole country. On some of the highest elevations in the area, the annual precipitation is over 1250 mm per year, while in the flat Polog valley the average is 800–900 mm. Almost one third of landslides in R. Macedonia have oc- curred in this region and many of them have caused damage to existing structures. In relation to lithological setting, 38% of the landslides in Macedonia have occurred in soil debris which cov- ers schistose or granitic bedrock, 11% in limestone (mostly rock- falls), and 31% in lacustrine sediments. Other units are less well represented. According to the landslide type, in the available re- ports, most of the landslides were defined as translational and rotational debris and earth slides and rock falls. Over 60% of the occurrences have been defined as intermediate and deep, and those mostly occurred in debris material overlying hard rock masses or in fluvioglacial and proluvial sediments. The depth of some of these landslides varies in the range from 5->25–30 me- tres. The remaining 40% of slides are shallow and mostly oc- curred in sandstones, lacustrine and flysch sediments. Rockfalls are usually classified as shallow slides. Almost 70% of the re- ported landslides in Macedonia have been caused by heavy rain- fall (PESHEVSKI, 2013). Due to rock falls and debris flows, there have even been several fatalities (PESHEVSKI et al., 2015b; HAQUE et al., 2016; PESHEVSKI et al., 2017). 3.2. Landslide inventory A landslide inventory map of the region was prepared using available geological maps (1:25.000 to 1:100.000), Google Earth, a Digital Elevation Model with resolution of 5 m, GIS–portal of the National Agency for Real Estate, ortho-photo images at a scale of 1:5.000, topographic (1:25.000) and geological maps (1:1.000 – 1:5.000) of past landslides. A total of 1172 example of slope instability were identified, most of which with undifferen- tiated types of slide mechanism. The size of landslides varies, but it can be stated that all of them had a clearly visible outline on some of the aforementioned data sources. For some of the landslides, published newspaper data was checked and data on their precise location was then confirmed by local residents. Un- fortunately, the exact date of occurrence is known for only a very small number of landslides. 3.3. Rating of units inside landslide conditioning parameter maps Rating of units inside landslide conditioning parameter maps were developed based on the available sources, which in this case included: the basic geological map of the R. Macedonia at a scale of 1:100000 (PETKOVSKI & POPOVSKI, 1982; MENKOVIĆ et al. 1982; PETKOVSKI & IVANOVSKI, 1979); digital eleva- tion model of R. Macedonia with an initial resolution of 5 m (AGENCY FOR REAL ESTATE CADASTRE, 2006); seismic intensity maps of the country according to the MSK scale for a return period of 100 and 500 years (SEISMOLOGICAL ASSO- CIATION OF SFR YUGOSLAVIA, 1978), the 1:1.000.000, map of average annual rainfall of the country (LAZAREVSKI, 1993), 1:1.000.000, land use map according to CORINE CLC2006 at 100 m resolution (EUROPEAN ENVIRONMENTAL AGENCY, 2009). The maps and their corresponding legends are presented in Fig. 3. Figure 2. a) Location of the Republic of Macedonia in Europe; b) Location of the case study region of Polog. G eo lo gi a C ro at ic a Geologia Croatica 72/186 3.4. Basic regional landslide susceptibility maps Two basic landslide susceptibility maps were produced, for the return period of maximum expected seismic intensity of 100 years (TLSR-100) and 500 years (TLSR-500). The results of the verifi- cation are presented in Tables. 1 & 2, as well as Fig. 4. Since more than 74% of the landslides belong in areas of high and very high susceptibility, it can be considered that both TLSR models provide a relatively good reflection of the landslide sus- ceptibility of the region. 3.5. Weighted landslide susceptibility maps Eight additional WTLSR maps are prepared with the application of sensitivity coefficients. The results are presented in Table 3 and Fig. 5. Figure 3. Rating of units within landslide conditioning parameters maps; a lithological type; b slope inclination; c seismic intensity, for return period 100 and; d for 500 years; e average annual rainfalls; f land use. Table 1. Landslide susceptibility classes for a return period of maximum expected seismic intensity of 100 years (TLSR-100). Susceptibility class TLSR-100 Value Area (km2) Area (%) Landslides Landslides (%) 1 – very low 0.709-3.596 420.35 17.44 16 1.37 2 – low 3.596-4.813 459.04 19.05 126 10.75 3 – medium 4.813-6.204 487.67 20.24 162 13.82 4 – high 6.204-7.630 448.26 18.60 267 22.78 5 – very high 7.630-9.577 594.29 24.66 601 51.28 Total: 2409.61 100 1172 100.00 G eologia C roatica Peshevski et al.: Preliminary regional landslide susceptibility assessment using limited data. 87 Table 2. Landslide susceptibility classes for a return period of maximum expected seismic intensity of 500 years (TLSR-500). Susceptibility class TLSR-500 Value Area (km2) Area (%) Landslides Landslides (%) 1 – very low 1.009-3.673 407.52 16.91 17 1.45 2 – low 3.673-4.905 478.76 19.87 130 11.09 3 – medium 4.905-6.304 485.68 20.16 148 12.63 4 – high 6.304-7.735 486.45 20.19 296 25.26 5 – very high 7.735-9.500 551.20 22.88 581 49.57 Total 2409.61 100 1172 100.00 Figure 4. Total Landslide Susceptibility Rating map for return period of maximum expected seismic intensity for a) 100 years; b) 500 years Table 3. Results from verification of the Weighted – Total Landslide Susceptibility Rating models (WTLSR). Model Susceptibility class Value of W-TLSR Area (km2) Area (%) No. of landslides Landslides in the zone (%) Difference from the initial model (%) WLT-TLSR 100 1– very low 0.124 – 0.698 429.74 17.83 25 2.13 0.76 2 – low 0.698 – 0.902 463.67 19.24 123 10.49 0.26 3 – medium 0.902 – 1.258 433.37 17.99 136 11.6 6.57 4 – high 1.258 – 1.643 414.89 17.22 239 20.39 2.39 5 – very high 1.643 – 2.051 667.94 27.72 649 55.38 4.1 Total 2409.61 100 1172 100.00 WLT-TLSR 500 1– very low 0.177 – 0.746 529.24 2.56 42 3.58 2.13 2 – low 0.746 – 0.957 471.05 32.45 132 11.26 0.17 3 – medium 0.957 – 1.286 333.62 20.4 112 9.56 3.07 4 – high 1.286 – 1.651 403.18 16.72 235 20.05 5.21 5 – very high 1.651 – 2.037 672.52 27.87 651 55.55 5.98 WSI-TLSR 100 1– very low 0.150 – 0.768 385.12 15.98 11 0.94 0.43 2 – low 0.768 – 1.110 404.76 16.8 65 5.55 5.2 3 – medium 1.110 – 1.408 618.88 25.68 230 19.62 5.8 4 – high 1.408 – 1.676 430.94 17.88 280 23.89 1.11 5 – very high 1.676 – 2.048 569.91 23.65 586 50.0 0.43 WSI-TLSR 500 1– very low 0.192 – 0.786 383.66 15.92 12 1.02 0.43 2 – low 0.786 – 1.126 406 16.85 70 5.97 5.12 3 – medium 1.126 – 1.415 614.57 25.5 219 18.69 6.06 4 – high 1.415 – 1.683 423.28 17.57 275 23.46 1.8 5 – very high 1.683 – 2.037 582.1 24.16 596 50.85 1.28 G eo lo gi a C ro at ic a Geologia Croatica 72/188 3.6. Analysis of the results From Tables 1, 2 & 3, it can be seen that the landslide suscepti- bility model with the highest percentage of landslides in very high and high susceptibility classes is the WLT-TLSR-100. According to this model, 55.06 % of the region’s territory belongs to zones with very low, low or medium landslide susceptibility, 24.22 % of landslides have been mapped in this area. Alternatively, ter- rains with high and very high susceptibility cover 44.94 % of the territory, with 75.77 % of the total landslides mapped. A landslide susceptibility model with the lowest percent of landslides in classes with very high and high susceptibility is the WAP-TLSR-100. According to this model 58.75% of the region’s Model Susceptibility class Value of W-TLSR Area (km2) Area (%) No. of landslides Landslides in the zone (%) Difference from the initial model (%) WAP-TLSR 100 1– very low 0.137 – 0.719 412.05 17.1 10 0.85 0.52 2 – low 0.719 – 1.014 472.92 19.63 141 12.03 1.28 3 – medium 1.014 – 1.287 530.51 22.02 183 15.61 1.79 4 – high 1.287 – 1.547 396.62 16.46 235 20.05 2.73 5 – very high 1.547 – 1.926 597.51 24.8 603 51.45 0.17 WAP–TLSR 500 1– very low 0.189 – 0.736 405.12 16.81 11 0.94 0.51 2 – low 0.736 – 1.027 481.14 19.97 140 11.95 0.86 3 – medium 1.027 – 1.297 529.64 21.98 175 14.93 2.3 4 – high 1.297 – 1.561 407.31 16.9 257 21.93 3.33 5 – very high 1.561 – 1.912 586.4 24.34 589 50.26 0.69 WLU-TLSR 100 1– very low 0.137 – 0.672 425.08 17.64 19 1.62 0.25 2 – low 0.672 – 0.868 451.39 18.73 123 10.49 0.26 3 – medium 0.868 – 1.109 488.93 20.29 160 13.65 0.17 4 – high 1.109 – 1.364 503.36 20.89 304 25.94 3.16 5 – very high 1.364 – 1.801 540.85 22.45 566 48.29 2.99 WLU-TLSR 500 1– very low 0.189 – 0.691 424.97 17.64 19 1.62 0.17 2 – low 0.691 – 0.885 460.33 19.1 126 10.75 0.34 3 – medium 0.885 – 1.123 482.96 20.04 151 12.88 0.25 4 – high 1.123 – 1.367 477.96 19.84 286 24.4 0.86 5 – very high 1.367 – 1.787 563.39 23.38 590 50.34 0.77 Figure 5. Weighted - Total Landslide Susceptibility Rating maps: a WLT-TLSR100; b WLT-TLSR500; c WSI-TLSR-100; d WSI-TLSR-500; e WAP-TLSR-100; f WAP-TLSR-500; g WLU-TLSR-100; h WLU-TLSR-500. G eologia C roatica Peshevski et al.: Preliminary regional landslide susceptibility assessment using limited data. 89 territory belongs to zones with very low, low and medium land- slide susceptibility, with the presence of 28.49% of the total num- ber of landslides. Terrains with high and very high susceptibility cover 41.26 % of the territory, with 71.50 % of the total landslides. In a comparison with other models, the WLT-TLSR model shows the highest prediction ability for the analyzed study area, as long as the size of the susceptibility zones is not considered. However, this does not necessarily mean that one of these models should be adopted as the most suitable for the study area, which is dis- cussed later. The sensitivity does not show a straightforward case if we consider the landslide percentage per susceptibility class as a ref- erence, and compare the TLSR with the WTLRS models. The model with the highest performance is not necessarily the most sensitive one, and in addition, it can be more sensitive to the change of the parameter that is not weighted with the highest value (0.3). The sensitivity varies across different landslide sus- ceptibility classes. In most cases, very low and low susceptibility classes show the least sensitivity to any conditioning factor change, except for slope inclination, because low susceptibility classes show high change (around 5.2% in both variants, WSI- TLSR-100 and WSI-TLSR-500). Thus, low and very low suscep- tibility is sensitive to slope inclination parameter change and in- sensitive to any other alteration. Medium susceptibility class is predominantly sensitive to lithological type and slope inclination change (around 5.5% in both variants, WLT/WSI-TLSR-100 and WLT/WSI-TLSR-500). High susceptibility is particularly sensi- tive to the change of lithological type and annual precipitation (5.2% for WLT-TLSR-500 and 3.3% WAP-TLSR-500 models), as well as to a change of land use weighting, but to a lesser extent (3.16% for the WLU-TLSR-100 model). Very high susceptibility is the most sensitive to change in lithological type, with rates of 4.1–5.98%. The total sensitivity (the sum of the difference across all susceptibility classes in all models) indicates that the most dramatic changes in landslide distribution are evident for litho- logical type (14.08% and 16.6%, for WLT-TLSR-100 and WLT- TLSR-500 variants, respectively), followed by slope inclination (12.97% and 14.69%, for WSI-TLSR-100 and WSI-TLSR-500 var- iants, respectively). The smallest changes in landslide distribution are noted for land use (6.83% for WLU-TLSR-100 and 2.39% WLU-TLSR-500). This supports earlier claims, i.e. the sensitiv- ity is not directly related to the performance. If we take into account the sensitivity of the different WTLSR models in relation to the change in size of the zones with high and very high susceptibility (in relation to the initial basic models TLSR100/500) and compare these with the change in number of landslides, then it can be seen that different models show changes from –1.83% to 1.68% for the size and –2.64% to 1.71% for the number of landslides (Tab. 4). Also, Table 4 shows that models WSI-TLSR100 and WSI- TLSR500 have the most significant change in size of the area un- der high and very high susceptibility, and the least change in number of landslides in these zones, both with negative trends. On the other hand, the models WAP-TLSR100 and WAP- TLSR500 show the most significant decrease in size of the area, but in this case also with a more significant decrease in the num- ber of landslides. The high and very high classes in models LU- TLSR are shown to be the least sensitive to the application of the sensitivity coefficients, while the WLT-TLSR models show cor- relation with a positive trend of increase both in the size and num- ber of landslides. This leads to the conclusion that the model WSI-TLSR100 is the most suitable to be applied for this particular case study re- gion, since the size of the area under high and very high suscep- tibility has reduced by 41.7 km2 from the initial model, and only 2 landslides had been moved to a lower susceptibility class zone. 4. DISSCUSSION AND CONCLUSIONS With the results obtained, it was confirmed that the selection of landslide conditioning parameters and the way in which they are being taken in the modelling have by far the greatest impact in any particular landslide susceptibility method. In practice vari- ous conditioning parameters are in use when performing land- slide susceptibility or hazard zoning: lithology, engineering-geo- logical characteristics of the rock masses, DEM derivatives, slope aspect, elevation, normalized difference vegetation index (NDVI) , and land use, landform classification, rainfall, runoff curve num- ber (CN), seismic influence, distance from structural elements (such as faults), distance to drainage systems, distance to linear infrastructure (roads & railways), anthropogenic factors, etc. The number of conditioning parameters in the models ranges from 3 (e.g. RAPOLLA et al., 2012) to >20 (e.g. LAGOMARSINO et al., 2014). All of these various susceptibility models have their advantages and disadvantages in relation to the reliability of data sources, questionable correlation with landslide development (for example the use of distance to road or distance to drainage fac- tors), subjectivity of the engineering-geological classification methods (rock mass classification), precision of measuring instru- ments, etc. As shown by MIHALIĆ ARBANAS et al., (2016), even very high resolution bare-earth DEM derivatives can have limitations in some aspects. With the methodology presented herein, it is shown that even in the case of very limited and in some cases poor amounts of data from various aspects (incom- plete landslide inventory, different scale landslide conditioning maps, missing very precise data on rainfall quantity, different periods of map production, etc.) it is still possible to produce rep- resentative landslide susceptibility zonation. Our modelling results show more than 74% of landslides fall- ing in zones of high and very high susceptibility classes, which is comparable with some more complex statistical studies. HE & BEIGHLEY, (2008) obtained 71% and KAVZOGLU et al., (2014) 69–77% precision by applying statistically based algo- rithms. Table 4. Change of the total size of zones with high and very high landslide susceptibility and corresponding change of number of landslides in these zones. TLSR-100 TLSR-500 WLT-TLSR 100 WLT-TLSR 500 WSI-TLSR 100 WSI-TLSR 500 WAP-TLSR 100 WAP-TLSR 500 WLU-TLSR 100 WLU-TLSR 500 Area (km2) 1042.55 1037.65 1082.83 1075.7 1000.85 1005.38 994.13 993.71 1044.21 1041.35 Change (km2) +40.28 +38.05 –41.7 –32.27 –48.42 –43.94 +1.66 +3.7 Percent of change +1,68 +1,57 –1,73 –1,34 –2 –1,83 +0,08 +0,15 No. of landslides 868 877 888 886 866 871 838 846 870 876 Change in No. of landslides +20 +9 –2 –6 –30 –31 +2 –1 Percent of change +1,71 +0,77 –0,17 –0,52 –2,56 –2,64 +0,17 –0,09 G eo lo gi a C ro at ic a Geologia Croatica 72/190 Based on results presented in our research, it is obvious that there is also the possibility to test the sensitivity of landslide sus- ceptibility models in a simple, yet acceptable manner. The opti- mization of preliminary prepared models should always be per- formed, resulting in a more robust model, with an adjusted, more realistic distribution of landslide susceptibility classes. It is hereby recommended that sensitivity testing be applied as a standard post-processing procedure in any regional heuristic landslide susceptibility analysis, and especially for cases with limited amounts of data for the landslide conditioning factors’ variability and properties of the landslides. In order to improve the sensitivity testing method presented herein, and test its reli- ability for use in more advanced landslide susceptibility ap- proaches, further studies are suggested: differentiation in respect to the landslide mechanism of landslides, use of other thematic maps for the conditional parameters (for example use of peak ground acceleration maps instead of seismic intensity maps, use of maximum rainfall intensity maps for 24, 48 and 72 hours in- stead of annual precipitation maps) using a larger scale of some of the parametric maps, etc. It should be noted once again that, even methods such as those presented in this paper can be used only as a first estima- tion of landslide susceptibility, and one can argue that they are based mostly on engineering judgment. They can however pro- vide indications of the aspects that deserve attention while allow- ing for the undertaking of more advanced studies in the future. For example, to point out zones of some regions where more de- tailed datasets need to be collected in order to apply more ad- vanced susceptibility assessment methods. These datasets will consist of a more extensive and detailed landslide inventory of the analyzed region, supported by field mapping/remapping of more significant landslides or zones with dense landslides, same scale thematic conditioning parameter maps, and updated data- sets on rainfall. Therefore, the presented approach is considered to be especially practical for large scale (regional) mapping in developing countries, where various obstacles from an organiza- tional, technical and economic nature exist. It is simple and makes reasonable the use of the limited but readily available data- sets in cases where greater detail is lacking, while still producing representative susceptibility zonation of a region. The possibility to vary ratings and weightings of sensitivity coefficients of the carefully selected conditioning parameters is considered to be another advantage of the approach. Due to the large scale, (the same as for other available regional susceptibility zonation meth- ods), the spatial precision should be confirmed by repeating the analysis over a period of 5 – 10 years. Such analysis should be performed using an updated landslide inventory map and, if available, refreshed datasets of the temporally dependent condi- tioning parameters. As for any other landslide susceptibility modelling approach, the one presented in this paper also has some limitations and re- strictions. These include the inability to use it for a coarser scale analysis where more exact correlation between the landslide con- ditioning factors and landslide occurrences needs to be supported by detailed statistical data, rainfall thresholds, and exact earth- quake intensities during sliding events. Another drawback is the different scale and age of production of the thematic maps used which range from the 5 m resolution DEM to the 1:1.000.000 seis- mic intensity maps and rainfall maps. This required certain ex- trapolations to be undertaken, which certainly influenced the out- come of the study. The extent of this inf luence for each extrapolated parameter should be further tested in follow-up studies. The landslide mapping technique which was applied is similar to that used in most studies at this scale, and is considered to be an error prone task as suggested by GALLI et al. (2008). Using this technique it is only possible to produce a landslide in- ventory with undifferentiated mechanisms of landslide move- ment. Therefore, to improve the quality of the inventory map, and consequently the quality of the validation processes, further land- slide susceptibility studies in the region should be supported by field checks. REFERENCES ABOLMASOV, B. & STOJKOV, K. (1994): The influence on the landslide in urban planning in Belgrade City.– In: Proceedings of VII IAEG Congress, Lisbon, Vol. III, p. 2161–2168. ABOLMASOV, B. & STOJKOV, K. (1995): Jedan primer izrade karte mogucnosti pojave i razvoja klizanja.– Drugi Simpozijum Istraživanje i sanacija klizišta, Don- ji Milanovac. AKGUN, A. (2012): A comparison of landslide susceptibility maps produced by logis- tic regression, multi-criteria decision, and likelihood ratio methods: a case study at İzmir, Turkey.– Landslides, 9/1, 93–106. doi: 10.1007/s10346-011-0283-7 ARDIZZONE, F., CARDINALI, M., CARRARA, A., GUZZETTI, F. & REICHEN- BACH, P. (2002): Impact of mapping errors on the reliability of landslide hazard maps.– Nat. Hazards Earth Syst. Sci., 2, 3–14. doi: 10.5194/nhess-2-3-2002 BARREDO, J.I., BENAVIDES, A., HERVÁS, J. & VAN WESTEN, C.J. (2000): Com- paring heuristic landslide hazard assessment techniques using GIS in the Tirajana basin, Gran Canaria Island, Spain.– International Journal of Applied Earth Observa- tion and Geoinformation, 2/1, 9–23. doi: 10.1016/S0303-2434(00)85022-9 BĂLTEANU, D., CHENDEŞ, V., SIMA, M. & ENCIU, P. (2010): A country-wide spatial assessment of landslide susceptibility in Romania.– Geomorphology, 124, 102–112. doi: 10.1016/j.geomorph.2010.03.005 CATANI, F., LAGOMARSINO, D., SEGONI, S., & TOFANI, V. (2013): Landslide susceptibility estimation by random forests technique: sensitivity and scaling is- sues.– Nat. Hazards Earth Syst. Sci., 13, 2815–2831. doi: 10.5194/nhess-13-2815- 2013 CHALKIAS, C., FERENTINOU, M. & POLYKRETIS, C. (2014): GIS-Based Land- slide Susceptibility Mapping on the Peloponnese Peninsula, Greece.– Geoscienc- es, 4, 176–190. doi:10.3390/geosciences4030176 CHANG, H. & KIM, N.K. (2004): The evaluation and the sensitivity analysis of GIS- based landslide susceptibility models.– Geosciences Journal, 8/4, 415–423. doi: 10.1007/BF02910477 CONSTANTIN, M., BEDNARIK, M., JURCHESCU, M.C. & VLAICU, M. (2011): Landslide susceptibility assessment using the bivariate statistical analysis and the index of entropy in the Sibiciu Basin (Romania).– Environmental Earth Sciences, 63/2, 397–406. doi: 10.1007/s12665-010-0724-y COROMINAS, J., VAN WESTEN, C., FRATTINI, P., CASCINI, L., MALET, J-P., FOTOPOULOU, S., CATANI, F., VAN DEN EECKHAUT, M., MAVROULI, O., AGLIARDI, F., PITILAKIS, K., WINTER, M.G., PASTOR, M., FERLISI, S., TOFANI, V., HERVAS, J. & SMITH, J. T. (2014): Recommendations for the quantitative analysis of landslide risk.– Bulletin of Engineering Geology and the Environment, 73, 209–263. doi: 10.1007/s10064-013-0538-8 DEMIR, G., AYTEKIN, M. & AKGUN, A. (2015): Landslide susceptibility mapping by frequency ratio and logistic regression methods: an example from Niksar–Re- sadiye (Tokat, Turkey).– Arabian Journal of Geosciences, 8/3, 1801–1812. doi: 10.1007/s12517-014-1332-z DEVKOTA, K.C., REGMI, A.D., POURGHASEMI, H.R., YOSHIDA, K., PRAD- HAN, B., RYU, I.C., DHITAL, M.R. & ALTHUWAYNEE, O.F. (2013): Land- slide susceptibility mapping using certainty factor, index of entropy and logistic regression models in GIS and their comparison at Mugling–Narayanghat road sec- tion in Nepal Himalaya.– Natural Hazards, 65/1, 135–165. doi: 10.1007/s11069- 012-0347-6 DRAGICEVIĆ, S., CAREVIĆ, I., KOSTADINOV, S., NOVKOVIĆ, I., ABOLMA- SOV, B., MILOJKOVIĆ, B. & SIMIĆ, D. (2012): Landslide susceptibility zo- nation in the Kolubara River Basin (western Serbia) – analysis of input data.– Car- pathian Journal of Earth and Environmental Sciences, 7/2, 37–47. EUROPEAN ENVIRONMENTAL AGENCY (2009): 2006 Corine land cover, 100 m – version 12/ 2009 ERENER, A. & DÜZGÜN, H.S.B. (2010): Improvement of statistical landslide sus- ceptibility mapping by using spatial and global regression methods in the case of More and Romsdal (Norway).– Landslides, 7, 55–68. doi: 10.1007/s10346-009- 0188-x ERCANOGLU, M., KASMER, O. & TEMIZ, N. (2008): Adaptation and comparison of expert opinion to analytical hierarchy process for landslide susceptibility map- ping.– Bulletin of Engineering Geology and the Environment, 67/4, 565–578. doi: 10.1007/s10064-008-0170-1 G eologia C roatica Peshevski et al.: Preliminary regional landslide susceptibility assessment using limited data. 91 FELL, R., COROMINAS, J., BONNARD, C., CASCINI, L., LEROI, L. & SAVAGE, W.Z. (2008): Guidelines for landslide susceptibility, hazard and risk zoning for land use planning.– Engineering Geology, 102/3–4, 85–98. doi: 10.1016/j.enggeo. 2008.03.022 FEIZIZADEH, B. & BLASCHKE, T. (2014): An uncertainty and sensitivity analysis approach for GIS-based multicriteria landslide susceptibility mapping.– Interna- tional Journal of Geographical Information Science, 28/3, 610–638. doi: 10.1080/ 13658816.2013.869821 FRATTINI, P., CROSTA, G. & CARRARA, A. (2010): Techniques for evaluating the performance of landslide susceptibility models. – Engineering Geology, 111/1–4, 62–72. doi: 10.1016/j.enggeo.2009.12.004 GALLI, M., ARDIZZONE, F., CARDINALI, M., GUZZETTI, F., REICHENBACH, P. (2008): Comparing landslide inventory maps.– Geomorphology, 94/3–4, 268– 289. doi: org/10.1016/j.geomorph.2006.09.023. GLADE, T., PETSCHKO, H., BELL, R., BAUER, C., GRANICA, K., HEISS, G., LEOPOLD, P., POMAROLI, G., PROSKE, H. & SCHWEIGL, J. (2012): Land- slide susceptibility maps for lower Austria, methods and challenges.– In: Proceed- ings of 12th Congress INTERPRAEVENT 2012, Grenoble, France, 497–508. GÜNTHER, A., EECKHAUT, M., REICHENBACH, P., HERVÁS, J., MALET, J.P., FOSTER, C. & GUZZETTI, F. (2013a): New developments in harmo- nized landslide susceptibility mapping over Europe in the framework of the European soil thematic strategy.− In: MARGOTTINI, C., CANUTI, P. & SASSA, K. (eds.): Landslides science and practice. Vol. 1., Landslide Inven- tory and Susceptibility and Hazard Zoning, Springer, Berlin, Heidelberg, 1, 297–301. doi: 10.1007/978-3-642-31325-7_39 GÜNTHER, A., REICHENBACH, R., MALET, J.P., EECKHAUT, M.V.D., HER- VÁS, J., DASHWOOD, C. & GUZZETTI, F. (2013b): Tier-based approaches for landslide susceptibility assessment in Europe.– Landslides, 10/5, 529–546. doi: 10.1007/s10346-012-0349-1 GUZZETTI, F., REICHENBACH, P., ARDIZZONE, F., CARDINALI, M. & GAL- LI, M. (2006): Estimating the quality of landslide susceptibility models.– Geomor- phology, 81/ 1–2, 166–184. doi: 10.1016/j.geomorph.2006.04.007 HAQUE, U., BLUM, P.,DA SILVA, P. F., ANDERSEN, P., PILZ, J., R. CHALOV, S., MALET, J.P., JEMEC AUFLIČ, M., ANDRES, N., POYIADJI, E., LA- MAS, P.C., ZHANG, W., PESHEVSKI, I., G. PÉTURSSON, H., KURT, T., DOBREV, N., GARCÍA-DAVALILLO, J.C., HALKIA, M., FERRI, S., GA- PRINDASHVILI, G., ENGSTRÖM, J. & KEELLINGS, D. (2016): Fatal land- slides in Europe.– Landslides,13/6, 1545–1554. doi: 10.1007/s10346-016-0689-3 HE, Y.P. & BEIGHLEY, R.E. (2008): GIS-Based Regional Landslide Susceptibility Mapping: A Case Study in Southern California.– Earth Surface Processes and Landforms, 33, 380–393. doi: 10.1002/esp.1562 HERVÁS, J. (2007): Guidelines for mapping areas at risk of landslides in Europe.– In: Proceedings of the experts meeting, JRC, Ispra, Italy, 23–24 Oct 2007. JRC Report EUR 23093 EN, Office for Official Publications of the European Communities, Luxembourg. p. 53. doi: 10.2788/63147 ILANLOO, M. (2011): A comparative study of fuzzy logic approach for landslide sus- ceptibility mapping using GIS: An experience of Karaj dam basin in Iran.– Proce- dia, Social and Behavioral Sciences, 19, 668–676. doi: 10.1016/j.sbspro.2011.05. 184 JOVANOVSKI, M., ABOLMASOV, B. & PESHEVSKI, I. (2013a): Analyses of Landslide Hazard Evaluation Factors Using Polynomial Interpolation.− In: MAR- GOTTINI, C., CANUTI, P. & SASSA, K. (eds.): Landslide Science and Practice. Vol.1., Landslide Inventory and Susceptibility and Hazard Zoning, Springer, Ber- lin, Heidelberg, 1, 561–566. doi: 10.1007/978-3-642-31325-7_73 JOVANOVSKI, M., MILEVSKI, I., PAPIC, J., PESEVSKI, I. & MARKOSKI, B. (2013b): Landslides in the Republic of Macedonia Triggered by Extreme Events in 2010.– In: LOCZY, D. (eds.): Geomorphological impacts of extreme weather. Springer Geography, Springer, Dordrecht, 265–279. doi: 10.1007/978-94-007- 6301-2_17 KAVZOGLU, T., SAHIN, E.K. & COLKESEN, I. (2014): Landslide susceptibility mapping using GIS-based multi-criteria decision analysis, support vector ma- chines, and logistic regression.– Landslides, 11/3, 425–439. doi: 10.1007/s10346- 013-0391-7 KOMAC, M. (2005): Probabilistic model of slope mass movement susceptibility – a case study of Bovec municipality, Slovenia.– Geologija, 48/2, 311–340. doi: 10. 5474/geologija.2005.026 KOMAC, M. & RIBICIC, M. (2006): Landslide susceptibility map of Slovenia at scale 1: 250,000.– Geologija, 49/2, 295–309. doi: 10.5474/geologija.2006.022 LAGOMARSINO, D., SEGONI, S., FANTI, R., CATANI, F. & CASAGLI, N. (2014): Regional Scale Landslide Susceptibility Mapping in Emilia Romagna (It- aly) as a Tool for Early Warning.– In: SASSA K., CANUTI, P. & YIN, Y. (eds.): Landslide Science for a Safer Geoenvironment. Vol. 2., Methods of Landslide Studies, 2, 443–449. doi: 10.1007/978-3-319-05050-8_69 LAZAREVSKI, A. (1993): Climate in Macedonia. Skopje (In Macedonian) LEE, S. & TALIB, J.A. (2005): Probabilistic landslide susceptibility and factor effect analysis.– Environmental Geology, 47/7, 982–990. doi: 10.1007/s00254-005- 1228-z LEE, M.J., CHOI, J.W., OH, H.J., WON, J.S., PARK, I. & LEE, S. (2012): Ensemble- based landslide susceptibility maps in Jinbu area, Korea.– Environmental Earth Sciences, 67/ 1, 23–37. doi: 10.1007/s12665-011-1477-y LIN, L., LIN, Q., & WANG, Y. (2017): Landslide susceptibility mapping on a global scale using the method of logistic regression.– Nat. Hazards Earth Syst. Sci., 17, 1411–1424. doi: 10.5194/nhess-17-1411-2017 LEOPOLD, P., HEISS, G., PETSCHKO, H., BELL, R. & GLADE, T. (2013): Sus- ceptibility Maps for Landslides Using Different Modelling Approaches.− In: MARGOTTINI, C., CANUTI, P. & SASSA, K. (eds.): Landslide Science and Practice. Vol. 1., Landslide Inventory and Susceptibility and Hazard Zoning, Springer, Berlin, Heidelberg, 1, 353–356. doi: 10.1007/978-3-642-31325-7_46 MAGLIULO, P., DI LISIO, A. & RUSSO, F. (2009): Comparison of GIS-based meth- odologies for the landslide susceptibility assessment.– GeoInformatica, 13/3, 253– 265. doi: 10.1007/s10707-008-0063-2 MARJANOVIĆ, M. (2010): Regional scale landslide susceptibility analysis using dif- ferent GIS-based approaches.– In: Geologically Active, WILLIAMS, A.L., PINCHES, G.M., CHIN, C.Y., MCMORRAN, T.J. & MASSEY, C.I. (eds.): Proceedings of 11th IAEG Congress, Auckland, New Zealand, 435–442. MARJANOVIĆ, M. (2014): Conventional and machine learning methods for landslide assessment in GIS.– 1st edition, Palacky University in Omolouc, Czech Republic, (monograph), 205 p. MENKOVIĆ, LJ., KAROVIĆ, J. & KOŠĆAL, M. (1982): Basic geological map of Socialist Federal Republic of Yugoslavia 1:100000, sheet Prizren. Prepared in pe- riod 1973–1979 by Geological Survey – Skopje. Printed by the Military Geograph- ical Institute of Yugoslavia. MIHALIĆ, S. (1998): Recommendations for Landslide Hazard and Risk Mapping in Croatia.– Geologia Croatica, 51/2, 195–204. MIHALIĆ ARBANAS, S., KRKAČ, M. & BERNAT, S. (2016): Application of in- novative technologies in landslide research in the area of the City of Zagreb (Cro- atia, Europe).– Geologia Croatica, 69/2, 231–243. doi:10.4154/gc.2016.18 NARANJO, J.L., WSTEN, C.J. AN & SETERS, R. (1994): Evaluating the use of train- ing areas in bivariate statistical landslide hazard analysis – a case study in Colom- bia.– ITC Journal (3): 292–300. OH, H.J., LEE, S., CHOTIKASATHIEN, W., KIM, C.H. & KWON, J.H. (2009): Predictive landslide susceptibility mapping using spatial information in the Pech- abun area of Thailand.– Environmental Geology, 57/3, 641–651. doi: 10.1007/ s00254-008-1342-9 PAMELA, P., SADISUN, I. & ARIFIANTI, Y. (2018): Weights of Evidence Method for Landslide Susceptibility Mapping in Takengon, Central Aceh, Indonesia.– IOP Confe-rence Series: Earth and Environmental Science, p. 118. doi: 10.1088/1755- 1315/118/1/012037 PESHEVSKI, I., JOVANOVSKI, M., MARKOSKI, B., PETRUSEVA, S. & SUSI- NOV, B. (2013): Landslide inventory map of the Republic of Macedonia, statistics and description of main historical landslide events.− In: MIHALIĆ ARBANAS, S. & ARBANAS, Ž. (eds.): Proceedings of the 1st Regional Symposium on Land- slides in the Adriatic-Balkan Region, Landslide and flood hazard assessment, Za- greb, Croatia, 207–212. PESHEVSKI, I. (2015a). Landslide susceptibility modelling using GIS technology, Doctoral thesis.− Faculty of Civil Engineering, Skopje, 181 p. PESHEVSKI, I., JOVANOVSKI, M., PAPIC, BR. J. & JOVICIC, L. (2015b): Prep- aration of landslide susceptibility maps using GIS, recent case study in R. Macedo- nia.– Scientific Journal of Civil Engineering, 4/1, 47–53. PESHEVSKI, I., PETERNEL, T. & JOVANOVSKI, M. (2017): Urgent need for ap- plication of integrated landslide risk management stretegies for the Polog region in R. of Macedonia.– In: MIKOŠ, M., VILÍMEK, V., YIN, Y. & SASSA, K. (eds.): Advancing Culture of living with landslides. Vol. 5., Landslides in different environments, Springer International Publishing, 5, 135–145. doi: 10.1007/978-3- 319-53483-1_43 PETKOVSKI, P. & IVANOVSKI, T. (1979): Basic Geological Map of Socialist Fed- eral Republic of Yugoslavia 1:100000, sheet Kicevo. Prepared in period 1963-1972 by Geological Survey – Skopje. Printed by the Military geographical institute of Yugoslavia. PETKOVSKI, P. & POPOVSKI, S. (1982): Basic geological map of Socialist Federal Republic of Yugoslavia 1: 100000, sheet Gostivar. Prepared in period 1967–1976 by Geological Survey – Skopje. Printed by the Military geographical institute of Yugoslavia. PETSCHKO, H., BRENNING, A., BELL, R., GOETZ, J. & GLADE, T. (2014): As- sessing the quality of landslide susceptibility maps – case study Lower Austria.– Natural Hazards and Earth System Sciences, 14, 95–118. doi: 10.5194/nhess-14- 95-2014 POSCH-TRÖZMÜLLER, G. & BÄK, R. (2013): Minimal Standards for Susceptibil- ity and Hazard Maps of Landslides and Rock Falls.− In: MARGOTTINI, C., CA- NUTI, P. & SASSA, K. (eds.): Landslide Science and Practice, Vol.1. Landslide Inventory and Susceptibility and Hazard Zoning, Springer, Berlin, Heidelberg, 1, 387–393. doi: 10.1007/978-3-642-31325-7_51 G eo lo gi a C ro at ic a Geologia Croatica 72/192 RAMESH, V. & ANBAZHAGAN, S. (2015): Landslide susceptibility mapping along Kolli hills Ghat road section (India) using frequency ratio, relative effect and fuzzy logic models.– Environmental Earth Sciences, 73/12, 8009–8021. doi: 10.1007/ s12665-014-3954-6 RAPOLLA, A., DI NOCERA, S., MATANO, F., PAOLETTI, V. & TARALLO, D. (2012): Susceptibility regional zonation of earthquake-induced landslides in Cam- pania, Southern Italy.– Natural Hazards, 61/1, 115–126. doi: 10.1007/s11069-011- 9790-z ROSSI, M., GUZZETTI, F., REICHENBACH, P., MONDINI, A.C. & PERUCCAC- CI, S. (2010): Optimal landslide susceptibility zonation based on multiple fore- casts.– Geomorphology, 114/3, 129–142. doi: 10.1016/j.geomorph.2009.06.020 ROY, C.J. & OBERKAMPF, W.L. (2011): A comprehensive framework for verifica- tioni, validation, and uncertainty quantification in scientific computing.– Compu- ter Methods in Applied Mechanics and Engineering, 200, 25–28, 2131–2144. doi: 10.1016/j.cma.2011.03.016 ROZOS, D., BATHRELLOS, G.D. & SKILLODIMOU, H.D. (2010): Comparison of the implementation of rock engineering system and analytic hierarchy process methods, upon landslide susceptibility mapping, using GIS: a case study from the Eastern Achaia County of Peloponnesus, Greece.– Environmental Earth Sciences, 63/1, 49–63. doi: 10.1007/s12665-010-0687-z SEISMOLOGICAL ASSOCIATION OF SFR YUGOSLAVIA (1987): Seismic Zon- ing Maps of SFRY 1:1,000,000) Degrees by MSK-64 Seismic Intensity Scale). For Return Periods of 50, 100, 200, 500, 1000 and 10000 years. SHAHABI, H., HASHIM, M. & AHMAD, B.B. (2015): Remote sensing and GIS- based landslide susceptibility mapping using frequency ratio, logistic regression, and fuzzy logic methods at the central Zab basin, Iran.– Environmental Earth Sci- ences, 73/12, 8647–8668. doi: 10.1007/s12665-015-4028-0 SHAHABI, H. & HASHIM, M. (2015): Landslide susceptibility mapping using GIS- based statistical models and Remote sensing data in tropical environment.– Scien- tific Reports, 5/3:15, 9899. doi: 10.1038/srep09899 SHIBIAO, B., QIANG, X., WANG, J. & PINGGEN, Z. (2013): Pre-conditioning Fac- tors and Susceptibility Assessments of Wenchuan Earthquake Landslide at the Zhouqu Segment of Bailongjiang Basin, China.– Journal of the Geological Soci- ety of India, 82, 575–582. doi: 10.1007/s12594-013-0190-8 SOETERS, R. & VAN WESTEN, C.J. (1996): Slope instability recognition, analysis, and zonation.– In: TURNER, A.K. & SCHUSTER, R.L. (eds.): Landslides, in- vestigation and mitigation. National Academy Press, Washington, 129–177. THIERY, Y., MALET, J.P., STERLACCHINI, S., PUISSANT, A. & MAQUAIRE, O. (2007): Landslide susceptibility assessment by bivariate methods at large scales: application to a complex mountainous environment. – Geomorphology, 92/1,2, 38–59. doi: /10.1016/j.geomorph.2007.02.020 TRIGILA, A., FRATTINI, P., CASAGLI, N., CATANI, F., CROSTA, G., ESPOSI- TO, C., IADANZA, C., LAGOMARSINO, D., MUGNOZZA, G., SEGONI, S., SPIZZICHINO, D., TOFANI, V. & LARI, S. (2013): Landslide susceptibi- lity mapping at national scale: The Italian case study.− In: MARGOTTINI, C., CANUTI, P. & SASSA, K. (eds.): Landslide Science and Practice, Vol. 1., Land- slide Inventory and Susceptibility and Hazard Zoning, Springer, Berlin, Heidel- berg, 1, 287–295. doi: 10.1007/978-3-642-31325-7_38 UMAR, Z., PRADHAN, B., AHMAD, A., JEBUR, M.N. & TEHRANY, M.S. (2014): Earthquake induced landslide susceptibility mapping using an integrated ensemble frequency ratio and logistic regression models in West Sumatera Province, Indo- nesia.– Catena, 118, 124-135. doi.org/10.1016/j.catena.2014.02.005 WANG, Q., LI, W., CHEN, W. & BAI, H. (2015): GIS-based assessment of landslide susceptibility using certainty factor and index of entropy models for the Qianyang County of Baoji city, China.– Journal of Earth System Science, 124/7, 1399–1415. doi: 10.1007/s12040-015-0624-3. WANG, L.J., GUO, M., SAWADA, K., LIN, J. & ZHANG, J. (2016): A comparative study of landslide susceptibility maps using logistic regression, frequency ratio, decision tree, weights of evidence and artificial neural network.– Geosciences Jour- nal, 20/1, 117–136. doi: 10.1007/s12303-015-0026-1 WANG, Q., WANG, D., HUANG, Y., WANG, Z., ZHANG, L., GUO, Q., CHEN, W., CHEN, W. & SANG, M. (2015): Landslide Susceptibility Mapping Based on Se- lected Optimal Combination of Landslide Predisposing Factors in a Large Catch- ment.– Sustainability, 7/12, 16653–16669. doi: 10.3390/su71215839 VAN WESTEN, C.J., & VAN ASCH, T.W.J. & SOETERS, R. (2006): Landslide ha- zard and risk zonation-why is it so difficult?– Bulletin of Engineering Geology and Environment, 65/2, 167–184. doi: 10.1007/s10064-005-0023-0 ZHU, A-X., WANG, R., QIAO J., QIN, C-Z., CHEN, Y., LIU, J., DU, F., LIN, Y. & ZHU, T. (2014): An expert knowledge-based approach to landslide susceptibility mapping using GIS and fuzzy logic.– Geomorphology, 214, 128–138. doi: 10. 1016/j.geomorph.2014.02.003