Landslide inventory mapping using LiDAR data in the City of Zagreb (Croatia)

Landslides in the Podsljeme area of the City of Zagreb cause significant economic losses, which have been increasing over the last several decades due to the urbanisation of hilly areas and the influence of climate changes. An airborne LiDAR digital terrain model (DTM) with a spatial resolution of 30 × 30 cm was used to prepare a landslide inventory map of the pilot area (21 km) with more than 700 identified landslides. The area of the smallest identified landslide is 43 m, while 90% of the landslides are between 100 and 2,000 m. The frequency–size distribution of all mapped landslides in the pilot area shows a very high level of landslide inventory completeness. Therefore, it is concluded that the LiDAR-based terrain model is a valuable tool for the preparation of detailed landslide inventories in heavily vegetated regions such as the hilly area of Medvednica Mt. ARTICLE HISTORY Received 17 June 2018 Revised 13 September 2019 Accepted 20 September 2019


Introduction
Landslide is the movement of a mass of rock, debris, or earth down a slope under the influence of gravity (Cruden & Varnes, 1996), while a landslide inventory map shows the location, spatial extent and type of landslides in a region (Guzzetti et al., 2012;Wieczorek, 1984). The spatial distribution of present and historical landslides in certain areas is a preliminary step towards landslide susceptibility, hazard, and risk mapping (Galli, Ardizzone, Cardinali, Guzzetti, & Reichenbach, 2008), while the final result of landslide mapping is the consideration of landslide phenomena as spatial constraints in land-use planning. Another purpose of landslide inventory maps is the information about landslide statistics, which includes the analyses of landslide area, shape, and location. The size of the smallest landslide (Guzzetti, Cardinali, Reichenbach, & Carrara, 2000) and frequency-size statistics (Malamud, Turcotte, Guzzetti, & Reichenbach, 2004) allow the verification of the reliability and the completeness of landslide inventory maps. Malamud et al. (2004) defined 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 influences, and vegetation growth (Bell, Petschko, Röhrs, & Dix, 2012;Malamud et al., 2004).
LiDAR (Light Detection and Ranging) is a remote sensing technique that uses pulses of light to gather information about the surface of the Earth below it. In this study, a high-resolution digital terrain model (DTM) derived from aerial LiDAR surveys serves as the base data set for creating derivative maps. Over the last few years, the LiDAR derivatives have been used to identify and map landslide morphology in areas that are cultivated (Ardizzone, Cardinali, Galli, Guzzetti, & Reichenbach, 2007) and partially or completely covered by dense vegetation (Petschko, Bell, & Glade, 2015;Razak, Straatsma, Van Westen, Malet, & de Jong, 2011;Van Den Eeckhaut et al., 2007). The landslide maps obtained through the visual analysis of LiDAR-derived DTMs had better statistics for the landslide size (area) than the inventories obtained through field mapping or the interpretation of aerial photographs (Ardizzone et al., 2007;Razak et al., 2011;Schulz, 2007;Van Den Eeckhaut et al., 2007).
Airborne laser scanning (ALS) was performed in December 2013 for the Podsljeme area of 180 km 2 after extreme weather periods in winter 2012/2013 when more than 50 landslides were (re)activated in the City of Zagreb. The main objective of this paper is to present a detailed landslide inventory map for the study area of 21 km 2 in the Podsljeme area based on visual interpretation of topographic derivatives of high-resolution LiDAR DTM. The other objective is the analyses of landslide statistics obtained from an inventory map and the analysis of the spatial distribution of landslides due to land use.

