Landslide inventory and characteristics, based on LiDAR scanning and optimised field investigations in the Kutina area, Croatia

This paper presents the preliminary results of analyses of landsliding processes derived from detailed LiDAR (Light Detection and Ranging) scans supported by field prospection on the south western slopes of Mt. Moslavačka gora, in the wider Kutina area. This area is known for frequent landslides, but dedicated regional landslide research has not been previously undertaken. High resolution LiDAR scanning and orthophoto imaging enabled the production of a reliable landslide inventory, but also enabled research on landslide properties and the morphology of the area. Field mapping and prospection, sampling and borehole coring assisted in the collection of information about the material characteristics and specific features of typical landslides. In the research area, which covers more than 71 km 2 , more than 1200 very small landslides were detected. The majority of landslides were discovered in just several geological units indicating their high susceptibility: Pleistocene silts and sands with clayey interlayers, followed by M 2 silty sands and gravels, and M 7 sands. Nearly half of the landslides are estimated to be of recent and younger age, while other landslides may be considered as being historical implying a “long tra-dition” of landslide events in the research area. Preliminary terrain surface roughness analysis also supported the conclusion that the inventory contains landslides of several historical generations which are still detectable. In addition to slides (1123), this research also discovered numerous earthflow processes (143), which are more frequent in the predominantly sandy units. The landslides in this area are largely located on the banks of the gullies and are directly related to the action of water. Regarding that situation and the engineering properties of the encountered geological units, four types of bank instabilities can be differentiated: slides on top of rock mas­ ses; slides in firm soil mixtures; landslides in sands; landslides in predominantly coherent soil complexes.


INTRODUCTION
Landslide is a general term used for the gravitational movements of the ground mass, whether it involves soil, rock or various combinations of both, as well as the landforms resulting from such movements (HIGHLAND & BOBROWSKY, 2008). Their occurrence is closely related to conditioning factors that, according to numerous existing examples of literature, frequently include: lithology, hydrological characteristics, the presence of older landslides, slope angle, vegetation and land use. Landslides are usually initiated by: ample precipitation, flash floods, human intervention, and/or seismic activity (ZÊZERE et al., 1999;MAHALINGAM et al., 2016). Therefore, identification of the specific factors that directly influence slope instability in certain regions is key for the prediction of such processes and should be used for landslide susceptibility mapping and assessment of potential hazards.
In order to detect conditions in which landslides are activated, it is essential to create a reliable landslide inventory. The conventional approach to build a landslide inventory comprises a combination of stereoscopic aerial photography and field prospection, which hasn't always proved successful, given its dependence on various factors such as the accessibility of the terrain or vegetation cover. However, the efficiency of the investigation of barely accessible areas has been enhanced with the rapid development of one remote sensing technique, Light Detection And Ranging land use risk through landslide susceptibility maps design), an Interreg IPA CBC project. The research area that covers 71 km 2 in the wider Kutina area is known for frequent landslides which locally endangers people, their properties and infrastructure. Previous research of landslides in this area were geotechnical investigations for landslide mitigation. Regional landslide studies have not previously been conducted, but few landslide locations were known from recent landslide reports (Report Landslide portal https://www.hgi-cgs.hr/prijava-klizista/) and historical landslides located on the existing Engineering Geological Map 1:500 000 (ČUBRILOVIĆ et al., 1967).
For all the aforementioned reasons, there are several objectives of this paper: to define and describe the workflow for the creation of a landslide inventory based on detailed LiDAR scanning; to present the landslide inventory and results of its study; to document the added value of complemented field and remote sensing research in landslide studies.
As for most of the Croatian examples, landslides in this area are very shallow and very small, according to FELL's (1994) classification (500-5000 m 3 ). In addition, most of the landslides are located in forested areas. The advantages of the use of airborne LiDAR data in creating a landslide inventory in forest areas are also commented on by others (GÖRÜM, 2019; Van den EECK-HAUT et al., 2007). Therefore, the LiDAR technique was used to enable the production of reliable HRDEM. The results of Li-DAR scanning facilitated the optimisation of field explorations that included: geological mapping, sampling, borehole coring, in-situ testing and landslide inventory verification.
With the use of raster maps derived from the initial HRDEM (hillshade, slope and contour map) a visual analysis was performed within Esri ArcGIS 10.0 software, for the development of the landslide inventory.
The landslide analysis has been performed on two levels. The first concerned a general analysis of landslide processes for the complete research area scanned by LiDAR, while the second entailed more detailed analysis and field exploration in two pilot areas, Kutina North and Kutina South. For the purpose of the de-tailed investigation of pilot areas, geological mapping at a scale of 1:5 000 and engineering-geological prospection were also performed.

