Kinematic model of the slow-moving Kostanjek landslide in Zagreb, Croatia

The interpretation of landslide kinematics provides important information for those responsible for the management of landslide risk. This paper presents an interpretation of the kinematics of the slow-moving Kostanjek landslide, located in the urbanized area of the city of Zagreb, Croatia. The sliding material (very weak to weak marls, often covered with clayey topsoil) exhibits plastic, rather than rigid behavior. Due to this reason, and low landslide velocities, landslide features, such as main scarps or lateral flanks, are barely noticeable or do not exist in most of the landslide area. The data used for the kinematic interpretation were obtained from 15 GNSS sensors, for the period of 2013-2019. The monitoring data revealed a different spatial and temporal distribution of landslide velocities, resulting as a consequence of geomorphological conditions and forces that govern the landslide movements. Temporally, eight periods of faster movements and seven periods of slower movements were determined. Spatially, velocities measured in the central part of the landslide were higher than on its boundaries. The interpretation of the surface (horizontal and vertical) displacements and the direction of movement reveal a new insight into the engineering geological model and provide important information for the management of the Kostanjek landslide risk.


Introduction
Landslides, generally defined as the movement of a mass of rock, debris, or earth down a slope (Cruden, 1991), play an important role in the evolution of a landscape (Crozier, 2010). Landslides include different materials, various types of movements and states of activity, as well as different velocities (Cruden and Varnes, 1996;Hungr et al., 2014). The combination of slow movement and plastic or rigid-plastic material behavior can form landslides without distinctive morphology, where landslide features are missing and landslide kinematics cannot be easily interpreted. Knowledge of kinematics helps to reveal temporally and spatially variable stresses acting within landslides, their boundary geometries, mechanical properties of materials composing landslides, external forcing conditions, and characteristics of future landslide movement (Schulz et al., 2017). Additionally, slow-moving landslides sometimes accelerate rapidly and fail catastrophically, causing widespread destruction and casualties (Lacroix et al., 2020), so knowledge about kinematics is important for those responsible for managing the risks posed by landslides (Glastonbury and Fell, 2008). Monitoring is one of the most important tools to understand landslide kinematics and dynamics (Angeli et al., 2000), while geomorphological mapping provides an additional methodology to evaluate and understand processes with surficial expressions that affect slope stability (Clayton, 2017).
The scope of this study is an interpretation of the Kostanjek landslide kinematics. Kostanjek is a slow-moving landslide, mainly composed of weak marls. It is located in the urbanized area of the city of Zagreb, the capital of the Republic of Croatia. The kinematic motion of the landslide was interpreted for the period of 2013-2019, mainly from continuous data gathered by the 15 GNSS sensors installed in the landslide area. GNSS provide high quality data related to precision and its temporal resolution, important for the interpretation of landslide dynamics and kinematics (Malet et al., 2002;Kavoura et al., 2020), as well as for the establishment of early warning systems (Wu et al., 2019). Additional data necessary for the kinematic interpretation of the Kostanjek landslide were obtained during detailed geomorphological mapping of sporadic morphological features appearing on the landslide. Directions and velocities of GNSS movements, with the interpretation of the morphological features, revealed new insight into the kinematics of the Kostanjek landslide, providing important information for the interpretation of the landslide

Methods
The Kostanjek landslide is the biggest landslide in the Republic of Croatia. Due to the importance of this landside (approximately 300 single-family houses and infrastructure networks are located on the moving landslide mass), a monitoring system was installed in the period of 2011-2014 (Krkač et al., 2017). The monitoring objective was landslide mitigation through the development of an early-warning system (Mihalić Arbanas et al., 2013).

