Veri cation of historical landslide inventory maps for the Podsljeme area in the City of Zagreb using LiDAR-based landslide inventory

Data on landslide locations in the Podsljeme area can be found in three historical landslide inventory maps. Inventory maps from 1979 and 2007 are made based on geomorphological eld mapping and historical records while an inventory map from 2017, compiled for the study area (21 km2), is based on LiDAR data acquired in December 2013. A comparison of three landslide inventory maps was performed with three tests based on landslide statistics, frequency-area distribution, geographical discrepancy of landslides and landslide density maps. The results show signi cant di erences in the number of identi ed landslides as well as in the size and distribution of landslides. A comparison of inventories also showed the unreliability of the existing historical inventories and usefulness of LiDAR data for preparation of complete landslide inventory maps. LiDAR-based inventory for the Podsljeme area could be a valuable tool for a wide range of users i.e., decision-makers, land developers and environmental and civil defense agencies. It is also necessary for landslide susceptibility and hazard mapping, which is a prerequisite for landslide risk reduction.


Introduction
Landslides are de ned as the movement of a mass of rock, debris, or earth down a slope (Cruden, 1991) and play an important role in the evolution of landscapes.The term "landslide" describes a wide variety of processes that result in the downward and outward movement of slope-forming materials including rock, soil, arti cial ll, or a combination of these.The materials may move by falling, toppling, sliding, spreading, or owing (Cruden and Varnes, 1996).Landslide inventory should give insight into the locations of landslide phenomena, the types, failure mechanisms, causal factors, date of landslide triggering, frequency of occurrence, volumes and the damage that has been caused.Systematic information on the type, abundance, and distribution of landslides in the form of inventory is the preliminary step toward landslide susceptibility, hazard, and risk assessment (Mihali Arbanas and Arbanas, 2014;Reichenbach et al., 2018).Another purpose of landslide inventory maps is the information about landslide statistics, which includes the analyses of the landslide area, shape, and location.The size of the smallest landslide (Guzzetti et al., 2000) and frequency-size statistics allow the veri cation of the reliability and the comple teness of landslide inventory maps.Malamud et al. (2004) de ned that a substantially complete landslide inventory map must include a substantial fraction of all landslides at all scales.Generally, historical landslide inventories do not include a substantial fraction of the smallest landslides because smaller landslides are often lost due to erosional processes, anthropic in uences, and vegetation growth (Malamud et al., 2004;Bell et al., 2012).
In the past, the most popular techniques for landslide mapping were aerial photo interpretation and geomorphological eld mapping (Guzzetti et al., 2012).When using aerial photos, the accuracy of an inventory depends on the scale, date and quality of the stereo pairs, vegetation cover density in the study area, on the type, quality, and characteristics of the stereoscopes, and also on the skills and experience of the interpreters.Geomorphological eld mapping provides very accurate landslide inventory if landslide boundaries are geo-referenced using GPS, but the most common problems dur- The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2019, pp.45-58, DOI: 10.17794/rgn.2019.1.5ing eld mapping are that landslides are often too large to be seen completely in the eld and old landslides are often partially or totally covered by vegetation.A common misconception among the researchers is that mapping landslides in the eld is more accurate than mapping landslides remotely, i.e., the perspective offered by a distant view of the landslide is preferable and can result in more accurate and more complete landslide maps (Guzzetti et al., 2012).
In the last two decades, new advanced methods and technologies are used for landslide inventory mapping over large areas (Guzzetti et al., 2012), including: (1) analysis of surface morphology using high resolution Li-DAR digital elevation models (DEMs), and (2) analysis and interpretation of satellite images, including panchromatic, multispectral, hyperspectral and synthetic aperture radar (SAR) images.Satellite images with spatial cell resolutions larger than 3 m (e.g., Spot, Landsat) and SAR images (e.g., Radarsat, ERS, Envisat) are very useful for the quanti cation of small displacements on single large landslides (Singhroy, 2005) In this research, we will focus on the hilly area of the southern foothills of Medvednica Mt., also known as the Podsljeme area, which is a very attractive residential part of the City of Zagreb.Although, the land use and development measures have been taken and unsafe construction practices was regulated, the relative proportion of human induced landslides in the Podsljeme area have been continuously increasing over the last 90 years (Jurak et al., 2008;Mihali Arbanas et al., 2014).Historical landslide data in the City of Zagreb can be found in landslide inventories from different periods, professional and scienti c papers, graduate and doctoral thesis, and in reports of geotechnical investigations (Mihali et al., 2012).In the framework of the Croatian-Japanese SATREPS FY2008 scienti c project (Mihali and Arbanas, 2013) the whole area of 180 km 2 was scanned by LiDAR in high resolution in 2013 aiming to make a landslide inventory map.Preliminary analyses of Li-DAR hillshade, slope and contour maps showed that visual interpretation of surface morphology could enable more reliable identi cation of small and shallow landslides, typical for the Podsljeme area (Mihali et al., 2013) than geomorphological led mapping and aerial photo interpretation.
This paper presents a comparison of two historical landslide inventories from the Podsljeme area (from 1979 and 2007) with LiDAR-based landslide inventory compiled for the study area (Bernat Gazibara et al., 2017b).Landslide inventories from 1979 and 2007 have a practical application in the City of Zagreb of ces for physical planning and construction, land-use planning and civil protection.The main purpose of the comparison was to evaluate the completeness of these inventories, as well as to evaluate the reliability of landslide data in both historical inventories.The quality and reliability analysis of landslide inventories in the City of Zagreb have not been performed so far and the novelty of this study is a comparison of landslide inventory maps prepared by eld mapping and interpretation of LiDAR data in the Podsljeme area.

