A case study in the research polygon in Glina and Dvor municipality, Croatia-landslide susceptibility assessment of geological units

In this paper, a preliminary analysis of the landslide inventory is presented for the wider area of the municipalities of Glina and Dvor, within Sisak-Moslavina County in Croatia, where LiDAR scanning for 45.85 km 2 was conducted. Landslide polygons were outlined based on the visual interpretation of HRDEM derivates. In total, 477 landslides were contoured with an average land - slide density of 9.85 per km 2 . Most of the landslides are characterised as moderate, shallow, and not recent. The spatial relationship between landslides and geological units is expressed with the landslide index. Subsequently, the geological units were grouped into four engineering geologi cal units representing different susceptibilities to landslides. The geological units most prone to landslides are the Eocene, Oligocene, Palaeocene and Jurassic sandstones. Even though all geological units were analysed here, the majority of landslides are within sandstones. A particu lar emphasis was on landslide occurrence in metamorphic and igneous rocks of the ophiolite sequence, a distinctive characteristic of the research area where less susceptibility to landslide processes was observed. Moreover, to further distinguish the differences between the units in the area a morphometric characteristic (relief) and drainage network was also analysed. The purpose of this analysis was to additionally confirm the landslide susceptibility assessment and the division of geological units into engineering geological units, which again implied the differ ent behaviours between landslides in igneous and metamorphic rocks compared to sandstones. Because the research area is poorly studied regarding landslide susceptibility, relief, and drain age networks, these findings will be a step forward in recognising the relationship between them and creating a base for the development of a landslide susceptibility map for this area. The most common definition of a landslide is movement of a mass of rock, debris, or earthdown a slope under gravily (1996)


INTRODUCTION
In the literature, the most used definition of a landslide is a mass of rock movement, of debris or earth down a slope under gravity (CRUDEN & VARNES, 1996). Intensive rainfall, earthquakes, and human activities are the main triggering factors of slope failure, whereas the occurrence of a landslide is controlled by different factors such as geology and geomorphology as the most important ones ( VAN WESTEN et al., 2008). Landslide processes can be better understood by preparing landslide inventory maps (SOETERS & VAN WESTEN, 1996;FELL et al., 2008;GHOSH et al., 2012). Landslide inventory maps document the extent of landslide phenomena in a region and show information about the distribution, types, pattern, recurrence, and statistics of slope failures (GUZZETTI, et al. 2012). It is important to prepare inventory maps for areas where landslides are frequent and abundant, and where slope failures are sparse or rare ( VAN DEN EECK-HAUT et al., 2007) so landslide susceptibility, hazard, and risk at the regional, national and continental scales could be determined more accurately (BRABB et al., 2000). An overview of traditional techniques and recent developments for landslide inventory mapping is given by VAN WESTEN et al. (2008) and GUZZETTI et al. (2012). One of the recently most critical techniques of landslide investigation includes the use of airborne light detection and ranging (LiDAR) method. LiDAR is used to create accurate and precise high-resolution digital elevation models (HR-

A case study in the research polygon in Glina and Dvor municipality, Croatia-landslide susceptibility assessment of geological units
Marina Filipović*, Ivan Mišur, Vlatko Gulam and Marija Horvat Croatian Geological Survey, Sachsova 2, HR-10000 Zagreb, Croatia; (*corresponding author: mfilipovic@hgi-cgs.hr) doi: 10.4154/gc.2022.04 Abstract In this paper, a preliminary analysis of the landslide inventory is presented for the wider area of the municipalities of Glina and Dvor, within Sisak-Moslavina County in Croatia, where LiDAR scanning for 45.85 km 2 was conducted. Landslide polygons were outlined based on the visual interpretation of HRDEM derivates. In total, 477 landslides were contoured with an average landslide density of 9.85 per km 2 . Most of the landslides are characterised as moderate, shallow, and not recent. The spatial relationship between landslides and geological units is expressed with the landslide index. Subsequently, the geological units were grouped into four engineering geological units representing different susceptibilities to landslides. The geological units most prone to landslides are the Eocene, Oligocene, Palaeocene and Jurassic sandstones. Even though all geological units were analysed here, the majority of landslides are within sandstones. A particular emphasis was on landslide occurrence in metamorphic and igneous rocks of the ophiolite sequence, a distinctive characteristic of the research area where less susceptibility to landslide processes was observed. Moreover, to further distinguish the differences between the units in the area a morphometric characteristic (relief) and drainage network was also analysed. The purpose of this analysis was to additionally confirm the landslide susceptibility assessment and the division of geological units into engineering geological units, which again implied the different behaviours between landslides in igneous and metamorphic rocks compared to sandstones. Because the research area is poorly studied regarding landslide susceptibility, relief, and drainage networks, these findings will be a step forward in recognising the relationship between them and creating a base for the development of a landslide susceptibility map for this area. The most common definition of a landslide is movement of a mass of rock, debris, or earthdown a slope under gravily (1996) DEM) in raster grid representations of the topography or true 3D point clouds with a high density of information (JABOYEDOFF et al., 2010) in the geographic information system (GIS) environment. As field reconnaissance is time-consuming, HRDEM derivates enable desktop mapping of landslides to be far more manageable, faster, and more precise ( VAN WESTEN et al., 2008;GUZZETTI et al. 2012). Moreover, advances in GIS technology have solved the problem of showing multiple forms of landslide information on one map (GUZZETTI et al., 2012). These new and emerging mapping methods have greatly facilitated the production and updating of landslide maps (GUZZETTI et al., 2012). This highly detailed topographic data enables interpretation of the morphological features, often associated with landslides and allows us to overcome a major limitation of mapping from aerial imagery since the morphology beneath the forest cover can be observed ( VAN DEN EECKHAUT et al., 2007).
All of this information can later serve as an input for the development of landslide susceptibility maps (LSM). Various methods and procedures have been described in the literature, explaining the basic principles and techniques to evaluate landslide susceptibility and evaluate the hazards posed therein (AHMED et al., 2014). It is rarely the case where all of the existing methods in the literature are applied during the creation of LSM because such analyses depend upon the quality and resolution of available data. Based on a review of published literature (HUTCHINSON, 1995;VAN WESTEN et al., 1997;GUZZETTI et al., 1999) LSM's categorised according to various criteria and always denoted as being quantitative or qualitative, direct or indirect, heuristic, probabilistic or deterministic, for differentiated or undifferentiated types of landslide mechanism (AHMED et al., 2014, PE-SHEVSKI et al., 2019. Within project safEarth (Transnational advanced management of land use risk through landslide susceptibility maps design), carried out from 2017 to 2019 under the Interreg IPA Cross-border Cooperation Programme Croatia-Bosnia and Herzegovina-Montenegro 2014-2020, six polygons covering an area of approximately 310 km 2 were scanned using LiDAR. In this paper, research results of one of the pilot areas in Sisak-Moslavina county covering 45.85 km 2 is presented. Bearing in mind that the research polygon is located in the municipalities of Dvor and Glina areas, where many people were displaced by the events of the Independence War in Croatia (1991Croatia ( -1995, large parts of the terrain are inaccessible and overgrown by vegetation, so the use of airborne laser scanning flights (LiDAR) was requisite. Furthermore, a scanned polygon in the study area was chosen due to characteristic geological features since a considerable part of the lithological units in the area is composed of metamorphic and igneous rocks, covering more than 20 % of the research area. The other parts of the terrain, mostly covered with sedimentary rocks, have also been a part of this study. Therefore, a landslide inventory by visually analysing HRDEM derivatives from LiDAR were made for this area for the first time. The main goal of this paper is to determine the influence of lithology, extracted from the geological map of Sisak-Moslavina county in the scale of 1:100 000 (HEĆIMOVIĆ & AVANIĆ, 2014), on the spatial distribution of landslides and to analyse geometric characteristics in different geological units (GU). Additionally, the relationship between geomorphological properties (i.e. relief), drainage density, and spring occurrence (extracted from the topographical map on the scale of 1:25 000 in areas where multiple springs appear) with the occurrence of landslides within different GU's was also examined. Also, a microscopic analysis was made to determine the major mineral composition of rocks and alteration processes caused by weathering, which consequently can lead to soil instabilities causing landslides. Furthermore, to deal with this great variety of sources prepared at different scales, a large part of the work included homogenising the data in a GIS environment (LUCIANETTI et al, 2019). As a result of these analyses, a landslide susceptibility assessment of GU's was made by grouping them into 4 geological engineering units (EGU).

