Spatial distribution and geometric characteristics of landslides with special reference to geological units in the area of Slavonski Brod, Croatia

A preliminary analysis of landslide spatial distribution and their geometric characteristics is presented for the area of Slavonski Brod, located in the northeastern part of Croatia and belonging to the Pannonian Basin System. A landslide inventory for the study area of 55.1 km2 is accomplished for the first time, based on the visual interpretation of a high resolution LiDAR digital terrain model. In total, 854 landslide polygons are delineated, corresponding to an average density of 15.5 landslides per square kilometre. The average landslide area is 839 m2, and most of the landslides can be classified as small landslides (76 %). The spatial relationship between landslides and geological units is analysed and expressed as a landslide index. The Late Pannonian sands with silts and gravel interlayers and Pliocene clay, sands, gravels, and coal are determined as the units that are most susceptible to landslide processes. The majority of landslides (85 %) are concentrated within these two units, for which a detailed analysis is performed, determining the morphometric parameters (slope and relief) and drainage network. The parameters’ classes that create favourable preconditions to slope instabilities are defined, based on the landslide density within individual classes. Besides, the geometric characteristics of landslides (size and shape) within these two units are compared. The results serve as the basis for further investigations. They help to foresee the area of future landslides through landslide susceptibility maps, and offer a better understanding of the influence of fluvial-denudation and slope processes on recent landscape evolution and form. pretation of stereoscopic aerial photographs (e.g., VAN DEN EECKHAUT et al., 2007; JABOYEDOFF et al., 2012; GUZZETTI et al, 2012; GÖRÜM, 2019). The application of visual interpretation techniques are tested for different landslide types and under a wide range of environmental conditions (e.g., ARDIZZONE et al., 2007; VAN DEN EECKHAUT et al., 2007; PETSCHKO et al., 2016; BERNAT GAZIBARA et al., 2019; PAWLUSZEK, 2019; GÖRÜM, 2019; JAGODNIK et al., 2020). Along with the increasing availability of HR digital elevation models, advances in Geographical Information System (GIS) have improved quantitative analysis of the earth’s surface, including landslide zonation studi es related to inventory, and susceptibility, hazard or risk assessment (VAN WESTEN, 2000; CHACON et al., 2006). Until now, in Croatia, LiDAR-based landslide inventories exist for smaller areas of the City of Zagreb and the Vinodol Valley. In the populated area of the City of Zagreb, 24 km2 (MIHALIĆ ARBANAS et al., 2016), and 21 km2 (BERNAT GAZIBARA et al., 2019) were examined. In the Vinodol Valley, 18.55 km2 were studied in the central part of the Dubračina river basin (TOŠEVSKI, 2013), and 0.48 km2 in the Slani Potok gully (JAGODNIK et al., 2020). The study area is one of the six pilot areas selected in the framework of the safEarth project (“Transnational advanced management of land use risk through landslide susceptibility maps design”), which was carried out from 2017 to 2019 under the Interreg IPA Cross-border Cooperation Programme CroatiaBosnia and Herzegovina-Montenegro 2014-2020. The main project goals were related to landslide susceptibility maps at different scales. For the project implementation, an ALS was performed, Article history: Manuscript received December 08, 2020 Revised manuscript accepted July 02, 2021 Available online February 28, 2022


INTRODUCTION
Landslides have been in the focus of many researchers in the field of engineering geology, geotechnics, and geomorphology for many years, especially because they present a geohazard event that can cause extensive damage and casualties (VARNES & IAEG, 1984;BELL, 2003;HAQUE et al., 2016;MIHALIĆ AR-BANAS et al., 2017), and significantly contribute to landform evolution (CROZIER, 2010;BELL et al., 2012). For landslide analysis related to both of these subjects, preparation of a landslide inventory map for a particular area is an initial step (GUZ-ZETTI et al., 2012), as it is believed that the locations of past and present landslides are the key to predicting future ones (VARNES & IAEG, 1984). A landslide inventory primarily records the landslide location and, when known, the type of movement, activity, date of occurrence, and other landslide characteristics (MALA-MUD et al., 2004;FELL et al., 2008).
Techniques of landslide mapping are constantly developing. In particular, detailed information about surface topography provided by airborne laser scanning (ALS), also known as airborne light detection and ranging (LiDAR), improved the investigation of the Earth's surface processes (JABOYEDOFF et al., 2012;TAROLLI, 2014). In landslide studies, the use of LiDAR has grown greatly during the last decade and has become a promising method for landslide mapping, monitoring, and modelling (DER- RON & JABOYEDOFF, 2010). Considering the LiDAR-based landslide inventories, many authors have highlighted the advantages of surface interpretation using bare earth high-resolution (HR) LiDAR digital terrain models (DTM) over conventional methods, i.e., geomorphological field mapping and visual inter-

