Preliminary analysis of a LiDAR-based landslide inventory in the area of Samobor, Croatia

The paper presents an analysis of the LiDAR-based landslide inventory for the area near Samobor, in northwestern Croatia with two main objectives: i) to define the geological units (obtained from Basic Geological Map of SFRY) most susceptible to landslides, and ii) to analyse the limitations of the Basic Geological Map and its applicability in landslide susceptibility map design. Within the study area of 63.8 km2, 874 landslide polygons were manually outlined, covering an area of 2.15 km2. The landslide outline confidence level, landslide index and the relief energy map were used to analyse the landslide susceptibility of a particular geological unit. By that, units in the same state of stress, i.e., in the same relief energy group were compared. This preliminary analysis has shown that the geological units Pl,Q, M3, and 1M3 are the most susceptible to landslides and that older geological units, Pc and K1,2, are also prone to landslides. Still, landslides within those older units can be considered as old and inactive. As for the limitations of the Basic Geological Map of SFRY, three things emerged, namely scale, the geological unit defining approach, and the neglect of regolith. Despite the limitations presented, the usability of the Basic Geological Map of SFRY in the development of small-scale landslide susceptibility maps is emphasized. However, instructions that should attribute engineering geological features to the geological units outlined in the Basic Geological Map should be prepared in the near future. landscaping, and construction codes; use of physical measures (drainage, slopegeometry modification, and structures) to pre­ vent or control landslides; and development of warning systems. Preventative directions focused on the first step aimed to de­ termine the spatial distribution of landslide-prone areas, assess the probability of landslide occurrence, and estimate the threat of elements within the area affected by landslides. This can be achieved by making several types of maps for the area of inter­ est, covering susceptibility, hazard, and risk (CHACÓN et al., 2006; FELL et al., 2008). However, in order to be able to produce such maps by accepting the general principle that “the past and the present are key to the future” (VARNES & IAEG, 1984), it is necessary to create detailed landslide inventories for a particu­ lar area (GUZZETTI et al., 2012). From that perspective, airborne light detection and ranging (LiDAR) data improved the creation of landslide inventories dur­ ing the last decade, especially when considering their precision and detail. In this study, a detailed LiDAR-based landslide inven­ tory is prepared for the Samobor area located in the northwestern part of Croatia. This area is characterized by high geological di­ versity which within the study area enables as many as 18 geo­ logical units (GUs) to be distinguished within Basic Geological map "of SFRY 1:100.000, Zagreb sheet (ŠIKIĆ et al., 1977), herein after referred to as BGM" (BGM). Based on the developed landslide inventory and available geological data, an extensive spatial analysis is performed, with the main purpose of gathering knowledge of landslide occurrence within each GU. It provides the basis for accomplishing the main objectives of this paper: i) to define the GUs within the study area most susceptible to landslides, and ii) to analyse the limitations of the available geological maps and their applicability in the landslide susceptibility map (LSM) design. Article history: Received October 01, 2021 Revised manuscript accepted November 01, 2021 Available online February 28, 2022


INTRODUCTION
Landslides are the most widespread geological event affecting 4.8 million people and caused more than 18.000 deaths worldwide in the period from 1998 to 2017 (CRED, 2018). In the same period, landslides, along with volcanoes and other mass movements caused damage of US $ 61 billion. Given that world population growth at the moment shows a linear trend (ROSER et al., 2013), an increase of negative effects related to landslides is to be ex pected and these expectations according to SCHUSTER (1996) are based on the following reasons: increased urbanization and development in landslide-prone areas, continued deforestation of landslide-prone areas, and increased regional precipitation caused by changing climatic patterns.
The systematic monitoring of landslide occurrences and the damage they cause is not established in Croatia on a national level. However, by inspecting the Register of Damages from Natural Disasters at the Ministry of Finance, it is possible to extract some data for the period from 2015 to 2019. According to the source, the total damage from landslides in that period was US $ 22 million, counting only reported landslides and their casualties in those municipalities or cities that have declared a natural di saster in the specified period. This data indicates that landslides in Croatia clearly present a harmful threat to humans, their properties, and the environment. Still, raising public awareness of landslides as geohazard events in Croatia remains a challenging task.
As the Earth's surface is the only human habitat, it is logical to preserve it, which can be achieved especially through preven tive activities, primarily through effective planning and manage ment (DAI et al., 2002). This approach includes restriction of de velopment in the landslide-prone area; use of excavation, grading,

Preliminary analysis of a LiDAR-based landslide inventory in the area of Samobor, Croatia
landscaping, and construction codes; use of physical measures (drainage, slope-geometry modification, and structures) to pre vent or control landslides; and development of warning systems.
Preventative directions focused on the first step aimed to de termine the spatial distribution of landslide-prone areas, assess the probability of landslide occurrence, and estimate the threat of elements within the area affected by landslides. This can be achieved by making several types of maps for the area of inter est, covering susceptibility, hazard, and risk (CHACÓN et al., 2006;FELL et al., 2008). However, in order to be able to produce such maps by accepting the general principle that "the past and the present are key to the future" (VARNES & IAEG, 1984), it is necessary to create detailed landslide inventories for a particu lar area (GUZZETTI et al., 2012).
From that perspective, airborne light detection and ranging (LiDAR) data improved the creation of landslide inventories dur ing the last decade, especially when considering their precision and detail. In this study, a detailed LiDAR-based landslide inven tory is prepared for the Samobor area located in the northwestern part of Croatia. This area is characterized by high geological di versity which within the study area enables as many as 18 geo logical units (GUs) to be distinguished within Basic Geological map "of SFRY 1:100.000, Zagreb sheet (ŠIKIĆ et al., 1977), herein after referred to as BGM" (BGM).
Based on the developed landslide inventory and available geological data, an extensive spatial analysis is performed, with the main purpose of gathering knowledge of landslide occurrence within each GU. It provides the basis for accomplishing the main objectives of this paper: i) to define the GUs within the study area most susceptible to landslides, and ii) to analyse the limitations of the available geological maps and their applicability in the landslide susceptibility map (LSM) design.