Study area
The investigated area is located in Sisak-Moslavina county, northwest of the city of Glina and southeast of the city of Dvor ( Figure 1). The area is positioned in the Zrinska gora Mt. surrounded by the Una river to the east, the Sava river to the south, and the Kupa river to the north, while in the west, it borders with the Bosnia and Herzegovina. The mountain has a very developed relief with steep slopes and deeply incised stream valleys. The study area comprises 48.45 km 2 with the lowest elevation point at 167 metres a.s.l. and the highest at 552 m.a.s.l. Around 25 % of the area is in an elevation category ranging from 201-300 m.a.s.l., whereas 40 % of the investigated area belongs to the elevation category ranging from 301-400 m.a.s.l., and 30 % of the area is in the category from 401-500 m.a.s.l. Only 4 % of the area is in elevation category 167-200 m.a.s.l., and 2 % in elevation category 501-552 m. According to slope classification by DEMEK (1972), more than half of the terrain (54 %) belongs to very sloping terrain and sloping terrain (24 %), whereas 13 % can be classified as very steep. 6 % of the area is slightly sloping terrain, and only 2 % of the area is plain.
The research area belongs to the geomorphological region of the Pannonian Basin, where hills/mountains of the broader research polygon belong to a micro-geomorphological region (BOGNAR, 2001). Neogene and Quaternary regional uplift elevated the area for approximately 500 m (ČUBRILOVIĆ et al., 1967) causing the formation of relief of high erosive energy. Therefore, the investigated polygon is characterised by a branched relief. Major watercourses of the Zrinska gora Mt. are the Sunja, Petrinjčica, Žirovac (Žirovnica), Maja, and Utinja streams, while in the research area, the primary streams are the Žirovac and its tributary Stupnica, which have their confluence with the Una river near the city of Dvor.
According to the KÖPPEN & WEGENER (1924) classification, ŠAPIĆ (2012) characterised the climate of this area as a moderately warm and wet. This type of climate is characterised by small amounts of precipitation during the winter, but without a dry period during the whole year. According to the THORN-THWAITHE (1931) classification, the investigated area belongs to the humid climate area. The lowest central air temperature is during January, while the highest is in July. Data taken from meteorological stations in the vicinity show that the average annual precipitation is approximately 1000 mm and is equally distributed in all stations. The driest month is February, while the highest rainfall period is during September and November (PINTARIĆ, 2020).
The predominant vegetation is beech, sessile oak, hornbeam, and chestnut forests, while grasslands, meadows, and pastures cover the remaining area.
Half of the research polygon belongs to the municipality of Dvor, and half to the municipality of Glina. In the nearby area, several smaller villages are still inhabited. As already mentioned in the previous section, unfortunately during the Croatian War of Independence, the area was deserted, and the population in these settlements has been significantly reduced. Most households are empty, and access roads are covered with vegetation making it difficult to move around. Despite the small population, several landslides have been reported in this area, especially during the 2014 and 2018 weather extremes. This was also one of the reasons why this research polygon was chosen for LiDAR scanning.

Geological setting
Zrinska gora Mt. represents an extension of the Sana-Una unit and belongs to the Inner Dinarides (NENADIĆ et al., 2010). The geology of the Zrinska gora Mt. is heterogeneous and complex. Most of the Zrinska gora Mt. is composed of Eocene flysch and the igneous-sedimentary complex of Jurassic and Lower Cretaceous age. The central part of Zrinska gora Mt. (Šamarica and the surrounding hills) comprises Palaeogene deposits. The categorisation is based on chronostratigraphic division and labelled as GU's, 1 -12, except for ββ which is lithostratigraphic category (Table 1, Figure 1). The most dominant lithological members of these units will be emphasised in section 3.2. In the map legend, units M 2 and K 2 3 do not have a number and are not a part of the analysis because they are outside the LiDAR scanned area.