Spatial distribution and geometric characteristics of landslides with special reference to geological units in the area of Slavonski Brod, Croatia
pretation of stereoscopic aerial photographs (e.g., VAN DEN EE-CKHAUT et al., 2007;JABOYEDOFF et al., 2012;GUZZETTI et al, 2012;GÖRÜM, 2019). The application of visual interpretation techniques are tested for different landslide types and under a wide range of environmental conditions (e.g., ARDIZZONE et al., 2007; VAN DEN EECKHAUT et al., 2007;PETSCHKO et al., 2016;BERNAT GAZIBARA et al., 2019;PAWLUSZEK, 2019;GÖRÜM, 2019;JAGODNIK et al., 2020). Along with the increasing availability of HR digital elevation models, advances in Geographical Information System (GIS) have improved quantitative analysis of the earth's surface, including landslide zonation studi es related to inventory, and susceptibility, hazard or risk assessment ( VAN WESTEN, 2000;CHACON et al., 2006).
Until now, in Croatia, LiDAR-based landslide inventories exist for smaller areas of the City of Zagreb and the Vinodol Valley. In the populated area of the City of Zagreb, 24 km 2 (MIHALIĆ ARBANAS et al., 2016), and 21 km 2 (BERNAT GAZIBARA et al., 2019) were examined. In the Vinodol Valley, 18.55 km 2 were studied in the central part of the Dubračina river basin (TOŠEVSKI, 2013), and 0.48 km 2 in the Slani Potok gully (JAG-ODNIK et al., 2020).
The study area is one of the six pilot areas selected in the framework of the safEarth project ("Transnational advanced management of land use risk through landslide susceptibility maps design"), which was carried out from 2017 to 2019 under the Interreg IPA Cross-border Cooperation Programme Croatia-Bosnia and Herzegovina-Montenegro 2014-2020. The main project goals were related to landslide susceptibility maps at different scales. For the project implementation, an ALS was performed, which resulted in HR DTMs for an area of the app. 300 km 2 in total. It gave the opportunity to accomplish landslide inventories over a wide area and increased the knowledge of landslide characteristics and their spatial distribution.
The pilot area presented in this paper (55.1 km 2 ) is located in the northeastern part of Croatia, in the hinterland of Slavonski Brod. It represents the southwestern part of the Pannonian Basin System (HORVÁTH & ROYDEN, 1981;ROYDEN, 1988) and belongs to the North Croatian Basin (PAVELIĆ & KOVAČIĆ, 2018). In this area, mostly Miocene, Pliocene, and Quaternary sediments are exposed on the surface (ŠPARICA et al., 1979a;ŠPARICA, 1986).
Most of the study area is covered by forest, and the most populated area, the centre of Slavonski Brod, is outside of the study area. However, the expansion of the settlements towards the north resulted in the highest city elevations of approx. 250 m. Landslide occurrence and damage to property have been recently documented after a period of heavy precipitation in March 2018. While severe damage occurred to private houses and local road infrastructure, a natural disaster was declared for the hilly areas of Slavonski Brod. This has led to increased interest in landslide processes and related problems by the local authorities.
Here, a LiDAR-based landslide inventory for the area of Slavonski Brod has been completed for the first time. The preliminary results presented aim to define the landslide spatial dis-tribution and geometric characteristics within different geological units (GUs). This was accomplished through several steps: (i) identifying and mapping landslide polygons using HR LIDAR DTM derivatives, (ii) correlating the landslide distribution with GUs identified within the study area, (iii) defining morphometric characteristics within the two most representative GUs, and (iv) examining in more detail the relationship between morphometric parameters and landslide occurrence, as well as differences in landslide geometric characteristics within these units.