GEOLOGICAL SETTING
The total research area covers more than 71 km 2 and is located on the southwest slopes of Mt. Moslavačka gora situated between the city of Kutina and the G. Jelenska settlement.
Considering the geological history, the whole area belongs to the Pannonian Basin System (PBS) which represents one of the Mediterranean back-arc basins the formation of which commenced in the Early and Middle Miocene, due to the collision of the African (Apulian) and European plates. The PBS is surrounded by the Alps, Carpathians, and Dinarides and includes several differently sized, deep sub-basins separated by a comparatively shallow complex of basement rocks (HORVÁTH & ROYDEN, 1981;ROYDEN, 1988) (Fig. 1).
Palaeogeographically, the PBS belongs to the area and bioprovince of Central Paratethys. During the Miocene a connection between Central Paratethys and the Mediterranean and the Indo-Pacific Ocean was repeatedly lost and re-established (STEIN-INGER et al., 1988, RÖGL;1996;HARZHAUSER & PILLER, 2007;KOVÁČ et al., 2018). The development of the basin took place in two phases. The first (syn-rift) phase was characterized by tectonic thinning of the crust and isostatic subsidence, while the second (post-rift) phase was marked by subsidence caused by the cooling of the lithosphere (HORVÁTH & ROYDEN, 1981;ROYDEN et al., 1983;ROYDEN, 1988). In the Croatian part of the PBS, the syn-rift phase lasted from the Ottnangian to the Middle Badenian, while the post-rift phase extended from the Late Badenian to the end of the Quaternary (PAVELIĆ & KOVAČIĆ, 2018 and references therein).
The existing geological map of the investigated area, at a scale of 1:100 000 (CRNKO, 2014), is based on chronostratigraphic criteria (Fig. 2). This map enabled understanding of the regional geology and the general engineering geological properties for the whole LiDAR scanned area of 71 km 2 . For the purpose Figure 1. A tectonic sketch of the Pannonian Basin System and its surroundings (after ROYDEN, 1988) showing the location of the research area. of this research (under the safEarth project), more detailed geological mapping at a scale of 1: 5000 was completed for two pilot areas (Kutina North and Kutina South; Fig. 3) in order to define the geological units based on lithofacies criteria. Besides the more precise outlining of geological boundaries, these maps allowed further differentiation of some lithostratigraphic units, which is needed to determine associations between material properties and landsliding processes at specific locations.
Because of the different scale and different attitude in these geological maps, some units have different symbols for the same Figure 2. The geological map of the research area, originally at a scale of 1:100.000 (according to: CRNKO, 2014). The map also displays the inventoried landslides and two pilot polygons. geological material (Figs. 2 and 3). This is the reason why some of the described units in the following text are represented with two symbols. The first symbol represents the chronostratigraphic unit from the 1:100.000 map and the second symbol represents the lithostratigraphic unit from the 1:5.000 map.
Mt. Moslavačka gora represents one of the main surface exposures of crystalline basement within the Neogene sediments of the PBS (PAMIĆ, 1990;PAMIĆ et al., 2002). According to KOROLIJA et al. (1986) and CRNKO & VRAGOVIĆ (2014) the Mt. Moslavačka gora is interpreted as a horst structure formed during the Miocene by vertical faulting. The largest part of the crystal basement consists of granite rocks, while metamorphic rocks of medium to a high degree of metamorphism are less represented. The granite-migmatite-metamorphic complex of Moslavačka gora is in tectonic contact with Neogene deposits, and only a small part of their transgressive contact has been preserved (CRNKO & VRAGOVIĆ, 2014;CRNKO, 2014). Granites of Mt. Moslavačka gora contain K-feldspar, oligoclase, quartz and biotite with secondary muscovite (PAMIĆ et al., 1984). Moslavačka gora granite, besides the usual granite ingredients, contains andalusite and sillimanite. According to PAMIĆ et al. (1984) on the north-western part of Mt. Moslavačka gora outcrops of fresh granite are rarely found. The structure of granites is finegrained (rarely medium to coarse-grained). Zones of grusification of granite appear on the surface and they are several metres thick. The metamorphic complex, which is preserved only in the form of smaller or larger enclaves within granite, consists of homogeneous and heterogeneous migmatites (Fig. 4), different types of amphibolites and amphibolite schists and paragneiss (CRNKO & VRAGOVIĆ, 2014) (Magmatic-Metamorphic Complex; MMK). According to PAMIĆ (1990) migmatites were developed gradually from the highest-grade amphibolite facies rocks and are represented by numerous varieties, but metatexites and stromatites are very common.
In the investigated area, Neogene deposits have a limited distribution on the surface. The oldest Miocene deposits are freshwater Ottnangian sediments (Daranovac fm. -M 2 ; Dar) discovered in the southeastern and southwestern part of Mt. Moslavačka gora.
The Daranovac formation. (M 2 ) -disconformably (transgressive or fault contact) overlies the crystalline basement. The lower part of the Daranovac fm. comprises a few metres of weakly lithified rock-fall breccia, conglomerates and gravels (Fig. 5a). At some places the sand directly covers the crystalline basement. The Middle and upper part of the Daranovac fm. is represented by coarse-grained deposits intercalated with sandy and structureless silty units that may contain thin conglomerate lenses. The coarse clastic parts of the sediment in the lower part of the formation attain the size of blocks and boulders, and originate from the rocks of the immediate basement (granite, gneiss, amphibolites, pegmatites, and rarely gabbro and quartzite) (CRNKO & VRAGOVIĆ, 2014). These deposits in the neighbouring PBS areas are interpreted by PAVELIĆ & KOVAČIĆ (1999) and PAVELIĆ (2001) as alluvial fan deposits.
Carpathian deposits (the Glavnica fm. -M 3 ) are mainly in tectonic contact with the Daranovac fm. and younger deposits. These deposits are poorly represented in the research area. They consist of sandy and silty marls, calcareous siltstones, marls and, rarely tuffitic sandstones and lenses of fine-grained conglomerates (CRNKO & VRAGOVIĆ, 2014).
Badenian deposits (Vra-Vej) lie transgressively over, or are in tectonic contact with a weathered crystalline basement or lower Miocene deposits. The main feature of Badenian sediments is a very dynamic vertical and lateral exchange of various lithified to unlithified, carbonate-free and calcareous clastic sediments. These are marine deposits consisting of conglomerates, weakly lithified, coarse-grained sandstones with locally developed bioaccumulated and clayey limestones (Vrapče fm. -M 4 ). Since volcanic activity continued from the Ottnangian to the Badenian, within the Vrapče fm. tuffites were also registered (CRNKO & VRAGOVIĆ, 2014). Structureless marls, that represent the Vejalnica fm. are the predominant lithological member of the middle and upper Badenian. This unit is mostly composed of thick marl beds with rare intercalations of coarse-grained clastics which were transported into the basin by gravitational flows, reaching the outer shelf (PAVELIĆ et al., 1998;ĆORIĆ et al., 2009).
Sarmatian sediments (the Dolje and Pećinka formations -M 5 ) occur on the southeastern and southwestern slopes of Mt. Moslavačka gora in several smaller separate areas (CRNKO & VRAGOVIĆ, 2014). The Sarmatian sediments are mainly in tectonic contact with older rocks. They consist of sandy and silty marls, sandstones, consolidated gravelly sandstones and limestones that were deposited in a marine environment with reduced salinity. The lower Pannonian deposits are present only on the southeastern and southern slopes of Mt. Moslavačka gora (CRNKO & VRAGOVIĆ, 2014). They are developed as thinbedded clayey limestones and marls (Croatica beds or Croatica fm. -M 6 1 ; Cro) and structureless marls to sandy marls of the Banatica beds (Medvedski Breg fm. -M 6 2 ; MeB) in the brackish, Lake Pannon. Unlike the Croatica fm., which indicates deposition in a shallow-water lake environment with low salinity, massive marls of the Medvedski Breg fm. suggest a deep-water brackish depositional environment (KOVAČIĆ & GRIZELJ, 2006). The Lower Panonian deposits are conformably overlain by the marly -sandy Abichi beds of the Upper Pannonian (Andraševec fm. -M 7 1 ). The Andraševec fm. was deposited in a prodelta -delta slope, lacustrine environment (KOVAČIĆ & GRIZELJ, 2006;KOVAČIĆ et al., 2004). Distribution of the Andraševec fm. in this region is limited to three smaller separate areas (CRNKO & VRAGOVIĆ, 2014). From the lower to the upper Pannonian, sedimentation was continuous, with more expressed lowering of certain parts of the horst structure of Mt. Moslavačka Gora. This explains why the Upper Pannonian deposits are more widespread on the surface, than lower units. The Upper Pannonian or "Rhomboidea beds" (Nova Gradiška fm. -M 7 2 ; NGr) are represented by multicoloured, coarse-grained to fine-grained sands, with a transition to silt (CRNKO & VRAGOVIĆ, 2014) (M 7 2 ; NGr) (Fig. 5b). The sands mostly contain micas (illite-muscovite) and in some cases clay minerals. According to their granulometric composition, the samples of the Nova Gradiška fm. belong to the wellsorted silty sands and sandy silts. Medium-grained to coarsegrained sands are rarely present. They are usually well sorted, and contain medium to well-rounded grains. In some cases, they contain significant quantity of clay minerals. According to KOVAČIĆ & GRIZELJ (2006), KOVAČIĆ et al. (2004) the deposits of the Nova Gradiška fm. have been deposited in a delta front, shallow brackish, lacustrine environment.
The sediments of the Middle and Upper Pliocene conformably overlie the Upper Pannonian deposits. They are known as the "Paludina beds" (Vrbova fm. -Pl; Vrv) and have the widest distribution of all Neogene formations of this area (CRNKO & VRAGOVIĆ, 2014). Fine-grained to coarse-grained sands with gravel lenses or intercalations are dominant (Vrv-s) (Fig. 5d). Within these sands, there are interlayers of silt and/or clayey silt several metres thick (Vrv-c) (Fig. 5c). Contacts between the sands and silts/clayey silts are gradual or sometimes sharp. The deposits are structureless with the very rare occurrence of layers. Stratification is manifested through sudden changes in the granulometric composition and especially lamination in the sands and silts. In the middle and upper part of the Vrbova fm., within the sands and silts, horizontal lamination is often well expressed (CRNKO & VRAGOVIĆ, 2014). Cross-and trough cross-bedding and lamination can also be observed. The general characteristics of stratification are the same for all levels of the Vrbova fm. According to KUREČIĆ (2017), Pliocene sediments were deposited in an alluvial environment or in a freshwater lake.
In the Pleistocene, loess or loess-like sediments (snl) covered all older formations. During the Holocene deluvial-proluvial sediments (dpr) formed as a result of surface leaching and torrents from the hills of Mt. Moslavačka gora. These are mostly unlithified or weakly lithified sediments such as sands, silts, marls and loess, which have been subjected to relatively easy surface leaching. Recent alluvial sediments have been formed in the valleys of the younger river and stream courses (a,ap) (CRNKO & VRAGOVIĆ, 2014) .

