Landslide susceptibility assessment of the City of Karlovac using the bivariate statistical analysis

A preliminary landslide susceptibility analysis on a regional scale of 1:100 000 using bivariate statistics was conducted for the City of Karlovac. The City administration compiled landslide inventory used in the analysis based on recorded landslides from 2014 to 2019 that caused significant damage to buildings or infrastructures. Analyses included 17 geofactors relevant to landslide occurrence and classified them into four groups: geomorphological (elevation, slope gradient, slope orientation, terrain curvature, terrain roughness), geological (lithology-rock type, proximity to geological contacts, proximity to faults), hydrological (proximity to drainage network, proximity to springs, proximity to temporary, permanent and to all streams, topographic wetness) and anthropogenic (proximity to traffic infrastructure, land cover using two classifications). Five scenarios were defined using a different combination of geofactors weighted by the Weights-ofEvidence (WoE) method, resulting in five different landslide susceptibility maps. The best landslide susceptibility map was selected upon the results of a ROC curve analysis, which was used to obtain success and prediction rates of each scenario. The novelty in the presented research is that a limited amount of thematic data and an incomplete landslide inventory map allows for the production of a preliminary landslide susceptibility map for usage in spatial planning. Also, this study provides a discussion regarding the used method, geofactors, defined scenarios and reliability of the results. The final preliminary landslide susceptibility map was derived using ten geofactors, which satisfied the pairwise CI test, and it is classified in four zones: low landslide susceptibility (57.05% of the area), medium landslide susceptibility (20.63% of the area), high landslide susceptibility (13.28% of the area), and very high landslide susceptibility (9.03% of the area), and has a success rate of 94% and a prediction rate of 93% making it a highly accurate source of preliminary information for the study area.


Introduction
Landslide susceptibility is defined as the spatial, time independent probability of landslides occurring in an area depending on the local terrain conditions (Guzzetti et al., 1999;Guzzetti, 2005). Landslide susceptibility can be obtained using different analytical approaches, i.e. heuristic, statistical and deterministic (Soeters and Van Westen, 1996). Statistical landslide hazard assessment has become very popular, especially with the use of Geographic Information Systems (GIS) and the possibility of applying data integration techniques that have been developed in other disciplines (Van Westen et al., 2005). Landslide susceptibility maps are based on the assumption that landslides are likely to occur under the same conditions as those under which they occurred in the recent past. Due to that reason, the preparation of landslide susceptibility maps requires a landslide inventory map used in combination with a series of environmental factors. All susceptibility maps depict a relative indication of the spatial probability of landslides in the form of zones. A landslide susceptibility map ranks the relative slope stability of an area into categories that range from stable to unstable. The identification and map portrayal of areas highly susceptible to damaging landslides are the first and necessary steps towards loss reduction (Mihalić Arbanas and Arbanas, 2015).
The first extensive papers on the use of spatial information in a digital context for landslide susceptibility mapping date back to the late seventies and early eighties of the last century. Among the pioneers in this field were Brabb et al. (1972) in California and Carrara et al. (1977) in Italy. All research on landslide susceptibility and hazard mapping uses digital tools for handling Rudarsko-geološko-naftni zbornik i autori (The Mining-Geology-Petroleum Engineering Bulletin and the authors) ©, 2021, pp. 149-170, DOI: 10.17794/rgn.2022.2.13 spatial data, i.e. Geographical Information Systems (GIS). Van Westen's dissertation (Van Westen, 1993) is the first comprehensive overview of the application of GIS technology in landslide hazard zonation, followed by the classification of landslide hazard analysis methods. A dissertation by Guzzetti (2006) provided numerous examples of analysis, assessment and zonation of landslide susceptibility, hazards and risks in Italy, which have been subjects of scientific research. An overview of the spatial data types required for landslide susceptibility assessments and the methods for obtaining these data were published in Van Westen et al. (2008). The methods and approaches proposed and tested to ascertain landslide susceptibility are shortly presented in Fell et al. 2008a,b. Despite numerous attempts and unquestionable progress, the general assumptions and the most popular methods and techniques used to assess landslide susceptibility have not changed significantly in the last few decades (Guzzetti, 2021). Reichenbach et al. (2018) critically reviewed the statistically-based landslide susceptibility assessment literature by systematically searching for and then compiling an extensive database of 565 peer-reviewed articles from 1983 to 2016. They have found that the Weight of Evidence (WoE) analysis is among the most common statistical method for landslide susceptibility modelling, together with logistic regression, neural network analysis, data-overlay, index-based and machine learning methods. WoE belongs to bivariate statistical analyses (BSA), one of the simplest statistical analysis methods and it is popular in numerous research fields, enabling environmental scientists to model various natural conditions. BSA techniques can be used as a simple geospatial analysis tool to determine the probabilistic correlation between dependent variables (produced using the inventory map of a landslide incidence) and independent variables (landslide causal factors) by computation of landslide densities and the significance of each factor. In BSA, the importance of each factor is investigated separately (Porwal et al., 2006;Guzzetti, 2021).
This paper describes landslide susceptibility modelling by adopting bivariate statistical analysis. The research was performed in the administrative area of the City of Karlovac (total area of 402 km 2 ), which possesses a database with 196 registered landslides in the form of a landslide inventory map. The study area is characterized by different geological and geomorphological conditions, including lowland and highland areas of the Pannonian Basin and Dinarides (see Figure 1a,b). An analyzed landslide inventory map (see Figure 1c) contains locations of all landslides registered by the City administration in the period from 2014-2019, and most of them were triggered by rainfall events during Cyclone Tamara (Mihalić Arbanas et al., 2017). Landslide and rainfall events are related to extreme weather conditions that are becoming more common in Croatia in recent years (Bernat Gazibara et al., 2018), causing material damage in a range of disasters. The application of the Weight of Evidence method, presented in this paper, using a series of landslide causal factors, resulted in the landslide susceptibility map on a small scale of 1:100 000.
The objective of this study is to derive a preliminary landslide susceptibility map on a small scale, using incomplete landslide inventory and limited input data about geofactors, for the application in spatial planning. The main tasks of the research were: (i) optimization of geofactors using a pairwise conditional independence (CI) test; (ii) landslide susceptibility analyses for five   scenarios using different combinations of landslide causal factor groups; (iii) verification of all susceptibility maps using AUC rates for defined scenarios. The research was organized to elaborate the assumption that a limited amount of input data can be used for a reliable prediction of "where" landslides are likely to occur.