Geographical and geomorphological settings
The study area is located in the northeastern continental part of Croatia (Fig. 1a) and covers an area of 55.1 km 2 . The area is situated northwest of the City of Slavonski Brod, the main centre of Brod-Posavina County (Fig. 1b), and covers the area of three municipalities, Slavonski Brod, Podcrkavlje, and Sibinj (Fig. 1c). Slavonski Brod is situated at the junction of two spatial units. The larger, southern part of the city area belongs to the lowland area along the Sava river, and the smaller, northern part, to the wider hilly belt, in the hinterland of which is Dilj Mt. with the highest peak of 461 m. Several suburban settlements are situated within the study area, where 11 % (app. 7,100 inhabitants) of the population of Slavonski Brod (BPŽ, 2019) reside in Podvinje and Brodsko Vinogorje. In these settlements, landslides were reported to local authorities in 2010 and 2018, when a large amount of precipitation was recorded.
According to the Climate Atlas of Croatia (ZANINOVIĆ et al., 2008), the climate of the Slavonski Brod area is characterized as a temperate continental climate, with harsh winters, short springs, and warm and humid summers. Based on a 30 year measurement period , and data recorded at the meteorological station "Slavonski Brod", in the study area the average annual air temperature is 10.7 °C. Average monthly temperatures rise until July when they reach a maximum of 21.0 °C. The avera ge annual amount of precipitation is 748.1 mm. Concerning the annual course of monthly precipitation, the main maximum occurs in July, and the main minimum in February. Still, certain years show a significant deviation in monthly amounts from the average precipitation conditions. Based on geomorphological regionalization (BOGNAR, 2001), the study area belongs to the mega-macro-geomorphological region of the Pannonian Basin, within which the hills of Dilj Mt. as a micro-geomorphological region stands out. In the study area, the elevation ranges from 94 m to 353 m, and most of the area (almost 73 %) is within the elevation class of 100-200 m (Fig.  1c). The elevation class of 300-350 m represents less than 1 % of the study area. The slopes vary from gentle to steep, up to 83°. The southern slopes of Dilj Mt. are very branched, expanding into a dense network of numerous gullies and valleys. Pediments and hills are formed by slope and fluvial-denudation processes (BOG-NAR, 1996). According to the area zonation by the degree of stability, on the Engineering Geological Map of SFRY in a scale of 1:500.000 (ČUBRILOVIĆ et al., 1967) a large part of the study area is defined as a "prevailing unstable area under natural conditions and under the man's activity mostly unstable".
GU8 -Alluvial sediments (Holocene) fill the stream valleys in the central and northeastern parts of the area and are composed of gravel, sand, and silt. These deposits originate from older sediments of the area, from the Badenian to the Pleistocene.
Based on the described GUs a simplified evolution of the study area can be proposed. During the Badenian in the deeper marine part (offshore) marls (Vejalnica Fm.), and in the shallow nearshore part lithothamnia limestones (Vrapče Fm.) were deposited. In the area of Dilj Mt. the Sarmatian is characterized by marine sedimentation in an environment with reduced salinity (Dolje Fm.). In contrast to the northern part of the Dilj Mt., in the study area, no sedimentation occurred during the Sarmatian. This is probably due to local tectonics during that time. Later, due to the Early Pannonian transgression, clayey limestones of the littoral environment begin to accumulate, followed by the Early and Late Pannonian marls of the Medvedski Breg Fm. The marls were deposited in a relatively deeper part of the newly formed brackish Pannonia lake (AVANIĆ et al., 2018a, b). During the Late Pannonian, marls with sandy interlayers accumulated in the prodelta environment (Andraševec Fm.), and sands with layers of silt and gravel of the delta front (Nova Gradiška Fm.), and during the Pliocene, sands, gravels, clays, and coals related to rivers and swamps. Sedimentation from the Lower to the Upper Pannonian has characteristics of a transgressive sequence, and from the Late Pannonian to the Pliocene of a regressive sequence. In the transgressive sequence, there was a deepening from the littoral to the lake basin, and in the regressive sequence, there was a shallowing or propagation of the river clastic system, when after the lake sediments, the deposits of prodelta, delta-front, alluvial, and swamp plains accumulated. The final terrestrial sedimentation occurred during the Quaternary and is represented mostly by loess, as well as by deluvial and alluvial sediments.

Input data
The HR DTM used in this study is a result of an ALS, which was performed in the early spring of 2018 by Slovenian company Fly-com Technologies d.o.o., fulfilling the requirements for at least 20 points per square meter and assuring the accuracy for each point of at most 10 cm in all directions. As a final product of raw LiDAR data post-processing, DTM with a 0.5 m resolution is obtained. At the same time, topographic scanning was performed resulting in HR orthophotos with a 10 cm resolution.
Based on geological reconnaissance and existing basic geological maps (ŠPARICA et al., 1979a;ŠPARICA, 1986), a geological map of the investigated area was compiled and individual GUs (FILJAK et al., 2016) were interpreted.