METHODS
The pilot areas were selected according to the lithological diversity in order to define conditions for the activation of landslides in various geological materials. Field research in these areas defined the engineering-geological properties and processes on the slopes of different geological units with the idea of applying all findings to the wider area. This also enabled a better understanding of the material properties, weathering zones and specific behaviour for particular lithostratigraphic units. Other than coring, sampling and performing in-situ measurements, field verification of the landslide inventory has been performed within the defined pilot areas. Analysis of the surface roughness for natural slopes and landslide areas was also performed to quantify the magnitude of disruption caused by landsliding. The roughness data also enabled discussion about relative landslide age estimation.

INVENTORY
Detailed aerial LiDAR scanning (20 dots per m 2 , and accuracy ±10 cm) of the research area, which covers more than 71 km 2 , enabled production of a high-resolution digital elevation model (HRDEM; 0,5x0,5 m). The analysis and visual inspection of HR-DEM derivates (hillshade, slope and contour map) and detailed orthophoto (10 cm resolution) enabled remote indication of the features corresponding to the mass movements, and supported creation of the landslide inventory containing more than 1200 objects (Fig. 6).
The effectiveness of the visual analysis of LiDAR DEM derivates in building landslide inventories is proved by numerous authors (PETSCHKO et al., 2014). Since not all landslides are recognised with the same level of certainty, in this research, landslides were classified according to the confidence in their identification. Allocation to the corresponding landslide identification confidence levels is based on the visibility and persistence of landslide boundaries, landslide body features and their magnitude, the freshness of terrain marks and deformations, and modifications of the original vegetation (Fig. 7). The procedure is highly subjective, but is based on simple and clear criteria as   shown in Figure 7. The degree of certainty of identification of the slide boundaries has been elaborated by WIECZOREK (1984). Similarly, Van Den EECKHAUT (2007) has differentiated three classes of slides according to freshness and preservation of the typical landslide features visible on LiDAR based derivates. Therefore, it is believed that modified criteria suggested in Figure  7 can attain a high degree of congruence between the different experts for future research projects.
While generating the inventory for the research area, different surficial processes were identified: sliding, flowing, creeping, rockfalling and excessive erosion. The processes were differentiated according to the morphological features, shape of the polygons, appearance of polygon edges and orthophoto images.
In summary, the landslide inventory provided, in addition to the precisely positioned contours, a database with classified information about: the type of process, visibility and expression of characteristic landslide features, relation to nearest streams or other watercourses and landslide identification confidence level. The combination of these data and that gathered by field explorations gave us the opportunity to gain a more complete understanding about landsliding processes and their relationship to geological and other influential factors (ZÊZERE et al., 1999).