Study area
The Kostanjek landslide is located in the city of Zagreb (see Figure 1), at the toe part of Medvednica Mountain. The landslide occupies an area of 1 km 2 , with an approximate volume of 32 × 10 6 m 3 (Stanić and Nonveiller, 1996). Its area is approximately 1,000 times bigger than the area of an average-sized landslide (730 m 2 ) in the city of Zagreb (Bernat Gazibara et al., 2019).
Kostanjek is a deep-seated translational landslide, with a sliding surface depth of up to 90 m (Ortolan, 1996). The landslide was caused by anthropogenic factors, i.e., excavation of marl at the toe part of the landslide used for cement production in the nearby cement factory (see Figure 1). In 1963, blasting techniques were introduced in the excavation of marl (Stanić and Nonveiller, 1996;Ortolan, 1996), and immediately after, damage to the factory buildings and other structures appeared (Ortolan and Pleško, 1992; Stanić and Nonveiller, 1996).
Excavation in the open-pit mine was stopped in 1988 after identifying mining activities as the main landslide triggering factor. A total volume of 5.3 × 10 6 m 3 of marl was excavated (Stanić and Nonveiller, 1996) in the toe part of the landslide. To this day, the surface displacements did not stop, and damage to the buildings and infrastructure progresses.
A summary of landslide investigations and data interpretations are presented in Ortolan (1996). The interpretation of the landslide model was performed based on the data collected by geodetic surveys, borehole data from 1931, 1972 and 1988 and data collected by geophysical surveys from 1989. Despite substantial damages in the landslide area, since its activation to 2012, there was no consistent geodetic monitoring of landslide movement. For the period of 1963-1988, historical data about landslide displacement were obtained from the interpretation of aerial photos (Ortolan and Pleško, 1992). The determined displacements were up to 6 metres. The difference in displacements (determined on approximately 110 points), i.e. the difference in velocities on the landslide surface, was interpreted as movement on three different sliding surfaces, the deepest at 90 m and two subparallel sliding surfaces at depths of 65 m and 50 m (Ortolan and Pleško, 1992). The positions of the sliding surfaces were defined on the basis of unfavorably oriented bedding planes in Miocene sediments. The displaced mass is composed of Sarmatian, Lower and Upper Pannonian sediments. The geological crosssection of the wider landslide area is presented in Figure 2. The dominating type of rock within the landslide area is marl, very weak to weak according to ISRM (2007), with low primary hydraulic permeability, from 10 -6 to 10 -8 cm/s. The bedding planes are dipping towards the SSE, approximately in the direction of the slope, with a dip angle up to 15° (the slopes within the landslide are generally gentle, except in the open marl pit in the central part, where steep slopes prevail). The Sarmatian period is characterized by regular rhythmical sedimentation of carbonate and clay-rich marl laminae. Marls sporadically contain thin layers of sandstone and limestone. The sliding surface of the landslide was partly formed in these deposits, along with the thin layer of clay (Krkač 2015), approximately 9 metres below the geological contact with overlaying Pannonian deposits (Ortolan 1996). The lower Pannonian deposits continuously overlay the Sarmatian deposits. The thickness of these deposits is approximately 30 metres. The main characteristic of this lithostratigraphic unit is a vertical alternation of weak marls and hard clayey limestones (Vrsaljko 2003). The upper Pannonian deposits consist of massive, often very weathered marls. This lithostratigraphic unit includes two lenses of coarse-grained sediments, which are chaotic, poorly cemented and sorted, with well-rounded grains. These two bodies probably originated as a consequence of slumps, i.e. synsedimentary submarine landslides (Vrsaljko et al. 2012).
Today, the most important triggering factors of the Kostanjek landslide reactivations are increased groundwater levels (GWL), and a consequent reduction of effective stresses. Maximal observed GWL depth (19 m), in the central part of the landslide, was measured in April 2014, while the minimal observed GWL (10.5 m) was measured in November 2013. According to Krkač et al.
(2019), the GWL depth is between 15 and 16 metres in the central part of the landslide, and triggers the movement of the landslide. The oscillations of the GWL near the landslide boundaries are significantly lower, up to 4 metres (Krkač 2015). The increased GWLs and landslide movements are the result of intensive precipita-tions and snow melt (Krkač et al., 2020). The amount of precipitation that causes GWL to rise depends on the soil moisture, above the water table, related to the wet and dry periods (Krkač, 2015). Generally, the climate of the city of Zagreb is characterized by warm summers and cold winters. The average annual precipitation measured in Zagreb (at the Zagreb-Grič meteorological station) is 889 mm (Krkač et al., 2020), with the maximal precipitations during the summer and autumn months. The minimal average precipitation occurs in February (45 mm), while the precipitation reaches its peak in June (96 mm). The second peak of monthly precipitation usually occurs in October (91 mm).