METHODS AND MATERIALS
Factors responsible for landslide occurrence can be divided into controlling and triggering factors. Slope failures can be triggered by rainfall, earthquakes, volcanic eruptions, and manmade activities. e.g. clear-cutting (PAPATHANASSIOU et al., 2021), whereas controlling factors can be represented by relevant thematic maps using GIS techniques (MORADI & REZAEI, 2014). Many landslide susceptibility studies focus on the identification of those parameters or variables that generally control the slope instability of a given area. Controlling variables can be identified by gathering the available geomorphological and environmental factor maps across the study area (AHMED et al., 2014). The first step for examining the correlation of the characteristics of the landslide occurrence, for example the landslide area with the controlling factors, is to develop an inventory map where all slope failures are reported (PAPATHANASSIOU et al., 2021). Therefore, the mein topic of this paper is analysis of LiDAR data and afterwards preparation of HRDEM derivates and finally, development of a landslide inventory. Further on, the research focuses on the correlation of a landslide inventory with: i) geology/lithology ii) relief and iii) drainage network. Spring analysis was also considered during the analyses but only within GU 12, where multiple springs appear. The data used to prepare these layers were derived from topographic maps, geological maps, and HR-DEM at different scales as already previously mentioned. Based on the analysis of these layers and field reconnaissance a categorization of lithological units into engineering geological units was made. The resultant statistical analysis and thematic maps represent a base for developing landslide susceptibility maps in the future.