Study area
The City of Zagreb is located in the western part of NW Croatia (Bognar, 2001) and it is part of the European Pannonian Basin (Malvi and Veli , 2011).Within the city limits, the relief changes in the north-south direction from lower mountain type (the ridge of Medvednica Mt.) and hilly type (the south-eastern foothills of Medvednica Mt., Podsljeme area), through lowland (the uvial oodplain of the Sava River) to the hilly type on the south (the hills of Vukomeri ke Gorice), in accordance with changes of geological settings (Mihali Arbanas et al., 2012).
The study area comprises 21 km 2 in the western part of the Podsljeme area (see Figure 1).This area is mostly urbanized and densely populated; the current land use includes 11.7 km 2 of arti cial surfaces, 4.7 km 2 of agricultural areas, and 4.6 km 2 of forests.The elevation ranges from 115 to 612 meters a.s.l., the prevailing slope angles (44 % of study area) range from 12° to 32° and 90 % slopes have slope angles >5°.The study area is comprised of Upper Miocene and Quaternary sediments (Šiki , 1995).The Upper Miocene deposits are strati ed marls, silts and sands with moderately to slightly-inclined bedding (bedding slope angle in the range of 10-20°).The top parts of Miocene deposits are ne-grained soils, mostly silts.The Quaternary deposits are heterogeneous mixtures of mostly impermeable clayey soils.The geological contact between the Miocene sandy silts and the Quaternary clayey soils is highly susceptible to sliding (Jurak et al., 2008).Preliminary analysis based on the spatial distribution of slope angle units (5-55°) and Miocene and Quaternary deposits showed that 76 % of the study area is susceptible to sliding.
In the Podsljeme area, active geomorphological processes are linear erosion as a result of ash oods from Medvednica Mt. and landslides (Jurak et al., 2008).Dominant landslide types in Podsljeme area are small and shallow landslides (Mihali Arbanas et al., 2016).Vrhovec landslide (see Figure 2) is a typical landslide in the Podsljeme area, formed along geological contact between Pleistocene (Q 1 ) ne-grained soils and Pontian (M 7 ) sandy-silty soils, initiated by human activities of uncontrolled disposal in the upper part of the slope.The preparatory causal factors of slope instabilities in the Podsljeme area depend on the geomehanical properties of soils, geomorphological processes and different types of man-made processes.Landslides are triggered primarily by rainfall (Bernat et al., 2014a) and man-made activities.For example, in the winter of 2012/2013, the prolonged heavy rainfall periods and the rapid melting of a thick snow cover (with a total of 378.7 mm from January to March) triggered more than 50 landslides in the Podsljeme area (Bernat et al., 2014b).The climate of the City of Zagreb is continental with a mild maritime in uence, with a mean annual precipitation (MAP) of 887 mm, while precipitation is mostly recorded in the period from May to November (Zaninovi et al., 2008).