STUDY AREA
The boundaries of the study area are defined within the safEarth project co-financed by the Interreg IPA CBC Programme Croatia -Bosnia and Herzegovina -Montenegro. One of the main goals of this project was to create LSMs at different scales and in dif ferent areas. In order to achieve this goal, an airborne LiDAR survey was performed on six pilot polygons, for which detailed landslide inventories were created. The selection of pilot polygons was based on the geological analysis, attempting to cover as wide a variability of geological environments as possible.
The special geological features of the Samobor polygon, which is presented within this paper, set it apart from other safEarth polygons exactly due to the large number of GUs and the fact that a part of the polygon is composed of carbonate rocks of Triassic age.
The study area covers approximately 63.8 square kilometres and is located 19 kilometres southwest of the centre of Zagreb. More precisely, it is situated between Sveta Nedelja, Samobor, and Rude in the north, and Drežnik Podokićki and Rakov potok in the south (Fig. 1). Most of the area belongs to the Sava river basin (38.38 km 2 ), and a smaller part to the Kupa river basin (25.46 km 2 ).

Geological settings
Landslides as exogenetic processes are the result of the interaction between several Earth systems; the lithosphere, atmosphere, and hydrosphere. Professionals and scientists often analyse different natural and human factors that can trigger the landslides (e.g., LI et al., 2012;ZEZERE et al., 1999;LUKIĆ et al., 2018). However, to bring the rock material into a state of motion, the rock material first must be in a state of labile equilibrium, i.e., a large percenta ge of the rock material strength must be "consumed".
This primarily depends on two factors, namely the strength of the rock material, and then the stress (i.e., potential difference) that the rock material suffers. As this paper for the first time com prehensively analyses the landslides within the study area, it is logical to present that area according to these two primary fac tors as follows: -rock material strength is described through geological characteristics; -potential difference is described through the endogenetic movements.

Geological units
The sheets of Basic Geological Map at scale 1:100.000 that cover the area of Croatia were prepared, which for the most part lasted through the second half of the 20 th century. It represents the most comprehensive overview of geology in the entire Croatian territory and as such represents an unavoidable basis in any type of research that requires knowledge of the geological conditions of a particular area.
According to BARNES & LISLE (2004), there are four main groups of geological maps, and the BGM falls into the category of Regional geological maps, and as such it is produced by col lecting field data together with the analysis of photogeological documentation. The geological settings of the study area ( Fig. 1) are extensively described within the BGM sheet of SFRY 1:100.000 (ŠIKIĆ et al., 1977) and an accompanying guide (ŠIKIĆ et al., 1979), according to which an overview of GUs pre sented within the Samobor area is given.
The study area is composed of Palaeozoic, Mesozoic-Palaeo gene, Neogene and Quaternary deposits.
The oldest formation of Permian age (P 2,3 ), is represented by brownish sandstones with interlayers of shale and siltstone with sporadic transitions to quartz conglomerate.
Sediments of the Lower Triassic (T 1 ) appear in the form of purple and reddish mica sandstone and siltstone and in a smaller amount of ooid limestone. In the upper part, ooid limestones with marl intercalations predominate.
The sediments of the Middle Triassic (T 2 ) are represented by the transition from diagenetic dolomitized stromatolite, dolo pelmicrite and fenestral dolopelmicrite. Also, sporadically, thicklayered and massive crystalline limestone and dark mudstones are observed.
The massive dolomite with weakly expressed stratification of the Upper Triassic (T 3 ) are characterized by frequent interca lation of dolopelmicrite, fenestral dolopelmicrite and characteristic thin laminated dolostromatolite.
The deposits of the Cretaceous period (K 1,2 ) are quite litho logically heterogeneous and are represented by mediumgrained greenish sandstone, with shale, radiolarite and chert alternation. Within the sedimentary rocks a larger mass of basic igneous rocks like diabase, spilitized diabase and spilite (ββ) are observed.
The Palaeocene (Pc) deposits lie transgressively on older Tri assic or Cretaceous deposits and they are mostly represented by polymictic conglomerate with imbrication of marl, siltite and sand stone, as well as bioclastic limestone and biolithite are observed.
The sediments of Ottnangian age ( 1 M 2 1 ) are lithologically heterogeneous with conglomerate, breccia, and gravel dominate while siltstones are subordinate.
The Badenian( 2 M 2 2 ) sediments are characterised by biocal carenite, conglomerate and biolithite limestone developed in the ridge and coastal facies, while silty marl and clayey limestone developed in deeper marine sedimentary environments.
Sarmatian deposits ( 1 M 3 1 ) disconformably overlie Badenian deposits in the form of biocalcarenite and sandstone deposits in the marine environment with reduced salinity, as well as charac teristic deep water laminated marl and silt with sporadic inter calations of sand.
The shallow lake deposits of the Lower Pannonian (M 3 1,2 ) better known as the "Croatica layers" are characterized by stratified clayey limestones of irregular strata that may be partly inter bedded with marl. In the upper Pannonian, marls are deposited in a deeper lacustrine environment with interlayers of sand and sandstone.
Concordant on the marl, are pro-delta and delta deposits of the Upper Pontian (Pl 1,2 ) in the form of sand, calcite-rich siltstone and subordinate marl.
Sediments of Pliocene-Pleistocene age (Pl, Q) were formed in the transitional, fluvial to lacustrine environment and are rep resented by interbeds of gravel with sand lenses, sand with gravel lenses, clayey to sandy silt and clay with sporadically strong li monitization.
Loess (l) and loess-like sediment (lb) are mostly silty clay and clayey silt, while proluvial deposits (pr) are in the form of coarse-grained, slightly rounded gravel mixed with sand and clay. Gravel, sand and clay are deposits of alluvium (a), and deposits of the second Sava terrace (a 2 ), gravel and sand, were formed by erosion and transportation of previously alluvial sediments.