Landslide inventory
The landslide inventory was derived from LiDAR post-processing HRDEM data with a spatial resolution of 0.5 x 0.5 m. Associated airborne laser scanning flights were undertaken in February and March 2018 by an external company fulfilling the requirements for at least 20 points per square metre and assuring the accuracy for each point of at most 10 cm in all directions. At the same time, topographic scanning was performed resulting in HR orthophotos with a 10 cm resolution. Morphological terrain characteristics were explored in the ArcGIS 10.2.1 environment.
LiDAR high-resolution digital terrain model (HRDEM) derivates for interpretation of the main landslide features included: i) curvature map, ii) hillshade maps calculated with sun azimuth angles of 315⁰ and 45⁰, and the sun elevation angles of 45⁰, iii) slope map displayed in green to red colour ramp, iv) orthophotos with a spatial resolution of 10 cm and v) contour line map with the elevation equidistance of 10 m. Furthermore, landslide polygons were outlined and the main landslide features were extracted by combining these derivates. Landslide inventory was made based on visual interpretation and by combining the above-mentioned HR-DEM derivates. Based on the experience in mapping landslides in other studies, a recommendation for using as many LiDAR derivates as possible, both singularly and in combinations, to recognise and precisely delineate the boundaries of individual landslide features is very important (JAGODNIK et al., 2020).
During this process, landslides were graded in a range from one to ten. The given numbers were based on visual interpretation of several main landslide characteristics (VARNES, 1978, VAN DEN EECKHAUT et al., 2012): i) the main scarp located in the upper part of a landslide as a most important feature in landslide recognition, ii) the main scarp often transforms into the flanks, which represent the landslide edges and they usually appear perpendicular and downslope to the main scarp, iii) radial and transverse cracks in a landslide body area which often contain some minor scarps and usually appear very rough, iv) the landslide toe as a curved margin of the displaced material. Therefore, number one represents a landslide with the least visible landslide features while number ten, is where almost all of the main landslide features are visible. Furthermore, based on these landslide rankings (LR) three groups were formed: low LR (grades 1-3), moderate LR (grades 4-6), and high LR (7-10). Additionally, because of the subdued morphological characteristics of some landslides, it is generally harder to consistently map landslides accurately. Typically, these subdued characteristics are associated with older landslides ( VAN DEN EECKHAUT et al., 2007, PETCHKO et al., 2016. Hence, low grade landslides could therefore represent older landslides.   Figure 3a. shows high LR where most of the main landslide features are visible such as the main scarp, hummocky topography etc. and landslide contours are easily drawn. In Figure 3b. moderate LR is presented, where the landslide main scarp and toe are clearly visible, but the borders of the landslide are not so clear and also the hummocky nature of the terrain is not so well emphasised. Finally, in Figure 3c. is an indication of the main scarp and toe, but it is difficult to follow the landslide contours. Also, it is hard to distinguish other landslide main characteristics.

Landslide geometric analysis
The geometric characteristics (shape, sizes) of landslides in this area within each GU were defined. Information on the GU's, with dominant lithology, of Sisak-Moslavina county at the scale of 1:100 000 were derived from the Croatian Geological Survey (HEĆIMOVIĆ & AVANIĆ, 2014). Landslide areas were derived automatically, whereas length and width ratios were extracted manually in the ArcGIS environment. Landslide areas in square metres are shown with Box and Whisker plots. Landslides are divided into five classes according to ŠESTANOVIĆ (2001): from very small (< 100 m), small (100 m -1,000 m), moderate (1,000 m -10,000 m), large (10,000 m -50,000 m), and very large (> 50,000 m) landslides. Basic statistical analysis of landslide cadastre data was made in ArcGIS 10.2.1 with Analysis Tools and Excel 2016. The landslide width to length ratio is expressed with the landslide aspect ratio (LA). According to the four-group classification of LA, landslides are categorised as transverse (LA≤0.8), isometric (0.8<LA≤1.2), longitudinal (1.2<LA≤3), and elongated (LA>3) (TIAN et al., 2017).
To compare the occurrence of landslides within different GU's a landslide index (LI) was calculated. The LI was used to quantitatively define the strength of the relationship between the occurrence of landslides in different GU's (LEE & TALIB, 2005). The landslide index is expressed as the ratio in percentage between the landslide area and the area of each GU (CONFORTI & IETTO, 2020). Consequently, the LI can be used to quantitatively assess the landslide susceptibility of each GU.

Relief analysis
Relief analysis is a numerical parameter determined by altitude difference between the highest and lowest points within the unit area and is calculated using the ArcGIS tool Focal statistics. The radius of a circle taken around each cell to calculate the elevation range is set to 240 m. The radius was chosen based on the approximate average slope length in the study area. In the local context, relief is conditioned by the specifics of the terrain and represents a parameter of the intensity of the development of exogenous processes (GERIĆ, 2019). Landslides, as any geological phenomena, are formed due to the simultaneous action of exogenous and endogenous forces (NIYAZOV & NURTAEV, 2017). In areas with higher relief, the intensity of erosion is higher, and in areas with lower relief, there is an increased accumulation of material. Therefore, maximum and minimum values represent the most unstable parts of the terrain, i.e. areas of the largest sediment concentrations. Therefore, this analysis was made because relief is considered one of the essential parameters influencing landslide occurrence. Regionally, relief reflects the youngest tectonic movements (LOZIĆ, 1995). A map of the relief provides an evaluation of tectonic uplift and fluvial erosive action (CONFORTI & IETTO, 2020). It is important to highlight that the interpretation of relief is always better when the influence of lithological composition on the erosion intensity and accumulation is determined (MARKOVIĆ, 1983). Relief is categorised based on the Natural breaks (Jenks) algorithm in 5 classes as follows: 13-66 m (very low relief), 67-91 m (low relief), 92-116 m (moderate relief), 107-147 m (high relief) and 147-215 m (very high relief). The Jenks natural breaks classification method was useful for comparing multiple maps created from different underlying information, giving the best group similar values and maximising the differences between classes (ZHOU et al, 2016).

Drainage network and spring features
Drainage provides a basis for understanding the initial gradient, variation in rock resistance, geological and geomorphologic his-tory of the drainage basin, or watershed (KUMAR et al., 2017). To improve the knowledge between slope and fluvial processes (NG, 2006), the drainage network was manually extracted using a topographic map of Croatia at the scale of 1:25 000. Moreover, some additional watercourses were added from HRDEM hillshade substrates. Drainage density is an important aspect of the drainage network composition, which measures the degree of drainage development within a region (HORTON, 1945). In areas close to streams, water erosion processes often cause slope undercutting, which often leads to landslides. Furthermore, debris and soil materials close to water bodies are prone to collapse during heavy rainfall (BATHRELLOS et al., 2009). A drainage density map was produced using the Line density tool in the ArcGIS spatial analyst tools to see how the formation of watercourses affects the occurrence of landslides in different GU's. Additionally, the distance to springs was also taken into consideration during the description of GU 12 behaviour, concerning the occurrence of landslides. Four categories representing low to high drainage density (I-IV) are extracted based on the Jenks Natural Breaks algorithm. Category 0 is also shown on the map and in the graphs but due to the great distance from landslides was not processed.

Field recognition and verification
Field campaigns included landslide recognition and verification of the landslide inventory in each geological unit in the research polygon ( Figure 4). Identification of landslides from LiDAR derivatives is verified with several field recognitions. Landslide inventory data can be seen in Figure 1. In October and November of 2019, the first two campaigns were focused on the processes in metamorphic and igneous rocks of the research area. During the remaining field campaigns (March and April of 2020), landslide verification was made for the rest of the lithological units of the research polygon. This research aimed to determine the depth of regolith, mineralogical and granulometric material composition using terrain methods. According to the Oxford Dictionary of Earth Sciences (ALLABY & ALLABY, 2003), regolith is a  layer of unconsolidated weathered material, including rock fragments, mineral grains, and all other superficial deposits on unaltered solid bedrock.

Analysis of the rock samples
During the field investigations, 34 rock samples were collected, from which 17 samples labelled MVIM (28 thin sections) were selected for micropetrographic analysis (Figures 1 and 5). Macroscopic and microscopic observations were performed at the Croatian Geological Survey to determine the textural and structural features and mineralogical composition of the rock samples and alteration processes. Photomicrographs were taken with a Zeiss Axio Lab. A1 HAL 35 and A5 petrographic microscopes equipped with a digital camera.

Landslide geometric analysis
The landslide inventory map was composed for an area of 48.45 km 2 in Sisak-Moslavina county. A total of 477 landslides were recognised, comprising an area of 2.32 km 2 (4,8%). Subsequently, the inventory analysis showed that the average number of land-slides per km 2 is 9.85 for the total area. Twelve different GU's are registered. In Table 1 a general overview of the main statistical parameters is given. It is evident that in GU's a, ap (GU 1); Pl, Q (GU 2), and T 2 (GU 10) landslides are not registered. In GU 1 alluvium deposits represent the youngest part of the terrain where sediments that have been eroded by water are redeposited in the lowest part of the terrain. Less than 1 % of the area within GU 2 and GU 10 is covered with LiDAR scan therefore, the absence of landslides was expected. Hence, some of the GU areas have less than 1 % of the area covered with landslides (C, D (GU 12); T 1 (GU 11), and K 1,2 (GU 5)). On the other hand, in some of the GU's (J 2,3_F (GU 8) and ββ (GU 6)), where more of the surface was detected with LiDAR scanning, landslide areas range from 1% to 2 %. GU, where the highest number of landslides were recorded, are J 2,3 (GU7); Pc (GU 4), and E, Ol (GU 3).
Looking at the size of the landslides polygons, four groups of similar size ranges prevail ( Figure 6). Average sizes range from 2700 m 2 to 3500 m 2 and belong to GU 5 and GU 8. The second groups are GU 3 and GU 4, where the average size is around 5700 m 2 . The third group consists of GU 4 and GU 7 with average sizes of around 6500 m 2 . The fourth group (GU 9 and GU 12)   Table 1.
has the highest average size for landslides ranging from 8000 m 2 to 9200 m 2 . According to ŠESTANOVIĆ (2001), all landslides in this polygon are considered moderate in size. It is also important to emphasise that maximum landslide sizes within all GU range from 13000 m 2 to 43000 m 2 . They are categorised as large, according to the scale mentioned above, except for the largest landslide (104440 m 2 ) in the polygon belonging to GU 3 which is categorised as very large. The smallest landslide (89 m 2 ) belongs to GU 7. Looking at Figure 6 it can be seen that only a few very small landslides (<100 m 2 ) are mapped. Moreover, landslide shapes and sizes (Figure 7) also show that most of them are longitudinal with an average length of 132 m and width of 51 m, which results in LA 2.58. Only in GU's 3 and 4 do average LA values exceed 3 and landslides are characterised as elongated. Nevertheless, the mean values within these units are slightly below 3. For all other GU units, values are between 2.5 and 2.8, so they are considered longitudinal. Field prospection, which also included landslide verification, gave insight into shallow soil slides prevalent in most GU's. Understanding this specific type of landslide mechanism could be dependent on the composition of land sliding material, which was identified during terrain observations as a mix of coarse-grained dust, sand, gravel, and often containing bedrock fragments. Moreover, underlining bedrock material is 0.5 to max 1.5 metres below the shallow regolith in all GU's of the research area.

Landslide index
The landslide index (LI), expressed as the ratio in percentage between the landslide area and the area of each lithology class in the researched polygon, shows values from a maximum of 10 %,   Table 1. 7 %, 5 % for the three GU's mainly consisted of sandstones (GU 3; GU 7; GU 4) (Figure 8). In GU 9 the LI is around 4 %. The dominant lithology in GU 9 is dolomite, so these values are quite unexpected. From the map (marked yellow polygon in Figure 1 with symbol A), it is notable that unit GU 9 is in contact with GU 5 where once again the dominant lithology is sandstone. Given the scale of the map used in this analysis, it is more likely that the geological boundary is not entirely accurate and that the landslides are developed in sandstone regolith that is hypsometrically above the dolomites. Lower values of LI around 1.5 % are recorded for igneous and metamorphic rocks (GU 8 and GU 6). The lowest values of LI are below 1 %, and they are observed for lime-stone, dolomite, sandstone and shale deposits (GU 10; GU 11 and GU 12).
From Figure 9 it can be seen that more than 60 % of all of the landslides in our research polygon are classified in the low LR range. GU's 1, 2 and 10 are omitted from the graph because no landslides were recorded there. The average grade for all of the GU's presented in this area is 2.7. According to the average grades, it could be concluded that the landslides on the investigated polygon are characterised by ambiguous landslide main features (VARNES, 1978), which may be an indication of older landslides. GU's mostly composed of sandstones, in this case, GU 3 and GU 7, have the most significant number of landslides in the third, high LR category, around 10 %. Also, moderate LR in these GU's is present (30%). The GU consisting primarily of metamorphic rocks (GU 8) has around 35 % of landslides in the second category and only 5 % of landslides are categorised in the high LR range. GU 6, GU 9, and GU 12 have 10 % of landslides in the moderate LR range and less than 5 % of landslides in the high LR range.

Relief analysis
The map of the morphometric parameter concerning relief is given in Figure 10. The general conclusion is that the research area is characterised by geomorphological heterogeneity. The first category has the lowest relief values ranging from 13-66 m which are mostly observed on floodplains and gently sloping terrain. The second category curve grows from 67-91 m, where gently sloping terrain is also observed. In the third category, values  range from 92-116 m. In this moderate relief category, slightly steeper terrain is observed. The fourth, high relief category ranging from 117-147 m, indicates steep terrain. The fifth and very high relief category (148-215 m) show a very high energy terrain, with steep slopes and deeply incised valleys. In the fourth and fifth categories, and also in part of the third relief category, higher values are an indication that the area is characterised by hard rocks and deep incised, narrow valleys.
From the values of relief categories in the various GU's, several differences could be distinguished. GU 1 and 2 are only represented in the first three categories with average relief values of around 80 m. GU 1 has a maximum value of 143 m even though its geological setting implies that the terrain should be plain, but due to the scale of the geological map, i.e. geological boundaries positioning, these values could be expected. Sandstones from GU  Figure 11a shows that all categories are present within GU's which cover the highest areas of the investigated polygon. Nevertheless, within GU's 6, 7, 8, 11, 12, higher categories dominate lower ones. This indicates the prevalence of harder rock materials and older geology. From Figure 11b it can be seen that in GU's 1, 2, more than 70 % of the surface is accounted for in the first Figure 11. Relief categories (I-very low relief, II-low relief, III-moderate relief, IV-high relief, and V-very high relief) within all GU's in square km (a) and percentages (b). Legend for each geological unit is presented in Table 1.

Figure 12.
Relief categories (I-very low relief, II-low relief, III-moderate relief, IV-high relief, and V-very high relief) within all GU's with mapped landslides in square km (a) and percentages (b). Legend for each geological unit is presented in Table 1. two categories. In contrast, more than 65 % of GU 3, 4, and 5 is in the second and third category, and less than 15 % is in the fourth category. For GU 6 and 11, 90 % of the area is in the last three categories. More than 65 % of GU's 7, 8, and 12 are in the second, third, and fifth categories. In GU 9, 85% of the area is in the first two categories. Figure 12 shows landslide areas within different relief categories among GU's. It is evident that in GU's where the largest number of landslides are recorded almost all relief categories are present (Figure 12a). Nevertheless, the predomination of categories I, II, and III exists. From Figure 12b it can be recognised that in GU 3, 28 % of the landslide area is in the first relief category and 33 and 36 % in the second and third categories respectively. For GU 7, 39 % is in category II, 37 % in category III, 14 % in category IV, and only 3 % in the last category. More than 80% of GU 8 and GU 12 are in categories II, III, and IV. In GU 5, 6 landslides are predominantly present in the last two categories (> 60 %). In GU 4, landslides are predominantly present in the third category. In GU 9 landslides are predominant in the first two categories, whereas in GU 11 one recorded landslide is in the fourth category.

Drainage network
Drainage network morphological analysis provided an insight into the spatial distribution of landslides in the context of drainage network development in the different GU's. Extracted drainage network order and stream density are shown in Figure 13. The investigation area has developed steeper and narrower incisions in places where the bedrock is closer to the surface. In places where sandstone rocks prevail, drainage networks are slightly wider. By visual checking, four individual basins ( Figure  13) are recognised in the research area, however complete catchments are not recorded with LIDAR, so it is assumed that this number is probably even higher.
In Figure 14, the spatial distribution of the drainage network categories within different GU's is shown. It is noticed that in   Table 1. Geologia Croatica 75/1 GU's 3, 7, 8, and 12, all categories are present (Figure 14a). From Figure 14b, it is observable that in GU 3 more than 50 % of the area is in the first two categories, whereas the second and third categories are equally distributed, only 11 % of the drainage area is in category IV. For GU 7 the first three categories are equally distributed throughout the whole area, and only 12 % of the area is in category IV. In GU 8 categories II, III, and IV are all almost equally distributed, whereas a slight domination of category III with 38 % is seen. In GU 12 the first three categories are equally distributed and there is also a slight domination of category II (40%) and the area in category IV accounts for 7 %. Given the small distribution of areas from other GU's, it is expected that the distribution of categories is not equal. In GU 5 all categories are visible, in GU 4, 6 and 9 first two categories are dominant and in GU 1 and 11 category II and III are dominant.
Looking at Figure 15a GU 3 and GU 7, where dominant lithological units are sandstones, drainage density is much higher compared to other GU's. In both GU 3 and 7, all categories are present but in different percentages (Figure 15b). For GU 3 the first three categories have similar percentages, whereas the last category has a value of 7 %. Geological unit 7 has similar values in the first two and last two categories, whereas the second category is slightly predominant with 38 %. In GUs 4, 9, and 12, landslides close to the drainage network are also present, but in much smaller quantities. In GU 4, only the first two categories are present with the predomination of category II at 71 %. A similar situation occurs with GU 9, the only difference being that the biggest area is category I with 93 %. In GU 12, 77 % is in category I. In GU 8 all categories are present with a predomination of category II at 49 %.

Microscopic analysis
From the microscopic analysis, two analysed samples defined as sandstones, six samples belong to the group of metasandstones, and metasiltstones and six are altered basalts. The other three samples are determined as diabase/dolerite, serpentinised peridotite, and radiolarian chert ( Table 2).
Metasediments of the Upper Palaeozoic complex (C, D), GU 12, are mainly composed of metasiltstone, and metasandstones with limestone interlayers. These rocks most probably are the product of hydrothermal alterations up to very low metamorphic grade. Alteration of the majority of the feldspars by processes of sericitisation and chloritization, is one of the major characteristics of these rocks. From these groups, 7 samples were analysed: MVIM 15,MVIM 16,MVIM 17,MVIM 18,MVIM 32,MVIM 33,and MVM 34 (Table 2). The main constituents of the mineral composition of the metasediments are quartz, plagioclase (often altered by sericite), white micas, and secondary chlorite, which is a product of the alteration. Lithoclasts and chert particles occur in smaller quantities.
The diabase/dolerite and altered basalts of the Middle Jurassic (ββ) (GU 6 and GU 8) are represented by analysed samples with glomeroporphyritic texture and porphyritic texture and amygdaloidal structure (MVIM 21, MVIM 23, and MVIM 24).  Table 1. Samples contain relicts of pyroxene affected by chloritization and uralitisation, resulting in secondary minerals of chlorite and actinolite. Plagioclase grains are albite, often affected by sericitisation and calcitization. Samples MVIM 2, MVIM 3, MVIM 5, and MVIM 26 have aphyric/vitrophiric and intersertal texture, where chlorite is a product of the devitrification of volcanic glass. The sample MVIM 11 is serpentinised peridotite. Olivine and orthopyroxene, under the influence of water, were turned into fibrous-leaf aggregates of serpentine (chrysotile/antigorite), located in olivine grain cracks. In addition to these minerals, talc, carbonates, and chlorite appear secondary. Iron from olivine most probably accumulates into a fine-grained aggregate of magnetite, and a part of the magnesium is bound into magnesite. Siltstones, sandstones, and cherts of the Middle-Upper Jurassic (GU 7), were represented with a sandstone sample (MVIM 12). The mineral composition of the sample is dominated by polygonal quartz grains, plagioclase, and white mica. The sample shows initial phase of chloritization and sericitisation processes rarely.
Sandstones, siltstone, and cherts from the Upper Cretaceous series (K 2 ) GU 5 are represented by diverse clastic rocks, marls, limestones, and cherts. The collected sample of sandstone (MVIM 7) is composed of polygonal quartz grains and numerous plagioclases and lithoclasts, mainly composed of quartz ag-gregates. The sample does not show foliation and acts as a slightly altered/weathered sandstone, with sparse opaque minerals probably resulting from limonitization processes.

DISCUSSION
The lithology of an area is a major factor in the formation of landslides in that area (REGMI et al., 2012). Thus, this study aimed to determine the influence of lithology, drainage, and relief on the distribution of landslides. From a lithological point of view the presented area shows a distinction between several GU's even though the geometric characteristics of the landslides are very similar. Moreover, some GU's are more prone to landslide processes than others (Figure 8). The performed analysis made it possible to group the area into four engineering geological units (EGU) as follows: 1) sandstones; 2) metamorphic and igneous rocks; 3) dolomites and metasandstones, and 4) sand and gravel ( Figure 16). Furthermore, these units should represent the basis for further analyses (e.g. landslide susceptibility map).
Landslide geometric analysis gave insight into the main landslide features. Firstly, the area was chosen because metamorphic and igneous rocks are present. The main objective was to determine the susceptibility to landslides of these units. Moreover, several units in which sandstone prevails were also covered with LiDAR scan and analysed herein. At first, it seemed that the Figure 16. EGU on the studied area.
whole area behaves as one geological unit (Figure 7) as the landslide geometric characteristics are very similar. Over the whole area, average landslides are in the moderate category according to ŠESTANOVIĆ (2001) (Figure 6). Also, from the length/width ratio it is evident that longitudinal and elongated landslides prevail (TIAN et al., 2017).
Terrain prospection also gave insight into the thickness of landslide material. Typical landslides are shallow to very shallow, not exceeding 1.5 metres. To understand which class in lithology has more influence on landslide formation, a landslide inventory map has to be combined with the lithological map of the study area (REGMI et al., 2012). The relationship between the GU's and landslide occurrence can be described with the LI. LI (Figure 8) shows a clear distinction between GU's.
Relief analysis confirmed that the terrain, in general, is characterised by heterogeneity due to the lithological and geomorphological characteristics of the investigated area (CONFORTI & IETTO, 2020). This was partly expected due to the geological characteristics of the mapped units of the area. The similarity of these units exists because older (Jurassic) geology, with slightly weathered rocks prevails. Hence, the rock mass is very near the surface and not many landslides, with low LR, are recorded (Figure 9). Most commonly, the landslide age is relative, and landslides can be defined as recent, old, or very old (MCCALPIN, 1984;ANTONINI et al., 1993). Here, high LR points to more recent, younger landslides appearing fresh on the LiDAR HRDEM derivates. Furthermore, within older landslides, usually with lower LR, more of their morphology is reshaped and the landslide features are less recognisable (MCCALPIN, 1984;KEATON & DEGRAFF, 1996;BELL et al., 2012;PETSCHKO et al., 2014). Moreover, from Figure 11 it can be seen that where more surface of a certain GU is present, more or less, all relief categories exist. Also, it can be noticed that in GU's where sandstone prevails, more surface has lower relief values whereas, in the area with metamorphic and igneous rocks, higher relief is prevalent. It is observed that in older GU's with igneous and metamorphic rocks there are somewhat fewer processes (e.g. erosion, landslides, etc.).
Drainage networks have a negative impact on landslide susceptibility due to undercutting slope toe and saturating slope material, thereby reducing the stability of the slope (DEMIR et al., 2013). The statistical analyses of drainage show a differentiation between several geological units. Firstly, the link between landslides and drainage network within some geological units can be explained due to the lithological composition. Lithology and rock structure play a vital role in the development of the drainage network in any drainage basin. The drainage patterns upon land surface develop as directed by the underlying lithology and rock structure (MUKHERJEE & JHA, 2011). In Figures 14 and 15 it is evident that in GU 3 and in GU's where sandstones prevail, the drainage network is much more dense. Drainage networks can be quantitatively described using parameters such as drainage density, which is the ratio of the total length of streams within a network to the surface area of the network. Secondly, drainage network categories are equally distributed in areas where more area of the GU was scanned with LiDAR and analysed here. Furthermore, in all other geological units, particularly the igneous and metamorphic ones it can be concluded that landslides are not dependent on the vicinity of and drainage network development.
From these results, several conclusions can be drawn. Given the heterogeneity of the presented GU's, it can be observed that some of them are more susceptible to landslides than others. Therefore, an engineering geological categorisation of the area into four zones with similar degrees of landslide susceptibility was made.
In the first EGU category (GU 3, 4, 5 and 7), where sandstones predominate, a slightly thicker regolith was recorded (> 1m) during field research. Regolith is a mixture of mud, sand, and gravel. In these units, the LI (Figure 8) ranges from 5 % to a maximum of 10 %. In the second EGU category (GU 6 and 8), where metamorphic and igneous rocks are observed, regolith values range from 0.5-1.0 m and it is also a mixture of sand and gravel with a very small amount of mud component. For this category the LI index is around 1.5 %. Regolith in the third category (GU 9, 10, and 12) has a thickness of less than 0.5 m and is a mixture of sand and gravel. The LI index shows values lower than 0.5 %. Only with GU 9, with a predominant dolomite lithology, the LI is around 4 %. This unit is classified as EGU 3 for several reasons. Firstly, during field prospection and cabinet analysis, it is observed that the boundaries between the geological units are not always precise, this is also attributed to the scale of the map used as a basis for these analyses (HEĆIMOVIĆ & AVANIĆ, 2014). Given that the landslides are present on the border of the previously mentioned dolomites and sandstones from GU 5, it is evident that the landslide started in GU 5, located hypsometrically above unit GU 9 (see the yellow polygon in Figure 1 with symbol A). Landslides are formed in a slightly thicker sandstone regolith, which subsequently slides on a solid impermeable dolomite surface. As for GU 12, it is categorised as EGU 3, even though sandstones and shale prevail which are placed in EGU 1. The decision was made because the unit is not so susceptible to landslides, as its landslide index is around 1 %. In the unit, landslides only appear near major important springs on the geological boundary with the sandstones of GU's 7 and 12 (see the green polygon in Figure 1 with symbol B). Also, the reason why the unit is not susceptible to landslide can be found in the composition of these sandstones. Microscopic analyses indicate that these sandstones are siliciclastic, and therefore during the formation of the regolith quartz, feldspars and mica are more abundant, whereas in other sandstones with slightly greater susceptibility to sliding, units may also contain clay minerals in the regolith. The fourth EGU group involves a mixture of sand, gravel, silt and clay from the Holocene (GU 1). In the aforementioned unit, there are no landslides, this is expected as the terrain is almost flat. It is important to emphasise that a mixture of gravel, sand, and clay from Quaternary deposits are not grouped into any EGU because a very small area is covered with LiDAR. Even though in the literature (JURAK et al, 1998) it is emphasised that the greatest number of landslides can be found in the area of Miocene (M 7 1-2 ) and Plioquaternary deposits (Pl,Q), if more surface would be covered, the unit would be classified as EGU 1.
Lithology is one of those parameters known to influence landslides in some regions because certain geological conditions accelerate weathering and prepare the rock for mass movements (GORETTI, 2010). Fresh rocks contain a large amount of quartz, with some calcite, sericite, and some feldspar, while the weathered rocks are rich in clay minerals and a few opaque minerals (REGMI et al., 2012). For this reason, microscopic analyses gave insight into the mineral composition of the GU's. GU 3 has not been investigated microscopically in more detail in this study as there are valid data available in the literature. Investigation of these sediments was done by TIBLJAŠ & ŠĆAVNIČAR (1985), TIBLJAŠ (1987), TIBLJAŠ et al. (2017a, TIBLJAŠ (2019) and RADIĆ (2019). They reported zeolites, stibnite and laumontite, in the form of cement and in veinlets. The authors, using X-ray diffraction on polycrystalline samples, additionally determined the Kübler (KI) and Árkaiev (AI) index as indicators of thermal changes and reported that clay minerals in the vicinity of Hrvatska Kostajnica are represented by smectites or vermiculites (swelling clay), corrensite and chlorite. These sediments correspond to our GU 3 deposits. At the same time, the presence of illite and chlorite-smectite interstratifications and minor swelling clay reported in the references of the same authors are linked to GU 5. Swelling clays are emphasised in GU 4 deposits. GU 5 clastic and metasedimentary rocks are difficult to distinguish from GU 3 units due to the absence of fossils and unclear geological relationships. It could be concluded that landslide occurrence and higher LI values in these deposits with prevalent siliciclastic sediments (GU 3, 4, and 5) can be explained by the presence of clay minerals. Clay minerals can take up large amounts of water in their interlayers, which raises pore pressure, consequently lowering the effective strength of the material. In the Middle to Upper Jurassic igneous rocks (GU 6 and 8) instabilities and sliding were not reported, therefore lower LI values were to be expected and sliding, therefore lower LI values were to be expected. In contrast, in the Middle-Upper Jurassic metamorphic rocks (GU 7) with prevalent sandstones, shallow landslides are formed. Alteration products contain clay minerals and iron oxides which can explain the higher number of landslides in this unit than units GU 6 and GU 8. Minerals found in these sandstones include quartz, white mica, feldspar which explains the lower LI values in GU 7 than the GU 3 sandstones where clay minerals are present.
The overall conclusion is that the area, in general, is characterised by very similar types of landslides. In all the lithological units, landslides are quite shallow, up to 1.5 metres, with bedrock material often very close to the surface. Moreover, landslides in these units are moderate to large, and sometimes even exceed a dozen square kilometres. Usually not all of the main landslide characteristics are visible and because of this, we suppose that most of the landslides in a studied polygon can be classified as old landslides (MCCALPIN, 1984;ANTONINI et al., 1993). Moreover, a very small number of landslides are highly ranked just for this reason (Figure 9). Still, several high graded landslides are registered near important springs in GU12, which typically is not prone to landslides, as stated above. Bearing in mind that the GU 12 on the boundary with the GU 3 with a predominant sandstone lithology, it can be assumed that the water infiltrates into the thicker sandstone regolith, and springs appear when it reaches the impermeable GU 12 metasandstone lithology (marked green polygon in Figure 1 with symbol B). Moreover, weather extremes very often affect the hydrological conditions, which can raise the pore water pressure and decrease the effective stress, which can consequently trigger landslides. Nevertheless, the extent of these phenomena also depends on the amount of precipitation, lithological composition, and terrain relief. However, morphological relief analysis showed that landslides are present on the borders of a steep and low relief where erosion processes are still active. Relief analysis also shows some differences between the analysed processes in GU's. A small number of landslides in GU 12 can be explained by the presence of the mineral quartz, which is physically and chemically resistant to weathering, white mica and plagioclase.