Study area
The study area is defined by the administrative border of the City of Karlovac located in the central part of Croatia, encompassing the area of contact between two geotectonic units of Internal and External Dinarides (Schmid et al., 2008;Tomljenović et al., 2008) as well as the Pannonian Basin. Geomorphologically, most of the study area belongs to the megageomorphological regions of the Pannonian Basin and, to a lesser extent, to the Dinaridic Mountain System (Bognar, 2001). Consequently, the relief generally changes in the northeastsouthwest direction from plain (the fluvial floodplain of the Kupa River) through lowland (Crna Mlaka Depression) to hilly type (the hills of Vukomeričke Gorice, Utinjsko-Tušilovički Hills and Vojnić Hills) according to changes of geological settings (Mihalić Arbanas et al., 2017). The Karlovac Depression is surrounded by the lowest slopes of Samobor Hills on the northwest, the Vukomeričke Gorice on the northeast and it is separated by the Korana River from Kordun Foothills.
The elevation ranges from 101 to 372 m a.s.l. and the altitude increases to the west, to the southwest and the south, in the direction of relief changes. The prevailing slope angles (70% of the study area) are <5° whereas only 20% of the study area is in a range from 5° to 10° and 9% is in a range from 10° to 20°. The study area comprises 402 km 2 , where the current land cover includes about 20 km 2 of artificial surfaces, about 200 km 2 of agricultural areas, and about 180 km 2 of forests. An area of about 7 km 2 belongs to water bodies of four large rivers (the Kupa, Korana, Dobra and Mrežnica rivers) and other types of inland waters. The City of Karlovac is the administrative centre of the Karlovac County and its largest urban area with a population of 59,016 residents (DZS 2011). The study area is comprised of Quaternary, Neogene and Pre-Neogene sediments and mostly sedimentary rocks. Figure 2 presents geological characteristics of the City area according to the Basic Geological Map on a scale of 1:100 000, compiled from the Karlovac Sheet ) and the Črnomelj Sheet (Bukovac et al., 1983). A detailed description of geological units is given in Madaš et al. (1989). Quaternary deposits prevail at the surface of the study area (about 85% of the area). The Holocene sediments are of fluvial origin, composed of sand, sandy clays, clays, gravels, silt and bog sediments. The Pleistocene deposits are loess and heterogeneous mixtures of mostly impermeable clayey soils. The Miocene and Pliocene deposits are present only at the surface in part of the hilly area (about 4% of the area) in the form of stratified sandstones, conglomerates, marls and limestones, as well as sands, gravels and clays. The Pre-Neogene rocks are composed of Mesozoic deposits of the Triassic, Jurassic and Cretaceous periods (about 4% of the area) and Upper Paleozoic deposits (about 7% of the area). Mesozoic deposits are composed of various types of chemogene, clastic, intrusive, effusive and metamorphic rocks due to complex geotectonic history, including flysh-type rocks. The Upper Paleozoic is composed of quartz conglomerates, sandstones, siltstones and shales.
The City of Karlovac is highly affected by landslide hazards. The climate of the City is continental with a mild maritime influence, with mean annual precipitation (MAP) of 1 100 mm and with precipitation mostly recorded in the period from May to December (Zaninović et al., 2008). However, landslides are also controlled by extreme hydrological events, i.e. floods which are caused by regional events, such as Cyclone Tamara from 2014, when the heavy rainfall triggered more than 100 landslides in the Karlovac area during one rainfall event. Due to heterogenic geological and geomorphological conditions, landslide types are of various sizes, from small to medium-large landslides (<10 hectares). The preparatory causal factors of slope instabilities in the study area depend on the geomechanical properties of rocks and soils, geomorphological conditions and processes, including river erosion and different types of man-made processes.

