Preliminary regional landslide susceptibility assessment using limited data

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 into account 5 landslide conditioning parameters: lithology, slope inclination, average annual rainfall, 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 landslide 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 dividing the region into five susceptibility zones according to Jenks natural breaks classification. Verification 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 additionally tested by multiplying the conditioning parameter’s rating points by sensitivity coefficients. In this way, additional landslide susceptibility models are obtained, named Weighted Total Landslide 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 methods. Additional WTLSR models were prepared, and their correlation identified the best model. The presented approach is considered to be very convenient for conducting preliminary regional landslide susceptibility assessments with the ability to fine-tune the results. Due to its simplicity, it can be applied to additional landslide conditioning parameters other than the one presented 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 restrictions 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; COROMINAS et al., 2014) have presented ideas of the minimal requirements 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 modelling is the qualitative experience (knowledge) based – heuristic approach, which is considered a rather controversial topic in landslide 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). Alternatively, according to DRAGICEVIĆ et al., (2012), heuristic methods 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


INTRODUCTION
Plenty of landslide susceptibility models have been developed by using GIS technology within the past three decades.These models 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 categorized according to various criteria, and always denoted as: quantitative or qualitative, direct or indirect, heuristic, probabilistic or deterministic, for differentiated or undifferentiated types of landslide mechanism.
One of the most applied methods in regional landslide modelling is the qualitative experience (knowledge) based -heuristic approach, which is considered a rather controversial topic in landslide 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).Alternatively, according to DRAGICEVIĆ et al., ( 2012), heuristic methods are empirical and subject to various levels of uncertainty, but 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 combine various conditioning factors and produce a susceptibility map at any scale, but despite this, it is rarely the case that a particular model can be applied that gives the same results for different study area.Also, using a large number of conditioning factors in a particular analysis is not a guarantee for successful susceptibility modelling.Factor selection will always be subjective, so the assessment of regional factors by the expert conducting 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 compared to heuristic approaches.In this regard, THIERY et al. (2007) reveal major limitations that concern all statistical methods.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 incorporate 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 susceptibility 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 available, 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 frequent 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 aspects of the model sensitivity (CATANI et al., 2013).It is suggested by ROSSI et al., (2010) that optimal susceptibility predictions 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 context, most recent research on landslide susceptibility model sensitivity 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 complex algorithms are being engaged to check the sensitivity of a landslide susceptibility model.Extraction of a particular landslide conditioning parameter and repetition of the validation process 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, reliability 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 always be performed before producing final susceptibility maps at all levels and scales.
When there are various constraints for some regions of interest such as: different scales of landslide conditioning parameters maps; limited amounts of geotechnical data of historical landslides; low numbers of registered landslides in databases; no access to multi temporal aircraft imagery; no resources to perform field surveys, etc., the heuristic approach is possibly the most suitable 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 Polog 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 susceptibility maps obtained from poor datasets is feasible and is a practical tool in landslide susceptibility zonation practice.

THE ARBITRARY -POLYNOMIAL METHOD
The approach used in this paper is named the arbitrary polynomial 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 assessment, and their selection is based on analysis of limited historical data of landslides that had a destructive effect on the infrastructure 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 susceptibility assessment.In our case, because of their importance for the analyzed region, they were considered as landslide conditioning parameters.For the same reasons, numerous landslide susceptibility 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-eters are taken into account when assessing global landslide susceptibility (e.g.LIN et al. 2017).
According to the method, a heuristic gradation is performed by assigning rating points to the landslide conditioning parameters.The value of the ratings is defined by expert assessment on the basis of particular conditioning factor importance for landslide development in a given region.
In order to investigate the sensitivity to changes of these ratings, 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 common 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 statistical exercise.

Lithology (LT-R)
Analysis of reports and thematic engineering-geological maps of historic landslides in the Republic of Macedonia showed that bedrock and superficial geology are among the most important factors 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 conditions 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 similar physical-mechanical behavior.The rating values were assigned to all lithological units, wherein higher values indicate a higher potential for landslide occurrence.

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, rating slope inclination is not a trivial problem, especially on regional 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 susceptibility 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 assigned even at the lowest value of 3° inclination (translational and rotational landslides in Pliocene sediments on relatively flat terrain), and for all terrain over 25° inclination (mostly rock falls and rockslides), the highest rating of 3 is assigned.This type of interpolation was used for the analyzed case study region.In more detailed statistical models, the slope inclination is treated separately according to the slope movement types (landslides, rockfalls, debris flows, etc).Since the presented approach is intended to give a general presentation of a region's susceptibility, the interpolation technique is considered as very practical.In cases where detailed databases on the slope movement type exist, different interpolation curves should be applied for each of the sliding mechanisms analysed.

