Regional-Scale Landslide Susceptibility Mapping Using Limited LiDAR-Based Landslide Inventories for Sisak-Moslavina County, Croatia

In this paper, for the first time, a regional-scale 1:100,000 landslide-susceptibility map (LSM) is presented for Sisak-Moslavina County in Croatia. The spatial relationship between landslide occurrence and landslide predictive factors (engineering geological units, relief, roughness, and distance to streams) is assessed using the integration of a statistically based frequency ratio (FR) into the analytical hierarchy process (AHP). Due to the lack of landslide inventory for the county, LiDAR-based inventories are completed for an area of 132 km2. From 1238 landslides, 549 are chosen to calculate the LSM and 689 for its verification. Additionally, landslides digitized from available geological maps and reported via the web portal “Report a landslide” are used for verification. The county is classified into four susceptibility classes, covering 36% with very-high and high and 64% with moderate and low susceptibility zones. The presented approach, using limited LiDAR data and the extrapolation of the correlation results to the entire county, is encouraging for primary regional-level studies, justifying the cost-benefit ratio. Still, the positioning of LiDAR polygons prerequires a basic statistical analysis of predictive factors.


Introduction
Over the last few decades, we have witnessed frequent occurrences of natural hazards, among which landslides are no exception. This trend is expected to continue in the future, according to [1], due to increased urbanization, continued deforestation, and increased regional precipitation caused by changing climate patterns. Keeping in mind that climate changes undoubtedly affect the stability of natural and engineered slopes [2,3], it is necessary to emphasize the need for future landslide risk reduction. Since landslides have a great impact on humans, their property, and the environment [4], knowledge improvement and raising awareness among the wider public should be highlighted. Considering landslides, one way that helps land-use planners and policy-makers to achieve sustainable land management is terrain zonation according to slope stability through susceptibility, hazard, and risk assessment [5,6].
Landslide susceptibility maps (LSM) show the spatial probability of landslide occurrence [7,8]. In areas where the possibility of their occurrence exists, the production of LSMs is an important step to define the geoenvironmental characteristics of the area in the aspect of spatial planning. The zonation of terrain according to the degree of landslide susceptibility is the foundation of rational land-use management, with an emphasis on safe and planned construction.
The basic concept of landslide-susceptibility assessment includes the spatial distribution of prevalent geoenvironmental factors to determine areas prone to landsliding without temporal implications [9]. LSM divides the area of interest into zones of the same region of the Pannonian Basin, which is subdivided into the macrogeomorphological regions of the northwestern Croatian valley and Slavonian massive mountains [23].
In Sisak-Moslavina County, 6% of the area is settled, 30% of the area is agricultural land, 40% is forested [24], and more than 18% of the area is classified as other agricultural and forest areas [25]. From the pedological perspective, gley, pseudogley, eluvial-illuvial, and cambic soils are dominant in this region [26].  [27].