Data and methods
This chapter describes input data for landslide susceptibility analyses, including a landslide inventory and spatial data divided into groups of geomorphological, geological, hydrological and anthropogenic geofactors. In this study, bivariate statistics was used to determine landslide susceptibility. The Weight of Evidence method was applied to calculate the weight factors of classes for each factor map and a pairwise test was used to determine factor maps' conditional independence.

Data
The input data for the landslide susceptibility analysis of the City of Karlovac area includes the landslide inventory and other spatial data used to derive factor maps. Landslide inventory was obtained from the City administration in vector point files indicating the locations of 196 landslides recorded in the period from 2014 to 2019. The landslide inventory was compiled based on information received from citizens or road patrol who informed the City administration responsible for landslide remediation or civil protection. Most reported landslides have damaged infrastructure, mainly roads and private properties. The landslide inventory of the City of Karlovac was split into two data sets by random selection in ArcGIS 10.8, a modelling set (140 landslides in the inventory) for modelling landslide susceptibility and a validation set (56 landslides in the inventory) (see Figure 3).
The landslide causal factors for the landslide susceptibility modelling in the City were selected expertly and based on data availability. The existing data for this study area includes topography in the form of a Digital Elevation Model (EU-DEM) with a resolution of 25 x 25 m downloaded from Copernicus Land Monitoring Service (EEA, 2018), digitalised geologic data from the Basic Geologic Map Bukovac et al., 1983) in a scale of 1:100 000, digitalised hydrological data (superficial streams and springs) from the topographic map in a scale 1:25 000, land cover data of the CO-RINE Land Cover dataset (NRC/LC, 2018) and traffic infrastructure from Open Street Maps (OSM, 2019). Most factor maps are derived by further processing (i.e. geomorphometric, topographic and distance analysis) of the original input data. Only land cover data and geological units are applied as original data in the model. A total of 17 landslide-conditioning factors, such as elevation, slope gradient, slope orientation, terrain curvature, terrain roughness, lithology (rock type), proximity to geological contact and faults, proximity to drainage network, springs, temporary streams, permanent streams and all streams, topographic wetness, land cover (a), land cover (b) and proximity to traffic infrastructure were used in the study. All analyses were performed in ArcGIS 10.0 software, where all 17 factor maps were converted to raster format with a spatial resolution of 25 m. The following section describes all derived factor maps.
Elevation is a very frequently used parameter in landslide susceptibility studies because landslides may form in specific relief ranges (Dai et al., 2001). The relief map of the study area generated from the EU-DEM is shown in Figure 4a. The range of elevations shown on the DEM indicates a difference in altitude of 271 m. A factor map of elevations was created by reclassifying the 25 m resolution DEM.
Slope gradient is often considered to be the most important morphometric parameter used to more effectively analyze and describe relief (Van Westen et al., 2008). A slope gradient factor map was created using the Spatial Analyst extension (Slope tool) in the ArcGIS 10.0 software, and it represents the spatial distribution of slope angle values in the range from 0 to 90 degrees. The slope values in the study area (see Figure 4b) range between 0° and 40°. However, the mean value of the slope is 3.97°, with a standard deviation of 4.04°.
The slope aspect identifies the downslope direction of the maximum rate of change in the elevation value from Surface curvature is the curvature of a line formed by intersecting a plane (in some chosen orientation) with the terrain surface. The curvature map was created with the Curvature tool (3D Analyst) in the ArcGIS 10.0 software. The curvature shapes define the possible processes on the slopes, whether the flow slows down or accelerates on the slope. The values for the study area vary between -6.84 and 5.64 (see Figure 4d). The mean value of the profile curvature data is 2.46, with a standard deviation of 0.96.
A terrain roughness map (see Figure 4e) was created using the Roughness tool in the Geomorphometry and Gradient Metrics Toolbox (Evans et al., 2014). Terrain Ruggedness Index (TRI) is defined according to Riley et al. (1999), and the algorithm calculates the sum change in elevation between a grid cell and its eight neighbour grid cells. Therefore, an increased TRI shows more significant local relief heterogeneity.
Geological data were obtained from the Basic Geological Map on a scale 1:100 000, the Karlovac Sheet ) and the Črnomelj Sheet (Bukovac et al., 1983). The digitalisation of geological maps and further spatial analysis resulted in three factor maps: lithological map, which shows chronostratigraphic units (see Figure 4f); proximity to geological contacts; and proximity to faults. Geological data obtained by digitizing the Basic Geological Map 1:100 000 are shown in Figure 2.
A drainage network is a set of all drainage systems in an area, i.e. a set of natural canals through which water constantly flows or temporary, and which are connected into a single stream and represent the smallest independent geomorphological component (Marković, 1983). In this paper, the drainage network was derived from a 25 m resolution DEM using several tools in the Spatial Analyst -Hydrology Toolbox in ArcGIS 10.0 and is presented as vector data in the shape of a line, as shown in Figure 4g.
The digitalisation of the topographic map with a scale 1:25 000 (TK25) available on the WMS server of the State Geodetic Administration resulted in vector data of superficial streams in the shape of a line and springs in the shape of a point (see Figure 4h). Streams were classified as permanent and temporary. Based on the digitalized input data, four factor maps were created: proximity to springs, proximity to (all) streams, proximity to temporary streams and proximity to permanent streams.
Topographic wetness is a steady state wetness index, and it is commonly used to quantify topographic control on hydrological processes. The topographic wetness map was derived from a DEM using the Compound Topographic Index tool in Geomorphometry and Gradient Metrics Toolbox (Evans et al., 2014) in ArcGIS 10.0 software, and it represents the Compound Topographic Index (CTI) shown in Figure 4i. Lower values of CTI indicate lower wetness (hill tops), and higher CTI values indicate higher terrain wetness (plain areas).
The traffic infrastructure network, including roads and railways, is obtained from the Open Street Map website in the form of vector data, as shown in Figure  4j. The factor map showing proximity to the traffic infrastructure was created based on the input data.
Information on the land cover was downloaded from the Copernicus Land Monitoring Service (Corine Land Cover) website as accessible raster files ready for use in GIS. Land cover, according to the Corine Land Cover (CLC), is classified into three hierarchically organised levels of detail. The first and second levels of classification were used in this paper, as shown in Table 1