Landslide inventory maps for Podsljeme area
In the Podsljeme area, data for landslide inventories was systematically collected four times in the last 50 years.The rst landslide inventory map for the hilly area of the City of Zagreb was prepared by the Croatian national geological survey in 1967 (Šiki , 1967).Landslide data from this inventory is not publicly available today.The second landslide inventory for the Podsljeme area is an analogue landslide inventory map from 1979

Historical landslide inventory map from 1979
The landslide inventory map from 1979 was made on the basis of geomorphological eld mapping and historical geotechnical documentation at a 1:10,000 scale (Polak et al., 1979).The landslide inventory covers an area of 105 km 2 and includes 812 landslide polygons, e.g., 406 active landslides, 294 inactive landslides and 112 creeps.The number of landslides is even bigger, because groups of small landslides are contoured with one larger landslide polygon.The total landslide area is 2.44 km 2 , i.e., 2.32 % of the inventory area, and landslide density is 9.7 landslides per square kilometer.Mapped landslides extend in size from 452 m 2 to 289,501 m 2 , while the most abundant landslides have an area in range from 5,000 to 6,000 m 2 .The landslide inventory from 1979 is an analogue map (see Figure 3a) and only particular landslides or groups of landslides are described in the accompanying report.A landslide susceptibility map depicts four zones of relative susceptibility derived from direct geomorphological mapping in the eld (see Figure 3b).Susceptibility zoning characterized that 10.5 % of the area is classi ed as stable terrain, 78 % as conditionally stable slopes, 11 % as conditionally unstable slopes and 0.5 % as unstable slopes.

Historical landslide inventory map from 2007
The landslide inventory from 2007 (see Figure 4) was made on the basis of geomorphological mapping, geotechnical reports and historical records at a 1:5,000 scale (Miklin et al. 2007).The landslide inventory covers an area of 175 km 2 and shows contours of 538 slides, 125 creeps, 15 ows, 14 rock falls, and 15 unclassi ed LiDAR system was used, and captured data at a pulse rate of 266 kHz with a surface point horizontal accuracy of 8 cm and a vertical accuracy of 4 cm (Bernat Gazibara et al., 2017b).The last return of LiDAR data, with an average density of 5 points per square meter, was used for creating a bare-earth DEM with a 30 cm resolution.The landslide inventory map has been produced on a small part (21 km 2 ) of the Podsljeme area.Topographic derivative datasets used for interpreting landslide morphology were hillshade maps, slope maps, contour lines, curvature and surface roughness maps.All landslides were mapped at a larger scale (1:100-1:500) to ensure correct delineation of the landslide boundaries.During landslide identi cation, an orthophoto from 2012 was used to check morphological forms along roads and houses, such as arti cial lls and cuts, which can have a similar appearance to landslides on DTM derivatives.It is possible that some very old or prehistorical slides were overlooked during the landslide mapping process because landslide morphology was not visible on derivative LiDAR maps.Furthermore, remediated landslides and very small landslides (<40 m 2 ) were not mapped because it was impossible to identify their typical landslide morphology from LiDAR DTM 0.3 x 0.3 m resolution.

LiDAR-based landslide inventory map from 2013
Airborne laser scanning (ALS) was undertaken for the Podsljeme area (total area of 180 km 2 ) in December 2013, which corresponds to winter leaf-off period in Croatia.Moreover, LiDAR data was acquired 8 months after extreme weather periods in the winter of 2012/2013 when more than 50 landslides were (re)activated in the City of Zagreb (Bernat et al., 2014a).In this study, the The landslide inventory map (see Figure 5) consists of 702 landslides, while the total landslide area is 0.5 km 2 which covers 2.43 % of the inventory area (21 km 2 ).Mean landslide density is 33.3 landslides per square kilometer and the smallest landslide identi ed in the pilot area has an area of 43 m 2 .Mapped landslides extend in size to a maximum of 8,064 m 2 , while the most frequent landslides in the inventory have an area of 400 m 2 .Most of the landslides can be dated as recently (re)activated due to sharp appearance or a high degree of preservation of the landslide morphology.Therefore, the landslide inventory map of the study area represents a combination of seasonal inventory (landslides (re)activated in the winter of 2012/2013) and historical landslide inventory (sum of landslides events activated over a period of tens or hundreds of years).