Landslide mapping
Noticeable morphological signs that remain after landslide occurrence can be recognized and mapped through the interpretation of digital topographic surface model (GUZZETTI et al., 2012). The methodological approach for landslide recognition is based on visual interpretation and geomorphic feature extraction (SCAIONI et al., 2014) using HR LiDAR DTM derivatives. For that purpose, several DTM derivatives most commonly used for landslide inventory preparation were extracted, using standard ArcGIS tools. These are: (i) hillshade maps (according to SCHULZ (2004) different sun azimuth angles of 315°, 135°, and 45°, and a constant altitude angle of 45° were used); (ii) slope map; (iii) profile curvature map; (iv) contour line map with elevation equidistance of 1 m (Fig. 3). DTM derivatives were used individually or in combination to delineate landslide polygons. Indi- Figure 2. Geological map of the study area modified according to 1:100,000 scale basic geological maps (ŠPARICA et al., 1979a;ŠPARICA, 1986), with landslide inventory.
vidual landslides are identified and digitized as polygons at a large scale close to 1:500, comprising the completely affected area from the landslide crown to the toe. When more different landslide generations were recognized, each landslide polygon was delineated separately whenever possible in order to obtain the real and representative landslide geometric characteristics. As most of the study area is covered with forest, the use of orthophotos was limited, although in some locations the landslide main characteristics were clearly recognized (e.g., main scarp, crown cracks, flank scarps) (Fig. 4).
All these morphological signs are generally sharper, better developed and more clearly visible for recently active landslides and are modified through time by different processes of weathering, erosion, and deposition (MCCALPIN, 1984). As the landslide ages, the more its morphology is reshaped and the less recognisable the landslide features become in the field (MCCALPIN, 1984;KEATON & DEGRAFF, 1996;BELL et al., 2012;PETSCHKO et al., 2014). Thus, the recognition of typical landslide characteristics is evaluated by several authors titled as landslide freshness ( VAN DEN EECKHAUT et al., 2007), confidence (BURNS & MADIN, 2009), certainty of mapping (PETSCHKO et al., 2016), or certainty of landslide identification (BERNAT GAZIBARA et al., 2019).
In this study, the landslide outline confidence level (LOCL) is presented as the degree of landslide contour certainty based on the visibility and expressiveness of the landslide main features and is expressed with scores from 1 to 10. Each of the following landslide features is rated with a confidence score from zero to two: (i) main scarp, (ii) lateral flanks, (iii) toe, (iv) internal landslide morphology (e.g. hummocky or undulating topography, internal cracks, minor scarps), and (v) the ability of discernment of depletion and accumulation areas. Based on the confidence scores, four LOCL classes are defined (Tab. 1).
A landslide inventory was produced as the basis for the analysis of landslide spatial distribution and geometric characteristics. All of the analyses were performed for landslides with a confidence score higher than two, i.e., for LOCL classes II, III, and IV. In addition, the minimum landslide area for mapping was set to   Landslides with high precise delineated boundaries 50 m 2 , given the scale of the input DTM. In the analysis performed, no distinction was made between the different landslide types, or according to their depth and relative ages. The LiDAR-based landslide inventory was verified during a field engineering geological inspection, conducted in March 2019.

Landslide inventory analysis
Preliminary analysis of the landslide inventory was performed, which aimed to define landslide size, shape, and spatial distribution in the study area.
Geometric parameters used to analyse landslide size and shape are presented for every landslide polygon with several features that were determined automatically (area), and manually (length, width) using standard ArcGIS tools, or calculated by equation (aspect ratio).
The relationship between GUs and landslide occurrence within the study area is evaluated and represented by the landslide index. The landslide index, expressed in percentages, is calculated as a ratio between landslide area in a given GU and the area of that unit. The landslide index of each GU is compared to the landslide index of the entire study area, also expressed in percentages and calculated as a ratio between the total landslide area and total study area (BARTELLETTI et al., 2017).
In terms of statistics, the correlation between landslides and morphometric parameters (slope, relief, distance to streams) of two major GUs (Fig. 2) is examined and expressed through landslide density for each parameter class. Individual landslide density for each parameter class is compared to the landslide density of each GU, which is calculated as a ratio between landslide number within an individual GU and an area of that unit. For that purpose, centroids of landslide polygons are created.
The slope angle was automatically derived from DTM resampled to 5 m resolution. In order to observe landslide distribution from the slope angle perspective, the terrain was classified according to DEMEK (1972), whose classification is based on the dominant processes that activate on a slope of a particular angle. Six slope classes (expressed in degrees) are used: 0-2 (plain terrain), 3-5 (slightly sloping terrain), 6-12 (sloping terrain), 13-32 (very sloping terrain), 33-55 (very steep terrain), and > 55 (cliff).
Relief is expressed as a difference between the maximum and minimum elevation in a certain area and is calculated using the ArcGIS tool Focal statistics. The radius of the circle taken around each cell to calculate the elevation range is set to 250 cells, i.e., 125 m. The six classes (expressed in metres) are defined automatically using natural breaks and are: 0-18 (flattened area), 19-34 (very low relief), 35-47 (low relief), 48-60 (moderate relief), 61-75 (high relief), and 76-109 (very high relief).
To calculate the distance to streams, a detailed stream network was automatically extracted from the DTM resampled to 5 m resolution, using the Arc Hydro tools (v2.0). The standard procedure included filling sinks, computing the flow direction, and flow accumulation. Stream definition is based on a flow accumulation and a specified threshold value, which represents the number of cells indicating the start of a stream. The different threshold areas were tested and a representative network was chosen based on visually cross-checking the main streams using a topographic map at a scale of 1:25,000. However, strong topographic diversity influenced the selection of different threshold areas that are associated with the various GUs. The stream order was automatically calculated using the method proposed by STRAHLER (1956).
The landslide inventory and spatial analysis were performed using ESRI ArcGIS 10.2. The official Croatian coordinate projection, HTRS96/TM, is used for all maps presented herein. The visualization of graphs was aided by Daniel's XL Toolbox add-in for Excel, version 7.3.2 (KRAUS, 2014).