Methods
Landslide susceptibility is a quantitative or qualitative assessment of the spatial probability of a particular type of landslide to occur in an area. An overview of landslide susceptibility assessment methods is given in many papers, e.g. in Corominas et al. (2013) and Reichenbach et al. (2018). Qualitative methods are inventory-based and knowledge-driven methods. Quantitative methods are based on data (data-driven, statistical methods) and on physically-based models. In this paper, a bivariate statistical method referred to as Weights-of-Evidence (Agterberg et al., 1990; Bonham-Carter et

al., 1989; Bonham-Carter, 1994) is applied in a GIS
environment to derive quantitative spatial information on the predisposition to landslides. Because of the size of the study area (402 km 2 ) and limited data availability, the WoE method is considered as the appropriate approach. General schematic representation of the processes of map development according to Van Westen et al. (2002) is shown in Figure 5. Each factor map overlaps with the landslide inventory map for the purpose of obtaining the frequency of landslides in each class of all factors. By comparing the landslide density in the factor classes with the landslide density in the entire study area, the relative influence of the observed landslide factor is determined. After the calculated densities, weighting factors can be determined using different methods (Coe et al., 2004). The disadvantage of the bivariate method is that it starts from the assumption that the sliding factors are mutually independent, which is not usually the case in natural environments (Van Westen et al., 2002). To mitigate the influence of dependent factors, i.e. violating the conditional independence, a pairwise test is applied to all factor maps and their classes. The Weight of Evidence method was developed by the Canadian Geological Survey (Agterberg et al., 1990;Bonham-Carter et al., 1989), and it was used to map mineral resources. Sabto (1991) applied the above method to landslide hazard analysis. This method is easy to use in the GIS interface, and Table 2 shows the variables used to calculate the weighting factors. (3)