Study area
The city of Zagreb is located in the western part of the mega-geomorphological region of the European Pannonian Basin (Bognar, 2001). The study area comprises 21 km 2 in the hilly area of the city of Zagreb (known as the Podsljeme area). The investigated slopes are located along the southeastern-facing slope of the Medvednica Mt (Figure 1). This area is mostly urbanised and densely populated; the current land use includes 11.7 km 2 of artificial surfaces, 4.7 km 2 of agricultural areas, and 4.6 km 2 of forests. Almost 90% of the studied slopes have slope angles >5°, while the prevailing slope angles have a range of 5°-12°(28% of pilot area) and 12°-32°(44% of pilot area), making this area potentially prone to sliding. The study area is composed of upper Miocene (55% of pilot area) and Quaternary (37% of pilot area) sediments (Šikić, 1995). The Upper Miocene deposits are stratified sands, silts and marls, with moderately to slightly inclined bedding and bedding slope angles that range from 10°to 20°with the direction of slope (Avanić et al., 2003;Vrsaljko et al., 2006). The Quaternary deposits are heterogeneous mixtures of unfoliated, mostly impermeable clayey-silty soils that are prone to sliding (Nonveiller, 1964;Polak, 1978). Based on previous landslide investigations, Jurak et al. (2008) concluded that the geological contact between the Miocene or Quaternary clayey silty soils and the Miocene sandy silts is highly susceptible to sliding due to different permeability of sediments and periodical water saturation. The preparatory causal factors of slope instabilities in Podsljeme area depend on geomehanical properties of soils, geomorphological processes and different types of man-made processes (Mihalić Arbanas, Krkač, & Bernat, 2016). Preliminary analysis based on two predominating landslide casual factors, deposit types prone to sliding (upper Miocene and Quaternary sediments) and slopes inclination of 5°-55°, showed that 76% of the pilot area can be identified as landslide-prone area.
The climate of the City of Zagreb is continental with a mild maritime influence, with a mean annual precipitation (MAP) of 887 mm, and precipitation falls mostly in the period from May to November. Landslides are caused primarily by rainfall (Bernat, Mihalić Arbanas, & Krkač, 2014a) and human activities. In winter 2012/2013, the prolonged heavy rainfall periods and the rapid melting of a thick snow cover caused abundant landslide events in the hilly area of the city of Zagreb (Bernat, Mihalić Arbanas, & Krkač, 2014b). According to the precipitation percentile calculated based on data from 1961 to 1990, winter 2012/2013. In NW Croatia can be described as extremely rainy due to seasonal precipitation percentile above the 98th percentile (DHMZ, 2013). Monthly precipitations in the same period of 2013 were two to three times higher than the average monthly values of the same period from 1862 to 2012 recorded at the Zagreb-Grič weather station. Analysis of the 3-month period from January to March showed that cumulative precipitation for the analysed period in 2013 has the highest value (378.7 mm) measured in the last 150 years. Information regarding the (re)activation of more than 50 landslides was conducted by the city administration responsible for landslide remediation or civil protection. All reported landslides caused damage to houses and roads, which implies that the (re)activation of landslides in forests were not reported (Bernat et al., 2014b). Moreover, the local authorities do not have insight into the landslide location, so the main problem with the current practice of landslide risk management in the City of Zagreb is the lack of a suitable landslide inventory and geotechnical report archive (Mihalić Arbanas et al., 2016).

Data and methods
LiDAR is an optical remote-sensing technique as used to measure range to and reflectance of objects on the earth surface, producing highly accurate point cloud datasets (Wehr & Luhr, 1999). Laser pulses emitted from a LiDAR system reflect from objects both on and above the ground surface e.g. vegetation and buildings. The first returned laser pulse is associated with the highest feature in the landscape like a treetop, the top of a building, or the ground, in which case only one return will be detected by the LiDAR system. Multiple returns are capable of detecting the elevations of several objects within the laser footprint of an outgoing laser pulse and the last return for bare-earth terrain models. Airborne Lidar technology is well suited for the generation of high resolution DTMs (Ackermann, 1999).
LiDAR data for the Podsljeme area was acquired in the framework of the Croatian-Japanese project 'Risk Identification and Land-Use Planning for Disaster Mitigation of Landslides and Floods in Croatia' . Airborne laser scanning (ALS) was undertaken in December 2013, which corresponds to the winter leaf-off period in Croatia. The LiDAR system used in this study 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. The last return LiDAR data, with an average density of five points per square metre, were used for creating a bare-earth DEM with a 0.3 m resolution.
The topographic derivative datasets used for interpreting the landslide morphology were hillshade maps, slope maps, contour lines, and curvature and surface roughness maps (Figure 2). Derivatives of the LiDAR-based DTM were computed in ArcGIS 10.0 from the LiDAR DTM, while for the calculation of curvature and surface roughness, an additional ArcGIS Toolbox Geomorphometry and Gradient Metrics (Evans, Oakleaf, Cushman, & Theobald, 2014) was used. Landslide identification on the derivatives of LiDAR DTM is based on the recognition of landslide features (e.g. concave main scarps, hummocky landslide bodies and convex landslide toes). Hillshade maps calculated with different azimuth angles (45°, 135°and 315°) were used to avoid shades covering landslide features hindering the delineation of the landslide boundary (Schulz, 2007). A slope map (Figure 2 (b)) was used to interpret steep areas, which may also indicate scarps, ridges, or landslide toes (Ardizzone et al., 2007). A curvature map, i.e. the second derivative of the surface (Figure 2(c)), and contour lines with a spacing of 0.5 m were helpful for the identification of concave and convex features such as landslide accumulation and depression zones, respectively. The terrain roughness was computed based on a method introduced by Riley, De Gloria, and Elliot (1999), and it represents the variability of elevation or slope in a local neighbourhood, which can be used for the identification of unstable slopes.
All landslide mapping was performed at a large scale (1:100-1:500) to ensure the correct delineation of the landslide boundaries. During landslide identification, an orthophoto map from 2012 was used to check the morphological forms along roads and houses, such as artificial fills and cuts, which can have a similar appearance to landslides on DTM derivatives. Each mapped landslide polygon was assigned with the landslide type, certainty of landslide identification and precision of mapping. The certainty of landslide identification was expressed as 'high certainty' where the landslide parts and morphology are clearly visible and easy to detect and 'low certainty' where the landslide parts are not clearly visible or are missing and the landslide morphology is unclear. The precision of mapping was expressed as 'high precision' where the landslide boundaries (crown, flanks, and toe) are clearly visible and 'low precision' where the landslide boundaries are smooth and unclear, i.e. where landslides are assumed based on a zone of depression and accumulation. It is possible that some very old or prehistorical slides were excluded during the landslide mapping process because the landslide morphology was not visible on the derivative LiDAR maps. Additionally, remediated landslides and very small landslides (<40 m 2 ) were not mapped because from the LiDAR DTM 0.3 × 0.3 m resolution, it was impossible to identify the typical landslide morphology.
Visually interpreted landslides on LiDAR DTM derivative maps were additionally field checked during December 2018 and January 2019. Systematic field mapping of all landslides from the inventory could not be carried out due to a relatively large area of the entire pilot area (21 km 2 ) and impassable parts of a terrain (neglected and overgrown land parcels). For this reason, landslide verification was conducted on 10% of randomly selected landslides in the inventory.