Endogenetic movements within the study area
The ascent of internal energy originating in the Earth's core drives a complicated set of geological processes (HUGGETT, 2017). This energy induces endogenetic forces that continuously elevate or build up parts of the earth's surface and hence the exo genetic processes fail to evenout the relief variations of the sur face (UPPAL, 2006). In that way, they give rise to the fundamen tal relief forms and condition the manifestation of exogenetic processes (ERSHOV et al., 1988).
From the aforementioned fundamentals, it is clear that before any pioneering research of the exogenetic processes, it is neces sary to determine whether the studied area is in a state of uplift or subsidence. Thus, if the studied area is characterized by a pro nounced uplift concerning the surrounding areas, the resulting height difference can be considered as the accumulated energy of that area. Such accumulated energy will cause a whole range of exogenetic processes, e.g., erosion (initiation of erosive surface flows), gullying (a set of exogenetic processes that take place along the erosive concentrated surface water flow), landslides, and rockfalls.
Since the intensity of exogenetic processes depends, among other things, on the rate of uplift, it is necessary to quantify that rate. One of the first maps dealing with the neotectonic activity in Croatia is the neotectonic map at a scale of 1:1,500.000 (Fig. 2) produced by the Federal geological institute from Belgrade (ĆIRIĆ, 1967).
The map shows that the studied area is characterized by ne otectonic uplift, while a good part of the southwestern part is characterized by neotectonic subsidence. Thereby, the scale of the presented map should be kept in mind.
While analysing the endogenetic movements in the study area, a proposal for future research has arisen. The correlation between neotectonic activity and the intensity of exogenetic pro cesses would present a valuable study, the results of which could be used as one of the crucial data sets in landslide susceptibility assessment on a national level.

Remote sensing data sets
As predicted by JABOYEDOFF et al. (2012), the LiDAR sensor has become a standard tool for landslide analysis in this modern time. In the framework of the safEarth project, an airborne sur vey was performed by Flycom Technologies d.o.o. during the spring of 2018. A LiDAR survey, with requested point density of 20 points per m 2 and an accuracy of ±10 cm in each direction for each collected point, enabled the derivation of a high-resolution DTM with a cell size of 0.5 m. A topographic survey resulted with digital orthophotos with a 10 cm resolution.
The DTM was used to create several crucial derivatives (e.g., ARDIZZONE et al., 2007;VAN DEN EECKHAUT et al., 2007;JABOYEDOFF et al., 2012) for the production of landslide in ventories: -hillshade grid with a cell size of 0.5 m, derived by a sun azimuth angle of 315° and altitude angle of 45°; -slope grid with a cell size of 0.5 m; -contour lines with one m elevation equidistance. PIKE (1988) emphasises that a geometric signature is a subset of the geomorphic signatures, a broader and largely undeveloped concept that includes much more than topography. By studying the works that cite the mentioned author, it is possible to interpret the "signature" that exclusively defines the morphometry of a geomorphological phenomenon (GUZZETTI et al., 2012), e.g., a landslide. Even though technically this possible interpretation is not entirely correct, one can conclude that a distinct signature of landslide morphometry should for certain be noticed on high resolution elevation models and as such can be used for landslide investigations (JABOYEDOFF et al., 2012).

Landslide inventory creation
To create a landslide inventory, landslide polygons in this study were outlined manually within the GIS environment using visual analysis and interpretation of the topographic surface pro vided by high-resolution DTM derivatives. Primarily, the com bination of two data sets was used, slope grid as a top layer with 50 % transparency and the hillshade grid as a bottom layer. In some cases where landslide boundaries were not completely clear and sharp, the use of 3D visualisation software (Esri ArcSceene) was necessary to obtain the best possible results (Fig. 3).
Each polygon within the landslide inventory is accompanied by a landslide outline confidence level (LOCL) aimed to describe the level of visual clarity of the individual landslides. Specifically, the LOCL of each landslide was assessed according to how clearly the boundary of the landslide can be traced and features within the landslide body can be observed (e.g., cracks, surface depressions, hummocky topography). The LOCL is expressed by a grade from one to ten, where the highest grades are ascribed to those landslides with the most visible features. A very similar ap proach has been presented by BURNS & MADIN (2009). To conduct the analyses, the LOCL grades were clustered here into four groups (Fig. 3).

Field prospection
The field prospection was performed as a part of this work, pri marily to gain new knowledge regarding the engineering geo logical properties of materials present in the study area. A field prospection at 190 observation points gave a detailed overview of: -granulometric composition of Quaternary deposits which are outlined on the BGM; -lithological composition and durability of the bedrock ma terial; -types and characteristics of regolith, primarily thickness and granulometric composition.
It should be emphasized that the data from the field prospect ing are insufficient for extensive statistical analyses. However, the experience gained on that occasion greatly helped in the ar ticulation of certain, for this paper, important conclusions.

Landslide index
Landslide index is used to express the spatial relationship between delineated landslide polygons and GUs. According to BARTEL LETTI et al. (2017), it is calculated as the ratio between the area of landslides occurring in a given GU and the area of that particular GU. In the same way, the landslide index of the study area is cal culated as the ratio between the total landslide area and the entire study area. Since landslide indices are expressed in percentages, it was possible to compare the landslide index of each GU with the landslide index of the study area and draw the first conclusions about the susceptibility of each GU to landslide occurrence.