Where:
Npix 1 -number of pixels for present landslides and present potential landslide causal factor; Npix 2 -number of pixels for present landslides and absent potential landslide causal factor; Npix 3 -number of pixels for absent landslides and present potential landslide causal factor; Npix 4 -number of pixels for absent landslides and absent potential landslide causal factor; W i + -positive weights (indicating the importance of the presence of the factor); W i --negative weights (indicating the importance of the absence of the factor); W map -weight of evidence value (weight factor). As stated by Bonham-Carter (1994), the pairwise conditional independence (CI) test should be run on all combinations of factor maps used in the analyses using bivariate statistics, i.e. the Weight of Evidence method. The pairwise CI test can reveal the magnitude of violating the assumption of factor maps independence and identify which maps are causing the most dependence. Using Equations (4), (5) and (6), the value is calculated for each pair of the factor maps (a and b) and later compared to the tabulated values of . The tabulated values depend on the significance level p selected expertly and the degree of freedom d calculated using Equation (7). Furthermore, for corrections on the factor maps showing dependence, a Yates correction (Walker and Lev, 1953) was done according to Equation (8). (4) Where: a, b -studied factor map a and factor map b, respectively; k -multiplication of the number of classes in the two studied factor maps i.e.: (5) where: x, y -the number of classes in factor map a and factor map b, respectively; f i(o) -number of observed landslides in the overlap area of two classes of two studied factor maps; f i(e) -expected number of landslides in the overlap area of two classes of two studied factor maps i.e.: where: N -total amount of observed landslides.
(8) The variables shown are four possible combinations obtained after overlapping a landslide map (binary landslide map) with a weight map (binary variable map). After defining the variables, Equations (1) and (2) are used to calculate the positive and negative weights, respectively. The final weight of evidence value, used for landslide susceptibility analyses, is acquired using Equation

Results
The GIS-based landslide susceptibility assessment for the City of Karlovac comprises the following main processing steps: (i) calculation of weights for factor maps using the modelling set of the landslides by applying the Weight of Evidence method; (ii) testing the CI using pairwise test; (iii) calculation of the posterior probability map (i.e. combination of the controlling factors to predict po-tential landslide occurrences); (iv) developing four different scenarios to determine landslide causal factor group importance; (v) modelling validation using the validation set of landslides from the inventory, and (vi) classification of the landslide susceptibility map.