Landslide inventory
In the study area, a total of 854 landslide polygons were delineated, corresponding to an average landslide density of 15.5 landslides per square kilometre. The total mapped landslide area is 0.72 km 2 (Tab. 2), which represents 1.31 % of the total area. The boundaries for three landslides exceed the study area, but the polygons are completely delineated within the DTM buffer zone, and their actual sizes are used for landslide area statistics. Thus, the total study area affected by landslides is, therefore, smaller by 0.004 km 2 .
The landslide areas range from 50.43 m 2 to 17,855.32 m 2 , with an average value of 838.93 m 2 (Tab. 2). The distribution of landslide area is presented in more detail via the frequency density and a cumulative number of landslides (Fig. 5). It is shown that 3.5 % of landslides are smaller than 100 m 2 , and 20.5 % are larger than 1,000 m 2 . Most landslides (76 %) range between 100 m 2 and 1,000 m 2 , of which 54.4 % are in the range of 100 m 2 to 500 m 2 . Within this range, the highest frequency density of 1.48 m 2 is presented for landslide areas ranging from 200 m 2 to 300 m 2 . Only five landslides are larger than 10,000 m 2 .
Considering different classes of LOCL, the number of landslides decreases as LOCL increases (Tab. 2). Thus, almost 60 % of landslides mapped are described with the LOCL class of II, and a little over 40 % of landslides with LOCL class of either III or IV. As LOCL gives an overview of the expression of landslide features, these results show that in the study area, all typical land- slide features are not completely visible for most of the landslides (60 %) and the boundaries are particularly hard to determine unambiguously. It can be interpreted as either landslide activation occurred a long time previously explaining why their sharp boundaries and morphological signs have been masked over time, or the landslides are younger but their main features are lost due to fast and strong erosion processes. In addition, the average landslide area increases as LOCL increases (Tab. 2), which is especially a result of a great variability of the smallest landslides within each LOCL class (the smallest landslide for LOCL IV is almost 7.5 times larger than the smallest landslide for LOCL II).
The landslide inventory is verified in the field by checking delineated polygons and identifying the morphological features typical of landslides. The verification is conducted for landslides with different confidence scores (from 3 to 8), and a range of areas between approximately 70 m 2 and 7,600 m 2 , covering the area within three different GUs (GU3, GU4, and GU5). From 97 landslides visited, 94 were verified, which represents 11 % of the Li-DAR-mapped landslides. The remaining three landslides are all located in the forest and are marked with LOCL of II, i.e. confidence score of 3, but in the field, typical landslide morphology was not confirmed.

Landslides and geological units
It is recognized in many landslide studies that the geological setting greatly influences landslide occurrence. As one of the main predisposing environmental factors, geological data are used as one of the main input layers in landslide susceptibility assessment ( VAN WESTEN et al., 2008). Thus, for the landslide spatial distribution analysis, first, the relationship between landslide occurrence and GUs is investigated. The relationship is expressed with the landslide index (Fig. 6), and the landslide statistics for each GU are presented (Tab. 3). The landslide index of each GU was compared to the landslide index of the entire study area (BAR-TELLETI et al., 2017), enabling determination of which units are more correlated to landslide source areas.
Considering the landslide area distribution for different GUs, it can be seen that landslides are mainly concentrated within GU4 and GU5, where 85 % of landslides are present (Fig. 6). Within GU4, 423 landslide polygons are delineated, and within GU5, 303 landslide polygons, which corresponds to landslide densities of 31.66 and 18.73 landslides per km 2 , respectively. These GUs refer to Late Pannonian sands with silts and gravel interlayers (GU4) and Pliocene clay, sands, gravels, and coal (GU5). The landslide index in these units is higher than the landslide index of the study  ; GU3 -Marls with sandy interlayers (Late Pannonian); GU4 -Sands with silts and gravel interlayers (Late Pannonian); GU5 -Clay, sands, gravels, and coal (Pliocene); GU6 -Loes (Pleistocene); GU7 -Deluvial sediments (Holocene); GU8 -Alluvial gravels, sands, and clay Holocene). The landslide index for each geological unit and for the study area is presented.
area (1.30 %) at 2.03 % and 2.11 %, respectively. At the same time, as already previously stated, these units are most represented in the study area, covering almost 30 km 2 (54 % of the study area). However, the highest landslide index (4,13 %) is calculated for GU1 (Badenian limestones and marls), representing only 21 landslides. Consequently, the area covered by GU1 is only 0.6 % of the study area. For all other GUs, the landslide index is less than 1, i.e., below the landslide index of the study area (Tab. 3).
According to the average landslide area, the largest landslides occur in GU5, and the smallest in GU4 (Tab. 3). The differences in landslide average width are not so noticeable, with values ranging from 20.05 m to 25.61 m. The values of landslide average length, compared to others, for GU5 and GU6 stands out, with values of 51.55 m and 51.86 m, respectively.