Landslide inventory map
The landslide inventory map was produced for a small part (21 km 2 ) of the Podsljeme area because of the time-consuming procedure of the computation of derivative maps and the manual delineation of landslide boundaries. In total, the landslide inventory consists of 702 landslides with a non-uniform distribution across the pilot area (Main Map). The prevailing types are shallow soil slides, and only 10 landslides were classified as earth flows. Nearly 65% of all the identified landslides on the LiDAR DTM were evaluated as 'high confidence' due to the clearly visible landslide features on the LiDAR derivatives, and nearly 60% of all mapped landslides were evaluated as 'high precision' due to the fresh and clearly visible landslide boundaries on the LiDAR derivatives.
Based on the landslide inventory, the total landslide area is 0.5 km 2 or 2.43% of the pilot area (21 km 2 ). The mean landslide density is 33.3 slope failures per square kilometres, and the smallest landslide identified in the pilot area has a planimetric area (A L ) of 43 m 2 (Bernat Gazibara, Krkač, Sečanj, & Mihalić Arbanas, 2017). The mapped landslides extend in size to a maximum of 8,064 m 2 , while the mean landslide area is approx. 700 m 2 (mean = 704 m 2 , median = 411 m 2 , std. dev. = 921 m 2 ). The most frequent landslides in the inventory have an area of approx. 400 m 2 , and 90% of the landslide bodies showed a size between 100 and 2,000 m 2 . The frequency-size distribution of all mapped landslides in the pilot area ( Figure 3) shows two scaling regimes: a positive power-law scaling for small landslides and a negative power-law scaling for medium and large landslides. The distribution obtained for the mapped landslides can be explained with the 'universal distribution' described by Malamud et al. (2004), which implies that the occurrence of a landslide in the pilot area is a result of either natural triggering factors (e.g. rainfall or rapid snowmelt) or anthropogenic interventions with regional impacts (e.g. widespread deforestation and urbanisation). The transition between the positive and the negative power-law relations can be used to distinguish between small and medium landslides. Based on the rollover at approximately A L = 400 m 2 , 48% of the mapped landslides are small (A L <400 m 2 ) and 52% are moderate and large (A L >400 m 2 ) in size.
Verification of LiDAR-based landslide inventory was performed for the 70 randomly selected landslides in the range from 68 to 5862 m 2 . Based on the certainty of visual landslide identification on LiDAR DTM, 39 randomly selected landslides (55%) were evaluated with 'high certainty' and 31 landslides (45%) with 'low certainty'. Field checking resulted with 15 confirmed landslides (21%) and 14 assumed landslides (20%) based on hummocky morphology, while 41 landslides (59%) were not accessible due to neglected and overgrown land parcels. Furthermore, results from field checking were compared with the certainty of visual landslide identification on LiDAR DTM and all 15 confirmed landslides were evaluated with 'high certainty'. In total, 24 landslides assumed and not accessible by field checking were evaluated with 'high certainty' while 28 inaccessible landslides and three assumed landslides were evaluated with 'low certainty' on LiDAR DTM.
The landslide density per land use unit (Table 1) shows that more than 80% of the mapped landslides appear in forests, 10% in agricultural areas, and only 7% on artificial surfaces. The frequency-size distributions of mapped landslides show that there is no difference in the landslide size among the land use units, but the mean landslide density in forests (84 ls/km 2 ) is six times higher than that in agricultural areas and on artificial surfaces (14 ls/km 2 ). One of the possible explanations is that LiDAR-based landslide inventories can often be incomplete on settlements and arable lands due to frequent anthropogenic influences (Bell et al., 2012;Petschko et al., 2015;Steger, Brenning, Bell, Petschko, & Glade, 2016). More probable is that the higher landslide density in forests is associated with prevailing steep slopes and forest gullies. In contrast, the urban areas in the pilot area are unlikely to be associated with the most favourable conditions for landslide initiation, as the majority of housing and infrastructure is constructed on relatively flat terrain, i.e. the prevailing slope angles (62%) in urban areas range from 0°to 12°. Furthermore, agricultural lands in the Podsljeme area are rarely cultivated, and landslide features tend to be preserved longer in the landscape. All reported landslides from winter 2012/ 2103 endangered residential buildings and roads, but most of the landslides appeared in forests and agricultural areas (Bernat et al., 2014b). Analysis of the landslide distribution in the pilot area shows that more than 50% of the mapped landslides are located within a distance less than 100 m from roads, and 43% of the mapped landslides are located within 50 m of buildings and residential houses. Based on the population density grid (1 × 1 km) for the Podsljeme area, approximately 43,000 people live within a zone of 50 m from 418  mapped landslides. It can be concluded that the landslides in the pilot area may affect further urbanisation and that elements at risk are highly exposed to landslide hazards.