Landslide density
To get a better insight into the spatial distribution of landslides over the entire investigated area, a landslide density map was created. Landslide inventory is transformed to a grid format, where the value "one" is added to each cell that contains a land slide, and value "zero" to each cell without a landslide. Using the Focal statistics ArcGIS tool, a landslide density is calculated for each cell within a circle neighbourhood with a radius of 150 m. A landslide density map is used to delineate the zones for detailed analysis, in order to analyse the spatial distribution of landslides within specific GUs.

Relief energy
A relief energy map is derived from the high-resolution DTM (resolution 0.5 cm). Using the Focal statistics ArcGIS Tool, the elevation range is calculated by subtracting the maximum and minimum elevation within the circle neighbourhood with a radius of 150 m. The relief energy map is used as a representation of the stress state within the study area.

Statistical sample analysis
Statistical sample analysis is performed to see if the area of the individual GU is statistically significant within the study area. To determine the statistical significance, a simple histogram chart is analysed, comparing the area of each GU with the average area of GUs, which is taken as a reference line (Fig. 4). The presented histogram shows that GUs a 2 , 1 M 2 1 , and ββ do not present a sta tistically significant sample, since the histogram bars of those units are well below the reference line. Accordingly, the area of the GUs pr, Pl 1 2 , and Pc have questionable statistical significance with histogram bars a little below the reference line. This should be kept in mind when discussing those GUs, hence the results reached could be misinterpreted. For the other 12 GUs, it can be pointed out that the area they cover determines a statistically sig nificant sample, i.e., the conclusions derived from the analysis performed are based on relevant and sufficient data.

Landslide inventory map
A total of 874 landslides were outlined within the investigated area. The delineated landslides cover an area of approximately 2.15 square kilometres, which is 3.4 % of the total study area. The average values of the landslide area, length, and width are 2.462 m 2 , 70 m, and 37 m, respectively.
The landslide inventory map shows that the vast majority of landslide polygons are concentrated in the central part of the study area (Fig. 5).

Landslide index
After the visual analysis of the landslide inventory, the landslide share for the individual GUs was calculated. Some basic statistics related to the GUs and accompanying the landslide inventory are summarised in Table 1. The landslide share analysis was con ducted for all 18 GUs.
In the analysis of landslide number and landslide index dis tribution within GUs (Fig. 6), all landslides (LOCL grades from 1 to 10) were taken into account.
Considering the landslide number (Fig. 6a), it is shown that the bars for the GUs Pl 1 2 , M 3 1,2 , 1 M 3 1 , 1 M 2 1 , Pc, K 1,2 and T 1 rise above the reference line which should classify them as suscepti ble or prone to landsliding. This analysis served as a starting point in the study aimed at identifying the "problematic" GUs, regard ing the statistical sample size (Fig. 4). The problems are related to the GUs Pl,Q and 1 M 2 1 , as there are 293 delineated landslide polygons within the unit Pl,Q and only 1 within the unite 1 M 2 1 (Table 1). From the above, it is obvious that the representation of an individual landslide as a point loses the data of its surface area, which equalizes the importance of small and large landslides, masking landslide susceptibility characteristics for a particular GU. However, GUs with large areas and a high number of land slides may have a landslide density per km 2 below the average of the whole study polygon (e.g., Pl,Q).
The analysis that involves the surface area of the landslides within the study polygon is presented through landslide indices (Fig. 6b). This analysis shows somewhat different results to the previous one, as the bar for the GU Pl,Q stands above the refer ence line, and for GU 1 M 2 1 it is well below the reference line. The reason for that could be large landslide polygons delineated within that GU and this assumption is confirmed in Table 1, where the GU Pl,Q is presented as the second one in the rank list of average landslide area for all GUs.  Other GUs are in a similar relationship to the reference line as in the previous analysis. Grey bars of the units Pl,Q, Pl 1 2 , M 3 1,2 , 1 M 3 1 , Pc, K 1,2 , and T 1 (Fig. 6b) are above the reference line which characterises the listed GUs as susceptible to landsliding. Such results did not entirely match the engineering impression that was formed during the field prospection within the study polygon and during the creation of the landslide inventory. This refers to the older GUs Pc, K 1,2 , and T 1 , and especially to unit Pc. Although landslides can be observed within the mentioned GUs, they share somewhat different characteristics than landslides within the units Pl,Q, Pl 1 2 , M 3 1,2 , and 1 M 3 1 .
Such observations encourage further landslide index analy sis in which landslides of certain LOCL were taken into the ac count. Landslides were grouped into four groups, based on the following LOCL grades: 1) from one to three (Fig. 7a), 2) four to five (Fig. 7b), 3) six to seven (Fig. 7c), and 4) eight to ten (Fig. 7d). Part of the LOCL grade can very likely be attributed to the age of a landslide (MCCALPIN, 1984;BELL et al., 2012) since the grade depends only on the clarity of landslide feature expression, which gradually vanishes on the terrain surface as the landslide ages.
In the group of younger GUs (Pl,Q, Pl 1 2 , M 3 1,2 , and 1 M 3 1 ), there is no obvious trend as the bars sporadically stand above or below the reference line in all the presented charts (Fig. 7). Still, all of these units have been represented with bars that are well above or very close below the reference line. According to the LOCL, grades attributed to the landslides, both, old and recent landslides are present within those units. That could imply that these young GUs are as prone to landslides nowadays as they were in recent geological history. On the other hand, the bars for older GUs (Pc, K 1,2 , and T 1 ) show some trend from the first (Fig.  7a) to the last chart (Fig. 7d). It can be observed that bars related to older GUs are mostly above the reference line in the charts in Figure 7a and 7b, and well below the reference line in Figure 7c and 7d. This means that landslides within the older GUs are mostly rated with lower LOCL grades which could led to the con clusion that landslides are characterized by "blurry" features and are probably older than the landslides within the younger GUs.
According to the accompanying lithological characteristics (a mixture of gravel, sand, silt, and clay), the alluvial unit should be classified as a unit very susceptible to landsliding. Still, be cause these are the youngest deposits (Quaternary age) the endo genetic movements have not had the opportunity to bring them into a state of labile equilibrium. Therefore, a high landslide share within this unit ( Fig. 6a and 6b; Fig. 7a) is unexpected and is a result of the scale of the geological map used, which is elaborated within section 5.1.1.