CONCLUSION
This study represents the first analyses of the LiDAR HRDEM derivates in identifying and mapping landslides in a 48.45 km 2 polygon in Sisak Moslavina county. A total of 477 landslides are mapped with a landslide density of 9.85 landslides per square km. Landslide sizes range from 89 m 2 to 104440 m 2 with average sizes of 4835 m 2 and a median of 2430 m 2 which can be characterised as moderate landslides. Further analyses of the specific geometric properties of these landslides can be used as a starting point in understanding the underlying structural processes controlling the movement of the landslides.
The LI was calculated to quantitatively define the strength of the relationship between the occurrence of landslides in different geological units. Consequently, LI was used as one of the main factors to quantitatively assess the landslide susceptibility of each GU. With the help of LI, it was possible to group the GUs into four EGU of rocks with similar properties. The first group of rocks, with the highest LI values, belongs to the Eocene-Oligocene (GU 3); Palaeocene (GU 4), Lower-Upper Cretaceous (GU 5), and Middle-Upper Jurassic (GU 7) where sandstone lithology prevails. Moreover, the Lower Triassic (GU 11) is also classified as EGU 1 even though the LI is very small. Therefore, grouping was made considering the prevailing sandstone lithology in the unit and that only an area of 1.35 km 2 is covered with LiDAR scanning. The second EGU group involves Middle-Upper Jurassic metamorphic and igneous rocks, of ββ (GU 6) and J 2,3 (GU 8). The third EGU groups involve Upper and Middle Triassic dolomites (GU 9 and GU 10) and Carboniferous and Devonian (GU 12) metasandstones. The fourth EGU group involves a mixture of sand, gravel, silt, and clay from the Holocene (GU 1). It is important to emphasise that a mixture of gravel, sand, and clay from Quaternary deposits are not grouped into EGU because a very small area covered with LiDAR.
Relief and drainage network parameters in regard to landslide occurrence were also analysed. The purpose of this analysis was to additionally confirm the division of EGU units based on geometric parameters and LI index. Firstly, analysis of relief showed that in GU's where sandstone prevails more surface has lower relief values whereas where metamorphic and igneous rocks are predominant higher relief is prevalent. Moreover, drainage network analysis indicate that in formations where sandstone lithology prevails, the drainage network is more dense and consequently causes weakening of slope stability, which is often the cause of landslides.
Microscopic analyses gave insight into the mineral composition of the GU's. Landslide occurrence in deposits with general sandstone lithology and higher LI values can be explained by the presence of clay minerals that can take up large amounts of water in their interlayers, which raises pore pressure, consequently lowering the effective strength of the material leading to landslide. Lower LI values and shallow landslides in the Middle-Upper Jurassic igneous rocks could be explained by the presence of more resistant minerals. Middle-Upper Jurassic metamorphic rocks with prevalent sandstones are associated with a slightly thicker regolith mainly composed of silt, sand, and gravel forming shallow landslides. This could be explained by the products of alteration which can result in the formation of clay minerals. Very low LI values can be attributed to more resistant minerals such as quartz.
The LiDAR HRDEM was advantageous for identifying and mapping landslides in this inaccessible and forested area for which a historical landslide inventory was made. Results of landslide spatial distribution in comparison to lithology, relief and drainage serve as a basis for the determination of landslide susceptibility of this area. Given the different scales of the maps used during the interpretation of the data it would be advantageous to include more detailed, larger scale, geological maps during the interpretation of LiDAR data in the future. Furthermore, LiDAR data enables the production of statistically based landslide susceptibility maps. This research is also a basis for further work on landslide susceptibility in areas with similar geological conditions and it offers a basis for a more detailed geomorphological analysis with regard to landslide susceptibility.

ACKNOWLEDGMENT
This study was carried out within the safEarth Interreg IPA Croatia-Bosnia and Herzegovina and Montenegro 2014-2020 project -Transnational advanced management of land use risk through landslide susceptibility maps design. The authors are thankful for the provision of LiDAR HRDEM data obtained during this project. The analysis of the six samples from the MVIM-21 sampling point was studied through the bachelor thesis of Matteo Blatančić. (University of Zagreb Faculty of Mining, Geology and Petroleum Engineering). Additionally, the authors appreciate Koraljka BAKRAČ (Head of Department of Geology, Croatian Geological Survey) for permission to use the Geological map of Sisak-Moslavina county 1:100.000 (HEĆIMOVIĆ & AVANIĆ, 2014) for the purpose of this scientific publication.