Research methodology
For the purpose of data analysis, landslides from the 1979 analogue inventory map and 2007 on-line inventory map were digitized to obtain digital landslide boundaries for the study area (21 km 2 ) located in the northern part of the Podsljeme area (see Figure 1).A comparison of three landslide inventory maps was performed using the following analyses: frequency-area distribution, geographical discrepancy of landslides and landslide density maps.
The rst test compares the frequency-area statistics of landslides obtained from three different landslide inventories.Frequency-size statistics (Malamud et al., 2004) allows veri cation of the reliability and the completeness of landslide inventory maps.Landslide frequencysize distributions are characterized by a positive powerlaw scaling for the small landslides and a negative power-law scaling for moderate and large landslides, separated from each other by a rollover.Generally, historical landslides inventories do not include a substantial fraction of the smallest landslides because geomorphological features of smaller landslides are often lost due to erosional processes, anthropic in uences, and vegetation growth (Bell et al., 2012).
The second test aims at evaluating the degree of cartographic matching, i.e., geographical discrepancy between the three historical landslide inventories.The method proposed by Carrara et al. (1992) was used to compute the overall error index (E) expressed by the Equation 1: (1) Where: E -error index; A 1 -total landslide area in the rst landslide inventory (km 2 ); A 2 -total landslide area in the second inventory (km 2 ).
From Equation 1, the degree of matching (M) between two inventory maps is expressed by Equation 2: (2) Where: E -error index; M -matching index.
If two inventory maps show exactly the same landslides in the same positions (a rather improbable situation) matching is perfect (M=1) and the error is nil (E=0).If two landslide maps do not match completely, cartographic matching is nil (M=0) and the error is maximum (E=1).
The third analysis compares the abundance of slope failures in the three landslide inventories.To obtain this comparison, the study area was subdivided into slope units, i.e., distinctly delimited portions of the terrain that contain a set of conditions that differ from the adjacent units across de nable boundaries (Guzzetti et al., 2000).Terrain subdivision was done on 5x5 m DTM, using ArcGIS tools Flow accumulation and Basin.The results are slope units de ned with the crest line and the drainage line.Based on three landslide inventory maps, the percentage of landslide area, i.e. landslide density, was computed for each slope unit.The result is a landslide density map for each analyzed landslide inventory, which provides insight into the expected occurrence of landslides in any part of the investigated area without leaving unclassi ed areas (Guzzetti et al., 2000).Furthermore, Bulut et al. (2000) proposed that landslide density maps can be used as a weak proxy of landslide susceptibility.7a, b), and can be described as sustainably incomplete.In the 2007 landslide inventory, small landslides are not present, while large landslides are numerous (28 landslide phenomena with an area larger than 30,000 m 2 ).

A comparison based on landslide statistics and frequency-size distribution
Contours of these large landslides are more probably borders of areas with many small instabilities, than contours of single landslides.The transition between the positive and the negative power-law relations can be used to distinguish between small and medium landslides and to identify the size of the most abundant (i.e., most frequent) landslide in the inventory.The size of the most abundant (i.e., most frequent) landslide in the LiDAR-based inventory is approx.400 m 2 , while in the landslide inventory from 1979, it is approx.2400 m 2 and in the landslide inventory from 2007, it is approx.12,000 m 2 .Based on the size of the most abundant landslide in the three inventory maps, it can be concluded that the area of the most abundant landslides decreases with an increase in the completeness of the inventories (Galli et al., 2008).

Comparison based on geographical position
Geographical discrepancy between the three landslide inventories was analysed based on the method proposed by Carrara et al. (1992).In the ArcGIS 10.0, pair-wise geographical union and intersection of the three inven-tory maps were performed to compute the error, E (Eqs. 1) and matching, M indexes (Eqs.2).The results of geographical comparison of landslide inventory maps in the study area are summarized in Table 2.The resulting mapping error is very high (ranging between 0.94 and 0.96) which corresponds to a very low degree of matching (ranging between 0.06 and 0.04).Landslide surface common to the two inventory maps from 1979 and 2007 is 1 %, the overlap between landslide inventory from 1979 and 2013 is only 0.22 % and the overlap between landslide inventory from 2007 and 2013 is 0.9 %.
The landslide inventory maps from 1979 and 2007 were made on the basis of geomorphological eld mapping at a 1:10,000 and 1:5,000 scale.The purpose of the following analysis was to identify the impact of the location, mapping and digitizing errors associated with the eld mapping (mapping without GPS technology in 1979) on mismatch, due to different techniques used for landslide mapping.A series of mapping errors, E (Eqs. 1) and matchings, M indexes (Eqs.2), were computed for systematically enlarged landslide areas (buffers of 1 m, 2 m, 5 m, 10 m and 20 m) for all three landslide inventories.The results of the analysis are shown in Figure 8.With an increasing buffer size, mapping error decreases gradually for buffers < 2 m and then decreases rapidly, reaching its minimum value for the largest buffer of 20 m.Conversely, map matching shows a maximum value for buffer of 20 m.Mapping match is approximately 4-15 %, while the reaming mismatch (85-96 %) can be assigned to different geomorphological interpretations (1979 vs. 2007) and different techniques used for landslide mapping (1979 vs. 2013 and 2007 vs. 2013).