Landslide density map
The landslide density map (Fig. 8), as a derivative of the landslide inventory, highlights very similar zones as the landslide inven tory map (Fig. 5). Yet the landslide density map more clearly sin gles out subzones of similar characteristics in terms of landslide density, around which a more detailed discussion is presented in section 5.
The subzones on the landslide density map marked as poly gons of detailed analyses are separated in the area of two GUs, namely Pl,Q and 2 M 2 2 . The polygons one, two, and three are placed within GU Pl,Q and were outlined using a landslide den sity map and an energy relief map (presented within the next chapter). Although located in the same GU, these polygons indi Figure 7. Landslide index distribution for the geological units and landslide index for the study area, considering: a) landslides with LOCL graded between one and three; b) landslides with LOCL graded between four and five; c) landslides with LOCL graded between six and seven; d) landslides with LOCL graded between eight and ten.
cate very contrasting environments in terms of landslide density and relief energy, discussed in more detail in section 5.2.2.
Polygons four, five, six, seven, and eight presented on the landslide density map are taken from the BGM and represent en tirely GU 2 M 2 2 . During the analyses, an interesting difference in the share of landslides was noticed within those polygons, which is also further discussed in section 5.2.2.

Relief energy map
In general, a relief energy map can be used to determine the zones in which rocks are found in a similar state of stress. Just by visual observation (Fig. 9), there is an obvious trend going from the northwest part of the study polygon (the highest values of relief energy) to the southeast part of the study polygon (the lowest values of relief energy).
In addition to the presented map, the analysis of relief energy for GUs was performed and the histogram chart was made (Fig.  10). This simple chart in a way follows the mentioned relief energy zones as the oldest geological units are located in the northwest part of the polygon with the highest average relief energy values.

Geological units prone to landslides
Numerous studies in the professional and scientific literature dis cuss the conditions under which landslides may occur. In fact, the complete procedures for obtaining landslide susceptibility are based on the analyses of various conditioning factors, i.e., geoenvironmental variables. Still, there are no established rules on which conditioning factors to select within the study, thus those selections are mostly based on expert opinion (JEBUR et al., 2014). According to REICHENBACH et al. (2018), based on re view of 565 articles, geo-environmental variables are classified into 23 classes within five clusters: i) morphological (slope, as pect, curvature, elevation, other morphometric, geomorphologi cal), ii) geological (geo-lithological, distance to fault, geo-struc tural), iii) land cover (land use/cover, soil, distance to road, tree, other anthropic), iv) hydrological (river/catchment, distance to river), and v) other variables (precipitation, earth observation, geotechnical, seismic, other, landslide related, other climatic). All of the factors mentioned above certainly affect the land slide processes, and some of them are relatively easy to measure and the related data are often available. Therefore, their use in the process of creating the susceptibility, hazard, and risk maps is logical and above all justified. However, when assessing the sus ceptibility of geological or EG units, or better yet when compared to each other, it seems necessary from an engineering point of view to return to the foundations of mechanical settings. In other words, to consider two factors that are crucial for landslide oc currence which are the stress conditions within certain parts of the relief, and the strength (internal resistance force) of the ma terial within a particular GU.
This means that conclusions about whether one GU is more susceptible to landslides than another can be drawn only when the analysed units are in a state of equal stress state, i.e., when their strength is consumed in equal or similar amounts. For such an estimate of landslide susceptibility characteristics for a par ticular GU, it is logical to use the combination of the two parameters presented in this paper, namely landslide index and relief energy.
Relief energy is defined as the height difference within the particular terrain surface, and as such, it can indicate a certain state of stress within the particular surface area. By accepting this approach, every energy relief group (Fig. 10) could be con sidered as a group that suffers similar stress. Such an approach may also imply that a comparison of landslide susceptibility levels of particular GUs can be performed only within a single energy relief group.
For estimation of the strength as a starting point, the land slide index can be used. However, for deeper analyses of landslide susceptibility within the one relief energy group, one should cer tainly consider the geological composition of the particular GU.

Relief energy group A
This group consists of Quaternary deposits the genesis of which is related to water and wind erosion. Namely, these young alluvial, proluvial, and loess sediments have been deposited on floodplains in a wider area and as such are certainly prone to landsliding. Still, endogenetic movements have not yet had the opportunity to uplift them. Consequently, this group is characterized by the lowest level of relief energy, in a range from one to fifteen m (Fig.  10). The low levels of the landslide index are therefore logical and expected. It is clear that in all the graphs (Fig. 6 and Fig. 7) the landslide index in the GUs of the relief energy group A is much lower than the average level of the landslide index for the entire study area. This observation is somewhat disturbed by a slightly larger landslide index for the GU a (alluvium - Fig. 6a and 6b;  Fig. 7a). This can be related to the scale of the geological map used, which is elaborated on in more detail in section 5.2.1.