Weight maps
In this study, 17 factor maps were created for the landslide susceptibility analysis. Derived factor maps can be divided into four groups of landslide causal factors: (i) geomorphological factors, including an elevation map, a slope gradient map, a slope aspect map, a terrain curvature map, and a terrain roughness map; (ii) geological factors, including a lithology map, proximity to the geological contacts and faults; (iii) hydrological factors, including proximity to the drainage network, springs, temporary streams, permanent streams and all streams and a topographic wetness map; and (iv) anthropogenic factors, including proximity to traffic infrastructure, two land cover maps derived based on the first and second level of hierarchical division. Classes of stretched uncategorical factor maps were defined expertly, and they are listed in the following paragraphs together with the most important relations to landslides. The slope gradient map is divided into five classes: 0-5°; 5-10°; 10-15°; 15-20°; and 20-40°. The greatest influence on the occurrence of landslides is given by the slope gradient classes 10-15°; and 20-40°, while in the 0-5° class, fewer landslides were observed than expected.
The slope orientation map is divided into nine classes: flat terrain, north (N), northeast (NE), east (E), southeast (SE), south (S), southwest (SW), west (W) and northwest (NW). Landslides were not recorded on flat terrains, while most of them were recorded on slopes oriented to the northeast and south, and it is concluded that these slopes have the greatest impact on the occurrence of landslides.
The terrain roughness map is divided into four classes: smooth terrain (0-1); slightly rough terrain (1-2); rough terrain (2-3); and very rough terrain (3-5). There are only 12 landslides in the class of smooth terrain and very rough terrain (significantly less than the expected number of landslides), while the classes of slightly rough terrain and rough terrain have a significant impact on the occurrence of landslides.
Lithology (rock type) is divided into 29 classes according to the chronostratigraphic units defined on the Basic Geological Map on a scale of 1:100 000 (see Figure 2). The class with the most significant influence on the occurrence of landslides are Pliocene and Pleistocene deposits (Pl, Q): sands, gravels, clays, sandstones and conglomerates. A significant number of landslides were recorded in the class of alluvium (aQ 2 ): sands, sandy clays, clays, gravels, silt, bog sediments.
The proximity to the geological contacts map is divided into five classes: 0 -50 m; 50 -100 m; 100 -200 m; 200 -400 m; 400 -800 m; and 800 -3800 m. In class 0 -50 m, 49 landslides were recorded, and in class 50 -100 m, 36 landslides were recorded. Only 4.38 landslides were expected in the mentioned classes, so it is concluded that most landslides occur near the geological contact (from 0 to 100 m). In the class proximity greater than 800 m to the geological contact, 105 landslides were expected, and only 17 were recorded.
The The proximity to springs map is divided into five classes: 0 -100 m; 100 -250 m; 250 -500 m; 500 -1000 m; and 1000 -12 000 m. Class 0 -100 m has a very small area with only one recorded landslide. In the classes 100 -250 m, 250 -500 m and 500 -1000 m, significantly more landslides were recorded than was expected.
The proximity to temporary streams is divided into six classes: 0 -100 m; 100 -200 m; 200 -300 m; 300 -400 m; 400 -500 m; and 500 -2,300 m. Given the difference between the number of mapped and expected landslides, all classes except 500 -2,300 m indicate the possibility of landslides.
The proximity to permanent streams map is divided into six classes: 0 -100 m; 100 -200 m; 200 -300 m; 300 -400 m; 400 -500 m; and 500 -3,600 m. The greatest impact is given by the class 0 -100 m because a relatively large number of landslides were recorded near permanent streams.
The proximity to all streams map is divided into six classes: 0 -100 m; 100 -200 m; 200 -300 m; 300 -400 m; 400 -500 m; and 500 -1,800 m. Similar to the proximity to temporary streams map, all classes except the 500 -1,800 class indicate the possibility of landslides, i.e. have more observed than expected landslides.
The proximity to the traffic infrastructure map is di-   pected and mapped landslides, and it is concluded that the roads contribute to the instability of the slopes. However, it should be considered that most of the landslides were recorded by the municipal warden in charge of road maintenance. Land cover map A, derived from the first level of hierarchical classification, has four classes: artificial surfaces, agricultural areas; forests and semi-natural areas; and water bodies. In the class of agricultural areas, the expected number of landslides is 68.93, and 106 land-slides were recorded. From the above, it is concluded that the most significant impact on the occurrence of landslides has the class of agricultural areas.
The land cover map B, derived from the second level of hierarchical classification, is divided into nine classes: urban fabric; industrial, commercial and transport units; artificial, non-agricultural vegetated areas; arable land; pastures; heterogeneous agricultural areas; forests; scrub and/or herbaceous vegetation associations; and inland waters. Almost all mapped landslides are in the

Pairwise CI test
Pairwise conditional independence (CI) test between all classes of all causal factors was tested using the pairwise CI test before deriving landslide susceptibility maps. Table 4 shows χ 2 values calculated for each combination of weight maps using Equations 4, 5 and 6. Tabulated χ 2 values are acquired by calculating the degree of freedom (see Equation 7) for the significance level of 1%.
Conditional independence or conditional dependence for each pair of weight maps is determined by comparing the two χ 2 values. Greater calculated χ 2 value than the tabulated shows that the pair has conditional dependence, i.e. the null hypothesis of conditional independence is rejected.  Firstly, weight maps that contain similar information are reduced to one map. For example, between four maps proximity to drainage network, proximity to temporary streams, proximity to permanent streams and proximity to all streams, only one map is chosen as optimal weight map, proximity to drainage network because it was derived using a DEM thus including both geomorphological and hydrological information. Land cover (A) is rejected as it showed a more conditional dependence with other weight maps compared to the land cover (B). Besides, weight maps that showed significant conditional dependence with the remaining weight maps were excluded. This filtering resulted in removing the elevation, the slope gradient and the terrain curvature weight map, as they showed conditional dependence with several remaining weight maps. Finally, for the remaining conditional dependent weight map pairs, the Yates correction (Equation 8) is applied, resulting in less CI violation. The Yates correction is necessary to apply for small expected frequencies (Bonham-Carter, 1994), however, in this paper, it was applied for the mitigation of CI limitations, reducing the amount of dependence pairs. Table 5 shows the final selection of weight maps after the three-step filtering used for landslide susceptibility modelling.