Comparison based on landslide density maps
The study area (21 km 2 ) was divided into 70 slope units de ned by Carrara et al. (1991).The prevailing orientation of slope units is north-south and the most abundant size of slope is 18 ha (min.= 1 ha, max.= 194 ha, mean = 30 ha, median = 18 ha, std.dev.= 35 ha).Based on three landslide inventories, the percentage of landslide area was computed for each slope unit and the results are shown in Figure 9.Moreover, spatial distribution of landslide density varies signi cantly for the three inventory maps, which was expected due to differences in landslide sizes and the mismatch between their inventories (see Figure 5).Guzzetti et al. (2000) suggest that slope units with a small portion of landslide area (< 3 %), e.g., landslide deposit crossing a stream line or a few small, super cial landslides, can be considered as a stable area.Based on a visual comparison of the slope units and the inventory maps (see Figure 9), and relative small landslides in the study area, the cutoff value was set at 1 % of landslide area, i.e., slope units with a proportion of landslide area less than 1 % can be considered as stable.Table 3 shows a number of stable and unstable slope units in each landslide density map obtained from three landslide inventories.The density    The degree of matching between the three density maps can be visually estimated from Figure 9.The least overlapping is between the density maps obtained from 2007 and the LiDAR-based landslide inventories.The best matching is between the density maps obtained from 1979 and the LiDAR-based landslide inventories, because both maps show small to medium density of landslides in slope units which corresponds to similar total landslide areas (0.61 km 2 in 1979 inventory and 0.51 km 2 in LiDAR-based inventory).The density map obtained from the 2007 inventory is substantially different from the other two density maps and shows signicantly higher landslide densities.In the central part of the study area, landslide densities reached values of 30-50 % which is primarily affected by large and very large landslide polygons representing mapped phenomena in the inventory.