Morphometric analysis
The diversity of GU4 and GU5 environments is represented by their morphometric features, including slope angle, relief, and drainage network (Fig. 7).
Although the range and average values of the slope angles within GU4 and GU5 are similar, the slope angle class distribution shows differences (Fig. 7a). Almost 90 % of the area of GU4 is within three classes in the range of 6° to 55°, with the most representative class of 13-32° (54 %). Most of the area of GU5, almost 95 %, is within three classes in the range of 3° to 32°, with the most representative class of 6-12° (49 %). Slope angle class with values above 55° is the least represented, with an area below 1 % for both GUs.
The values of relief for GU4 and GU5 range from 5 m to 109 m, and from 9 m to 100 m, respectively (Tab. 4). For GU4, 81 % of the area is within three relief classes in the range of 48 m to 109 m, with the most representative class of 61-75 m (34 %). For GU5, 87 % of the area is within three relief classes in the range of 19 m to 60 m, with the most representative class of 35-47 m (44 %). For both GUs the relief class with values below 19 m is least represented covering approx. 1 % in both GUs (Fig. 7b).
Concerning the classes defined to analyse the distance to streams, 65 % of the area of GU4 is within the first three classes, with distances to streams less than 60 m. Within this unit, most of the area (25 %) is in the stream buffer zone of 0-20 m. For GU5, all classes are almost equally represented, with significant representation by class (43 %) where the distance to stream is more than 100 m (Fig. 7c). These results indirectly indicate that the stream network within GU4 is more branched and denser than within GU5. This is also shown by the results of the analysis of stream length for different stream orders (Tab. 4). Since the LiDAR polygon does not include the entire basin, only streams that are completely developed (streams of a certain order and all lower orders) in a single GU, are considered. There are half as many 1 st order streams within GU4 than GU5, although the total length of 1 st order streams does not differ much between these two units. Consequently, there is a big difference in their average stream lengths, 87.70 m within GU4 and 164.06 m within GU5. Thus, GU4 can be described as a unit with a short but dense series of 1 st order streams, and GU5 as a unit where 1 st order streams are rarer but much longer. The same trend is presented for 2 nd order streams. Furthermore, the difference in the stream number and stream lengths decreases for 3 rd order streams, and is even less for 4 th order streams.

Landslide spatial distribution
In two zones characterized by different geological and geomorphological conditions, landslide spatial distribution is analysed. Classes of parameters that create favourable preconditions to slope instability, presented in the previous sub-chapter, are defined based on landslide density (Fig. 8).
Considering the slope angle, landslides within GU4 are mainly concentrated within classes of 13-32° (52 % landslides) and 33-55° (45 % landslides). The highest landslide density is also related to these classes, with the highest peak of 115 landslides per km 2 within the class of 33-55° (Fig. 8a). Within GU5, a large number of landslides is also related to classes of 13-32° (64 % landslides) and 33-55° (12 % landslides). Still, landslides are also concentrated on somewhat milder slopes, within a class of 6-12° (23 % landslides). The highest landslide density is related to the class of 33-55°, with 79 landslides per km 2 (Fig. 8a), same as for GU4. The differences in the relationship between slope angle and  landslide occurrence can be expressed by the average values of slope angle within landslide polygons and the rest of the area. Thus, for GU4 the average slope angle within areas affected by landslides is 27°, and for GU5 is 18°. For the area not affected by landslides, the average slope angle for GU4 is 17° and GU5 is 12°. In a view of relief, 96 % of landslides within GU4 are concentrated within classes with relief values higher than 48 m, with the highest concentration within the class of 61-75 m (45 % landslides). The highest landslide density is also related to these classes (Fig. 8b). Within GU5, the highest landslide concentration is related to classes with relief values from 35 m to 75 m (64 % landslides). The highest landslide density is, however related to the class of >76 m (Fig. 8b). The differences in the relationship between relief and landslide occurrence expressed by the average values of relief within landslide polygons and the rest of the area are more significant compared to the slope angle. For GU4 the average relief within areas affected by landslides is 65 m, and for GU5 it is 51 m. For the area not affected by landslides, the average relief for GU4 is 61 m and for GU5 it is 45 m.
Considering the distance to streams, 92 % of landslides within GU4 are concentrated within the classes 0-20 m (62 % landslides) and 20-40 m (30 % landslides). The highest landslide density is also related to these classes, with the highest peak of 78 landslides per km 2 within the class of 0-20 m (Fig. 8c). Within GU5, 81 % of landslides are concentrated within classes with a distance below 60 m. Still, the highest landslide density is related to the class of 0-20 m, with a value of 52.6 landslides per km 2 (Fig. 8c).

