Use of a LiDAR-derived landslide inventory map in assessing Influencing factors for landslide susceptibility of geological units in the Petrinja area (Croatia)

A landslide inventory was created for an area of 22.6 km2 near Petrinja city in northern Croatia, based on the high-resolution LiDAR data complemented by orthophoto maps. A total of 216 landslides were identified, covering 2.91 % of that area. Landslide polygons were overlain on geological units based on the Basic Geological Map of SFRY at a scale of 1:100.000 that is the largest scale geological map available for the whole of Croatia. The relationship between landslides and geological units was expressed as a landslide index. Three geological units displayed increased landslide susceptibility. A Pliocene unit clearly had the largest susceptibility, followed by a Palaeocene-Eocene unit, and finally a Badenian unit. Landslide density was analyzed within these geological units to identify influencing factors for landslide initiation. Each geological unit revealed different influencing factors. The Pliocene unit is mostly influenced by bedding plane orientation and local relief. Heterogeneousness lithology is the dominant factor in the PaleoceneEocene unit, while the Badenian unit demonstrated the least certain interpretation as there are multiple factors involved. The forest road is presumed to be crucial, followed by spring occurrences and proximity to the tectonic boundary. The Basic Geological Map of SFRY proved to be a viable source of geological information for the creation of landslide susceptibility maps at a scale of up to 1:100.000, but with limitations in the case of lithologically heterogeneous geological units. Larger scale maps require more detailed research as landslide susceptibility factors vary in each geological unit. tions, thus avoiding or at least significantly lowering the direct consequences of landslides (CASCINI, 2008). In that regard, landslide susceptibility maps are the simplest solution. By adding more information, landslide susceptibility maps can be upgraded to hazard or risk maps (HERVÁS & BOBROWSKY, 2009). The first step in that process is development of a landslide inventory (WIECZOREK, 1983). Landslide susceptibility maps should cover large areas in order to facilitate urban planning. For that purpose, the cost of producing such maps should be kept low. To achieve that, already existing data should be used as much as possible. That can be problematic, as it often implies combining data produced at different scales and for different purposes. As compromises must be made in that process, it is important to critically evaluate the given outputs. Their final scale, reliability, purpose, and limitations should be clearly defined. In order to test the possibilities of assessing landslide susceptibility in Croatia, influencing factors were assessed for a pilot area near the city of Petrinja. It is one of six pilot areas in the safEarth project (“Transnational advanced management of land use risk through landslide susceptibility maps design”), which was carried out from 2017 to 2019 under the Interreg IPA Crossborder Cooperation Programme Croatia-Bosnia and Herzegovina-Montenegro 2014–2020. As there is no landslide inventory for the study area, one had to be created. In the scope of the safEarth project, airborne laser scanning (ALS), also known as airborne light detection and scanning (LiDAR) was performed for an area of approximately 300 km2, including the study area Article history: Received October 07, 2021 Revised manuscript accepted December 24, 2021 Available online February 28, 2022