Relief energy group B
Relief energy group B has a moderate level of energy in relation to other groups formed within this paper with the range of aver age relief energy (Fig. 10) from 25 m (GU Pl,Q) to 43 m (GU Pc). Values of the landslide indices are very scattered within relief energy group B (Fig. 6 and Fig. 7), and that is especially empha sized in Fig. 7. This scattering indicates that the GUs are very diverse in terms of susceptibility to landsliding, but it can also indicate different bedrock material, different regolith material, and finally different types and mechanisms of landsliding.
The following brief analysis of the landslide index related to the GUs of this group can serve as a starting point for future rec ommendations associated with the reclassification of the GUs defined on the BGM into the EG units. This can be especially important for the creation of the EG maps at a small scale as a crucial input data in LSM design.
Without going into deeper analyses, one could conclude that within the relief energy group B according to the landslide index (Fig. 6b), there are two groups of GUs. Namely, those with bars standing above the reference line (landslide susceptible units -Pl,Q, Pl 1 2 , M 3 1,2 , 1 M 3 1 , Pc, and K 1,2 ) and those with bars standing below (landslide unsusceptible units -2 M 2 2 , 1 M 2 1 , and ββ). Such conclusions would be very hasty for several reasons, which will be briefly elaborated below.

Statistical sample
The low landslide index for the 1 M 2 1 and ββ units presented in the scope of this research can be attributed to an insufficient statistical sample (Fig. 4). Although from an engineering point of view it would be justified to classify these GUs as unsuscep tible to landsliding according to their lithology, it would be wise to support such conclusions in future research by analysing the landslide inventory in areas where these GUs are better repre sented. Although the high landslide index is defined for GUs Pl 1 2 and Pc (Fig. 6b), it should be repeated that the statistical sample analysis for those units confirmed a questionable statistical sig nificance (Fig. 4).
GU Pl,Q covers a very large part of the study area (35.9 %) and according to the landslide index (Fig. 6b) could be labelled as a high landslide susceptible unit. Similarly, GU 2 M 2 2 covers a statistically significant area (6.7 %) and according to the landslide index (Fig. 6b) could be labelled as a low landslide susceptible unit. Still, the lithology of these units is very diverse so the land slide susceptibility level should be defined with caution. The litho logical heterogeneity of the parent rock which strongly influ ences the landslide density within those units is discussed as a part of one of the BGM limitations in section 5.2.2.
Landslide outline confidence level (LOCL) As already described, the LOCL generally defines the visibility and sharpness of a landslide on DTM derivatives used in land slide inventory creation. By accepting the assumption that the LOCL grade can be assigned to the landslide age, it is possible to conclude that GUs M 3 1,2 and 1 M 3 1 are still somewhat more sus ceptible to landslides in recent climatic conditions than GU K 1,2 , due to the higher LOCL grades. The answer to the question why unit K 1,2 is less susceptible to landslides in recent climatic condi tions than units M 3 1,2 and 1 M 3 1 , should be sought in the EG char acteristics of the near-surface parts of the above mentioned GUs.
GUs M 3 1,2 , 1 M 3 1 , and K 1,2 are characterized by sliding sur faces located in the regolith formed above the bedrock. Thus, several important features of both regolith and bedrock need to be considered to assess their susceptibility to landslides, and those are: thickness and granulometric composition of the rego lith; susceptibility of the bedrock to the mechanical weathering, i.e. formation of the regolith; and the permeability of the bedrock.
The granulometric composition, especially the proportion of the coarse-grained fraction, strongly influences the physico-me chanical properties of the regolith-like materials (PARK et al., 2018;INDRAWAN et al., 2006;SHAKOOR & COOK, 1990). With the increase of the coarse-grained fraction, the water per meability also increases, which prevents the growth of pore pres sure which strongly affects the slope stability. Similarly, the wa ter permeability of the bedrock strongly contributes to the oscillation of the pore pressure at the contact with the upper rego lith. Therefore, if the bedrock is characterised by low permeability (e.g., marl) it is very likely that the high pore pressure zone will be formed at the contact surface of the regolith and bedrock, presenting a potential sliding surface.
By taking all this into account in combination with the knowledge gained during the fieldwork and by studying the literature, it is possible to distinguish two different models.
For GUs M 3 1,2 and 1 M 3 1 the regolith can be defined as 1 to 5 m thick, with a very small proportion of a coarse-grained frac tion, and low permeability. The bedrock is susceptible to the for mation of a thicker regolith and is poorly water-permeable. With such characteristics, this EG model is defined as a high landslide susceptible one, due to the possible rapid disturbance of the me chanical stability on slopes.
For GU K 1,2, the regolith can be defined as 1 to 3 m thick, with a variable proportion of the coarse-grained fraction, and variable levels of water permeability. The bedrock is susceptible to the formation of a mediumthick regolith and with variable level of water permeability. With such characteristics this EG model is defined as being of low landslide susceptibility espe cially in the current climatic conditions.
The approach presented seems to be a good direction that can be used in the reclassification of the GUs defined within the BGM into the EGUs required for LSM development. This ap proach has already shown good results in landslide susceptibility mapping of the Sisak-Moslavina County (BOSTJANČIĆ et al., 2021), which certainly encourages future more serious attempts in a similar direction. Figure 11. a) Relief energy for the alluvial GU within the boundary outlined on the map at a scale of 1:100.000; b) Relief energy for the alluvial GU within the boundary outlined on the map at a scale of 1:5,000.

Relief energy group C
This group has the highest average relief energy within the study polygon, ranging from 53 m to 64 m, which may be a conse quence not only of more pronounced endogenetic uplift, but also the strength of the bedrock material. Landslide indices of GUs within this group indicate extremely low (T 2 and T 3 ) and low (T 1 and P 2,3 ) landslide susceptibilities (Fig. 6b). Slightly higher land slide indices for GUs T 1 and P 2,3 are very likely due to a thicker regolith superficial layer with characteristics very similar to GU K 1,2 . Geological units T 2 and T 3 are characterised by a very thin and permeable regolith superficial layer (up to 1 m), and perme able bedrock (dolomite and limestone) which is why they can be classified as having an extremely low susceptibility to landslides.