Landslide geometric characteristics
Several landslide characteristics are considered, including area, length, width, and aspect ratio (Table 5), in order to compare landslide size and shape within GU4 and GU5.
Comparing landslide sizes, the average landslide area in GU4 (642.24 m 2 ) is almost half the size of the average landslide area in GU5 (1,127.97 m 2 ). The range of landslide areas within GU4 and GU5 is quite different (see box plots, Fig. 9a) though the outliers are not represented, thus the maximum landslide areas are not shown.
The difference in landslide width is not as significant as landslide length, as was the case with the analysis of landslide inventory for each GU (Tab. 3). The average landslide lengths for GU4 Figure 8. Landslides and landslide density (number of landslides per area of each parameter class) for geological units 4 and 5, concerning: slope angle (a); relative relief (b); and distance to streams (c). Landslide density for each geological unit (number of landslides per area of the geological unit) is also presented. and GU5 are 34.10 m and 51.55 m, respectively (Tab. 5). This difference influences not only the landslide area but also the landslide shape. The average aspect ratio is thus slightly larger for landslides within GU5. The average aspect ratio for GU4 and GU5 is 1.88 and 2.27, respectively (Tab. 5). The aspect ratio classes indicate that for both GUs most landslides occur within two classes and landslides are mainly longitudinal. Further, it can be stated that an isometric shape is more associated with landslides within GU4 and an elongated shape with landslides within GU5 (Fig. 9b).

DISCUSSION
Landslide inventories are the basis for assessing landslide susceptibility, hazard, and risk assessment (SOETERS & VAN WESTEN, 1996;VAN WESTEN et al., 2008). Without knowledge about landslide locations, the methods for susceptibility models that predict landslides based on past conditions remain inapplicable. The landslide inventory accomplished for the study area includes landslides that occur over time (periods of tens, hundreds, or even more years), and are not associated with a trigger event. Thus, this inventory can be interpreted as a historical geomorphological inventory (MALAMUD et al., 2004;GUZ-ZETTI et al., 2012).
The landslide inventories are widely used to analyse landslide statistics, including their size, shape, and location (MAL-AMUD et al., 2004). Considering the landslide size and shape, geometricl parameters and their ratios can be further analysed to characterize landslides and assess their volumes, motion mechanisms, runout distances, and hazard (TIAN et al., 2017). Considering the landslide distribution, it has been proven that landslides mostly occur at preferred locations. Among others, geologicallithological and morphometric parameters have been used as the main environmental predisposing factors to analyse the effects on landslide spatial distribution ( VAN WESTEN et al., 2008;REICHENBACH et al., 2018). The distribution of landslides concerning the main factors controlling shallow landslides (lithotype, regolith thickness, slope angle, and land use) is discussed for the Picentino river basin in southern Italy (LUPIANO et al., 2019). CONFORTI & IETTO (2020) highlighted that the lithology, fault density, slope gradient, and relief have an important role in determining landslide occurrence. GUZZETTI et al. (2008) examined the relationship between landslide distribution and local morphological and lithological settings, in the Upper Tiber River basin in central Italy. Their analysis confirmed that geological conditions influence the abundance of landslides in the catchment. JACOBS et al. (2016) also concluded that landslide processes largely depend on the prevailing lithology. In their study, slope angle is distinguished as the main controlling topographic factor for landslides. BARTELLETTI et al. (2017) investigated the influence of geological-morphological and land use settings on the occurrence of shallow landslides. They used the landslide index to identify significant parameters correlated with landslide source areas. In lithologically homogeneous areas, (as each of the presented Gus may be considered), COCO & BUC-COLINI (2015) indicated that the slope morphometry, especially angle, height, and form were the most influential factors on landslide occurrence. They also confirm that landslides are coupled with a drainage network.
In this study, the geological setting is proven to be the significant influencing factor that controls landslide occurrence and abundance. Besides analysing the number of landslides within individual GUs, the correlation between GUs and landslide occurrence is estimated and expressed with the landslide index. The highest landslide index is defined for the unit represented by the Badenian limestones and marls (GU1). Since the area covered by GU1 represents only 0.6 % of the study area where only 21 landslides are mapped, further investigations are needed in areas where this unit is much more representative. The next two units for which the landslide index is above the landslide index of the study area refer to Late Pannonian sands with silts and gravel interlayers (GU4) and Pliocene clay, sands, gravels, and coal (GU5). As 85 % of the mapped landslides are concentrated within these units, they can be considered to be the most susceptible to landslide processes in the study area.
When taken into consideration the uniform geographical, climatic, and tectonic settings for GU4 and GU5, it can be reasonably assumed that the geological setting is the prevailing factor in the evolution of this area. Significant differences in morphometric features are strongly correlated to the geological  setting, i.e. lithological features and material characteristics, which are clearly visible on Fig. 7. Still, it is necessary to emphasize the scale of the input of geological data (Basic geological map in scale 1:100,000), which may affect the presented results in the boundary areas of GUs. This can also explain the fact that 2.69 % of landslides is placed within GU8 (alluvial gravel, sands, and clay), as presented in Table 3.
The unified results of detailed analysis performed for GU4 and GU5, elaborating the relationship between morphometric features and landslide source areas, and geometrical characteristics of landslides, are presented in Table 6.
Observing the classes of morphometric parameters analysed, it is evident that the predominant classes for slope and relief within these two units differ for one to two class categories (Tab. 6). Considering the distance to streams, the dominant classes differ even more. Significant morphometric characteristics within GU4 refer to high relief, steep slopes, the dense stream network, and branched gullies. Dominant distance from stream classes are less than 60 m. Morphometric characteristics within GU5 refer to low relief, less steep slopes, and a less dense stream network, when compared to GU4, with dominant distance from stream classes greater than 100 m.
Regardless of these differences, parameters' classes defined as the most influenced on landslide spatial distribution are almost equal for both of these units. The classes for which the landslide density of parameter class is higher or very close to the landslide density of the entire GU are singled out. Based on these results, landslides are predominantly occurring on very steep and very sloping terrain, with very high, high, and moderate relief, and are strongly coupled to streams (Tab. 6). The fluvial-denudation processes, together with slope processes have the dominant influence on the development of local landscape. The geomechanical properties of sands are a dominant factor for current landscape formation in GU4. Their drainage capabilities and low cohesion causes the development of a dense gully network due to strong erosion processes. Consequently, the majority of the landslides are associated with the steep gully banks. Typical landslides are small and shallow, often without clearly visible landslide elements in the foot area where the landslide material is removed by permanent or intermittent streams. In contrast to GU4, low drainage properties of GU5 govern landscape shapes in their respective area, where the soil complex comprises impermeable clay interlayers or lenses. Landslides are predominantly triggered by high groundwater pore pressures. Landslides are slightly larger and more elongated, when compared to landslides within GU4, which implies the flow of material when the landslide is activated. The differences in landslide geometric characteristics within these units implies that geological characteristics and consequent morphometric features influence not only the landslide density and spatial distribution but also the landslide size and shape.