Geological Setting and Description of Geological Units
The description of the geological units and development of the study area is mainly based on the BGM guidelines and the Mining-Geological Study of Sisak-Moslavina County [28]. As a part of the study, a 1:100,000-scale geological map of the county was produced [27], based on nine BGM sheets.
The oldest rocks, Paleozoic granites ( The oldest rocks, Paleozoic granites (┌), and various older metamorphic rocks (Mi), are found in the mountains of Moslavačka Gora and Papuk. Devonian and Carboniferous marine sediments in the Trgovska Gora Mountain are deposited in the flysch-turbiditic environment, represented by shales, sandstones, and limestone interlayers (C, D). Turbidite sedimentation continued during the Permian (P), forming interlayers of clays, siltites, and sandstones in lower horizons and gradually, coarser sediments like sandstones and conglomerates in the Late Permian (Petrova Gora Mountain).
Early Triassic shallow marine sediments (sandstones, siltites, shales, marls, and limestones-T1) are followed by deeper marine limestones, marls, cherts, and even pyroclastites (T2). In the Late Triassic, a shallow marine environment dominated, producing dolomites (T3). During the Middle and Late Jurassic, the spreading of the lithosphere in the Banija region caused uplifting of the ultramafic rocks and the formation of ophiolitic complexes. They are composed of various magmatic (diabase, spilites-ββ), metamorphic (shales, ), and various older metamorphic rocks (Mi), are found in the mountains of Moslavačka Gora and Papuk. Devonian and Carboniferous marine sediments in the Trgovska Gora Mountain are deposited in the flysch-turbiditic environment, represented by shales, sandstones, and limestone interlayers (C, D). Turbidite sedimentation continued during the Permian (P), forming interlayers of clays, siltites, and sandstones in lower horizons and gradually, coarser sediments like sandstones and conglomerates in the Late Permian (Petrova Gora Mountain).
Flysch sediments were continuously deposited at the edge of the former Pannonian basin from the Late Cretaceous to the Paleogene. The structural characteristics of 100 to 150 m thick Senonian conglomerates, sandstones, and marls (K 2 3 ) indicate the distal type of turbidite sequences [29]. Paleocene flysch (Pc) is dominantly characterized by sandstones, shales, and marls, with the thickness of the turbiditic sequences diminishing from approximately one to two m to several decimeters in distal zones [30]. Paleocene-Eocene flysch (Pc, E) is characterized by conglomerates, sandstones, and shales, rarely with marls and limestones, with an approximate thickness of 200 m.
After a period of flysch sedimentation, the shallowing of the area is documented with freshwater sedimentation during the Eocene-Oligocene, when the river, lacustrine, and marsh sediments were deposited with local occurrences of coals (conglomerates, sandstones, marls, clays, and coal-E, Ol). The thickness of the whole series is 150 to 200 m, with 50 m of brown coal that was experimentally exploited. During the Oligocene and Early Miocene, the North Croatian Basin is characterized by long emersion and strong tectonics, which formed new pre-Neogene relief [31]. In the Ottnangian, the terrain subsided, which enabled the deposition of gravels, sands, and clays with coals (M 2 ), of which the estimated thickness is 250 to 350 m. The short terrigenous phase during the Karpatian was followed by marine transgression in the Early Badenian. Depending on preset terrain morphology, either freshwater, brackish or marine clastites, or evaporites and coals deposited (M 4 ). Firstly, at the base of the Badenian sequence, coarser clastic sediments deposited (conglomerates, gravels, and sandstones). The carbonate composite rises with the progress of marine transgression (up to > 90%) during the Late Badenian, which enabled the deposition of reef limestones and, in the deeper sea, marls [32]. The thickness of Badenian sediments varies greatly, with a maximum of 300 m.
A similar depositional environment continued during the Sarmatian (M 5 ). In coastal and shallow zones of the basin, coarser-grained clastics composed of polymictic conglomerates, gravels, sandstones, and sands were formed. Bioclastic limestones are characteristic of the shallow marine environment and of thinly plated marls with rare interbeds of turbiditic sands and sandstones deposited in deeper seas. The thickness of the Sarmatian deposits is, at maximum, 100 m.
Sedimentation continued from the Sarmatian to the Pannonian, in a high-energy kaspibrackish environment with sands and sandy marls [29] or low energy shoreline with clayey limestones. The laminated and platy limestones with few-centimeter-thick beddings are very abundant (Croatica Fm.) in the Lower Pannonian. In the Middle Pannonian, laminated to thin-layered sands and sandy marls alternate with clayey limestones, and poorly lithified sandstones. During Late Pannonian subsidence, massive to weakly layered marls and sporadic clayey limestones and clays were deposited in the more distal part of the lacustrine basin (Medvedski breg Fm.; [33]). The thickness of the Pannonian deposits (M 6 ) is between 120 to 200 m.
During the Pontian, the region was gradually uplifting with constant reduction of the sedimentation area. Clayey and silty marls deposited in the Early Pontian were gradually replaced by silts and sands in the Late Pontian (M 7 ) [29]. According to the latest lithostratigraphic nomenclature, Pontian deposits are subdivided into the Andraševec and Hum Zabočki Fm [34,35]. The Andraševec Fm is composed of the alternation of sands, silts, and marls [36] and is continuously overlain by the Hum Zabočki Fm, which is mostly composed of fine-to medium-grained sands, silty marls, clays, and coal seams. The thickness of the Pontian deposits varies between 150 to 450 m.
After the Pontian, this region was tectonically uplifted, creating the protruded relief of the Trgovska and Petrova gora mountains. For the Pliocene, clays, sands, and gravels deposits (Pl) with a thickness of 200 to 400 m [37], associated with lacustrine and river environments, are characteristic. Pliocene sediments are equivalent to the Vrbova Fm. of the Sava basin [33,38]. During this time, the foresaid mountains became the main source of molasse gravels, sands, and clays that were deposited during the Plio-Quaternary. Plio-Quaternary clastics (Pl, Q), with variable thickness in a range of 50 to 100m, unconformably cover various older deposits from the Jurassic to the Pontian ages [29,37]. Sands and gravels alternate, irregularly forming interlayers and lens-shaped bodies. Clays are deposited as centimeter-to decimeter-thick interlayers and lens-shaped bodies within the sands. Quaternary sediments occupy almost 50% of Sisak-Moslavina County. At the beginning of the Pleistocene, tectonic processes reactivated existing faults and river valleys were formed. That intensified erosion, and sedimentation started, mostly in the Sava but also other river valleys.
In the Early Pleistocene, aeolian clayey-sandy silts (loess-l) were deposited at the edges of Sava valley, with a thickness up to 30 m.
During the Holocene, sedimentation was mostly related to river valleys and flatlands and is characterized by alluvial sands, gravels, silts, and clays (a, ap). The content and geometry of sediment bodies depend on facies: terrace, still water, flooded areas, and recent flows. Oxbow lake sediments form characteristically arched and elongated depressions, filled with mud, silt, and sand deposited during floods (jb). In the morphologically lowest parts of the Sava and Kupa River valleys, swamp areas are found. They mostly deposit dark clays and clayey silts (b) rich in organic matter [37], one to four m thick. Locally, the products of weathering, leaching, and torrential flows on slopes form poorly sorted and chaotic sediments on top of bedrock. Diluvial and proluvial sediments (dpr) are predominantly represented by sands and silts, with variable thickness that locally and rarely exceeds ten m.

Datasets
One of the principal assumptions employed in statistically based susceptibility mapping is that future landslides will more likely occur under similar geological and geomorphological conditions which led to their activation in the past [39]. Thus, landslide occurrence can be determined by the evaluation of inventory against different conditioning factors to define their influence on landslide events [40,41]. It is crucial to have an accurate landslide inventory for the area studied, based on which suitable conditions for their activation can be defined. Spatial association analysis in this study is performed using LiDAR-based inventories and four environmental precondition factors. The landslide causative factors (triggers) that are more associated to hazard assessment were excluded from the analysis, as was the influence of anthropogenic activities.

Landslide Inventories
Within this paper, several landslide inventories are used. The inventories differ in data-collection methodology, the scale of the presentation, and the period in which data were collected ( Table 1).
The  The identification and delineation of landslide polygons is based on visual surface interpretation techniques (e.g., [22,[44][45][46]). Thereby, HR DTM derivates (slope, contour lines, hillshade, and curvature) and orthophotos with a 10 cm resolution are analyzed using ArcGIS 10.2.1. During desktop landslide mapping, they were classified according to the landslide identification confidence (LIC), using scores from one to ten. LIC is presented as the degree of visibility; freshness; and the recognition of typical landslide features such as (i) main scarp, (ii) lateral flanks, (iii) toe, and (iv) internal morphology (e.g., minor scarps, internal cracks, hummocky, or rolling topography). All these features are soothed over time and are less recognizable in the field and on the digital elevation models (DEM) [47,48]. The LIC, in that sense, can serve as an indicator of relative landslide age. Thus, all landslides with LICs of one or two are considered inactive and old and are excluded from the analysis.
In total, 1238 landslide polygons with LIC higher than two are considered in the calculation, corresponding to an area of 2.25 km 2 (1.7% of the LiDAR polygons). The landslide area ranges from 6.88 m 2 to 104,440.27 m 2 , with an average value of 1832.84 m 2 ( Table 2). It can be seen that the majority of polygons have low LIC values (3)(4)(5)(6), indicating that active and younger landslides are still in the minority among the research polygons. Furthermore, a field prospection was performed, during which the depth of regolith and granulometric composition was estimated. A regolith is, according to the Dictionary of Earth Sciences [49], a layer of unconsolidated weathered material, including rock fragments, mineral grains, and all other superficial deposits, that rests on unaltered, solid bedrock. Field prospection discovered that mostly shallow and smaller landslides prevail in these terrains.
A total of 549 landslide polygons, with an LIC greater than five, were used in statistical analysis to estimate susceptibility. The remaining 689 polygons with an LIC of three or four were used for the verification of the derived LSM.
In addition to LiDAR-based inventories, which are considered to be recent, historical inventories are digitized from available geological maps (Table 1). These landslides are presented as points, and their spatial distribution is given in Figure 2. In total, 152 landslides are mapped, 53 digitized from EGM, 34 digitized from BGMs, and 65 attributed to working field maps. The last presented inventory is taken from the publicly available web portal "Report a landslide". This portal was established during safEarth project implementation in 2018 and is still active. Up to now, 57 landslides within Sisak-Moslavina County have been reported, mainly by the directorate of civil protection and local communities.

Landslide Predictive Factors
In the present study, four frequently used landslide predisposing factors were analyzed to map the susceptible areas: (i) engineering geological (EG) units, (ii) relief, (iii) roughness, and (iv) distance to streams, for all of which theoretical reasons exist that justify their use in susceptibility assessment [12]. The raster maps of all predictive factors were prepared in a GIS environment with a 20 m spatial resolution, based on the 1:100,000 scale geological map of Sisak-Moslavina County [27] for EG units and Croatian DEM with 20 m resolution for all other factor maps.

EG Units
Due to the lack of 1:100,000-scale EG maps for Croatian territory, geological units established through the available geological map (presented in Section 2) were used to define the EG units. The reclassification is based on pure EG logic and experience gained during the field prospection. Therefore, the simplified EG rock classification is used (Table 3), according to which two main groups are determined: soil and rock. Regarding landslide occurrence, the main difference between these groups is in the position of the failure surface against the bedrock material. Landslides activated in the soil-group have deeper failure surfaces, which generally affect the parent soil material. Their activation is strongly influenced by the characteristics of the soil material. The homogeneity of soil material also has an impact on landslide characteristics and their occurrence. Thus, from the aspect of EG, it is very important to know if a particular geological unit consists of a homogeneous coherent material (e.g., mud and clay), a homogeneous noncoherent material (e.g., gravel and sand), or a heterogeneous mixture of those two. Accordingly, the soil-group is divided into three corresponding subgroups (Table 3).
Landslides activated in rock-groups have a shallow failure surface, which is usually positioned at the contact of the bedrock and the superficial material, the regolith. The failure surface, in weak rocks, can be situated in the bedrock, but these cases are rare and isolated in the study area. Landslide properties are primarily controlled by the characteristics of superficial material, which are not described within the BGM's guides. According to engineering experience, regolith characteristics (thickness, granulometric, and mineral composition) are primarily governed by the characteristics and/or type of the bedrock. Additionally, landslide activation can be strongly influenced by the impermeability of the bedrock, resulting in high pore pressure at the contact of the parent material and regolith.
The rock-group is divided into two subgroups, weak and durable rocks. Weak rocks, such as marls, siltstone, and shales, are very susceptible to weathering processes, the result of which is a thick (more than two m) regolith, characterized mainly by a clayey-silt composition. Durable rocks are characterized by a thinner regolith with a much higher portion of a coarse-grained fraction. This subgroup includes classes of clastic sedimentary, carbonate, igneous, and metamorphic rocks.
Finally, for the study area, 29 geological units were reclassified into seven EG units (Table 4, Figure 3a).

Relief
The relief represents a local altitude difference between the highest and lowest height within the unit area and gives an indication of the potential energy for soil erosion and landsliding [50]. In areas with higher relief, higher runoff and lower infiltration are expected [51]. The relief map is computed from the DEM using focal statistics, calculating the elevation range for the circle area with a radius of five cells around each cell.

Roughness
Parameters correlated to surface roughness are quite commonly used to describe landslide activity in quantitative geomorphology [52]. Roughness is calculated from the slope map, which is derived directly from the DEM. The standard deviation of the slope is applied as one of the methods that have been used for studies at different scales [53]. Surface roughness was calculated using focal statistics. The radius taken around each cell to calculate the standard deviation was set to five cells. A 100 m window was chosen according to the DEM resolution and landslide dimensions. Roughness is classified based on the natural breaks (Jenks) algorithm into six classes (Figure 3c): 1 (0-1 • ), 2 (2-3 • ), 3 (4-6 • ), 4 (7-9 • ), 5 (10-12 • ), and 6 (13-21 • ).

Distance to Streams
Distance to streams is also one of the common conditioning factors used in landslidesusceptibility assessment [12]. Streams may negatively affect slope stability by eroding the toe of the slope or saturating the slope, which consequently leads to landsliding [54]. The stream network was automatically extracted from the DEM using Arc Hydro Tools v2.0. The derived drainage network was visually crosschecked on a Croatian topographic map at a scale of 1:25,000. Additionally, some minor rivers were manually extracted from the topographic map and major watercourses were removed, in order to retain lowerorder streams for the calculation; these are designated as headwater streams and medium streams [55]. The six buffer zones around these streams are defined for factor classification (Figure 3d), and the remaining areas are labeled as "no data", to which the smallest final weight is associated. Thus, seven classes are defined as follows: 1 (no data), 2 (0-100 m), 3 (101-200 m), 4 (201-300 m), 5 (301-400 m), 6 (401-500 m), and 7 (> 500 m).

Analytical Methods
To derivate the LSM for Sisak-Moslavina County, the bivariate statistical FR method and its integration in the AHP is used, following several steps: (i) weighting the classes within each predictive factor; (ii) weighting each predictive factor; (iii) calculating the landslide susceptibility index (LSI), (iv) reclassifying the LSI values, and (v) verifying the LSM.
FR is based on the relationship between the spatial distribution of landslides and various classes of predictive factors [56]. The FR of each factor class is calculated separately for each predictive factor using Equation (1).
where A Li is the area of landslides in the i-th factor class, A i is the area of the i-th factor class, A LT is the total area of landslides, and A T is the total area studied.
To enable the comparability between the FR of different predictive factors, FR values are normalized (FRn) into a range from zero to one using unity-based normalization, according to Equation (2).
where FRn represents normalized FR value, FR is the calculated value before normalization, and FR min and FR max are the minimum and maximum FR values of classes within the corresponding predictive factor. Besides the different impacts of factor classes, the influence of individual factors on landslide occurrence is also defined. To quantify the impact of each factor, the AHP method was used [57], of which the basic principle is to construct a comparison matrix expressing the relative importance of a set of attributes. Using the procedure described by [58] and [41], first, the predictor rate (PR) was determined for each predictive factor based on the degree of spatial association with landslides expressed through FR values (Equation (3)), instead of subjectively expressed by expert knowledge.
where the numerator represents the absolute difference between the maximum and minimum FR within each factor and the denominator represents the lowest absolute difference of all factors. Second, a comparison matrix is obtained, showing the relative importance of each predictive factor, and used for the estimation of eigenvectors by normalizing the pairwise results for each column. Third, fractional weights (FW) are calculated by averaging the eigenvectors, performed by dividing the eigenvectors' row sum by the number of input predictive factors. Last, the FWs were converted into integer weights (IW) by dividing each of the FWs by the smallest FW of all factors and multiplying with the number 10. The final factor class weights (FW) combine the factor classes' weights and factors' weights and are calculated by multiplying the normalized FR values (FRn) and IW values. Using the ArcGIS tool Raster calculator, the final weights of all overlapped factor maps are summed up to derive the LSI.
In order to present the studied area with zones of similar susceptibility to landslide occurrence, the LSI is reclassified into four classes based on its relationship with the cumulative number of landslides digitized from geological maps.
The verification of the LSM was performed using three different landslide datasets: (i) LiDAR-based landslide polygons with LICs of three or four, (ii) landslides digitized from different geological maps, and (iii) landslide locations collected via the web portal "Report a landslide". Landslide-susceptibility zones were verified with relative landslide density (RLD), calculated as a ratio between the landslide area or landslide number within each susceptibility zone and the area of corresponding susceptibility zone, all expressed in percentages.

The Weighting of Predictive Factors
Using the aforementioned approach, the first results for the LSM at a regional scale for Sisak-Moslavina County are presented. The results of the FR method, used to rank quantitatively the classes of each predictive factor to landslide susceptibility, are presented in Table 5. Factor classes that create the most favorable preconditions for landslide processes are identified as those with FR values greater than 1. Considering the EG units, landslides are mainly concentrated within EG unit 1, where 42% of landslide area is present. It can be considered to be the most critical landslide-prone unit, based on having the highest FR value of 1.75. A significant portion of landslide area is placed within EG units 4 and 5, comprising 16 and 13% of the landslide area, respectively. For those two units, FR values are greater than one. Although quite a bit of the landslide area, 18%, is associated with EG unit 3, due to the relatively large area of this unit, the FR is below one. Within relief, classes 3 and 4 have the highest FR values, of 1.79 and 1.35, respectively. Within these units, 79% of the landslide area is present. The FR value for class 5 is close to 1 (0.95), and this class contains 14% of the landslide area. Observing the FR results regarding roughness, classes 3 and 4 stand out with the highest FR values, of 1.39% and 1.17. 67% of the landslide area is placed within those classes. Analyzing the distance to streams, it is shown that FR values enlarge as distance increases. Thus, the highest FR, 1.54, is related to class 7, within which 60% of the landslide area is present.
Integrating the FR weights in the AHP, the prediction rate (PR) and a pairwise comparison matrix of factors' relative importance to landslide susceptibility are determined ( Table 6). Fractional weights (FW) and integer weights (IW), calculated over the estimated eigenvectors, are presented in Table 7. IWs are also expressed in percentages to get a better overview of factors' contribution to landslide susceptibility. Therefore, the influential role of EG units, relief, roughness, and distance to streams in landslide-susceptibility is expressed as 27%, 31%, 21%, and 20%, respectively (Table 7).

Landslide Susceptibility Map and Verification
Combining the normalized FR and IW, the final weights of each factor class are determined (Table 5) and used to calculate the LSI. The calculated LSI values vary within a range of 0 to 47 (Figures 4 and 5a), wherein higher values indicate greater susceptibility to landslide occurrence and vice versa. The derived map in its original form presents the LSI distribution over the study area (Figure 5a) with a cell size of 20 × 20 m. It is the resolution of all factor maps overlapped. The generated LSM is in the final step adjusted for regional scale (Figure 5b) by: (i) classifying the LSI, (ii) averaging the LSI, and (iii) resizing the cell size to 50 × 50 m.
To describe the different levels of landslide susceptibility, LSI values were classified into four susceptibility classes: low, moderate, high, and very high. Within this study, classification is based on the relation between LSI values and the percentage of the cumulative number of landslides, considering 152 landslide locations digitized from geological maps. Low and moderate classes are defined as areas where a maximum of 20% of cumulative landslides are present, and high and very high classes are defined as areas where at least 20% of cumulative landslides are present. The other two threshold values are set at 5 and 55% of the cumulative number of landslides. Based on those criteria, threshold LSI values were interpreted as 20, 26, and 32 ( Figure 4).
The LSM is adjusted to a suitable scale of 1:100,000, at which the raster can be reliably used, by averaging the LSI value within each 100 m circle neighborhood with the mean statistic and setting the cell size to 50 × 50 m, according to the mathematical relationship between the map scale and the image resolution [59].
The representation of each landslide-susceptibility zone in Sisak-Moslavina County is shown in Table 8. According to the results, very high and high landslide-susceptibility zones cover 36% of the county area, while moderate and low susceptibility zones cover 64%. * A Li is the area of landslides in the i-th class of the factor; the total landslide area is 1.08 km 2 . ** A i is the area of the i-th factor class; the total study area is 132.43 km 2 .    1st set-landslides digitized from maps; 2nd set-LiDAR-based landslide polygons (3,4); 3rd set-web portal "Report a landslide"; 4th set-LiDAR-based landslide polygons (>5).
LSM is verified using three landslide datasets, as presented in Table 8, in order to ascertain their reliability for the intended purpose. The overall verification results show that, for all verification sets, at least 70% of the landslide locations (points or areas) match the most susceptible zones, i.e., areas of very-high and high susceptibility. Furthermore, at most, 30% of the landslide locations (points or areas) match least susceptible zones, i.e., areas of moderate and low susceptibility. Among verification datasets, the results obtained for the first and second set are very similar and demonstrate the stronger spatial relationship between landslide locations and zones of higher susceptibility when compared to the third set.
The RLD indicates the general quality of the LSM [18]. Considering all three verification datasets, it is evident that the RLD slightly increases, moving from the lowest to the highest susceptibility class ( Table 8). The RLD within the low susceptibility zone is below 0.1 for the first and second verification set, and below 0.32 for the third set. Further, the RLD increases up to a maximum of 2.33 within a very high susceptibility class for the second verification set. Overall, it can be concluded that the RLD within both low and moderate susceptibility zones is below 0.5, and that within the high and very-high susceptibility zones, it is around two. For comparison, the fourth set is added, which relates to the LiDAR-based landslide polygons with an LIC greater than five. There, the RLD ranges from 0.02 (low susceptibility zone) to 3.26 (very high susceptibility zone).

Limitations
LSM techniques have evolved greatly over the last few decades and numerous techniques have been proposed and applied by many researchers. Although the performance of some methods has been proven to be better than others, there is no single method superior in all environments, and each technique has its unique advantages and limitations [12,60]. Therefore, optimal susceptibility zonation by combining several methods is recommended [12,61]. In this study, the integration of FR into AHP is used, to avoid subjectivity in defining the relative importance of predictive factors.
A well-known constraint during the derivation of LSM for larger regions is the relative scarcity of available data, relating to both the landslide inventory and predictive factor maps. To overcome the limitation of the nonexistence of a complete inventory for the entire county, limited LiDAR polygons are used to produce landslide inventories. The polygons are carefully selected to incorporate not only known landslide-susceptible geological units, but also other units and environments for which the landsliding potential is not known. LiDAR polygons represent the testing area where the spatial correlation between landslide occurrence and predictive factors was estimated. Correlation results are extrapolated to the entire county. The biggest limitation during the LSM model development was related to the size of the area covered with LiDAR polygons, which is, in the present study, limited to approximately three percent of the county area. The results of weighting factors used to define the LSI are strongly influenced by the representation of all predictive factor classes within LiDAR polygons. It is concluded that, for such large areas and regions, the positioning of LiDAR polygons prerequires a basic statistical analysis of predictive factors.
If LiDAR polygons do not cover all the classes of predictive factors proportionally, further ALS should be directed and/or redistributed to capture the missing data.
The proportions of each class of predictive factors within LiDAR polygons are presented in Figure 6. It is shown that the area of EG units within LiDAR polygons is not equally distributed (Figure 6a). Still, EG unit 1, which is, from engineering geological experience, known as susceptible to landsliding, has the largest FR, although only 1.18% of this unit is covered with LiDAR polygons. The classes for other predictive factors are evenly covered with LiDAR polygons, except for the first classes (Figure 6b-d). Those classes are related to lowland areas (Sava and Kupa valleys) and gentle slopes within the county where landslides are not so frequent. Thus, those areas were not considered to be priorities and were not in focus during the first ALS scanning of this region. Our future investigations will be focused on covering at least five percent (red line on the graphs in Figure 6, right vertical axis) of each class with HR LiDAR data. The reliability of the results is also influenced by the accuracy of the data. LiDARbased inventories are precise inventories, considering the HR DTMs used for landslide delineation. Still, visual DTM interpretation involves a certain degree of subjectivity. To maintain the same criteria during the landslide delineation by three different experts and enable the comparability of inventories, the LIC for each landslide is defined. The geological setting for the study area has been defined by available geological maps in a scale of 1:100,000, which is appropriate for regional-scale LSM. Other factor maps were calculated from available DEMs with 20 m resolution. Therefore, the pixel size of the DEM controls the resolution of the susceptibility estimation, but a resolution of 20 m is considered to be a suitable resolution for regional-level analysis.
Considering LSM verification, several limitations could also be pointed out. First, landslide locations used for verification, which are digitized from available geological maps, are presented as points. Their spatial accuracy, with respect to the scale of the input maps, can be questionable and can strongly affect verification results, especially when the landslide point is close to the boundary of two different susceptibility zones. Besides limited historical landslide locations, which is almost impossible to present on small-scale maps, weaknesses can be also related to the age of the maps, mapping methodology, georeferencing procedure, and a lack of data in some parts of the study area. Second, not all data collected from the web portal "Report a landslide" can be considered completely accurate until all locations are checked by experts in the field. Nevertheless, most of this data can be considered reliable, as a large number of landslides are reported from the Directorate of Civil Protection, which is involved in most crisis situations related to hazardous events. Some of the landslides related to the latter verification set are located within the green zones (low susceptibility) of LSM, which are associated with the banks of the Kupa River. Considering the resolution of input DEMs, it was difficult to highlight the local geomorphological properties within these areas. Thus, these local instabilities have an impact on the verification results. The main constraint arises from limited area and inventory coverage, which is already explained where LiDAR polygons are discussed.
Keeping in mind the limitations caused by showing landslides as points, a part of LiDAR-based inventories with LICs of three or four, as a verification dataset, is also used.
It is also important to point out that the presented LSM is based on landslides, which are mainly initiated by long-lasting/intensive precipitation. Seismic activity and climate change extremes that could initiate landslides in "unusual" locations in the future are not considered here.

Purpose and Method of Use
The LSM does not predict exactly where landslides will occur. Therefore, it is important to emphasize that red zones on the map (zones of very high susceptibility) do not indicate explicitly landslide locations. The susceptibility zones represent areas of different landslide spatial probabilities, i.e., the likelihood of their occurrence due to local terrain conditions [7], neglecting the factor of time.
For the LSM interpretation, it is very important to be aware of the purpose for which the map is made, which greatly influences the corresponding map scale and the methodology used. Besides, the final classification of LSI into susceptibility zones directly influences the final appearance of the map, which is essential for its practical application. For proper LSM use, it is important to prepare a guide for potential end-users, which should contain the description of the methodology used, instructions for map interpretation, and a review of its advantages and limitations. Besides, the LSM presents the summary of the results based on currently available data and the current state of the terrain, which is subject to change over time. Thus, it is recommended that the LSM be updated every few years, incorporating new data, knowledge, and methodologies.
Considering all the limitations presented, the use of the presented 1:100,000 scale LSM is limited to regional-scale studies. The derived LSM could be used for: (i) development of county spatial plans, (ii) identification of landslide problem areas at the regional level, (iii) selection of locations where it is necessary to conduct more detailed EG surveys and create more detailed landslide-susceptibility maps on a local scale (e.g., 1:25,000), (iv) elaboration of disaster risk-management strategy, (v) planning of regional development projects, (vi) setting engineering constraints during the development of large projects, and (vii) informing the local communities and general public.
From all of the above, it can be concluded that the main purpose of the regional-scale LSM is to point out the most susceptible zones in a wider area where detailed research should be conducted and large-scale LSMs produced. In this context, regional LSMs cannot replace detailed engineering geological surveys and should not be used for design purposes. In such a way, the LSM derived for Sisak-Moslavina County can be used to rank the municipalities and cities based on priority for further investigation (Figure 7).
The research results indicate municipalities and cities where >60% of the area is determined to be highly susceptible to landsliding. For those, such as Gvozd, Hrvatska Kostajnica, Glina, Dvor, Topusko, and Majur, detailed research should have priority. They should be followed by Donji Kukuruzari and Petrinja, where >40% of the area is associated with high susceptibility zones. Among others, within Novska and Kutina, more than ten percent of the area belongs to a very-high susceptibility zone. Jasenovac and Martinska Ves are completely within the zone of least susceptibility, implying the possibility of rare landslide occurrences related to extreme natural events or induced by human interventions. Moreover, within municipalities where low susceptibility zones prevail, such as Sisak and Popovača, detailed investigation should be focused only on local zones determined to be vulnerable.

Conclusions
In this paper, a regional-scale 1:100,000 LSM is presented for the 4,466 km 2 area of Sisak-Moslavina County in Croatia. To derive the LSM, a statistically based FR method is used for weighting the classes of each predictive factor, and the results were integrated using the AHP method to classify the importance of each factor in landslide-susceptibility mapping. The results of this approach are presented by the spatial distribution of LSI, whose values range from 0 to 47. After LSI reclassification, the terrain was divided into four landslide-susceptibility zones: low, moderate, high, and very high, covering 50%, 14%, 17%, and 19% of the area, respectively.
The most important concern that defined the goal of this paper was whether the use of a limited amount of LiDAR data over such a large area would yield a sufficiently good and usable landslide-susceptibility map. Based on the verification results, it can be concluded that the approach presented in this paper is encouraging for primary regional-level studies, justifying the cost-benefit ratio.
However, it should be emphasized that such a limited amount of key and expensive data in the process of LSM production, i.e., LiDAR data, would not be convenient enough if the analysis of the EG units and environmental parameters were not performed before the definition of LiDAR polygon positioning. Such an approach resulted in a statistically balanced landslide inventory sample and is certainly recommended for similar research in the future.
The results achieved in the form of an LSM at a scale of 1:100,000 can be used for different purposes on a regional and local level, with awareness of its advantages and limitations emphasized in the paper. The regional-scale LSM is a good reference point to detect areas for which detailed research is necessary. The more detailed the research is, the more expansive investigation will be. Thus, for the most financially demanding phase, the main objective is to reduce the investigation area and to direct resources and detailed research to the most vulnerable sites in the region.
As part of a comprehensive landslide risk-management system, the development of LSMs is the basis for further hazard and risk assessment. Therefore, the systematic monitoring of landslide occurrences and the damage they cause is crucial and is possible only with the cooperation of engineering geologists, local communities, and civil protection.
Author Contributions: All authors conceptualized the main research purpose of the paper. I.B. and M.F. performed the analysis, prepared the data, tables, and figures, and wrote the majority of the first draft of the paper. V.G. and D.P. contributed with paper review and editing. All authors have read and agreed to the published version of the manuscript.