FIELD INVESTIGATION
Field investigations in two selected pilot areas comprised mapping, sampling, borehole coring, in-situ testing and LiDAR based landslide inventory verification.
Field mapping included both geological and engineering geological mapping at a scale of 1:5000 during which characteristic landslides and soil profiles were recognised in more detail, and their data classified for databases. The rarely exposed outcrops were used for profiling, sampling and probing with a pocket pen-etrometer. Another purpose of the field prospection was the verification of LiDAR based landslide inventory.
The pocket penetrometer test method is used to rapidly evaluate the consistency and approximate unconfined compressive strength of soils in-situ. An ELE International pocket penetrometer with a 6 mm piston was used, where readings are directly expressed as kg/cm 2 . As there is no standard on how to use a pocket penetrometer or evaluate the data, a special procedure was applied here which enables the relative comparison of most of the materials of concern. To enable the relative comparison of 'harder' soils as well, readings slightly over 4.5 were registered as >4.5 kg/cm 2 , while the cases where no penetration was achieved as >>4.5 kg/cm 2 . This enabled the identification of harder sections in soil mixtures.
Characteristic soil profiles and samples of each unit were also acquired by shallow core drilling (Fig. 8).

ROUGHNESS
The terrain surface roughness analysis is based on a LiDAR-derived DEM with a resolution of 0,5x0,5 m. For the preliminary roughness research, a simple calculation of standard deviation of slope (SDS) was used. This method is universally applicable for geomorphological and landslide analysis because it correctly identifies the local variability of the surface, breaks-of-slope and smooth sloping areas (GROHMAN et al., 2011;GARRISS, 2019). The standard deviation of the slope is also successfully used for research in different scales (GROHMAN et al., 2011;ZHANG et al., 1999).
The SDS calculation was performed in a 3x3 m moving window, using the Focal Statistics tool in ArcGIS. Window size is dictated primarily by the intention to quantify the variability of the terrain surface in very small landslides, but also to identify characteristic morphology of the unaffected terrain. The calculation returns a dimensionless index that is representative of an intrinsic property of the surface, invariant with respect to rotation or translation (BERTI et al., 2013).

RESULTS
The geological map, originally created at a scale of 1:100 000 (CRNKO, 2014), displays the spatial distribution of geological units which were previously described in Materials and Methods (Fig. 2). It is obvious that the northern segment of the research area is composed of igneous and metamorphic rocks, mostly migmatites and granites. Younger, Neogene sediments of various origin and facies surround crystalline basement and constitute most of the central part of the research area. Southern and western edges of the research area are dominated by Pliocene sands or silty clays and deluvial sediments. Loess sediments irregularly cover older bases.
In order to learn about the engineering behaviour of the major geological constituents in the wider area, two pilot areas for field explorations at a scale of 1:5000 were located (Fig. 3). The field explorations enabled insight into many important additional parameters and features which can't be perceived from remote sensing, and aren't known from previous surveys. The major results of such a combined approach, analysing remote and in-situ data in the pilot areas are presented below. It is important to point out that outside the pilot areas, only remote sensing and previous geological data were analysed. Figure 8. In-situ coring with the Atlas Copco Cobra system and corresponding core. The recovered cores are 100 mm in diameter.