Conclusions
For the first time, LiDAR DTM maps were used for detailed landslide mapping in the hilly area of the City of Zagreb, i.e. a highly urbanised area with a large number of known landslides and a non-existing landslide inventory map. Visual interpretation of the LiDAR DTM derivative maps resulted in a landslide inventory map of the pilot area (21 km 2 ) in the Podsljeme area with 702 landslides and a mean landslide density of 33.3 ls/km 2 . Based on the certainty of landslide identification and field verification, it can be concluded that the derivatives of the LiDAR-based DTM are very useful tools for mapping landslides in the Podsljeme area, despite the characteristics of the landslides (small and shallow landslide bodies), the high level of urbanisation (modification of natural landscape) and the dense vegetation cover (forested areas). The frequency-size distribution shows that the landslide inventory of the pilot area is substantially complete, in order that can be used as a representative landslide dataset for landslide susceptibility assessment. Although, the resulting LiDAR-based landslide inventory does not contain very small landslides (<40 m 2 ), remediated landslides and old, historic landslides lost due to erosional processes and land-use practices. Most of the mapped landslides can be dated as recently (re)activated due to the appearance or degree of preservation of the landslide morphology and the presence of smaller landslides in the inventory. Therefore, the landslide inventory map of the pilot area represents a combination of seasonal inventory (landslides (re)activated in winter 2012/2013) and historical landslide inventory (sum of landslide events activated over a period of tens or hundreds of years). A seasonal landslide inventory could not be completed because the landslides (re)activated in winter 2012/2013 were not visible on historical Google Earth satellite images from September 2012 to August 2013. Most of the mapped landslides are located under the forest, so one of the advantages of using LiDAR in the Podsljeme area is being able to map potential landslides under dense vegetation cover. Despite the small areas, landslides in the Podsljeme area cause significant damages to buildings, infrastructure, and crops because of the high landslide density and close proximity to roads, buildings and residential houses. Detailed and complete landslide inventory maps are directly useable by planners and decision makers as a basis for spatial development and land-use planning and as a requirement for site-specific investigations prior to building design and construction projects.

Software
The derivative maps from the DTM and the landslide inventory map were originally constructed using Arc-Map 10.0, while the end map was constructed in Corel Draw X6.