Discussion
Three tests designed to assess speci c aspects of the quality of a landslide map were performed to compare three landslide inventory maps in the Podsljeme area.Based on landslide statistics for the study area, the size of the most abundant landslide in the inventory from 1979 is 5 times larger (2,407 m 2 ), and in the landslide inventory from 2007, even twenty-eight times larger (11,981 m 2 ) than the most frequent landslide size in the LiDAR-based inventory (427 m 2 ).It can be concluded that the size of landslides, and even landslide zones, in both historical inventories from 1979 and 2007 are unrealistic.According to previous scienti c investigations (Mihali Arbanas et al., 2016), dominant landslide types in the Podsljeme area are small, super cial to moderately shallow landslides (landslide volume <10 5 m 3 ; landslide depth <20 m).In fact, analysis of landslide sizes from preliminary or detailed investigation reports performed in the period from 1968 to 2008 showed that most of the investigated landslides in the Podsljeme area cover less than 8,000 m 2 (Mihali et al., 2012).The exception is the large, deep-seated Kostanjek landslide in the western part of the Podsljeme area with a landslide area of approx. 1 km 2 (Krka et al., 2017).
Based on the frequency-area distribution of the three landslide inventories, it can be concluded that historical landslide inventories from 1979 and 2007 are generally incomplete, because small landslides are missing, while the LiDAR-based inventory is sustainably complete.In the historical landslide inventories from 1979 and 2007, small landslides are not present, which is an important drawback re ecting the signi cant limitation of practical application of both inventories.In the last few years, initiations or reactivations of small and super cial landslides in the Podsljeme area were recorded two or three times per year (Bernat Gazibara et al., 2017a), usually during wet periods (November-April), causing damage on private houses and roads.Inventory compiled on the basis of LiDAR DTM acquired in December 2013 presents all landslides with more realistic landslide contours and gives information about the exact position of landslides in a whole range of landslide sizes (small, shallow and large, deep-seated landslides).
The three landslide inventories were used to calculate the percentage of landslide area in slope units, which can be used for a proxy estimation of landslide susceptibility (Bulut et al., 2000) in the study area.Since the LiDARbased inventory is proved to be the most complete landslide inventory map, a derived density map can be considered a more reliable proxy estimation of landslide susceptibility, comparing to the other two analysed inventories.The derived density map shows that 45 % of slope units (32 % of the study area) can be described as stable, and 55 % of slope units (68 % of the study area) can be considered as potentially unstable or landslide prone areas.The historical inventory from 2007 shows much higher landslide densities per slope units and gives an unrealistic picture of the spatial distribution of landslide susceptibility as well as locations of assumed landslide prone areas in the study area.Based on the quality, reliability and completeness of the two historical landslide inventories, it can be concluded that these inventories cannot be used for a susceptibility and hazard assessment in the study area.
A very high mismatch (E=85-96 %) between the three landslide inventory maps could be due to the application of different techniques used for landslide identi cation and mapping.Geomorphological eld mapping in the study area is not an adequate method for landslide inventory mapping because most of the landslides (more than 80%) are located in forested and semi-natural areas (Bernat Gazibara et al., 2017b).The ability to follow a landslide boundary accurately in the forest is limited by the reduced visibility of the slope, a failure due to dense vegetation cover, the local perspective, the size of the landslide, and the fact that the landslide boundary is often indistinct (Santangelo et al., 2010).LiDAR DTM enables a distant view of the slope and landslide, which is preferable during landslide inventory mapping to reach better results, and a more accurate and more complete landslide inventory map (Guzzetti et al., 2012).
The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2019, pp.45-58, DOI: 10.17794/rgn.2019.1.5 The most important bene t of LiDAR is that it enables landslide mapping in forested terrain because laser pulses can penetrate the canopy, unlike optical aerial or satellite images that do not penetrate the canopy.Another advantage of using LiDAR data is that it enables a (semi) automated detection of landslides (Van Den Eeckhaut et al. 2012) based on morphological parameters.Comparisons over the past few years indicated that landslide recognition using LiDAR DTMs is up to 5-10 times quicker than traditional photogrammetric techniques in the same landscape (Pirasteh & Li, 2017).

Conclusions
The historical landslide inventories from 1979 and 2007 and the detailed LiDAR-based landslide inventory from 2013 provide very different information on the size, distribution, and density of landslides in the study area (21 km 2 ).Historical landslide inventory maps from 1979 in a 1:10,000 scale and from 2007 in a 1:5,000 scale are the result of geomorphological eld mapping.Based on eld mapping, the size of a single landslide or landslide zone (area with many small landslides) is in a range from 467 to 43,228 m 2 in the 1979 inventory and from 282 to 317,262 m 2 in the 2007 inventory.The average size of landslides in the LiDAR-based inventory is 730 m 2 , and 75 % of the landslide bodies showed a size between 160 and 2,000 m 2 .Besides different average sizes of landslides in inventories, landslide inventory maps also differ in the spatial distribution of landslides.Cartographic error between analyzed landslide inventory maps is approximately 85-96 % and it can be assigned to different geomorphological interpretations (1979 vs. 2007) and different techniques used for landslide mapping (1979 vs. 2013 and 2007 vs. 2013).Differences in landslide size and spatial distribution also result in large differences in landslide density maps obtained for the 1979, 2017 and 2013 landslide inventory maps.The highest percentage of unstable area was obtained with the 2007 landslide inventory, while the density map obtained from the 1979 landslide inventory shows that only 58.7 % of the study area is unstable.The landslide density map obtained from the LiDAR-based inventory shows that 68 % of the study area (14.3 km 2 ) is unstable or landslide prone, and based on the highest degree of landslide inventory completeness, it can be concluded that this density map gives the best proxy estimation of landslide susceptibility in the study area.
A comparison of the three inventories at the study area makes it possible to draw conclusions for the whole Podsljeme area (180 km 2 ).The LiDAR DTM acquired in December 2013 enables preparation of a complete landslide inventory map with a whole range of sizes (small, shallow and large, deep-seated landslides) starting with a landslide area of 40 m 2 .
The LiDAR-based landslide inventory map, in contrast to the existing historical inventory maps, is suitable for practical application, planning, development or civil protection purposes in the Podsljeme area, because of the high reliability of its landslide contours.Landslide inventory from 2013 contains numerous small landslides which are easier to remediate, in case that they present a risk for people, material properties or the environment.The spatial distribution of landslides in the study area shows different patterns which implies the usefulness of a susceptibility assessment in the Podsljeme zone, as preliminary information for measures aimed at risk prevention.
Although the LiDAR-based inventory provides realistically accurate landslide information, visual interpretation of LiDAR-based DTM derivatives is still a timeconsuming process and thus limited to large-scale studies.The advantage of using LiDAR data is that it enables a semi-automated detection of landslides based on morphological parameters of mapped landslides and further research will focus on the automatization of landslide mapping and the creation of landslide inventory and hazard maps for the entire Podsljeme area (180 km 2 ).