CONCLUSION
In this paper, preliminary analysis of the initial LIDAR-based landslide inventory for the area of Slavonski Brod is presented for the first time. Using HR LiDAR DTM derivatives (hillshade, slope, curvature, and contour line maps), the geomorphological historical landslide inventory for an area of 55.1 km 2 was completed. In total, 854 landslide polygons were delineated, corresponding to an average landslide density of 15.5 landslides per square kilometre. The landslide area ranges from 50.43 m 2 to 17,855.32 m 2 , with an average value of 838.93 m 2 , according to which they can be considered as small landslides.
Given the strong relationship between landslide occurrence and geological setting in many studies, first, the spatial landslide distribution was evaluated, according to eight GUs defined in the study area. The correlation between GUs and landslide occurrence is estimated and expressed with the landslide index. According to differences in the landslide index, the geological setting can be considered as a significant feature in landslide susceptibility assessment. Late Pannonian sands with silts and gravel interlayers (GU4) and Pliocene clay, sands, gravels, and coal (GU5) are those units where the majority of landslides are concentrated. 85 % of landslide polygons mapped are present within these units. Besides, these units are predominantly represented in the study area, covering an area of almost 30 km 2 . The landslide index for both of these units is above the landslide index of the study area. According to all the above, GU4 and GU5 were chosen for detailed investigation, in order to analyse the landslide spatial distribution and compare landslide geometrical characteristics.
Defining the relationship between landslide distribution and environmental parameters enables a better understanding of landslide characteristics and occurrence, and can help to predict the future landslides. Parameters analysed in this study (slope, relief, and distance to stream) are proven to influence landslide spatial distribution, and classes with a stronger correlation with landslide density are defined.  The results of landslide spatial distribution compared to geological features, slope morphometry, and stream network, provide useful information for further studies related to landslide susceptibility and hazard. In that sense, the basic information of significant predisposing factors controlling landslide occurrence in the study area is given. Together with the landslide inventory created, a large proportion of relevant input data is prepared. Those data are of great importance and can be used to produce statistically-based landslide susceptibility maps and to predict areas of higher landslide susceptibility. It would allow the local authorities to direct detailed engineering geological investigations and financial resources toward landslide-prone areas, highly relevant for the hilly area of Slavonski Brod.