Monitoring of the Kostanjek landslide
The monitoring system at the Kostanjek landslide (see Figure 3) consists of multiple sensor networks for continuous observations of (Krkač et al., 2019): (1) external triggers (rain gauge and accelerometers); (2) hydrological properties (pore pressure gauges and water level sensors in boreholes and domestic wells, water level sensors at outflow weirs); (3) movement/activity (GNSS sensors, extensometers, borehole extensometers and inclinometer).
The monitoring of the landslide movement is mainly based on the GNSS network for measuring continuous surface displacements. The GNSS locations (see Figure  3) can be grouped as follows : above the landslide crown, i.e. outside the landslide area to check the assumption that this area is stable (GNSS 01); in the hinterlands of the open mine pit slope cuts (GNSS 04, 08 and 11); inside of the abandoned open mine pit (GNSS 05, 09 and 12); along the western and northwestern landslide boundary (GNSS 06, 07, 10, 13, 14 and 15); and in the northern part of the landslide (GNSS 02 and 03).
The GNSS data used in this study are 24-h post-processed. The precision of GNSS measurements, calculated as the root mean square error on the 24-h post-processing position (at 2σ, 95% confidence), is 3.2-4.6 mm in planimetry and 6.1-10.5 mm in altimetry (Krkač et al. 2017). The total temporal data coverage of GNSS since its installation is higher than 95% for most receivers, except GNSS 05, which has been continuously suffering from power supply issues.

Results
During the monitoring period (2013-2019) all GNSS sensors showed a significant displacement, i.e. the measured displacements were greater than the measurement errors, except for GNSS 01 which is located outside the landslide (see Figure 3). Cumulative horizontal and vertical displacement data from all 15 permanent GNSS stations in the landslide area are presented in Figures 4 and 5. As can be seen from these figures, the largest horizontal (650 mm) and the largest vertical (+412 mm) displacements during the six year period were measured with GNSS 09, which is located in the central part of the abandoned marl pit (see Figure 3). The smallest horizontal displacement (94 mm) in the same period was measured near the SW landslide boundary with GNSS 15, while the largest subsidence (-266 mm) was measured by GNSS 04, in the eastern part of the landslide. The total horizontal displacements measured with the GNSS sensors during the period of 2013-2019, near the landslide boundaries, were between 94 and 529 mm, and the total horizontal displacements in the central part of the landslide were between 555 and 650 mm. The difference in the measured displacements makes a difference in calculated annual velocities, between different parts of the landslide, approximately 1.2 to 5.5 times. The average annual velocities of all the GNSS sensors, at the Kostanjek landslide, were calculated from the cumulative horizontal displacements, for the period of 2013-2019, and interpolated over the entire landside area (see Figure 6). Additional information for the interpolation of landslide velocities was gathered by numerous field reconnaissance mapping during the monitoring period. Apart from the damage on the structures and infrastructure, from the field investigations it was not possible to observe any significant displacements along the landslide boundaries. The lack of clearly visible displacements near the landslide boundaries, as well as different velocities measured by the GNSS sensors, implies that the landslide moves significantly slower near its boundaries. The direction of movements across almost the entire landslide were consistent, in the direction of S-SSW (179-200°), except for GNSS stations 11, 12 and 13, which were moving in the direction of SSE (153-159°).
The total vertical displacements over the period of 2013-2019 ranged from -266 mm (subsidence) in the eastern part of the landslide (GNSS 04) to 413 mm (uplift) in the central part of the abandoned marl pit (GNSS 09). The other GNSS sensors did not measure such significant elevation changes. The total vertical displacements of the other 13 sensors are between -108 mm (subsidence) in the northern part of the landslide (GNSS 02) and 50 mm (uplift) at the toe part of the landslide in the open marl pit area (GNSS 05). The data reveals that most of the area in the upper part of the landslide exhibits subsidence, while the lower part of the landslide exhibits uplift, especially in the area of the abandoned marl pit (see Figure 7). The interpretation of vertical displacements was additionally confirmed by the results of field mapping. From these field investigations, it was impossible to observe significant vertical changes along with most parts of the landslide boundary. Clearly expressed vertical movements (uplift) were mapped only near the southern landslide boundary, where the landslide boundary intersects the access road to the aban- doned marl pit, and along the western landslide boundary, near GNSS 14 (see Figure 8a). Moreover, the greatest subsidence area was interpreted from the GNSS 04 data and reconnaissance mapping of the surrounding area. Namely, GNSS 04 is located in a depression, formed by the subsidence of the terrain along the subvertical fracture distanced only a few metres east from the GNSS (see Figure 8b).
In general, the periods of faster movement cannot be recognized from the vertical movement patterns, except for GNSS sensors 04 and 09 characterized by a significant amount of vertical displacement. Some GNSS sensors, such as GNSS 05 and 12, during the monitoring period displayed seasonal cyclic displacements, characterized by repeating patterns of up-and down movements, with about ±50 mm/yr amplitude for GNSS 05 and ±30 mm/yr amplitude for GNSS 12. The peak of the upward movement occurred during the summer months, while the downward movement occurred during the winter months.