INTRODUCTION
Landslides are one of the natural hazards that has a major impact on human lives. They cause major economic losses, but also direct loss of life, and are present worldwide. Quantification of those losses is difficult and incomplete. According to the World Health Organization, between 1998 -2013, landslides have caused more than 18.000 deaths and affected 4.8 million people (CRED, 2018). The European Environment Agency (2010), reported 312 deaths due to 70 major landslides in the period of 1998in Europe. HAQUE et al. (2016 presented a spatiotemporal distribution of deadly landslides in 27 European countries over 20 years (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) with a total of 1370 deaths and an average economic loss of approximately 4.7 billion Euros per year. All these studies are based on reports from major landslides, and all agree that the numbers are certainly underestimated as the majority of landslides are unreported. This is especially true for countries outside "global landslide hotspots" (Himalayan arc, Central America, etc.) with numerous landslides with catastrophic consequences (KLOSE et al., 2015). According to KLOSE et al. (2015), most of the landslides in Germany are individually of smaller consequence, so are mostly unreported. However, when summarized and collated, the costs of those smaller events are higher than the combined major, 100-year events.
Once triggered, landslides can cause significant damage and loss of life. Remediation of the affected slopes and infrastructure also brings significant costs. The most cost-effective means of reducing those costs is prevention. By identifying areas susceptible to landslides, urbanization can avoid or adapt to the condi-

Use of a LiDAR-derived landslide inventory map in assessing Influencing factors for landslide susceptibility of geological units in the Petrinja area (Croatia)
Tihomir Frangen*, Mirja Pavić, Vlatko Gulam and Tomislav Kurečić Croatian Geological Survey, Sachsova 2, 10000 Zagreb, Croatia * corresponding author: tihomir.frangen@hgi-cgs.hr doi: 10.4154/gc.2022.10 Abstract A landslide inventory was created for an area of 22.6 km 2 near Petrinja city in northern Croatia, based on the high-resolution LiDAR data complemented by orthophoto maps. A total of 216 landslides were identified, covering 2.91 % of that area. Landslide polygons were overlain on geological units based on the Basic Geological Map of SFRY at a scale of 1:100.000 that is the largest scale geological map available for the whole of Croatia. The relationship between landslides and geological units was expressed as a landslide index. Three geological units displayed increased landslide susceptibility. A Pliocene unit clearly had the largest susceptibility, followed by a Palaeocene-Eocene unit, and finally a Badenian unit. Landslide density was analyzed within these geological units to identify influencing factors for landslide initiation. Each geological unit revealed different influencing factors. The Pliocene unit is mostly influenced by bedding plane orientation and local relief. Heterogeneousness lithology is the dominant factor in the Paleocene-Eocene unit, while the Badenian unit demonstrated the least certain interpretation as there are multiple factors involved. The forest road is presumed to be crucial, followed by spring occurrences and proximity to the tectonic boundary. The Basic Geological Map of SFRY proved to be a viable source of geological information for the creation of landslide susceptibility maps at a scale of up to 1:100.000, but with limitations in the case of lithologically heterogeneous geological units. Larger scale maps require more detailed research as landslide susceptibility factors vary in each geological unit. tions, thus avoiding or at least significantly lowering the direct consequences of landslides (CASCINI, 2008). In that regard, landslide susceptibility maps are the simplest solution. By adding more information, landslide susceptibility maps can be upgraded to hazard or risk maps (HERVÁS & BOBROWSKY, 2009). The first step in that process is development of a landslide inventory (WIECZOREK, 1983).
Landslide susceptibility maps should cover large areas in order to facilitate urban planning. For that purpose, the cost of producing such maps should be kept low. To achieve that, already existing data should be used as much as possible. That can be problematic, as it often implies combining data produced at different scales and for different purposes. As compromises must be made in that process, it is important to critically evaluate the given outputs. Their final scale, reliability, purpose, and limitations should be clearly defined.
In order to test the possibilities of assessing landslide susceptibility in Croatia, influencing factors were assessed for a pilot area near the city of Petrinja. It is one of six pilot areas in the safEarth project ("Transnational advanced management of land use risk through landslide susceptibility maps design"), which was carried out from 2017 to 2019 under the Interreg IPA Crossborder Cooperation Programme Croatia-Bosnia and Herzegovina-Montenegro 2014-2020. As there is no landslide inventory for the study area, one had to be created. In the scope of the safEarth project, airborne laser scanning (ALS), also known as airborne light detection and scanning (LiDAR) was performed for an area of approximately 300 km 2 , including the study area presented here. LiDAR is a well-known data source for derivation of a landslide inventory (JABOYEDOFF et al., 2012) and is especially well suited for forested areas (SLATTON et al., 2007), as in the study area, as it can penetrate the canopy and derive bare earth digital terrain models (DTM). Besides the landslide inventory, an accurate and precise DTM derived from LiDAR data facilitates preparation of other required data, including slope, aspect, and relief energy. Parallel to the LiDAR scanning, aerial photographs were taken that allowed construction of an orthophoto map which supplemented the LiDAR data in the analysis.
Besides the landslide inventory and LiDAR data derivatives, geology provides crucial data for assessing landslide susceptibili ty (VAN WESTEN et al., 2008). For this purpose, geological data used was from the Basic Geological Map of SFRY 1:100.000 (PIKIJA, 1987a), as it is the largest scale map available for the whole of Croatia and the study area. Its usefulness and limitations were tested, within the scope of geological units (GUs) present in the study area.
To this day, there is no universally defined method for deriving landslide susceptibility models (REICHENBACH et al., 2018). The reasons for this are many, from different geological, climatic, and physiographical settings to data availability and the purpose for deriving models. This study does not aim to change any of this, as it covers a small area with limited influencing factors. Despite that, the aim is to guide future landslide investigations in the region, or similar areas, in methods of data usage, their influence, and limitations.

STUDY AREA
The study area is located in the northern part of Croatia and is defined within the safEarth project co-financed by Interreg IPA CBC (Croatia -Bosnia and Herzegovina -Montenegro). It covers an area of 22.6 km 2 (Fig. 1). The study area belongs to Sisačko -moslavačka County, SW of the Petrinja City. It is bounded by the Kupa, Petrinjčica, and Sanja rivers, north, east, and west, respectively. The river Utinja flows in a SE-NW direction towards the Kupa river passing through the middle of the study area. The relief is characterized by smaller hills intersected by creek valleys, with an average height of 150-250 metres, with the highest peak being Kosinjak (324 m a.s.l.). The climate is temperate continental under the mild influence of the Mediterranean climate of the northern Adriatic. The atmospheric conditions are very variable and intense exchanges of weather situations occur during the year (ZANINOVIĆ et al., 2008). The average monthly temperature ranges from -0.2 ° C for the lowest average (January) and 21.5 ° C as the highest average (June) for the nearest city of Sisak for which data are available (for the period 1949 -2019). For the same period of climate monitoring, the Meteorological and Hydrological Service of Croatia (DHMZ) recorded the average precipitation amount being the lowest in March (55.1 mm) and highest in November (93.5 mm). The average annual temperature in the City of Petrinja is 10.9 ° C (period 2013 -2020, DHMZ), while the average annual rainfall is around 1088.85 mm with a maximum of 1400 mm, measured at the closest meteorological station in Sisak (over the period 1949-2019, DHMZ). The Sisak meteo station is located approx. 13 km in a direct line to the centre of the study area polygon.

GEOLOGICAL SETTINGS
The investigated area is located at the SW margin of the Pannonian Basin System (HORVÁTH et al., 2015;PAVELIĆ & KOVAČIĆ, 2018) along the border of the Sava depression (SAFTIĆ et al., 2003). Only the NE marginal part belongs to the Neogene-Quaternary tectonic unit of the Sava depression, while the majority of the area belongs to the relatively elevated area of the tectonic unit of the Inner Dinarides (PIKIJA, 1987a(PIKIJA, , 1987b. These large-scale structures are bordered by Dinaric NW-SEoriented fault systems. It is still a highly tectonically active area with numerous documented landslides triggered by earthquakes and neotectonic movements along the main faults (MARKUŠIĆ et al., 2021;POLLAK et al., 2021).
The area is almost completely covered with Neogene-Quaternary deposits with occasional occurrences of pre-Neogene basement Palaeogene rocks (PIKIJA, 1987a). Detailed lithological data concerning the deposits of the research area are scarce, due to the chronostratigraphic approach used when the newest geological map of the area was developed (PIKIJA, 1987a). However, the basic lithology of the mapped units is given (PIKIJA, 1987b) and described in the following text in chronostratigraphic order using original mapping nomenclature (Fig. 2).

Palaeocene-Eocene (Pc,E): Conglomerates, sandstones, siltstones, and sandy marls
These deposits are exposed in the central part of the investigated area (Fig. 2), along the hills surrounding the river Utinja. According to PIKIJA (1987b), the Palaeocene-Eocene complex in the wider area is represented by conglomerates intercalated with other lithological members, namely sandstones, siltstones, sandy marls, and limestones. The conglomerates are the most distributed lithology, they are polymictic with variable content of sandy matrix, and are mainly carbonate-free. The clasts originate from the Upper Cretaceous sediments, volcanic rocks, and pyroclastics from the vicinity. Therefore, the clast composition varies. Mafic effusive rocks, granite, gneiss, crystalline schists, metasandstones, chert, and sandstones were determined amongst others (PIKIJA, 1987b). Clasts are 1-20 cm in diameter, rarely 20-40 cm. An estimated thickness of the deposits is around 200 m.

Badenian (M 2 2 ): Conglomerates, sandstones, limestones, and marls
Badenian deposits are transgressive to older (Palaeogene-Eocene) sedimentary rocks in the southern parts of the investigated area (Fig. 2). To the north, the Badenian unit is in tectonic contact with surrounding GUs or is mostly covered with the Holocene alluvium of the river Utinja. In the basal part, the Badenian deposits are represented by an accumulation of coarse-grained sediments (conglomerates, gravel, sandstones, biocalcarenites). In the upper part of the sequence, fine-grained sediments predominate with increased carbonate content up to 98% (clays and marls, bioclastic and biogenic limestones) (PIKIJA, 1987b). Generally, the most prominent lithology within the Badenian deposits is represented by shallow-water deposits i.e., reef and near-reef facies deposits: lithothamnium, coralline-lithothamnium, and bryozoan-lithothamnium limestones, and forereef biocalcrudites and biocalcarenites. These deposits are characterized by high porosity (up to 45%). Coeval deposited marls, limestone-marls and micritic limestones can also be observed in the area. The thickness of Badenian sediments varies up to a maximum of 200 m.
Lower Sarmatian ( 1 M 3 1 ): Marls, limestones, conglomerates, and sands Sarmatian deposits occur mostly in the SW part of the investigated area. They are continuously deposited onto Badenian deposits (Fig. 2). Composed of various lithologies shallow and deepwater deposits can be observed. Occasionally coarsegrained shallow water clastic rocks (polymictic conglomerates, conglomerate-sandstones, and gravelly sands) crop out, followed by shallow-water ooid limestones, marl-limestones, and calcarenites with porosity up to 40% (PIKIJA, 1987b). In stagnant and deep-water environments, thin-plate and laminated deposits are developed, also with high lithological variability. Marl predominates, followed by clayey limestones, and limestones. Occasionally they are intercalated with sands and biocalcarenites. The thickness of the Sarmatian deposits in the wider area is estimated to be up to 70 m.

Lower Pannonian ( 1 M 3 1,2 ): Limestones and marls
Lower Pannonian deposits comprise a thin belt following the zone with Sarmatian deposits along the SW part of the investigated area ( Fig. 2). They are represented by facies of thin-bedded limestones, and marls. Deposited continuously on the Lower Sarmatian deposits, thin-bedded micritic limestones or silty limesto nes are homogenous with a carbonate content of between 85 and 95% (PIKIJA, 1987b). Marls and limestone-marls are partly laminated (resembling the Sarmatian deposits) with carbonate contents between 50 and 80%. The thickness of the Lower Pannonian deposits in the area could be estimated as between 20-50 m.

Upper Pannonian ( 2 M 3 1,2 ): Marls
Upper Pannonian deposits cover most of the SW part of the investigated area. Deposited conformably on the Lower Pannonian deposits (Fig. 2), they show lithological similarities with the underlying unit (PIKIJA, 1987b). A gradual decrease in the carbonate content is observable. The majority of the unit is composed of yellow to brown/gray massive marls (35-75% carbonates) or rarely clayey marls. The overall thickness of the Upper Pannonian deposits can reach 200 m.

Pliocene (Pl 2,3 ): Clay, silt, gravel, sand, sandstones, conglomerates, lignite
Pliocene deposits occupy the NE part of the investigated area ( Fig. 2). They are discordantly deposited on the basement composed of various units, such as Pc/E, and Badenian deposits, and are overlain by Quaternary dpr and alluvial sediments (PIKIJA, 1987a) covering the largest area of the investigated polygon. These deposits known as the Viviparus beds are widely distributed in the SW part of the Pannonian Basin MANDIC et al., 2015;NEUMAYR & PAUL, 1875;RUNDIĆ et al., 2016). Pliocene deposits are mostly clastic and heterogenous: sands, gravel, silts, and clay with a general coarsening upwards trend, and occasionally conglomerates and lignite. In the nearby area of the Vukomeričke gorice (toward the NW) Pliocene deposits can be divided into five basic lithofacies units: facies of massive and normally graduated gravels, facies of cross stratified sands and gravelly sands, facies of massive and laminated sands, facies of clay silts, and heterolithic facies (KUREČIĆ, 2017). Quartz is the most dominant constituent of the sand, lithic fragments are less common, ranging from 13-54%, while feldspar grains are the least represented . Sandstones (lithoarenites) are limonitized and form thin interlayers within the sand. Occasionally polymictic conglomerates can be observed. The well-rounded clasts are up to 1 cm in diameter, composed of radiolarian chert, silicified vitreous tuff, quartzite, and quartz schist. The thickness of the deposits in the area has been estimated as 200-400 m (PIKIJA, 1987b).

Plio-Quaternary (Pl, Q): Gravel, sands, and clay
Pliocene-Quaternary deposits occur only sporadically, lying discordantly over various older lithological units (Fig. 2). In the area, they are composed of sand, gravel, clay, and rarely sandstones and conglomerates (PIKIJA, 1987b). Sands and gravels predominate in the area. Finer-grained sediments are commonly moderately sorted, occasionally fine quartz sands occur (up to 96% SiO 2 , PIKIJA 1987b). There is great variability of sediment mixtures ranging from clayey sands to gravelly sands and gravels with rounded pebbles up to 3 cm in diameter. Gray to brown clayey sediments occur in form of lenses and thin layers within sands. Those sediments are mineralogically similar to the Pliocene sediments with the predominance of zircon, tourmaline, epidote, and rutile within the heavy mineral fraction. Pebbles in gravel consist of chert, pyroclastics, and sandstones. Lithified deposits occur only rarely. The estimated thickness of the unit is up to 100 m.

Holocene (dpr): Diluvial-proluvial silt, sand, and gravel
Those poorly sorted and chaotic sediments occur only in thin zones along the contact zone of valleys and hills in the NE part of the investigated area. The composition of sediments is closely related to the composition of the source rocks in the close vicinity. The thickness of the dpr sediments is estimated at 10 m (PIKIJA, 1987b).

Holocene (a, a 1 , ap): Alluvial sand, gravel, silt and clay
On the NE edge of the investigated area, there are occurrences of the River Kupa terrace sediments (a 1 ). Building up flat areas, these sediments consist mostly of silt, intertwined with sand and gravel layers (PIKIJA, 1987b). Sporadically, there are also occurrences of floodplain clayey-sandy silts (ap) with a maximum thickness of up to 5 m. Sediments of recent alluvium (a) are associated with the river Utinja. They show a variable granulometric composition in the range from silt to gravel, mineralogically reflecting the local drainage area.
According to the area zonation by the degree of stability, on the Engineering Geological Map of the SFRY at a scale of 1:500.000 (ČUBRILOVIĆ et al., 1967), the largest part of the study area is defined as a "prevailing unstable area under natural conditions and under the man's action mostly unstable". The lesser part is defined as "prevailing stable area under natural conditions while under the man's action they can be prevailingly unstable". Only flat, alluvium areas are defined as stable areas.

DATA AND METHODS
The basis of this study was LiDAR data that will be explained in detail in the next section. LiDAR, along with its derivatives and together with orthophotos, was used for creating the landslide inventory.

Remote sensing data sets
The basis of this study is LiDAR data acquired during the safEarth project. Data acquisition was performed using a heli-copter, during the spring of 2018, by Flycom Technologies d.o.o. LiDAR used for the acquisition was Riegl LMS-Q780 and for RGB photographs a Hasselblad H60 camera with 50 mm lens was used. Ground control points were utilised for assuring the specified accuracy.
Laser scanning produced a point cloud with a minimum point density of 20 points per m 2 . The accuracy of every point was in the range of ±10 cm in each direction. Full waveform data were recorded which facilitated the separation of the ground from vegetation (MALLET & BRETAR, 2009). The described point cloud was used for deriving a Digital Terrain Model (DTM) with a cell size of 0.5 x 0.5 m. The DTM represents bare earth, stripped of vegetation and all man-made objects, and is therefore suitable for analysis of landslides and other morphological processes (RAZAK et al., 2011).
Photographs that were taken with an RGB camera were corrected with elevation data derived from LiDAR. This produced orthophoto maps with 0.1 x 0.1 m pixel size.
LiDAR-derived DTM was used for the creation of several derivatives useful for terrain visualization, and hence, landslide mapping: 1. Hillshade grid -shaded relief in the same resolution as DTM. Two grids were produced with different positions of the light source. The position of the first light source was an azimuth of 315° and 45° elevation. The other had the same elevation, but the azimuth was 45°. Combining them enabled landslides to be spotted in all orientations. 2. Slope angle grid -a slope angle of the terrain was also produced at a resolution of 0.5 x 0.5 m. It was represented in a continuous colour scale as that enabled the best visual representation of the terrain. 3. Aspect grid -used for obtaining the exact orientation of every cell. Resolution was also 0.5 x 0.5 m. 4. Contour lines -contours with equidistance of 1 m were used for visual help in recognizing changes in the slope shape due to landslides.

Creation of a Landslide inventory
The landslide inventory was created for the area that was covered with LiDAR data, as that was the source of input data for landslide detection. All LiDAR derivatives described in the previous subchapter were used in 2D using ESRI ArcMap 10.2 software and in 3D using ESRI ArcScene software. Landslides were contoured subjectively, using expert judgment (GUZZETTI et al., 2012). To avoid differences in a subjective approach, the whole inventory was created by one expert. In this way, all inaccuracies due to bias were the same across the whole inventory. All visible morphological signs of a landslide were evaluated. For an area to qualify as a landslide, at least 50% of its boundary length had to be recognizable (main scarp, lateral flanks, toe). The internal body of a landslide has to show some signs of movement (undulated terrain, cracks, minor scarps, water pooling, etc.). The whole area with material displacement was contoured, regardless of whether the material was carried away or deposited (from the top of the main scarp to the bottom of the toe).
All landslides in the inventory were given a grade, on a scale of 1 to 10.A lower-grade indicates a landslide with less pronounced features and overall lower reliability. In contrast, a higher grade indicates more pronounced features with higher reliability, which also implies more recent activity. Grades 9 and 10 are reserved for landslides with clearly visible features on all datasets (LiDAR and orthophoto). Since the area is forested over a large part, most landslides are not visible on orthophotos, unless they are very recent. Lower grades of 1 and 2 are used for areas with subdued landslide features and a possibility that they are wrongly interpreted as landslides.
In the areas with overlapping landslide bodies, or multiple landslide generations, each polygon was delineated separately, where possible. Each landslide body is unique, and no overlapping is allowed.

Geological data
GUs used for analysis are almost exclusively from the Basic Geological Map of SFRY at a scale of 1:100.000 (PIKIJA, 1987a). That data was used because it is the largest scale geological map available for the whole of Croatia and the only one that covers the entire study area. Although there are some obvious inaccuracies in it, (inevitably due to the small scale and available technology at the time of production), there have been no corrections applied to it, as one of the goals of this study is an assessment of their use in landslide susceptibility assessment in the Croatian territory.

Fieldwork
Fieldwork was carried out on two occasions. The purpose of the first fieldwork was the verification of the landslide inventory made exclusively with LiDAR data. Verification has been done on the easternmost part of the area. A second field trip was targeted at specific areas, which had anomalies in the number of landslides that could not be resolved using only GIS analysis, or needed field verification of the presumed explanations.

GIS analysis
All of the afore-mentioned data were aggregated in GIS using ESRI ArcGIS 10.2 with Spatial Analyst software.
LiDAR-derived data were analyzed using standard ArcGIS tools. Depending on the aims and purpose, data were extracted for different target groups (GUs, landslides, landslides in specific GUs). Basic statistical analysis was undertaken on such extracted data sets.

Landslide index
The landslide index is calculated as a ratio between landslide area in a given GU and the area of that GU and is expressed as a in percentage. It shows relationships between the GUs and landslide occurrence within the study area. It was also calculated for the entire area for comparison (BARTELLETTI et al., 2017;CONFORTI & IETTO, 2020).

Local relief
The map of local relief represents the maximum difference in elevation for the unit area (BUCCI et al., 2016;CONFORTI & IETTO, 2020;DELLA SETA et al., 2004;MOLIN et al., 2012). Historically, the unit area used was 1 km 2 . In more recent times, the unit area varies depending on the intended usage and the author's preferences. In this study, a local relief map was produced using the Focal statistics tool in ArcGIS. With it, for each cell in a raster, the maximum range in elevation was determined in a circle with a radius of 250 m. The radius was chosen based on the approximate average slope length in the study area.

Aspect and slope
Aspect and slope angle were calculated for each landslide body as an average of all the cells in it. They were also calculated for the whole area of each GU, on the basis of 5 x 5 m cells. Additionally, for each 5 x 5 m cell, dip direction and dip angle were calculated and their poles were analyzed in Stereonet 1.1 software (ALLMENDINGER et al., 2011;CARDOZO & ALLMEND-INGER, 2013). Due to hardware processing limitations, a subset of 10000 random cells was extracted from a GU. Their poles were projected on a stereonet and their density was shown using the Kamb contouring method (KAMB, 1959). In addition, the mean vector was calculated to get the average orientation of the slopes in a particular GU.

Inventory analysis
The investigated polygon, covered by LiDAR around Petrinja, has an area of 22.57 km 2 . The total area of landslides in that poly-gon is 0.66 km 2 , or 2.91 % of the total area, spread over 216 landslides. The average landslide density is 9.6 landslides per km 2 .
The number of landslides based on their area is shown in Fig. 3.
The smallest landslide has an area of 59 m 2 and the largest 51,673 m 2 . They are the only landslides belonging to the very small and very large landslide classes, respectively. Small landslides represent 44.9 % of the total and the majority belong to the moderate area class at 48.1 %. In accordance with this, the avera ge and median landslide area also falls into the moderate size class, with 3,047 m 2 and 1,156 m 2 , respectively. Large area landslides represent only 6 % of all landslides.
Based on a landslide appearance in the lidar data, each landslide was given a grade in the range of one to ten. A higher grade indicates a more recent, active landslide while low grades imply old and uncertain landslides. None of the landslides were assigned to grade ten. Roughly half of the landslides belong to the first two categories. Generally, the number of landslides tapers down with each higher grade (Fig. 4).  The size of the individual landslides is quite uniform across the grades, although the largest landslides are mostly related to lower grades (Fig. 5).
The dominant aspect of the landslides is towards northeast, north, and northwest (Fig. 6).

Geological units
Of eleven GUs in the studied area, five of them do not have any recognizable landslides, and two of them have only one landslide. The majority of landslides are situated in Pl 2,3 GU, which is also the most widespread (36.9% of surface area) (Table 1) (Fig. 7).   Based on the results, only three GUs were used for further analysis. Those are the Pliocene unit (Pl 2,3 ) , Palaeocene-Eocene unit (Pc,E), and the Badenian unit (M 2 2 ). Although the Holocene Diluvial-proluvial unit (dpr) has a relatively high landslide index, it was omitted from further analysis as that index results from only one landslide which raised the landslide index as it covers a large area. However, it is situated on the boundary with the Pliocene unit, so there is a large possibility that it belongs to that unit. Verification was impossible since the area is suburban and covered with regolith and vegetation.
The average slope angle of landslides in the three separated GUs is presented in Fig. 8.

Pl 2,3 unit
The Pliocene unit has 182 registered landslides, which is 84% of all landslides in the study area, while the GU covers 36.85% of the studied area. With 6.7% of its surface area covered by landslides, it has the biggest landslide susceptibility of all GUs in the area.
Overlaying landslide density on a local relief map shows increased values of local relief in areas with landslides (Fig. 9). In Pliocene deposits, the average local relief value in the landslide area is 60.5 m, while it is only 37.4 m in the area without landslides.
Analysis of slope and landslide aspect shows the pronounced inclination of landslides towards the northeast (Fig. 10). At the same time, the slope aspect has a much more homogenous orientation, with a slight inclination towards the north and south.
Projecting poles of 5 x 5 m cells of all slopes in this GU, we get a fairly homogeneous spread over all orientations. This is further confirmed by the mean vector with the pole dip direction of 196.6° and dip angle of 88.9°, so there is only a slight tendency of terrain sloping towards the north -northeast (Fig. 11a).
In the case of slope cells within landslides, there is a pronounced inclination towards the northeast. The mean pole vector of all cells has a dip direction of 209.1° and a dip angle of 83° (Fig.  11b).
The same analysis has been undertaken with bedding planes orientation of the Pliocene strata. As the area is mostly covered by regolith and vegetation, bedding planes are hard to recognize, there are only eight measured layers in the literature in this area and four of them are horizontal (PIKIJA, 1983(PIKIJA, , 1987a. Therefore, these results should be taken with caution. However, they show    a clear general inclination of bedding planes towards the northeast, with a mean vector pole dip direction of 214° and dip angle of 83.8° (Fig. 11c)

Pc,E unit
Although the Palaeocene-Eocene unit contains only 13 landslides, it still has a high landslide index of 4.2 as it is represented by 1.23 km 2 (5.45% of the whole area). Furthermore, it is split into two distinct parts. One is in the eastern part of the investigated area, the other is in the western part. All but one landslide are situated in the eastern part (Fig. 9).
A field investigation determined that although both areas belong to the Pc,E chronostratigraphic unit, their lithology is different. The Eastern part is composed of sandstones and siltstones while the western part is dominated by conglomerates (Fig. 12).

M 2 2 unit
The Badenian unit is the second-ranked GU in regards to the number of landslides, containing 16 landslides. However, with a 3.12 km 2 surface area, the Badenian unit covers 13.82% of the investigated area, hence, its landslide index is relatively low at 0.9. That puts it below the landslide index of the whole investigated area, which is 2.9.
Looking at Fig. 9, it is clar that landslides are not uniformly spread in the Badenian unit. There is an obvious hotspot in the easternmost part of the investigated area.
While the Badenian unit is not lithologically homogenous, the observed difference in landslide density cannot be attributed to lithological variation. The area affected by landslides has a limestone lithology (Fig. 13a), which would imply lower landslide Figure 13. Scarp of the landslide in the regolith of the Badenian limestones geological unit (a) and a channel formed from the road runoff water that leads to a landslide body below the road (b).  density as limestones form thinner regolith, with less clay content and more sharp-edged limestone fragments.
After field investigation, several contributing factors were recognized. The dominant one is a forest road which is in direct contact with the majority of the landslides (Fig. 13). It disturbed the mass balance by cutting into the slope and affected the slope by concentrating the drainage of runoff water from the road surface (Fig. 13b). Another possible contributing factor is the geological boundary with Pc,E GU, together with spring occurrences alongside it (Fig. 14b).
There is no available data for the construction of the forest road. It cannot be seen on the orthophoto map from 1968. (Fig. 14c). There are also no visible landslides, but since the area is forested, it cannot be concluded that there are none present. Orthophoto map made during the year 2018 clearly shows the forest road and one landslide can be recognized (Fig. 14d).
Inventory in this study lacks information on the triggering time. It is a geomorphological historical inventory in which the age of the landslides is not known (GUZZETTI et al., 2012;MALAMUD et al., 2004). Landslides are given a grade, which is influenced by age, but also other factors including human impact, erosion intensity, type of movement, etc. (BELL et al., 2012;McCALPIN, 1984). Since environmental conditions change over time, older landslides might not be relevant to today's conditions. However, determining the landslide age and appropriate criteria for the exclusion of certain landslides is difficult to determine and should be a topic of further investigation. The grades given to landslides are not used for analysis, but are useful for evaluating confidence in the results. The main criteria for choosing relevant landslides were visibility of at least half of their boundary. However, it is hard to evaluate the validity of that condition. Although there are many slopes with signs of creeping movement, they are left out of the inventory and further analysis, as they do not have a defined boundary (CRUDEN & VARNES, 1996) and delineating creeping movements proved to be very subjective.
Landslides in the studied area are predominantly moderate and small in area. The largest landslides are mostly related to lower grades and that can have multiple explanations. One possibility is that the conditions in the past were different from today, so there was a tendency for bigger landslides. If that is the case, that landslides would be best left out of the landslide susceptibility analysis. However, it is also possible that they represent multiple smaller landslides, with mutual boundaries masked by erosion. Therefore, leaving them out of analysis would cause a significant underestimation of landslide susceptibility. Differentiation of these two cases is problematic. In the absence of conclusive evidence for attributing lower grade landslides to different conditions in the past, it is safer to include them in the analysis.

Geological units
GUs are a crucial element in landslide formation as they supply materials of different properties. However, determining their influence is often problematic. The most common geological information in landslide susceptibility models consists of chronostratigraphic units in the bedrock, as shown on standard geological maps (REICHENBACH et al., 2018). As stated in REICHEN-BACH et al. (2018), that approach may be problematic as the relationship between chronostratigraphic units and mechanical properties of materials involved in landsliding may be unclear. Chronostratigraphic units may contain a variety of lithological formations with different mechanical properties. This problem can be seen in this study, in the case of Pc,E GU, where there are two distinct lithologies with different landslide susceptibilities.
Since the majority of landslides in the inventory are shallow, translational earth slides (CRUDEN & VARNES, 1996), they are mostly originating in the regolith (ALLABY, 2020). Although regolith composition is dependent on the bedrock, it can deviate from it, especially near geological boundaries or in the case of rugged terrain.
There is also a scale discrepancy. Geological maps used for analysis are at the scale of 1:100.000, while lidar data is equivalent to scale 1:1.000 (TOBLER, 1987). Such a big difference inherently leads to errors, especially if the fact is ignored.
With all of these problems, it is justifiable to evaluate the worthiness of using the Basic Geological Maps of SFRY 1:100.000 for landslide susceptibility maps. Despite all the aforementioned problems, there are no other geological maps at a larger scale, which cover the entire region. Especially if we consider evaluating landslide susceptibility for the whole of Croatia. Although GUs in it are chronostratigraphic, lithological units are differentiated in some instances. In many chronostratigraphic units, the lithological composition is similar enough for combining their mechanical properties. While that is not an optimal source of geological data input for landslide susceptibility, it is viable if all the above problems are addressed. A preferred method for addressing such issues is fieldwork, as demonstrated in this study. As that is not always feasible, remote sensing methods could be applied in an attempt to differentiate lithologies or refine geological boundaries, but that has to be evaluated for each case. A possibility of other, more suitable geological maps for problematic areas should also be investigated.
In summary therefore, due caution is needed in interpreting the landslide susceptibility of GUs in the studied area.
In the studied area, the Pliocene GU has by far the most landslides and the most area coverage. The landslide index is also highest for the Pliocene unit. It can be concluded that the Pliocene GU is the unit most susceptible to landslides in the given area.
Judging by the landslide index, which takes into account the number of landslides in relation to the overall area, the Palaeocene-Eocene GU has the second-highest landslide susceptibility. The Badenian GU has more landslides, but over a larger area, hence its landslide index is lower. It has to be noted that the number of landslides and area of both of these GUs is low, hence the reliability of this interpretation is also low.
All other GUs in the study area cover a small area, or have only a few or no landslides. Therefore, no conclusions in regards to landslide susceptibility can be made for those units. One of those GUs is Pl,Q which is known for high landslide susceptibility (JURAK et al., 1998), but has no landslides in the study area. Since it is represented by an area of only 0.35 km 2 , it is recommended that it be disregarded from landslide susceptibility evaluation.
Furthermore, due caution is advised with the Miocene units: 2 M 3 1 , 1 M 3 1,2 , 1 M 3 1 . With zero, one, and three landslides per unit respectively, they leave an impression of a relatively stable area. It has to be noted that the average slope angle in those GUs is around 10° (Table 1). Comparing that to the GUs with a larger number of landslides, the Badenian and Palaeocene-Eocene units are much steeper with approximately 20°. Only the Pliocene unit has a similar average slope, and the most landslides in the study area, which confirms it as the most landslide susceptible GU. Comparing slope angle in landslide bodies, for GUs with more than ten landslides, the Badenian units are the steepest, Palaeocene-Eocene unit somewhat less steep, and the Pliocene substantially less (Fig. 8). This reflects the mechanical properties of materials in the respective GUs, which are weakest in the Palaeocene unit, as landslides are triggered in less steep slopes. Slope in landslide bodies also correlates with the landslide index for those GUs, with a higher landslide index indicating weaker mechanical properties.
Concentrating further on the three GUs with more than ten landslides (Pliocene, Palaeocene-Eocene, and Badenian), it can be seen that landslides are not evenly spread over their respective areas (Fig. 9). However, the causes of this are unique for each GU.

Pc,E unit
The simplest case is in the Palaeocene-Eocene unit. Although it is one chronostratigraphic unit, it consists of two lithological units with different mechanical properties. The western part has only one landslide as it is composed of conglomerates that have a very thin regolith. The Eastern part is composed of various clastic rocks such as sandstones and siltstones. They have a significantly lower mechanical strength and durability, and therefore are likely to develop a thicker regolith with less sharp-edged fragments. Differentiating those two lithologies is difficult without fieldwork. Determining landslide susceptibility for the whole GU, without considering different lithologies, could provide substantial inaccuracies. Since there is no better source for geological data available for the whole of Croatia, discarding the Basic Geological Map would not be rational. However, landslide susceptibility should be limited to a small scale, reflecting the quality of the input data ( VAN WESTEN et al., 2008). If a large scale is needed, a further investigation that will differentiate lithological units is mandated.

M 2 2 unit
The Badenian unit also has two distinct lithological units within it. While the initial assumption for the uneven spread of landslides was that same lithological differentiation, fieldwork disproved this. Although limestones have better mechanical properties and durability than marls, sandstones, and poorly cemented conglomerates, they produce enough regolith with sufficient thickness for landslide formation. Despite this, there are not many landslides developed in the whole Badenian unit, except an area in the easternmost part of the unit (Fig. 14). Field investigation revealed the forest road as the main cause. Most of the landslides are positioned near the road or in contact with it. There is a clear connection with landslides and cutting into the slope or associ-ated with runoff water from the road surface. There are no records for the time of activation for these landslides. Since other factors including the proximity of the tectonic boundary with the Palaeocene-Eocene unit and spring occurrences are present, no definitive trigger can be identified. It could be that the combination of all the aforementioned factors in the same area is responsible for a large number of landslides. Examples like this are hard to predict, especially on smaller scales. Roads can be a significant factor in landslide activation, but done with proper drainage, they can be stabilizing factors in the slope. Differentiating those two cases is difficult with the available data.
While occurrences of groundwater and springs in connection with landslides are well established (ANBALAGAN, 1992;GETACHEW & METEN, 2021;GÖKCEOGLU & AKSOY, 1996), their relationship is not straightforward. The majority of larger springs have no connections with landslides. The most common, landslide-inducing springs, are very small. As such, they are unmarked on topographical maps. Locating them is time-consuming and not practical for larger areas. Locations of the springs in this example are predetermined by geological setting. Specifically, contact of the permeable Badenian limestones with non-permeable Palaeocene-Eocene deposits. While determining groundwater level is crucial for any single landslide investigation, establishing it for a wide area is difficult and often unattainable.

Pl 2,3 unit
The Pliocene GU is the dominant unit in the study area, in both the number of landslides and area coverage. It also exhibits contrasting landslide densities (Fig. 9) and again, it is hard to resolve the primary reason for it. Figure 10. shows a tendency for a northeast landslide orientation. That corresponds to the general orientation of the slopes inclined to the Petrinjčica river valley that is the lowest in the area. In contrast, to the southwest side of the Pliocene unit, the Utinja stream flows which is approximately 55 m higher. Consequently, the northeastern slopes have much higher relief energy than the southwestern slopes.
Surprisingly, the aspect of all 5 x 5 m cells does not follow that trend. North and south aspects are prevalent over east and west orientations. Since the aspect disregards the angle of the slope, their orientation was also analyzed by projecting plane poles on a stereographic plot (Fig. 11). The main vector was calculated for all slopes and slopes within landslides. In the case of all slopes, this has confirmed the relatively even slope distribution, with only a slight tendency towards the north -northeast. In the case of landslides, the mean pole vector indicates the northeast plane aspect, although not as clear as in Fig. 10. Besides relief energy, the prevailing factor indicating the northeast landslide orientation could be due to the geological setting. The Pliocene unit tends to dip towards the northeast which could favour that direction for landslides. This is confirmed by only a small number of measured bedding planes as the area is covered, so there are only a few outcrops, and even then, bedding planes are not always visible. However, the whole regional setting of the Pliocene unit indicates the northeast aspect. Quantifying the influence of relief energy and bedding plane orientation is difficult, but they are the most probable culprits for uneven landslide distribution in the Pliocene unit.

Landslide susceptibility factors and scale
Landslide susceptibility can be used for different purposes and at different scales, and therefore many methods and factors have been developed for determining it (CASCINI, 2008;FELL et al., 2008;REICHENBACH et al., 2018;SALEEM et al., 2019;ZHOU et al., 2016). Despite all this effort, there is no universally accepted standard method for any scale or purpose. This study does not attempt to resolve this, as the study area is too insignificant in size or factor variability. Despite that, some conclusions can be made that may serve as guidance in a future attempt at standard method development, at least in the region where this study was conducted.
The results of this study imply a scale of 1:100.000 as the largest scale for a more simplistic approach. There is a Basic Geological Map available at that scale for the whole country of Croatia. As stated before, chronostratigraphic units used in that map are a problem, but would be acceptable in most cases. However, that should be tested for each case. GUs are one of the most important factors, as seen with the help of the landslide index (Fig. 7). The example of the dpr GU reminds us of a need for caution as there can be problems related to this data. Maps in that scale, or smaller, can be compiled in less time and with fewer financial resources.
Susceptibility maps at larger scales can lead to significant misconceptions if a more detailed approach is not undertaken. Examples from three GUs in this study demonstrate this. The Pliocene, Miocene, and Palaeocene-Eocene GUs have very contrasting areas within them regarding landslide densities. When the landslides are averaged over the whole area of the corresponding GU, each GU can be assessed in landslide susceptibility with enough precision for scales of 1:100.000 or smaller. In contrast, the same approach at larger scales would introduce a lot of inaccuracies. Furthermore, every GU had unique factors that caused drastically different landslide susceptibilities. Considering the modest area of this study and such heterogeneity, it is hard to hypothesize a universal method for determining landslide susceptibility at large scales, which would apply to large regions with a lot more GUs. This implies that scales larger than 1:100.000 require much more extensive data. In the region of this study, that involves extensive fieldwork as such data is not readily available remotely.

CONCLUSION
A landslide inventory was created for an area of 22.6 km 2 near Petrinja city. LiDAR scanning and orthophoto maps were used for producing the inventory. A total of 216 landslides were detected, which cover 0.66 square km or 2.91 % of the total area. LiDAR data proved an exceptional data source as it provides high precision and accuracy of the ground surface, even in forested areas. The most notable difficulties concerning creating the landslide inventory presented in this study are subjectivity in demarcating the landslide outlines and the unknown times of landslide activation. Grading landslides based on the discernibility of their features in the LiDAR and orthophoto data may indicate a relative landslide age, but it is too ambiguous for reliable conclusions to be drawn.
The geology of the study area was described based on a Basic Geological Map of SFRY at a scale of 1:100.000. On that foundation 11 GUs were identified in the area. Only three GUs had more than ten landslide bodies within them, so a more detailed analysis was conducted on them. The vast majority of landslides were discovered in the Pliocene GU with 182 landslides or 84%. It is followed by the Badenian unit with 16 landslides and the Palaeocene-Eocene unit with 13 landslides. Those numbers were normalized based on the total surface area of each GU and were expressed as a landslide index. The Pliocene GU has the largest landslide susceptibility with a landslide index of 6.7. It is followed by the Palaeocene-Eocene unit with 4.2 and finally the Badenian unit with 0.9. The average landslide index for the whole area was 2.9. The Holocene Diluvial-proluvial unit (dpr) has a large landslide index of 1.4, but that is based on only one large landslide situated on the border with the Pliocene GU, so it is dismissed as unreliable for further conclusions. Slope angle in landslide bo dies also confirms the Pliocene GU as the most susceptible, followed by the Palaeocene-Eocene unit and finally the Badenian unit. Other GUs present in the area should not be regarded as not susceptible to landslides, based only on this study. Their represented surface area and/or other risk factors are low (slope, local relief). That excludes them from any reliable conclusions on landslide susceptibility in other regions. Only Holocene alluvial sediments (a, a 1 , ap) could be regarded as not landslide susceptible as they are found only in flat areas. An exception to that rule would be landslides along the river banks that are not present in the study area.
All three investigated GUs have an uneven distribution of landslides in their respective areas. In all three cases, the determining factors are different. In the case of the Pliocene GU, there is an indication that the dominant factors for larger landslide susceptibility are higher local relief and unfavourable bedding plane orientation on affected slopes. Different lithologies is the reason in the Palaeocene-Eocene GU, with conglomerates having very low susceptibility and sandstones, siltstones, and sandy marls having much larger susceptibility. Human impact in the form of a road is the most probable factor in the Badeninan GU, but it cannot be proven without the known activation times of the landslides.
If the scale of 1:100.000 is considered, Basic Geological Maps of SFRY are a good enough source for geological data input for landslide susceptibility maps, but with certain limitations. In the case of lithologically heterogeneous GUs, a more detailed approach is needed. For more detailed scales, Basic Geological Maps of SFRY are not adequate by themselves.
Based on the results of this study, landslide susceptibility maps at scales larger than 1:100.000 demand more individual approaches in determining the most influencing factors.