Figure 1 :
Figure 1: Relief map of the City of Zagreb with the outline of three existing landslide inventory maps and study area.Histograms show the elevation and slope angle for the study area.

Figure 2 :
Figure 2: Composite display of the Vrhovec landslide on: (a) Orthophoto; (b) Hillshade map overlain by a transparent slope map.

Figure 5 :
Figure 5: Landslide inventory map based on LiDAR data acquired in 2013 (Bernat Gazibara et al., 2017b).Landslides are presented on a hillshade map overlain by a transparent slope map.

A
comparison of landslide contours from three analyzed inventories, Polak et al. (1979), Miklin et al. (2007) and LiDAR-based inventory (Bernat Gazibara et al., 2017b), is presented on Figure 6.Descriptive statistics are listed in Table 1.In the area of 21 km 2 , historical landslide inventories from 1979 and 2007 consist of approx.160 landslides, while LiDAR-based inventory has 702 landslides or approx.4.5 times more mapped landslides.The total landslide area in inventory from 2007 (3.73 km 2 ) is seven times larger than the landslide area in the inventory from 1979 (0.61 km 2 ) and in the LiDAR-based inventory from 2013 (0.51 km 2 ).The largest landslide in the inventory from 2007 extends over 317,262 m 2 , which is seven times larger than the size of the largest landslide in the inventory from 1979 (43,228 m 2 ) and almost forty times larger than the size of the largest landslide in the LiDAR based inventory (8,064 m 2 ).The smallest landslide in the LiDAR-based The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2019, pp.45-58, DOI: 10.17794/rgn.2019.1.5inventory (43 m 2 ) is ten times smaller than the smallest landslide in the landslide inventory from 1979 (467 m 2 ) and six times smaller than the smallest landslide in the landslide inventory from 2007 (282 m 2 ).Frequency-size distribution of landslides in the three inventory maps in the study area is shown in Figure 7. LiDAR-based landslide inventory shows 'universal distribution' (see Figure 6c) described by Malamud et al. (2004) and an inventory map can be described as sustainably complete.Historical landslide inventories from 1979 and 2007 do not show normal distribution (see Figure

Figure 8 :
Figure 8: Mapping error, E and mapping match, M showing degree of cartographic matching for three pair-wise combinations of landslide inventory maps in the study area (21 km 2 )

Figure 9 :
Figure 9: Landslide density maps for the study area (21km 2 ) based on three landslide inventory maps

Table 1 :
Characteristics of the three inventory maps available for the study area (21km 2 ) in the Podsljeme area

Table 2 :
Landslide inventory comparison in the study area (21 km 2 ), Podsljeme area

Table 3 :
Comparison of stable and unstable slope units based on landslide density computed for the three landslide inventory maps available for the study area

Landslide inventories No. of stable slopes Area (km 2 ) Area (%) No. of unstable slopes Area (km 2 ) Area (%)
The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2019, pp.45-58, DOI: 10.17794/rgn.2019.1.5map obtained from the 1979 landslide inventory shows 40 stable slope units (57 % of slope units), while the landslide density map made from the landslide inventory in 2007 has only 20 stable slope units (29 % of slope units) and the LiDAR-based density map has 32 stable slope units (45 % of slope units).