Rainfall (AP-R)
Rainfall is governed by many factors in a region, especially elevation, mountain range orientation, climate zone, etc.In this context, 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 lifethreatening 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 landslide movement rates, rating values for this conditional parameter can be assigned arbitrarily according to, for example, average annual rainfall records.

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 reported 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 seismic data are available for a particular region, it can be rated in a similar manner.

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 Environmental Agency, ver.2006 was used as reference.A landslide inventory 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.

Calculation of Total and Weighted-Total Landslide Susceptibility Ratings
Each specific combination of conditioning parameters is connected 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: 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 (earthquake) 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 sensitivity coefficients.In this way, a Weighted Total Landslide Susceptibility 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 proportions/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 coefficients 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 argued 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 production 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: 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.

WLT-TLSR
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 interest.

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 km 2 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 infrastructure 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 regional 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, serpentinites, 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, metaconglomerates, 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.
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, intensity 8).It is worth mentioning that many other earthquakes ranging 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 characterized 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 precipitation 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 occurred 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 covers schistose or granitic bedrock, 11% in limestone (mostly rockfalls), and 31% in lacustrine sediments.Other units are less well represented.According to the landslide type, in the available reports, 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 metres.The remaining 40% of slides are shallow and mostly occurred 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 rainfall (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).

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 undifferentiated 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.Unfortunately, the exact date of occurrence is known for only a very small number of landslides.

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 elevation 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.

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 verification 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 susceptibility of the region.

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.

Analysis of the results
From Tables 1, 2 & 3, it can be seen that the landslide susceptibility 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, terrains 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  territory belongs to zones with very low, low and medium landslide susceptibility, with the presence of 28.49% of the total number 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 discussed later.The sensitivity does not show a straightforward case if we consider the landslide percentage per susceptibility class as a reference, 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 susceptibility 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 susceptibility is sensitive to slope inclination parameter change and insensitive 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 sensitive 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 lithological 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 variants, 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 sensitivity 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 number 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 correlation with a positive trend of increase both in the size and number of landslides.
This leads to the conclusion that the model WSI-TLSR100 is the most suitable to be applied for this particular case study region, since the size of the area under high and very high susceptibility has reduced by 41.7 km 2 from the initial model, and only 2 landslides had been moved to a lower susceptibility class zone.

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 various conditioning parameters are in use when performing landslide susceptibility or hazard zoning: lithology, engineering-geological characteristics of the rock masses, DEM derivatives, slope aspect, elevation, normalized difference vegetation index (NDVI) , and land use, landform classification, rainfall, runoff curve number (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 factors), subjectivity of the engineering-geological classification methods (rock mass classification), precision of measuring instruments, 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 (incomplete 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 representative landslide susceptibility zonation.
Our modelling results show more than 74% of landslides falling in zones of high and very high susceptibility classes, which is comparable with some more complex statistical studies.HE &BEIGHLEY, (2008) obtained 71% andKAVZOGLU et al., (2014) 69-77% precision by applying statistically based algorithms.Based on results presented in our research, it is obvious that there is also the possibility to test the sensitivity of landslide susceptibility models in a simple, yet acceptable manner.The optimization of preliminary prepared models should always be performed, 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 reliability for use in more advanced landslide susceptibility approaches, 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 instead 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 estimation of landslide susceptibility, and one can argue that they are based mostly on engineering judgment.They can however provide indications of the aspects that deserve attention while allowing for the undertaking of more advanced studies in the future.For example, to point out zones of some regions where more detailed datasets need to be collected in order to apply more advanced 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 datasets 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 organizational, technical and economic nature exist.It is simple and makes reasonable the use of the limited but readily available datasets 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 methods), 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 conditioning parameters.
As for any other landslide susceptibility modelling approach, the one presented in this paper also has some limitations and restrictions.These include the inability to use it for a coarser scale analysis where more exact correlation between the landslide conditioning factors and landslide occurrences needs to be supported by detailed statistical data, rainfall thresholds, and exact earthquake 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.000seismic intensity maps and rainfall maps.This required certain extrapolations to be undertaken, which certainly influenced the outcome of the study.The extent of this influence 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 inventory with undifferentiated mechanisms of landslide movement.Therefore, to improve the quality of the inventory map, and consequently the quality of the validation processes, further landslide susceptibility studies in the region should be supported by field checks.

Figure 1 .
Figure 1.Flow-chart of the arbitrary polynomial method.

Figure 2 .
Figure 2. a) Location of the Republic of Macedonia in Europe; b) Location of the case study region of Polog.

Figure 3 .
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.

Figure 4 .
Figure 4. Total Landslide Susceptibility Rating map for return period of maximum expected seismic intensity for a) 100 years; b) 500 years

Table 1 .
Landslide susceptibility classes for a return period of maximum expected seismic intensity of 100 years (TLSR-100).

Table 2 .
Landslide susceptibility classes for a return period of maximum expected seismic intensity of 500 years (TLSR-500).

Table 3 .
Results from verification of the Weighted -Total Landslide Susceptibility Rating models (WTLSR).

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.