Limitations of the BGM in LSM creation
One of the pre-defined goals of this paper was to test the usability of the BGM in the process of LSM creation. Interesting observations elaborated later in this section try to present the wellknown limitations of this map in a simple and illustrative way. However, here, based on logical and numerically measurable details, closely related to the occurrence of landslides, an attempt was made to show the level of the limitations of the geological maps that cover the whole Croatian territory.

Scale
The scale of any type of map affects the level of detail that can be shown on the map and the spatial accuracy of the cartographic elements (PERKINS & PARRY, 1990), e.g., geological bounda ries and faults. How the scale affects the precision of the line ele ment (in this case the geological boundary) is shown by a simple analysis of the boundary of the alluvial deposits.
During the analysis of the Quaternary geological group, for an alluvial GU, a slightly higher landslide index was observed ( Fig. 6 and Fig. 7a). Besides, a somewhat higher value of the average relief energy was also noticed (Fig. 10). As these Quaternary sediments are deposited on floodplains, such values of landslide index and average relief energy are not logical. This inaccuracy is attributed to the scale of the BGM, which was used to deline ate the GUs polygons. In order to confirm the stated claim, anal ysis of the polygon outlining part of the alluvial deposits in the study area is presented.
In Figure 11 two different geological boundaries are shown. The first is taken from the BGM at a scale of 1:100,000 that was used in this study (Fig. 11a) and the second one is taken from the geological map at a scale of 1:5,000 which was made for the part of study area within the safEarth project (Fig. 11b). Presented polygons were used to clip the relief energy map. The cells that were left within the polygons were analysed and compared.
The polygon outlined from the small-scale map contains parts of high energy relief, while the polygon outlined from the large-scale map contains mild energy relief. The basic statistics of these differences are given in Table 2.
The age and the genesis of the alluvial material imply that these deposits are placed in plains where endogenetic movements have not yet uplifted the terrain. In that sense, the low average relief energy is largely to be expected. As the lower value of this parameter is obtained using the large-scale geological map, it cer tainly seems more accurate and thus more acceptable.
Accordingly, the BGM should be used with caution in the process of LSM creation, i.e., it can be used for the production of LSM at the same or similar scale, and in no circumstances at scales larger than 1:25,000.

Geological units defining approach
Instructions for BGM design define the approach for outlining the GUs within the different types of rocks (igneous, metamorphic and sedimentary; SAVEZNI GEOLOŠKI ZAVOD, 1962). Thus, for sedimentary rocks, the first criterion is age, while for igneous and metamorphic rocks it is lithological composition. The pre sented approach in separating sedimentary GUs from the aspect of engineering geology causes quite a problem, as lithology largely conditions the engineering geological properties. With such an approach, GUs with very heterogeneous lithological com position have been outlined, which makes the assessment of their susceptibility to landslides very difficult, and in some cases even impossible. This problem is described through the analysis of two GUs outlined within the study area, Pl,Q and 2 M 2 2 .
Pl,Q This unit consists of gravelly, sandy, silty and clayey sediments, mostly deposited in the form of lenses. Solid rocks are present in the form of breccias and conglomerates but their prevalence is almost negligible.
Only a few scientific papers analyse Pl,Q deposits (NOVO-SEL-ŠKORIĆ et al., 1986;ŠIKIĆ, 1995;HEĆIMOVIĆ, 2009), but unfortunately, none of them seem to analyse in more detail the lithological heterogeneity of these deposits. In terms of land slide activation, the importance of detecting the spatial distribu tion of individual fine and coarse-grained lithological members within Pl,Q deposits is presented. This GU is characterised by landslides with a deeper failure surface and the main triggering factor is the increase in pore pres sure due to the existence of lenses of coarse fraction material in the fine-grained soil mass.
It should certainly be emphasized that the explanation given below regarding the landslide triggering mechanism within the GU Pl,Q came from insufficiently detailed research, which sug gests that the given conclusions should be confirmed in the future by conducting research that is even more detailed.
Three zones within GU Pl,Q based on different landslide density were separated (Fig. 8), in which the main criterion was a contrasting spatial distribution of landslides. Since all three polygons were isolated in the GU Pl,Q with the primary assump tion that they share similar lithological characteristics (with the same resistance force, i.e., strength), the cause of the different landslides density had to be sought in the state of stress within each outlined polygon. Thus, the polygons were added to the en ergy relief map (Fig. 9) to determine their average relief energy. Finally, for all three zones, the landslide share and average relief energy were overlapped and analysed (Fig. 12).
Zone one, which has the highest landslide share, somewhat unexpectedly does not have the highest average relief energy. Zone two with the lowest landslide share, unexpectedly, has the highest average relief energy. Even the data related to zone three have somewhat unexpected results because that zone has the low est average relief energy and a relatively high landslide share. The analysis presented led to re-examination of the first as sumption, which is that the entire GU Pl,Q is characterized by a similar lithology. During preliminary field profiling in all three zones, the field determination of soil samples was conducted. It suggested that zones one and three are characterized by a very heterogeneous lithological environment, i.e., an environment of very diverse granulometric composition in which the fine-grained fraction is somewhat predominant. Zone 2 is more likely charac terized by coarse-grained material (sand and gravel).
This conclusion made it possible to design the two charac teristic EG models for the GU Pl,Q (Fig. 13). Model A (Fig. 13a) is related to zones one and three (Fig. 12), and model B (Fig. 13b) to zone two (Fig. 12). Model A describes landslide activation caused by increased pore pressure in the fine-grained material and erosion of the landslide toe. In model B, landslide activation in the current state of stress is not possible because the present lithology drains very easily, i.e., an increase in pore pressure is not possible. Landslide foot erosion is also not possible because surface erosive flows cannot develop due to the very rapid infil tration of surface water into deeper horizons The GU Pl,Q has been indicated in many previous works as a unit susceptible to landslides (JURAK et al., 1998;MIHALIĆ et al., 2008). The graphs in Fig. 7 certainly confirm these claims, however, when designing future LSMs (especially at a smaller scale for which new LiDAR substrates will not be made), this unit should be considered as highly heterogeneous. The recommen dations (guidelines) that should follow such maps must indicate such properties of these deposits, which would certainly require very detailed research necessary for the production of LSMs at a larger scale.