Discussion
The Kostanjek landslide kinematics were interpreted for the period of 2013-2019, based on the reliable data gathered by the continuous measurements of 15 GNSS sensors. During the monitoring period, occasional reconnaissance field mappings were performed to confirm changes in landslide features in some parts of the area. GNSS data undoubtedly revealed the highest landslide velocities in the central part of the landslide and significantly smaller velocities near the landslide boundaries. This is mostly in agreement with the results of the photogrammetric survey of aerial photos conducted by Ortolan and Pleško (1992). They concluded that differences in velocity are the result of a mass movement on three different sliding surfaces, subparallel to the bedding planes. According to this interpretation, a sliding mass in the central part of the Kostanjek landslide, located above all three sliding surfaces, exhibits the highest velocities, while a sliding mass near the landslide boundaries, located above one (the deepest) sliding surface, exhibits the lowest velocities. However, the displacement vectors, from the open marl pit were not available in the study by Ortolan and Pleško (1992), so this area was interpreted as the area with the smallest velocities, which slides on the deepest sliding surface. GNSS monitoring from the period of 2013-2019 presented in this paper reveals the opposite, i.e. the area of the highest velocities is located in the abandoned marl pit. Different velocities within the landslide were observed on some other landslides, for example on the mudslide in Wealden Beds (Allison and Brunsden, 1990) or the translational landslide Utiku (Massey et al.,  2013). The probable reason for the different spatial distribution of the movement, between the landslide center and its flanks, is increased friction along the flanks, between the moving mass and the stable terrain, resulting in reduced velocities. The result of the friction between the moving mass and the stable terrain is additionally expressed as bulging in the western landside flank area, near GNSS 14, observed during field mapping. This area is outside of the bulging zone and theoretically should subside, but during the period of 2013-2019, it eventually uplifted 5 mm. The reduction in the landslide velocity near the main scarp probably occurs due to an extension of the sliding mass in this area as the landslide moves down the slope. The extension phenomenon also involves thinning of the sliding mass, subsequent subsidence, and sliding mass disintegration. The reduction of the landslide velocity in the landslide toe area can be explained by the compression of the moving mass in collision with the stable terrain. This compression results in an uplift of the terrain and with the disintegration of the landslide mass. Except for the difference in velocities due to the geomorphological conditions, the described spatial distribution of movements within the landslide can be additionally explained by the fact that groundwa- ter level changes in the central part of the landslide govern the Kostanjek landslide movement (Krkač, 2015). Due to this reason, the velocities are highest in the central part of the landslide, while much lower and relatively constant near the boundaries, especially near the western flank (GNSS 13 and 15), i.e. slower movements along the landslide boundaries are governed by the faster moving blocks in the central part of the landslide. The highest uplift occurs in the abandoned open marl pit area where data from GNSS 09 showed an uplift of 413 mm, almost ten times higher than GNSS 05, also located near the toe part of the landslide. The probable reason for this phenomenon is the movement of two different rigid-plastic blocks in the direction of the open marl pit. These movement directions are visible from horizontal vectors of displacement (GNSS 08 and 04 towards SSW and GNSS 11 and 12 towards SSE). GNSS 11 and 12 movement directions suggest that this part of the landslide mass, west from the open marl pit, is detached from the main body and is moving towards the SSE. The detachment of the block occurs along the fracture, in the hinterland of the western slopes of the open marl pit, expressed locally in the form of subsidence and damaged structures. The other block, east from the open marl pit, moves towards the SSW, as most of the sliding mass, but is detached from the main body by the distinctive subvertical fracture. West from the fracture, a significant amount of subsidence was measured (-266 mm) by GNSS 04, in the period of 2013-2019. Two blocks moving towards the open marl pit generate a zone of intensive compression, resulting in significant uplift and highest velocities in an open marl pit. In this area, the sliding surface is interpreted at a relatively small depth (Ortolan, 1996), so the load of the material does not provide sufficient resistance to active forces caused by two blocks from the east and the west.
Deformations in the landslide boundary areas express plastic, rather than rigid, behavior of the Miocene very weak to weak marls, covered by silty-clayey superficial deposits of various depths. Plastic behavior can explain a relatively small number of clearly expressed fractures in the landslide area, such as the subvertical fracture in the eastern part of the landslide, located near GNSS 04. Due to this reason, it is not possible to determine the boundaries of the two blocks moving towards the open marl pit, resulting with a zone of intensive compression. Most of the surface manifestations of landslide activity are expressed on the artificial structures, such as infrastructures and buildings, in the form of cracks, joints and inclination of structures, while on natural terrain, the deformations are mostly expressed as uplift and subsidence (Mihalić Arbanas et al., 2016). Thus, the landslide boundaries of the Kostanjek landslide cannot be mapped accurately based on the morphology of landslide features despite almost 60 years of its continuous activity. Additional reasons for difficulties in observing the landslide features are relatively low velocities, continuous reparation of damage on private houses and other infrastructures, and artificial changes of the morphology of the terrain.
However, when clearly expressed landslide features are missing at the surface, point data such as GNSS cannot provide a full picture of landslide kinematics. Additional information about surface movement on this kind of landslide, with minimal deformation features, can be obtained from remote sensing, such as Synthetic Aperture Radar interferometry (InSAR) or the Unmanned Aerial Vehicle (UAV) based photogrammetry method. These methods are often used to investigate and interpret landslides (e.g. Ohki et al. 2020; Godone et al., 2020). The main advantage of remote sensing is the ability to acquire spatially continuous data, even with centimetre precision, which can be very useful when they have to be integrated with conventional ground-based techniques (Tofani et al. 2013). Thus, remote sensing data (although not continuous) would provide valuable