Landslide susceptibility modelling, verification and zonation
For each class of the selected factor maps listed in Table 5, a weight value was defined applying the WoE method resulting in the W map value (see Table 3). This study created the landslide susceptibility maps by summing (overlapping) the factor maps according to the assigned weight values based on five scenarios. The first, Scenario 0, is used to create the final landslide susceptibility map as it contains all the weight maps which passed the three-step filtering. The other four scenarios are defined to determine which landslide causal factor group has the most significance to landslide susceptibility modelling. Therefore, in each of the four scenarios (Scenario I-IV) different group of landslide causal factors is used, as shown in Table 6.
The success and prediction rate of all scenarios was determined for all five landslide susceptibility maps, and the results are shown in Table 6. Furthermore, ROC curves for Scenario 0 used for the final susceptibility map are shown in Figure 6, presenting the success and prediction rates. Scenario 0 was selected as the optimal in order to use all eight weight maps which passed the three-step filtering for the final landslide susceptibility map. How  ever, four additional scenarios were defined, proving that only the anthropogenic group of landslide causal factors causes significant changes in the AUC values. Scenarios I, II and III resulted in almost equal AUC values compared to the optimal Scenario 0. In contrast, Scenario IV showed a success rate of only 79% and a prediction rate of 74%, which is significantly lower than the other scenarios. The landslide susceptibility map created based on Scenario 0 using the Weight of Evidence method was classified into four classes (low, medium, high, and very high) using the cut-off values defined according to Bernat Gazibara (2019) (see Figure 7). The classification was done using the ROC curve derived from all landslides, both the modelling set and the validation set.

Discussion
The proposed methodology to assess landslide susceptibility at a 1:100 000 scale is based on a bivariate method calibrated on a modelling set of landslides and tested by a validation set of landslides. The limitation of the used landslide inventory of the City of Karlovac is that it was compiled based on information collected by the municipal warden or citizens. It mainly contains reported landslides that caused damage to either buildings or infrastructures. This means that numerous landslides in green areas (forests, pastures, etc.) are missing from the inventory. The mean landslide density of 0.48 landslides/km 2 at the area of the City is significantly lower than expected in terrains of similar geoenvironmental conditions (Bernat Gazibara et al., 2019). However, since landslides are distributed across the whole study area, for the purpose of this research, which is a preliminary landslide susceptibility map, the inventory is proven to be representative.
Modelling was started with 17 weight maps, followed by a pairwise CI test which resulted in a dependence of numerous combinations between weight map pairs. The main reason for this is that several weight maps included similar information and needed to be filtered out and excluded from the modelling. Further filtering and corrections lowered the CI violation to its minimum, ending up with only 3 out of 45 possible combinations showing dependence. The proposed method resulted in satisfactory AUC values for both success and prediction rates even though the modelling was done without a slope gradient map, commonly used in landslide susceptibility modelling. This proves that reducing the CI violating, i.e. decreasing the distortion of the results in a bivariate method, can be successfully performed by using a pairwise CI test without losing the performance of the model measured in AUC values.
Therefore, the final landslide susceptibility map was selected based on a maximum weight map, which contained minimum conditional dependence. Usage of 10 unique landslide causal factors with a significant number of classes in each provided that the final zonation has a detailed range of information. Final landslide susceptibility map zonation greatly influences its application, because the distribution of the zones could have an effect on spatial planning and land use management. The western and southern parts of the study area present a good transition from low to very high landslide susceptibility zones, as well as equal zone distribution. On the other hand, the north-eastern and central parts of the study area are generally classified as a low landslide susceptibility zone, which is expected considering the landslide causal factors at the area. The mentioned areas are predominately plain terrains where landslides are unlikely to occur. However, even in such areas, in a much lesser surface area, high and very high landslide susceptibility zones are distinguishable.
In this paper, map zonation was done using the ROC curve cut-off values defined by Bernat Gazibara (2019), which resulted in a highly accurate and reliable classified landslide susceptibility map. Only 9% of the study area classified as very high landslide susceptibility contains 85% of mapped landslides, making it an efficient ratio pointing out the most susceptible areas. Therefore, additional zonation methods were not applied in this study since the mentioned method provided satisfactory information about specific areas in the City of Karlovac regarding landslide susceptibility. However, it is important to point out the importance of the zonation method since different methods can result in very different zone distributions, leading to different map interpretation and its application. Such is the case with the proximity to traffic infrastructure, where the high influence of a single class shown in the Weight of Evidence analysis is well preserved even after the final overlapping of all weight maps proving that even a single factor can greatly impact the final zonation map.
Four additional scenarios were defined to test the influence of each group of landslide causal factors on the model performance. The results showed a direct link between the anthropogenic group of landslide causal factors and the AUC values. It is the only scenario that resulted in significantly different results compared to the final Scenario 0. Removing proximity to traffic and land use (B) resulted in far worse results in the AUC values, which is expected since 50 to 100 m proximity to traffic class in the proximity to traffic weight map has a significant W map value. Since the landslide inventory mainly consists of landslides near the traffic infrastructure, the AUC values are directly linked to the landslide inventory. This observation proves how certain classes of weight maps can directly influence the quality of landslide susceptibility modelling because of their direct correlation to the landslide inventory, i.e. proving the significance of using a representative landslide inventory. Also, using a representative proximity to traffic weight map, which proved to be an important precondition for the occurrence of landslides, should be completed by roads from satellite images on a larger scale in order to deliver representative W map values. The importance of a representative landslide inventory and representative weight maps is crucial to receive detailed results and high precision conclusions regarding certain geofactors and their relevance to landslide susceptibility. Also, using contrast and studentized contrast in the classification process would allow for the definition of different breakpoints when creating factor maps, however, in this paper, classes were defined expertly. Different breakpoints, i.e. the domain of each class in the weight maps leads to different landslide susceptibility modelling, which was not studied in this paper. According to Neuhäuser et al. (2011), contrast and studentized contrast allow the classification to reflect the original spatial association of landslide causal factors and the landslides.