M 2
The Badenian geological complex stands out with a lower land slide index compared to other GUs within the relief energy group b ( Fig. 6 and Fig. 7). Although the BGM (ŠIKIĆ et al., 1977) con sider this complex as a single GU, from the lithological point of view (VRSALJKO et al., 2006;VRSALJKO et al., 2007) these  are very diverse deposits, which can be divided into at least two EG units: 1) EG unit 1 -biocalcarenite, conglomerate, and biolithite limestone; 2) EG unit 2 -silty marl and clayey limestone. These lithological units differ in water permeability and level of weathering durability. As such, they create regolith complexes of different thicknesses and granulometric features with very di verse water-permeable values (Table 3).
From all the above, it is obvious that the units described in this way also differ in their landslide susceptibility level, which was attempted to visualize through the following very simple analysis.
The polygons of Badenian deposits ( 2 M 2 2 ) from the BGM (Fig. 1) have been grouped according to the visual criterion of landslide density (Fig. 8). The polygons are labelled on the map from four to eight. For each polygon, the corresponding landslide index was calculated and the obtained results are presented in Figure. 14.
As expected, the outlined polygons have very different land slide indices, which very likely confirms the thesis of diverse li thology. The assumption, which has been verified to some extent by fieldwork, suggests that polygons five, six, and eight are part of EG unit 1 due to their very low landslide index. EG units 1 and 2 equally build polygon seven, due to the landslide index close to the landslide index of the entire study area. Polygon four is probably related to EG unit 2 due to the high landslide index. The pres ence of limestone within polygon four can also be confirmed by the large number of sinkholes which are characteristic of this li thology.

Neglecting the regolith
Regolith is a general term for the layer of unconsolidated or uncemented, weathered material which overlies bedrock over much of the land surface. It includes rock fragments, mineral grains, and all other superficial deposits such as unlithified in situ sap rolite, ash, colluvium, alluvium or drift (KEAREY, 2001;ALLABY & ALLABY, 2003). Superficial deposits that can be characterized as regolith, according to the presented definition, are mostly of Quaternary age, and as such, following the BGM instructions are ignored in most cases (e.g., they are not outlined if their thickness is less than 5 m; SAVEZNI GEOLOŠKI ZA VOD, 1962). During the fieldwork in the study area, it has been confirmed for older GUs (M 3 1,2 ; 1 M 3 1 ; 2 M 2 2 ; 1 M 2 1 ; Pc; K 1,2 , ββ, T 3 , T 2 , T 1, and P 2,3 ) that landslides occur almost exclusively in the regolith layer. Thereby, it can be concluded that for research re lated to landslides, this is the most important limitation of the BGM.

Conclusion
Here, a LiDAR-based landslide inventory is presented and ana lysed for the first time, for 63.8 km 2 of the Samobor area, (located in the northwestern part of Croatia). Using high-resolution DTM LiDAR derivatives, 874 landslide polygons were outlined, cover ing an area of 2.15 km2 (3.4 %).
According to the BGM at a scale of 1:100.000, the study area is characterised by strong geological diversity, where 18 GUs were distinguished. This enabled extensive spatial analysis to correlate the GUs and landslide occurrence and to define the units most susceptible to landslides. Thereby, two factors crucial for landslide occurrence were considered. Namely, the stress state expressed as relief energy, and the strength of the material described by the landslide index. The analysis is based on the approach that the same average relief energy implies the same average stress state. For this reason, landslide indices were compared for GUs within the same relief energy group. The results of statistical sample analysis were also taken into account, which clearly showed the representation of each unit in the observed area.
As a final result, GUs Pl,Q (gravel, sand, clay, and silt), M 3 1,2 (clayey limestone and -marl), and 1 M 3 1 (biocalcarenite, sand stone, marl, silt with intercalation of sand) stood out as being the most susceptible to landslides. The analysis has also shown that older GUs, Pc (conglomerate, marl, siltite, and sandstone) and K 1,2 (sandstone, shale, radiolarite and chert), are prone to landslides. Still, according to lower LOCL grades, which are given to each landslide based on the visibility and sharpness of landslide fea tures, the landslides within those units can be considered as old and inactive.
Analysis of landslide inventory together with the field prospection revealed some significant limitations of the BGM when considering their usage in LSM design. In the paper, spe cial attention was given to the limitations related to the scale, geological unit defining approach, and neglecting the regolith. However, the analysis also confirmed that the BGM has an im portant application in smallscale landslide susceptibility map ping, especially as a valuable set of maps that cover the whole territory of Croatia based on equal methodology. It is very diffi cult to expect that in the near future a superficial map for Croatia will be developed. Therefore, to try to reduce the negative im pacts of the BGM limitations, it seems logical to propose the in structions to attribute EG features to GUs outlined on the BGM. Using such instructions in the process of small-scale LSM deri vation will increase the accuracy of assessed landslide suscepti bility, which will encourage the use of LSMs, especially in the spatial planning sector.  Figure 14. Landslide indices for the Badenian polygons ( 2 M 2 2 ) shown within Fig. 8.