Conclusions
The Kostanjek landslide is a slow-moving, large, deep seated, translational landslide, with very poorly expressed landslide features that prevent the interpretation of the landslide geometry and its kinematics based on detailed mapping. For this reason, the landslide was equipped with a network of 15 GNSS sensors to monitor surface movement, aimed at long-term reliable and precise monitoring. The data time series gathered by GNSS sensors, supplemented with occasional field checking of landslide features, enabled the interpretation of landslide kinematics during the period of 2013-2019 and evaluation of the landslide model. During the monitoring period, all 14 GNSS sensors within the known landslide boundary, measured statistically significant displacements. The highest measured horizontal displacement, up to 650 mm, was measured in the central part of the landslide, i.e. in the open marl pit area, while the lowest velocities were measured near the landslide flanks. The greatest measured uplift, up to 413 mm, was also measured in the open marl pit area, while the greatest subsidence was measured in the eastern part of the landslide, near the subvertical fracture located above the open marl pit. During the monitoring period, eight periods of faster movements and seven periods of slower movements occurred. The highest velocity (4.5 mm/day) was measured during the first week of April 2013, within the first period of faster movement.
The difference in velocities on the Kostanjek landslide is the result of geomorphological characteristics and forces that govern the landslide movements. Lower velocities along the landslide flanks occurred due to the friction between the moving mass and the stable terrain; due to an extension near the main scarp, where the moving mass gets thin and disintegrates; and due to compression in the toe part, where the moving mass uplifts and disintegrates. The highest velocities occur in the central part of the landslide where the groundwater level changes govern the Kostanjek landslide movement. The information about the directions of movement, uplift in the open marl pit area and subsidence in its hinterlands, revealed new insight into the Kostanjek landslide kinematics, important for understanding the landslide dynamics and engineering geological model. Knowledge of the kinematic model is also crucial for landslide management, planning an early warning system, and other elements, such as evacuation routes in case of risky landslide movements.