Conclusions
This paper has demonstrated the necessity of using specific and adapted procedures for indirect landslide susceptibility assessment by bivariate methods, especially at a 1:100 000 scale, for complex environments with some uncertainty in landslide inventory data. The proposed procedure, based on a limited amount of thematic data and a landslide inventory map, consists of: data collection, spatial data processing in GIS, development of weight maps and defining of geofactor classes, analysis of the impact of individual classes of geofactors on landslides, CI analysis and weight map filtering, development of landslide susceptibility maps, verification of the susceptibility maps, and reclassification of the final susceptibility map into four classes of landslide susceptibility. The novelty in the presented research is that a limited amount of thematic data and an incomplete landslide inventory map allows for the production of a preliminary landslide susceptibility map for usage in spatial planning.
A reliable landslide susceptibility map was derived based on the analysis of ten geofactors: slope aspect, terrain roughness, lithology (rock type), proximity to geological contacts, proximity to faults, topographic wetness, proximity to springs, proximity to drainage network, proximity to traffic infrastructure and land cover (B classification). Pairwise CI test proved to be a successful method as it excluded seven out of ten geofactors, greatly influencing the overall performance of the model. The final preliminary landslide susceptibility map has a success rate of 94% and a prediction rate of 93%. The results proved a good correlation between the used geofactors and the occurrence of landslides. Three out of four additional scenarios resulted in similar success and prediction rates demonstrating that even with fewer geofactor maps, i.e. having even more limited data the reliability of susceptibility modelling would remain satisfactory. An exception being the anthropogenic geofactor group closely related to the landslide inventory map, as explained in the previous section.
Susceptibility modelling was done with easily accessible and high cost-benefit ratio input data, defining a robust and reproducible procedure for preliminary landslide susceptibility assessments. Our work shows that applying the simple Weight of Evidence method in the statistical model can satisfactorily recognise landslideprone areas in a complex environment on a small scale of 1:100 000. In case of more complete and more representative landslide data, as well as DEM of higher resolution, it would be possible to create a susceptibility map on a medium or large scale. Moreover, further analyses would be necessary to collect more detailed input data for all factor maps to create a more reliable landslide susceptibility map.
Considering the availability and resolution of spatial data used to analyse landslide susceptibility in the City of Karlovac, there are also limitations for applying this map. Namely, the derived preliminary landslide susceptibility map can be used as a basis for future research to narrow down the study area to a larger scale because it provides preliminary information about the study area on a small scale. Taking into consideration that 85% of the landslides are zoned in only 9.03% of the study area provides relevant information regarding the priorities for the City's administration management. Overlapping the zonation done in this paper with settlement boundaries can also provide information on a larger scale, point-ing our settlements in high or very high risk zones of landslide hazard. Most of the City's urban area is in the high and very high landslide susceptibility zone, which implies the necessity of further landslide risk management. The derived landslide susceptibility map can be used as a guideline for selecting areas for landslide susceptibility, hazard and risk zonation on a larger scale. Further landslide research in the City of Karlovac will lead to a better understanding of landslide prone areas, which is crucial for improving land use policy and reducing damage to buildings and infrastructures.
Analyses performed in this study resulted in a final version of the preliminary landslide susceptibility map that can be used for informative purposes as a basis for select areas in conducting further spatial analyses of landslide hazards.