SOIL PROPERTIES
During field mapping and borehole core determination, characteristic shallow profiles for each unit were defined. The field soil classification for each layer enabled insight into the predominant granulometry for a particular unit (Fig. 9). Figure 9 documents the overall predominance of silty and sandy materials. Silty materials primarily constitute the weathering zones of marly beds (Vra-Vej) or loess (snl), but are also important constituents of soil mixtures (al, Vrv, Dar). Sands are dominant in the NGr sands and Vrv-s, Dar soil mixtures. Although clayey material is present in most of the units, the clay fraction is dominant only in some layers of the dpr, Vrv-c and Vrv-s. Gravelly interbeds are noted in the MMK, Dar, Vra-Vej but only very rarely in the Vrv-s.
According to GRIZELJ et al. (2017), carbonate minerals, clay minerals and quartz are the main constituents of all Miocene pelitic sedimentary rocks of the Croatian part of the PBS. Exceptions include the analysed Badenian and late Pannonian sediments from Moslavačka gora Mt. which do not contain carbonate minerals. Among clay minerals, swelling minerals (smectite or illite-smectite) are the main constituents, while detrital illite and kaolinite are present in small quantities.
Field consistency estimation defines most of the coherent soils from the different units as firm (Dar, dpr, MMK, snl), while the weathered zone above the marls and carbonate bedrock Vra-Vej is classified from firm to stiff (according to ISO 14688-1:2002ISO 14688-1: , 2002. Vrv units have an even wider consistency span, from soft to stiff (Fig. 10), which is a consequence of different granulometry of the interlayers.
The penetrometer testing was performed on both the exposed outcrops and the borehole cores. Altogether 620 readings were collected for most of the units encountered during the field research (Tab. 1). Average readings for all units are over 1 kg/ cm 2 . The coherent materials can be considered as medium stiff to very stiff, or even hard (Tab. 1; according to the Texas Department of Transportation, 1999). A lower consistency is recognised for the Vrv-c, Vrv-s and MMK units but these units also have a higher percentage (> 20 %) of readings outwith the instrument scale (Tab. 1). This is congruent with the field estimations because it indicates the heterogeneity of the material and documents the presence of both very stiff and incoherent interlayers. The Al, dpr and snl units are considered as firm but do not have a lot of 'out of scale' readings which could be the outcome of their material and genesis resulting in relatively homogenous properties. Dar sediments are also firm, but have a lot of sections where the instrument could not penetrate, indicating very stiff and coarse material. NGr sands are stiff, which was expected since they lack coherent components. Segments of heterogenous weathered zone of Vra-Vej are also stiff, but more than 65 % of readings are performed in very stiff material.
The testing depth is mostly between 30 cm and 2.5 m, rarely up to 6 m. According to penetrometer probing at various depths, two general trends can be recognised. NGr, snl, al and partly Dar have a clear improvement of material mechanical properties with increasing depth of the sampling. This could be a very important factor influencing the depth of the landslide plane. In other units there is no clear trend of material 'improvement' related to increase in depth, therefore the presented averages could be considered as representative of the whole characteristic profile.

LANDSLIDES
The landslide inventory for the whole research area contains 1267 entries of which 359 were validated during fieldwork. Among all the registered polygons, only one is classified as rockfall. Therefore this study is fully dedicated to the investigation of slides (1123) and earthflows (143), which are considered as the dominant slope gravitational processes for the area.
The results of landslide cover identification show that the majority of field verified landslides are found in forests and woodland (62,11 %), with a lesser percentage in shrubs and meadows (33,68 %), while the lowest incidence has been associated with arable land (4,21 %). Given the fact that this data is obtained for less than 10 % of the landslide inventory, estimation of the correlation of landslide occurrence with land use was not performed.
The average slide area is 995 m 2 , which is considered as very small regarding both affected area and depleted mass volume (500-5000 m 3 ; FELL, 1994). The average earthflow area is considerably larger at 2 866 m 2 . The estimated depth of the landslide slip surface is mostly 1-5 m for the Dar and Vrv sediments, but is even shallower for the MMK unit (Fig. 11). Dar sediments also have a larger proportion of very shallow landslides, which is not the case for the Vrv formation.

Landslide susceptibility
The proportion of particular units subjected to landsliding can indicate landslide susceptibility. As for some other parameters, analysis of the unit's susceptibility is done separately for the pilot (detailed mapping) and total research area.
The analysis of the entire inventory suggests that Pl sediments (in pilot areas Vrv), are by far the most susceptible to sliding and flowing processes (Tab. 2). Namely, more than 5 % of the area corresponding to the Pl unit is affected by slides and more than 3 % by earthflows. It can be noted that numerous earthflow processes have been identified in the Pl unit. During detailed geological mapping at a scale 1:5000, the same unit has been divided into areas with a predominant sandstone content (Vrv-s) and the clayey layers or lenses (Vrv-c). Comparison of the two Pl components shows that the Vrv-c, or clayey layers, are more susceptible to flowing processes (Tab. 2). Therefore, distinction of two lithostratigraphic segments within the Pl unit enabled a better understanding of the behaviour of this complex deposit. Particular Miocene deposits, M 2 (Dar) and M 7 (NGr), can also be considered as susceptible to landsliding given that more than 1.5% of their area is subjected to sliding. The same units also have a noticeable percentage of earthflow processes (around 0.5 %) which is expected with regard to the predomination of sandy fractions.
The susceptibility to landsliding of other units is mostly much lower. Units A, b, M 5 and v are poorly represented in the research area and, throughout this investigation, no evidence of landsliding processes in these units was detected. This however doesn't mean that these units are not susceptible to landsliding.

Landslide field verification and identification confidence
As presented above, the landslide inventory is based on the visual inspection of remote sensing data.
Field verification in the pilot area confirmed almost 300 landslides and rejected nearly 50. Most of the rejected 'landslides' were initially detected with very low confidence, but after field prospection were disregarded, due to a lack of typical landslide elements. Some of the rejected 'landslides' are the consequence of excessive erosion processes on the steep gully banks, which have a similar morphological appearance to landslides. Field inspection also allowed recognition of landslide generations or initiated merging of previously detached landslide polygons. The corrected and improved inventory is subjected to analysis in this paper. Figure 11. The estimated depths of 71 landslides from the pilot areas. The methodology for landslide identification confidence is defined as described above with results of the evaluation presented in Figure 12. It is evident that most of the inventoried landslides in the research area are identified with low confidence. The reason for lower visibility and clarity of the landslide elements in a number of cases can be twofold: 'older' landslides have masked contours and settled body elements due to the effects of time, or very shallow and slow landslides with gentle and small defects. Well-developed and easily recognisable landslides are found in units that also are highly susceptible to landsliding (Fig. 12), including the Pl (Vrv), M 7 (NGr) and M 2 (Dar) deposits. Although there are significantly less landslides in the weathering zones of loess (l) and metamorphic bedrocks (Mi), these units also have "fresh" landslides of high or medium confidence levels.

Landslide spatial distribution
It is noticed that the majority of the inventoried landslides are located above streams or other water courses and usually form grapelike clusters (Fig. 13). Namely, the research area is characterised with a dense network of gullies, channels and valleys (Fig.  3, Fig. 13). The banks are usually the steepest segments of the terrain. The active incision of the water course into the soil material starts sapping the banks which initiates landsliding.
Such slides are dominant over the whole research area, and are most frequent in all units, excluding the Dar deposits (Fig.  14). This could indicate that the main triggering factor in almost all lithographic units is the action of water. Higher precipitation causes an increase in water flow and even flash floods at the base of the gullies, which induces erosion of the foot of the slope and initiates landsliding. Also, more permeable sandy materials which are frequent in the research area (Pl (Vrv); M 2 (Dar); M 7 (NGr)) are more easily saturated with water, which may initiate earthflows.

DISCUSSION
The presented results demonstrate that the predominantly clayey, silty or sandy lithofacies units (Vrv, Dar, NGr) are prone to landsliding (Tab. 2). Weathering zones of hard bedrocks (igneous/ metamorphic rocks, marls or limestones) may have similar granulometry, but are less susceptible to sliding (MMK, MeB, snl, Vra-Vej) due to their different mineralogical content and limited thickness of the weathering zones. Alluvial materials also have unfavourable granulometry and mineral content, however these sediments are mostly found on mild slopes or planes and therefore have a lower number of active processes. The susceptibility of other units which are poorly represented in the research area is unknown (A, b, M 5 , v).

LANDSLIDE TYPE DIFFERENTIATION
Two types of landslides are differentiated according to the failure and transport mechanism which influence their present morphological features. Sliding could be described as a mass movement with a distinct zone of weakness (failure surface), that separates the displaced mass from more stable underlying material (VAR-NES, 1978). The displaced material in slides may move along a failure surface with little internal deformation. The flow is de-  fined as a continuous movement in which failure surfaces are usually not preserved (CRUDEN & VARNES, 1996). The displaced mass undergoes liquefaction and runs or flows as a viscose liquid. The gradation from slides to flows depend on water content and material mobility.
During the investigation of the area, the slides and earthflows were differentiated according to specific morphological features and shapes of the polygons. The slides are usually recognised as relatively compact, with more regular contours than earthflows. From the orthogonal viewpoint, the main body of a slide is usually oval or slightly elongated, but also irregular or fun-shaped (Fig. 3, Fig. 15). The majority of the slides are symmetrical, with a characteristic arc-shaped main scarp and toe. In some cases, secondary scarps and other cracks are also visible.
Slides mainly occur on the banks of the gullies and the triggered material moves perpendicular to the direction of the stream (Fig. 15). Because of that, the slide foots are frequently eroded by the streams. The earthflows usually follow water channels, gullies or valleys (Fig. 15). Therefore, transportation of depleted ma-terial is parallel to the local depression and stream, sometimes causing skewed or even cornered contours (Fig. 15).
Earthflows are usually elongated and have a characteristic "hourglass" shape (Fig. 3). Their main scarp contour is frequently irregular or jagged (Fig. 15). In rare cases the main scarp is arcshaped similar to the vast majority of the slides. In an ideal case, where there is enough room in the foot plane for colluvium material, the earthflow toe has a fan-like shape. In other cases, the toe geometry is dictated by the local topography in the zone of accumulation.

LANDSLIDING IN GULLY BANKS
The research area is characterised by a relatively dense network of gullies. Over geological time, the drainage canals are cut into the bedrock forming various types of gullies and valleys. Deepening and widening of the gullies together with the action of water causes high pore pressure, liquefaction and seepage forces, all of which can initiate landsliding in their banks (RECKENDORF, 2009). Naturally, these processes are accelerated by high precipitation or extremes because of changes at the watershed and/or reach level.
The bank instabilities are also influenced by the size and morphology of gullies and physical and mechanical properties of the rock material forming the banks. Considering these facts, four different models of landsliding in gully banks can be differentiated: 1. Slides on top of rock masses -extremely small and shallow slides in the base of steep and narrow gorges (Fig. 16a) The depleted material represents the regolith layer on top of the rock mass. The size of the instability is dictated by the thickness of the soil cover on top of the bedrock. In cases of thin cover material, landslides are very shallow and rare.
It is visible that harder bedrock material, with colluvium cover, (Mi, MMK) forms steep and narrow gorges. In such settings, landslides are infrequent, but several extremely small and shallow (< 2 m) slides are recognised (Fig. 16a).
2. Slides in firm soil mixtures -extremely small to very small and shallow slides on steep banks (Fig. 16b) The soil material is a poorly sorted mixture of the different fractions (from silt to gravel), and is dominantly incoherent, firm  to stiff, susceptible to water erosion and slightly permeable. For those reasons, flows are rare but if steep banks are formed slides may be very frequent.
Such processes are characteristic for Dar sediments which form characteristic wide, deep ravines with steep banks in which sliding is frequent. The slides can involve the top or bottom of the banks and are frequently very shallow (<2-3 m) and small (Fig. 16b). In these materials, flows are rarely evidenced.
3. Landslides in sands -very small and shallow slides and flows on moderate slopes of ravines (Fig. 16c) Pure NGr sands have a specific behaviour because of the very lightly cemented or compacted grains. The sands demonstrate slight cohesion and are stiff, but are highly permeable and susceptible to water erosion. This material forms very wide and deep ravines with moderate or even steep slopes. The slides and flows are frequent, cover a bigger area and are up to 5 m deep (Fig. 16c). The streambank slides are frequently initiated by the rapid incision of the stream and exceeded critical bank height. Typical flows are also very common.
4. Landslides in dominantly coherent soil complexes -very small to small and shallow landslides on mild slopes of ravines (Fig. 16d) The soil complex is composed of permeable (silty sands) and impermeable (clays and silts) interbeds or lenses. This material forms very wide ravines with mildly inclined slopes (Fig. 16d).
Ravines with relatively mild banks are characteristic of the Pliocene sediments. Specific lithology identifies this unit as the most susceptible to landsliding. Both slides and flows are very frequent. In comparison to other units studied, these landslides are the largest with estimated depth to 5 metres. Irregularly shaped slides and flows are frequent. Several generations of some slides also occur. This unit also has the highest percentage of highly rated landslides in concordance with recent landsliding which has been reported in the past decade by local authorities.

RELATIVE LANDSLIDE AGE AND ROUGHNESS
It is known that landslide features are modified over time by natural processes, therefore, in this research, attention is also fo-cused on the magnitude of deformation and its relationship to the relative age of the landslides.
Although criteria used here for relative age classification are similar to the recommendations of other authors (McCALPIN, 1984;KEATON & DeGRAFF, 1996;BELL et al., 2012;PETSCHKO et al., 2014), the absolute landslide ages remain unknown and may differ significantly to the McCALPIN (1984) classification (which covers landslides from less than hundred to more thousand years). It's mainly because other studies predominantly discuss much larger landslides of different types and geology, and in different climate conditions, all of which influence the erosion processes and time until they become unidentifiable.
As already explained, during visual inspection of the LiDAR derivate maps, landslides were differentiated according to the magnitudes of deformations and clarity of the specific landslide features producing landslide identification confidence (Fig. 7, Fig.  12). These data were supplemented during field prospection with observations about the freshness of landslide scarps, contours, fractures and other features. Similar complementation of LiDAR and field data was used to map landslides in forest areas of Flemish Ardennes in Belgium (Van Den EECKHAUT, 2007). The field verification also enabled the identification of vegetation cover for altogether 95 objects, together with gathering the information on tilted and bent threes, as well as generations of woods. In cases where landslides affect infrastructure or households, the approximate date of landslide activation is determined from local residents or existing landslide reports. If landslides were (re) activated between airborne LiDAR scanning and field prospection, they were considered as recent (< 1 year). From all of the aforementioned attributes, four classes of landslide relative ages were designated: recent, young, mature, old. Similarly to the Mc-Calpin classification (1984), young, mature and old landslides can be considered as inactive, but recent landslides can be active or inactive.
The relative ages of more than 120 landslides is estimated during engineering geological mapping (Fig. 17). According to the evaluation, Vrv-S sands have the highest percentage (20%) of recent landslides, followed by the MMK weathered zone (15%) and Dar deposits (12%). It is important to draw attention to the relatively low number of landslides prospected in the MMK unit that may influence the relevance of this conclusion.
Although it is shown that estimation of relative landslide age is possible, there are insufficient data to correlate these classes to an absolute scale. Consequently, absolute landslide ages remain unknown.
Another problem of relative age estimations are slow movement slides or creeps (BELL et al., 2012). That type of movement produces a very mildly undulated surface without landslide contours which may be similar to old landslides.
The overall average roughness of the terrain surface in the research area is 1.76. The roughness of natural slopes for particular units differ according to lithology and topographic position (Fig. 18). It is clear from the presented results, that hard bedrocks have higher average roughness (Fig. 18 -blue bar; Mi, Γ) with mean SDS values around 2. Higher roughness is also noted for the Pl clastic sediments which are very susceptible to water erosion and have a dense network of gullies and water incised channels, but are also frequently cultivated (discussed further below). Gabbro and amphibolite also have slightly higher SDS values (> 1.7) than all the other sedimentary deposits. Namely, all Miocene clastic sediments, loess, deluvial and alluvial Quaternary sediments have SDS values below 1.7.
Landslides have much higher average SDS values than intact areas: 2.95 for slides and 3.00 for earthflows. In all the affected units, the average landslide SDS values are much higher than for natural slopes (Fig. 18). The rise in SDS values is from around 25 % to over 100 %.
Most of the lithological units have average landslide roughness SDS > 2.5 (Fig. 18). The landslide freshness or even visibility can be disrupted by human activity, but anthropogenic "masking" of landslide features can be considered here as minimal since just several registered landslides are found around households or on cultivated land.
If it is supposed that the same lithology in similar natural environment generate analogous landslides with comparable magnitudes of features which are then equally weathered/eroded during time, the present roughness may correspond to landslide age. Therefore, the roughness quantification could help in the determination of relative landslide ages. For this purpose, the roughness for each particular landslide is averaged (Fig. 19). The presented example demonstrates how a subjective visual impression of landslide roughness can be quantified, providing quantitative roughness data. It is evident that landslide no. 159 has the highest roughness with an average SDS 159 = 3.32 as opposed to landslide no. 323 with an average SDS 323 = 1.87. Somewhere in between is landslide no. 321 with average SDS 321 = 2.54. Since all of the landslides are formed in comparable morphological and geological environments, the field estimations of their relative ages are compatible   Accordingly, the comparison of all relative age estimations and corresponding roughness values displays a clear trend of diminishing landslide roughness with time (Fig. 20). That means that roughness might be used to quantify relative landslide ages. This should be done with great care because roughness is also influenced by other factors including: lithology and material properties, landsliding type, volume and depth of the depleted material, and the period of the dormant stage.
Therefore, the determination of more details for each particular landslide is needed for further landslide roughness studies. More exact research of roughness and landslide age correlations should involve determination of their absolute age (LAHUSEN et al., 2016).
It's important to point out that in this analysis, the human impact on the average roughness of slopes was not evaluated. Namely, it is apparent that a great resolution of the DEM and slope data detects edges of household parcels, road edges, cuttings and plowing marks, all of which raise the roughness of the natural terrain. Alternatively, flattened terrain at infrastructure locations smooths natural terrain irregularities. Therefore, in further studies anthropogenic interventions, natural ravines and gullies should be filtered from the detailed analysis of roughness of natural slopes.

CONCLUSIONS
The paper presents the results of LiDAR scanning and field research of landslides in the wider Kutina area. High resolution Li-DAR scanning of the research area enabled the production of a very detailed terrain model supporting development of the landslide inventory which, up to now, contains 1267 records. Registered landslides are very small and therefore are not catastrophic, but at some locations landslides caused material damage and still represent a threat. Also, landslides affected more than 2 % of the research area and collectively influenced the local geomorphological setting.
The landslide identification was based on LiDAR derivates and orthophoto images. Since not all of the characteristic features are always recognised, the inventoried landslides were differentiated according to the identification confidence levels into three classes: low, medium and high. The confidence levels are based on the visibility and persistence of landslide boundaries and the magnitude of landslide body features. During field prospection, these data were supplemented with observations about the freshness of landslide features and changes or defects in vegetation, altogether enabling estimation of relative landslide age. Although it is shown that estimation of relative landslide age is possible, there are not enough data to correlate these classes to an absolute scale.
In addition, roughness analysis based on the standard deviation of slope (SDS) was performed. The analysis displays significantly higher average roughness of landslide areas (SDS slides =2.95; SDS earthflow =3.00) than the surface roughness of unaffected slopes (SDS natural =1.76). It is presented that the roughness of the landslides depend on several parameters: lithology and material properties, landsliding type, size and depth of the depleted material, and on the period of the dormant stage. Therefore, we can conclude that the roughness of a particular landslide is also a consequence of its age. That trend is also evident in a direct comparison of relative landslide ages with their roughness. Despite this, further research is needed to test the usability of the roughness data for estimation of relative landslide age on larger areas.
In the research area, different geological materials are found: igneous-metamorphic complex, marls, loess, different clastic sediments, deluvial and alluvial materials. The study detected the highest percentage of recent landslides in geological units which are very susceptible to landsliding. For that reason, the majority of the future processes can be expected in the following geological materials. The Pleistocene silts and sands (Vrv-s) with clayey interlayers (Vrv-c) are by far the most susceptible to landsliding with more than 5% of the area affected by slides and more than 3 % of the area by earthflows. M 2 (Dar) silty sands and gravels, and M 7 (NGr) sands can be considered as very susceptible to landsliding because of > 1.5 % area subjected to landsliding. The same units have a noticeable percentage of earthflow processes (around 0.5 %).
The research area is characterised by a dense network of gullies, channels and valleys. Their banks are usually the steepest segments of the terrain which stimulates landsliding. Slides on the gully banks are predominant over the whole research area, and are the most frequent in all units, excluding the Dar deposits. As a consequence, landslide polygons are arranged in grapelike clusters.
Hillshade and slope raster maps also indicate that various lithologies affected the morphology of the gullies. This also has significant impact on instability processes. For example, hard rock masses with a relatively thin soil cover form steep gully banks and extremely small and shallow slides at their bases. On the other end, coherent soils form significantly milder valley flanks with very small to small and deeper landslides. In addition, incoherent materials, dominantly sands or sandy soil mixtures are also prone to flowing processes forming earthflows. Following all relevant results of this research, four different models of landsliding related to bank instabilities are differentiated and described: slides on top of rock masses; slides in firm soil mixtures; slides in sands; slides in dominantly coherent soil complexes.
This study also documents and describes differences in morphological features and contours between slides and earthflows.
The presented results and conclusions ascertain the advantages of and need for complementary field investigations and Li-DAR scanning for landslide processes exploration. Since very small and shallow landslides predominate, high precision and resolution for LiDAR scanning are vital. Detailed LiDAR scanning revealed many unknown and important facts about the geomorphology and landslide processes of the research area, while also directing and optimising field mapping and sampling.

ACKNOWLEDGEMENT
This paper is the result of research conducted in the framework of safEarth (Transnational advanced management of land use risk through landslide susceptibility maps design), co-financed by Interreg IPA Cross-border Cooperation Programme Croatia -Bosnia and Herzegovina -Montenegro 2014-2020. The paper is also partially result of discussions and trainings performed within GeoTwinn project, that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 809943.
We would like to thank our enthusiastic and dynamic research group members: Laszlo Podolszki, Vlatko Gulam, Iris Bostjančić, Mario Dolić, Marina Filipović and Tihomir Frangen, for conducting a field investigation, as well as Radovan Filjak, Radovan Avanić and Tomislav Kurečić for detailed geological mapping of the two pilot areas and discussions on specific properties of each lithostratigraphic units. We also thank the reviewers for their helpful comments and advices.