Integrated approach for landslide hazard assessment in the High City of Antananarivo, Madagascar (UNESCO tentative site)

The High City of Antananarivo is one of the most important cultural heritage sites of Madagascar, on the UNESCO Tentative List since 2016. Built on the hilltop of the Analamanga Hill, a granite ridge overlooking the Ikopa River valley, it is renowned for its baroque-style palaces, such as the Rova royal complex, and neo-Gothic cathedrals dating back to the nineteenth century. During the winter of 2015, the twin cyclones Bansi and Chedza hit the urban area of Antananarivo, triggering floods and shallow landslides, as well as causing thousands of evacuees and many casualties. Between 2018 and 2019 several rockfalls occurred from the rock cliffs of the Analamanga hills, destroying housings and killing over 30 people. Both events showed that landslides can pose a high risk to the safety of the inhabitants, the infrastructure, and the cultural heritage of the High City. To assess the landslide hazard in the Analamanga Hill area, an integrated approach was adopted by means of the following actions: (i) creation of a multitemporal detailed scale landslide map; (ii) geotechnical characterization of the involved materials; (iii) analysis of landslide susceptibility in soils/loose deposits; (iv) runout analysis of debris flows channeling within large creek gullies; (v) landslide kinematic analysis of the rockmass; (vi) simulation of rockfall trajectories; (vii) analysis of rainfall data. The results show that the main factors affecting landslides are slope, lithology, creek-gully erosion, and anthropization, while most of the landslide events are clearly triggered by heavy rainfall. The landslide-prone areas (the phenomena include shallow landslides, rock falls, and debris flows) are located primarily along the cliff bounding the western hill slope, the southeastern sector (where abandoned quarries form large slope cuts), and subordinately in the steep creek catchment just east of the Rova. The thematic maps produced represent fundamental land use management tools to be used in Geo Disaster Risk Reduction (GDRR) by scientists, practitioners and the decision-makers involved in the High City protection and conservation. The study conducted represents an important contribution for improving the knowledge on landslide processes in an area with limited data such as Madagascar, and may be reproduced in cultural heritage sites characterized by similar geomorphological and urban scenarios.


Introduction
Disaster risk due to natural events is often much higher in developing countries, where phenomena such as fires, landslides, earthquakes, and floods usually lead to serious consequences with a tragic periodicity, due to the lack of resources and investments in understanding the associated hazards and risks (Petley 2012;Andriamamonjisoa and Hubert-Ferrari 2019). This is mainly due to the fact that developing countries are often located to a great extent in tropical zones, largely affected by extreme climatic events (such as cyclones, typhoons), and in active volcanic-seismic areas, which both in turns can trigger catastrophic flooding and landsliding (Alcántara-Ayala 2002). These countries are often affected by economic, political, social, and cultural issues that increase their susceptibility to natural hazards, causing casualties, loss of fertile soils, and pastureland, as well as severe damage of property and infrastructure, including tangible cultural heritage (Devoli et al. 2007). During the last decades, the population growth in Africa has nearly tripled (UNDESA/PD 2012), leading to the abandonment of rural areas and the uncontrolled urbanization of urban and suburban areas (Cobbinah et al. 2015). This is a challenge particularly in Madagascar, where the thick lateritic soil cover (regolith) resulting from a deep chemical weathering of the crystalline bedrock and the high rainfall intensity in the rainy season often result in flooding, landslides, gullying, and soil erosion (Thomas 1994;Lateef et al. 2010;Cox et al. 2010;Rakotobe et al. 2016). In this framework, the highest price for these disasters is mainly paid by the population belonging to the lower classes, especially those who are often forced to live in dilapidated buildings and in areas very exposed to natural disasters, also due to anthropological interventions, faulty or even non-existent urban development plans, and little investment for improving landslide risk perception among the local population (Petley 2012;Rakototsimba and Bettencourt 2010;Andriamamonjisoa and Hubert-Ferrari 2019). Even the cultural heritage of Madagascar is very often subject to the danger of natural events, also due to lack of prevention plans and mitigation interventions.
The scientific community can effectively support authorities and decision makers by producing landslide hazard assessment, to be used in development plans and/or for the implementation of the appropriate risk mitigation policies (Corominas et al. 2014). In recent years, there has been a constant growth in the awareness of landslide-related risks, which produced interesting results in the field of landslide detection and mapping, and in the methods for hazard assessment (Catani et al. 2005(Catani et al. , 2013. The literature background on the topic is extensive and hard to summarize; exhaustive overviews are provided by Brabb (1984); Varnes and IAEG (1984); Hansen (1984); Gupta and Joshi (1990); Carrara et al. (1991); Fell (1994); Soeters and van Westen (1996); Turner and Schuster (1996); Original Paper Aleotti and Chowdhury (1999); Dai et al. (2002) and van Westen (2004) ;Fell et al. 2008;Rossi et al. 2010;Catani et al. 2013;Fressard et al. (2014); Reichenbach et al. (2018); Canavesi et al. (2020) and Ozturk et al. (2021). All these authors agree that the key for a proper landslide risk assessment is knowing past events, to be achieved by direct or indirect surveys; in this framework the first step is to implement a landslide inventory map (Varnes 1984). According to Varnes and IAEG (1984), the standard procedures for landslide risk assessment are based on the assessment of hazard, vulnerability, and exposure. Hazard in particular, is defined as the probability of occurrence in a given area and time period of a slope instability phenomenon of a given intensity (Catani et al. 2005). Nevertheless, predictions based solely on spatial probability of occurrence are very common, since it is relatively simpler to achieve them: in such cases the term "landslide susceptibility" is adopted (Dai et al. 2002).
The High City of Antananarivo developed in the sixteenth century as the first core of the city from the top of the highest of the rocky ridges dotting the Ikopa River alluvial plain: the Analamanga Hill. Hosting several important churches, cathedrals, and palaces (among which the Rova royal complex is the most important from a cultural-religious point of view), it represents one of the principal cultural heritage sites of Madagascar. In the winter of 2015, the twin cyclones cyclone Bansi and Chedza hit the ar ea of Antananarivo, triggering floods and landslides (BNGRC 2015;Frodella et al. 2020) causing 80 fatalities and about $40 million of damage, equal to the 0.4% of the GDP (Andriamamonjisoa and Hubert-Ferrari 2019). In the following years, other landslide phenomena such as debris flows and especially rockfalls caused over 15 casualties (BNGRC 2019). These events highlighted the need to recognize the different geomorphological processes acting in the High City to assess landslide hazard and identify the areas more subject to risk. In this work, as a first step towards landslide hazard mitigation and a disaster risk reduction strategy for building resilience of the community and implement a sustainable management of the High City and related structures (roads and pathways), an integrated methodology was implemented by means of the following actions: (i) creation of a multitemporal detailed scale landslide map by means of field and remote sensing data; (ii) geotechnical characterization of the involved materials; (iii) analysis of landslide susceptibility on soil-loose deposits using the Random Forest algorithm (Breiman 2001); (iv) runout analysis of debris flows channeling within large creek gullies by using the SAGA-GIS Gravitational Process Path (GPP) model (Wichmann 2017); (v) landslide kinematic analysis of the rockmass using the Diana-k tool (Gigli and Casagli 2011;Gigli et al. 2012); (vi) simulation of rockfall trajectories using the ISGeoMassi code (CDM 2019); (vii) assessment of landslide temporal frequency by integrating the past landslides occurrence and the rainfall data analysis.

Study area
The High City of Antananarivo (18.55′ South; 47.32′ East) is situated in the Central Highlands of Madagascar on flat top of the Analamanga Hill, a granitic ridge elevating about 200 m above the Ikopa River plain (length of 3 km and a maximum height of 1436 m a.s.l.; Ciampalini et al. 2019) (Fig. 1). The site was chosen as the first capital of the Kingdom of Madagascar in the eighteenth century and therefore represents the first original core of the city. In the nineteenth century, during the French colonial period, the High City expanded from the hilltop to the hillslopes (Middle City), reaching the alluvial plain (Low City). New stone brick palaces, high dignitaries' buildings, churches, and chapels were built, such as the simil-baroque Andafiavaratra palace, and the gothic Cathedrals of Andohalo and Ambohipotsy, and especially the new Royal Palace (the Rova) with its royal chapel and tombs (Frodella et al. 2020) (Fig. 2). Therefore, the High City of Antananarivo is enrolled in the UNESCO World Heritage tentative list since 2016 ("La Haute Ville d'Antananarivo"; https:// whc. unesco. org/ en/ tenta tivel ists/ 6078). http:// whc. Unesco. org/ en/ tenta tivel ists). Its core zone encompasses the High City, while the buffer zone includes the Analamanga Hill slopes and the surrounding terraced rice fields (Fig. 1c). Nowadays, the Antananarivo urban area spreads over a surface of more than 80 km 2 , exceeding 1.7 million inhabitants, its urbanization having deeply influenced the natural landscape around the Analamanga Hill. The geomorphological setting of the Analamanga Hill is characterized by asymmetric slopes (Ciampalini et al. 2019;Frodella et al. 2020;Fig. 3): (i) a steeper western hillside, showing a subvertical granite rock cliff in the shoulder-backslope sector, while the colluvial footslope is represented by low energy surfaces formed by residual lateritic soils; (ii) an eastern hillside with a general minor slope angle (mean values around 30°) locally affected by linear creek erosion and sheet-overland flow erosion; (iii) a southern sector characterized by slopes around 40°, where diffuse and intense quarrying activity has carved several large niches; (iv) a northern urbanized slope gently linking with the river plain (except for a tunnel road cut). These western colluvial footslopes are severely affected by "lavakas", a large gully erosional feature typical of the Central Highlands of Madagascar (Ramifehiarivo et al. 2017). The natural hydrographic pattern of the hill shows small creeks with a radial-like pattern, often characterized by linear erosion (Fig. 3).
Geologically, the Analamanga Hill is formed by a migmatiticgranitoid batholite with its weathering products (Lee et al. 1989;Kroener et al. 2000;Dewandel et al. 2006) (Fig. 4). According to Frodella et al. 2020, these are formed by the following sequence, from bottom to top: (i) granite fresh rock ("F"), cropping out in the steep rock walls of the western backslopes, and in the dismissed quarries of south-eastern footslopes; (ii) laminated and slightlymoderately weathered granite (= saprock; "SW-MW"), overlying the bedrock especially on top of the rock cliffs; (iii) from highly to completely weathered granite (= saprolite; "HW-CW"), discontinuously covering the hilltop and the eastern hillside (according to boreholes and geophysical surveys they show a thickness from 5 to over 20 m in the center of the hilltop; Hyvert 1994; Andriamamonjisoa and Hubert-Ferrari 2019; GASC'ART Forage 2019); (iv) completely decomposed residual soils (= laterite; "RS"), encircling both the hill's footslopes; (v) loose sandy eluvial/colluvial cover (= "SC"; slightly pebbly sands), irregularly cropping out especially in the hilltop and the eastern hillside. The western and eastern slopes are covered by fan-like shapes of slope deposits (= "SD"), formed by cobbles, blocks and scattered boulders in a coarse sandy matrix (Fig. 4a). Loose heterogeneous anthropic detrital deposits (= "AD"), cover the above-mentioned sequence in the more intensively urbanized areas. In the river plain organic-rich clays and silts outcrop (alluvial deposits = AL). The granite bedrock is affected by three main discontinuity sets: J1 (185/80) and J2 (120/80), representing sub-vertical discontinuities with WNW-ESE and NE-SW trends (both sub-parallel with the direction of the main creeks); J3 (280/40), representing granite exfoliation discontinuities sub-parallel to the slope.

Materials and methods
A multi-temporal detailed scale landslide map was obtained using performed a heuristic approach (Fig. 5). This involved (i) geological/geomorphological field surveys (Frodella et al. 2020); (ii) geotechnical laboratory analyses on soil and rock samples collected during the surveys; (iii) interpretation of remote sensing data, among which very high resolution (VHR) multispectral Pleiades images (0.5 m of spatial resolution, from 2015 to 2019), and two VHR DEMs: a 2 m spatial resolution acquired in February 2015, and a second 1-m Pleiade DSM acquired on October 2016); (iv) homogenization in GIS environment and interpretation based on expert knowledge. On this basis, a landslide hazard assessment was performed in the complex geomorphological setting of the Analamanga Hill by integrating several analyses: (i) shallow landslide susceptibility; (ii) earth-debris flow runout simulations; (iii) rockmass 3D kinematic analysis; (iv) rock fall runout simulations. Finally, rainfall data from 2015 to 2019 and a landslide temporal inventory, obtained by means of surveys with the local population, were analyzed, so as to identify the rainfall conditions associated with the triggering of slope instabilities. According to the field surveys, the landslide processes on the Analamanga Hill developing in the soil and loose deposits involve only a few meters of thickness; therefore, these phenomena will also be referred as "shallow landslides".

Geotechnical analysis
To perform a preliminary geotechnical analysis of the bedrock and deposits cropping out on the Analamanga Hill, 10 granite rock and 10 soil samples were collected from the site and analyzed in the laboratories. In detail, for the rock samples, the following activities were performed according to ISRM (1985) recommendations: (i) classification of the samples; (ii) unit weight determination; (iii) point load test; (iv) tilt test. Point load tests were performed to examine the uniaxial compressive strength (UCS) and the uniaxial traction strength (UTS) of the analyzed samples (ISRM 1985). The rock samples were sawn to obtain 10 specimens with a rectangular parallelepiped-like shape (maximum and minimum mean width values of 81 and 35 mm, respectively, while maximum and minimum thickness of 31 and 46 mm, respectively). The granite rock shows a phaneritic equi-granular crystalline texture, with brownish oxidation coatings; therefore, value of mi (Hoek-Brown rock mass constant) = 29 was selected for the granite intact rock (RocData 2020). The basic friction angle was also determined by means of the tilt tests executed on artificially sawn discontinuity surfaces on the rock samples. A total of 3 samples were analyzed. For each sample from 3 to 9 tests were performed and the lowest value was considered to represent the analyzed rock basic friction angle. Regarding the soil samples, the following tests were performed following the ASTM International Geotechnical Engineering Standards: (i) classification of the samples; (ii) water content and unit weight determination; (iii) grain size analyses (wet sieving/sedimentation); (iv) shear test. The evaluation of the shear strength under drained and consolidated conditions was carried out using Casagrande's direct cut device. The consolidation pressures used for the analyses The adopted integrated approach for landslide hazard assessment in the study area (50 kN, 100 kN, and 150 kN) were chosen in relation to the on-site slope instability problems. During the test, in addition to the shear strength, the horizontal, and vertical deformation of the specimens were also measured.

Multitemporal detailed scale landslide map
The detailed scale landslide map landslide inventory was organized following the approach suggested by Soeters and van Westen (1996): mapping from the available DEMs and satellite imaging and performing a field survey and validation. Considering the landslide shallow depth verified in the field, the volume was considered equal to the surface (Catani et al. 2005). On these bases, a multitemporal detailed scale landslide map was obtained from 2015 to 2020, involving the winter 2015 shallow landslide event and other more recent single events (with special regard to the rockfalls that occurred in 2018 and 2019). The shallow landslide boundaries were digitized in a GIS environment by combining DEM and satellite image interpretation with two field GPS campaigns carried out in 2017 and 2019. During the latter, blocks belonging to the 2018-2019 rock fall events were also mapped. While the landslide type was described using the internationally accepted classification by Cruden and Varnes (1996), the DEM allowed measurement of both landslide magnitude (perimeter, area, volume) and morphometric basic parameters (maximum length, width, top of the headscarp, tip of the toe). In 2020, a new smaller landslide inventory, represented by a point-type shapefile, was carried out by integrating field surveys, GPS-RTK mapping, interviews with local population and checking possible news of the event on the web. Finally, the outputs of both databases were confronted and integrated with two technical reports of the local Civil Protection (BNGRC 2015(BNGRC , 2019, showing maps, location, description, and helicopter-based pictures of the landslides in the period 2015-2020.

Shallow landslide susceptibility
The Random Forest model (Breiman 2001), a machine-learning algorithm for nonparametric multivariate classification, was used to create the shallow landslide susceptibility map. Among its advantages, this method allows the adoption of categorical/numerical variables and the exploration of many explanatory variables, excluding assumptions regarding the data distribution (Catani et al. 2013). The machine-learning algorithm was therefore fed with many input parameters: elevation, slope, aspect, flow direction, flow accumulation, profile curvature, planar curvature, topographic wetness, and stream power indexes (obtained from a 1-m resolution Pleiade-1 DEM), lithology and NDVI map (which was used to discriminate forested area from bare soil/rock/urbanized areas). Landslide and geological data were obtained from Ciampalini et al. (2019) and Frodella et al. (2020). Besides the aforementioned parameters, a random values (0-1 interval) map was created and used as input to verify both the quality of the trained model and the presence of overfitting issues. Random parameter is not linked in any way to the presence/absence of landslides, so in a good model, this parameter must be the least important among the input data. To map landslide susceptibility, a starting database, made both of landslide and non-landslide pixels, was created. To train the Random Forest model, not all the available pixels can be used, since this can lead to overfitting issues (a model trained in such a way will be able to identify only known landslide); for this reason, the database of pixels was divided into 2 subsets, both containing landslides and non-landslide pixels. One subset was then divided into 2 parts: one containing the 70% of the pixels of the subset (used to train the model) and one containing 30% of the pixels of the subset (to test the trained model) and each part was defined to have landslide and non-landslide pixels equally distributed (50-50%). The trained model was then applied to all the data.

Earth-debris flows runout simulations
The intensity of potential earth-debris flows was assessed using the GPP model, an open-source code working in a SAGA-GIS environment (Wichmann 2017). It is used to simulate gravitational processes through the possibility of estimating the path of the material, the source, and the accumulation areas, respectively. The GPP model is a conceptual model which integrates components for process path determination, run-out, sink filling and material deposition, making it configurable for different processes like rockfalls and avalanches, and debris-mud flows. In detail, the model consists of stochastic models (Random walk, Markov chain, Monte Carlo simulation), physically based models, and empirical approaches aimed at simulating the path of moving masses on a DEM (Wichmann 2017). The model works from the regional to the local scale and requires different input data (DEM, friction angle, source areas) to assess the intensity of the modelled phenomena, combining both modelled flow thickness and velocity. For each component, various modeling approaches can be chosen by the user, such as geometric gradient/Fahrboeschung (Heim 1932) and shadow angle (Hungr and Evans 1988); mono-bi parametric model (Scheiddeger 1975;Perla et al. 1980); Sink filling (Gamma 2000); on stop (Wichmann 2017); slope, and-or velocity based (Gamma 2000). Based on field evidences the GPP model was applied in correspondence to the large gullies where earth-debris flows may occur, which were considered as source areas, while the outcomes of the laboratory analysis were used as geotechnical input data.

Rockmass 3D landslide kinematic global analysis
The kinematic analysis allows to highlight the rock wall sectors which are more prone to landsliding, by defining when a specific instability mechanism is kinematically feasible, given the slope geometry and the orientation of the rock mass discontinuities (Gigli et al. 2012). In hazard studies, this approach can be adopted to locate the landslide-prone areas of the rock mass . This approach analyses the instability mechanisms such as plane failure (pf) (Hoek and Bray 1981); wedge failure (wf) (Hoek and Bray 1981); block toppling (bt) (Goodman and Bray 1976); flexural toppling (ft) (Goodman and Bray 1976;Hudson and Harrison 1997). Casagli and Pini (1993) created a "kinematic hazard index" for each instability mechanism, by calculating the poles and intersection of the discontinuities located in critical areas of the stereographic projection. The kinematical hazard index is calculated as follows: where N and I are the total number of poles and intersections, respectively.
The use of specific software can allow for the analysis of a great number of discontinuities characterized by different friction angles. Based on the friction angles of the intersecting planes and the shape of the wedges, the lines of intersection and the equivalent friction angle are calculated automatically. Discontinuity dip and dip direction measures collected during the field surveys by means of geological compass were analyzed with RocScience Dips © software (Rocscience 2020). These were used as input data together with a 3D mesh obtained from the Pleiade DEM and (Fig. 4e). The discontinuity shear strength was considered purely frictional, and the mean peak friction angle was precautionally set to φ = 37°, based on lithological considerations (RocData 2020). The analyses of the discontinuity orientations were performed by employing the DiAna-k software (Gigli et al. 2012(Gigli et al. , 2014. The code expands the principles of kinematic analysis to overhanging slope sectors, and the susceptibility with respect to the different instability phenomena. The Global Kinematic Index (GKI) obtained comprises the spatial probability of occurrence of all the mechanisms combined, and therefore represents the areas with higher probability of rock block detachments.

3D rock fall analysis
In the study area, the rock fall analyses were be performed by using a 3D rock fall code named ISGeoMassi (CDM 2019). The code uses a hybrid lumped mass method (see Dorren 2003 andLi andLan 2015) for an exhaustive review of 3D probabilistic rockfall models) associated with a statistical analysis, allowing to model the variability of the input parameters for simulating the block impact behavior (here considered as a material point). The flight phase is regulated by the laws of dynamics (neglecting friction with the air), while the impact and the rotational slip are schematized with reference to the energy restitution (normal = Kn, and tangential = Kt), and the friction coefficients between the blocks and the slope (D°). These values regulate the interaction between the blocks and the ground, and need to be related to different ground cover polygons, in which the slope was subdivided in an ArcGIS environment using as a reference the Pleiade-1 image ( Table 1). The software simultaneously calculates many hypothetical trajectories (that will not affect each other) of the same isolated block descending along the slope. The slope geometry is represented by a mesh of triangles, and the rebound analysis and rotational slip on each triangle is performed with reference to the plane that contains it. The analysis of the various phases of the block motion (flight/rebound/rotational slip) continues until the block stops. This occurs when the energy and/ or the translational speed falls below a certain value of threshold, customizable by the user. In the need to outline a complex phenomenon, considering the uncertainty of the parameters that govern the analysis (topography, block-terrain interaction, initial conditions, etc.), the code uses a statistical model: some parameters including block size, initial speed, return coefficients, friction angle, roughness can be associated with a normal distribution, defined by mean and standard deviation values.

Rainfall data analysis
The approaches to evaluate the rain condition associated with the initiation of landslides can be separated into two main categories: deterministic and statistic (or empirical) approaches. Deterministic approaches are highly reliable but require a large amount of data and can be applied only to single slopes or small basins, although in recent years some attempt at developing spatially distributed models have been made (Salvatici et al. 2018). Statistical and empirical approaches usually require less data, as for instance rain data and landslide triggering location and date, and are mainly based on the search for a direct relationship between rainfall and landslides, without directly considering other parameters (Segoni et al. 2018). Among the statistical approaches, those based on the rainfall thresholds definition are the most widely used. In detail these are mainly based on the identification of the rain conditions associated or not with landslides, considering the duration (D) the cumulated rainfall (E) and the mean intensity (I) of a rainfall event, according to the model first proposed by Caine (1980). Some authors tried to use other parameters, such as soil moisture (Abraham et al. 2021), antecedent rainfall (Aleotti 2004), a combination of several rainfall scenarios (Segoni et al. 2018) or developing 3D thresholds, made by the combination of short-and long-term rainfall (Rosi et al. 2021). To evaluate the probability of a landslide event, the analysis of the rainfall conditions associated to their triggering was performed: information on landslide events from 2015 to 2020 was collected as well as the daily rainfall data of the same period (Climate Engine 2015). In this work, only 12 dated landslide events between 2015 and 2020 were available from the landslide inventory, so the rainfall conditions associated to their initiation have been analyzed by using a 2-steps approach: at first the cumulated rainfall and the mean rainfall intensity over different rainfall duration prior each landslide have been analyzed to identify the most influent rain condition on landslide triggering; then a I-D threshold, based on daily rainfall data, was defined, according to the procedure described in Abraham et al. (2021).

Geotechnical analysis
The results of the point load test performed on the rock samples (Table 2) show a mean value of UCS and UTS of 161.17 and 5.37 Mpa, respectively, falling within medium-high values according to literature (Hoek and Brown 1980;Hoek et al. 1992). The results of the tilt test report an average value of basic friction angle φ b of 31°, well matching the values of moderately rough granitoid rocks (Buocz et al. 2010;Alejano et al. 2012). Regarding the soil samples, a classification was carried out using the USCS (Unified Soil Classification System) system. From a geotechnical point of view the eluvial-colluvial deposits (SC) are formed by silty sands to well sorted sands with medium-high porosity (40%) and 37.6° of friction angle, while the laterite deposits (RS) can be classified as silts and fine sands to well sorted sands with low porosity (13.6%) and 33.9° of friction angle ( Fig. 6; Tables 3 and 4).

Multitemporal detailed scale landslide map
The obtained detailed scale landslide map (Fig. 7) is formed by two products: (i) a first inventory representing the shallow landslides triggered by the winter 2015 event, consisting of 78 shapefile polygons with attribute tables showing location (latitude/longitude of the centroid, perimeter, area, maximum length/width, crown/toe elevation a.sl.); (ii) a second point-inventory of 35 rockfalls/slides and debris flow sources are reported, including a description on location, type, volume, date of occurrence, town quarter, n° casualties/injured, description of the involved damage (Table 5). Four main landslide phenomena were detected (Frodella et al. 2020), according to which the involved materials are grouped as follows: (i) rock (F-SW-MW); (ii) debris (HW-CW-AD); (iii) earth (RS): 1) Unchanneled shallow earth-debris rotational/translational slides (Fig. 7b). These phenomena have small volumes and reduced mobility. These landslides are mainly situated both in the western slope of the ridge, in two creek catchments of the central-southern sector, and in the eastern slope in its centralsouthern sector, at the apex of the main gullies (Fig. 7). Since by field evidence these involve the first one-two meters of weathered bedrock, these landslides are represented by shallow phenomena.
2) Channeled earth-debris rotational/translational slides. These phenomena affect the laterite soils mainly in the main gullies of the western hillside (Fig. 7c), and occur following intense rainfall events during the cyclones may evolve into channeled earthdebris flows.
3) Rock falls and slides (planar and wedge mechanisms), from decimetric to decametric dimension (Fig. 7d). These form on the rock cliffs, both in the western and southern hillside. 4) Complex rotational/translational slide phenomena, affecting rock and debris (Fig. 7e). One of these phenomena, occurred in January 2019 in the Andohalo quarter, involving a pathway and the underlying earth/rock deposits, evolving as a debris fall-avalanche which destroyed a few housings and causing several casualties.

Shallow landslide susceptibility
Once the susceptibility model was successfully trained and tested (AUC = 0.9), it was applied to the whole study area. According to the results, the most important parameters to identify landslide susceptible areas were the NDVI, slope orientation (aspect), and slope curvature (both planar and profile curvature), while the less important were the flow accumulation and, as expected, the control parameter made of random values. The final output is a raster with a 2-m cell size, where each pixel has a percentage value expressing the landslide spatial probability of occurrence. The susceptibility values span from 0 to 1, and are grouped in 4 classes, according to the ArcMap "natural breaks" method: low, moderate, high, and very high (Fig. 8a). Roadways, pathways, together with the cultural heritage buildings shown in Fig. 2 were also reported. The obtained map shows how most of the study area is characterized by moderate to very high landslide susceptibility. The shallow landslide-prone sectors (high-very high susceptibility) are located along the Middle City quarters, particularly at the foot of the central sector of the western hillslope (from north to south: Mahamasina, Ankadilalana, Tzimbazaza quarters), where shallow debris rotational/translational slides involve mainly the eluvial-colluvial and the lateritic deposits (Fig. 8b). Here, the large gullies are rapidly expanding, as indicated by their retreating scarps damaging the road pavement (Fig. 8c), and the presence of several potential slope instabilities at the headcut of their apex (Fig. 8d). The eastern slope is characterized by a lower landslide susceptibility, except for two creek basins in its middle-southern sector: the first is located east of the Rova in the Manjakamiadana quarter, while the second is located east of the Ambohipotsy quarter (Fig. 8e). The northerm hill slope sector shows a general stability due to the widespread urban cover and the lower slope angles, except for the area around the entrance of the roadway tunnel through the hill, in the Ambohijatovo quarter. The southern sector reveals some potentially unstable areas at the top of the dismissed quarry slope cuts in the Anjahana quarter.

Earth-debris flow assessment
As for landslide susceptibility, four hazard classes were obtained based on the number of paths on each cell of the Lidar, which range from 1 to 250,779. Each cell value was then rescaled to 100 using the ArcMap tool "raster calculator". Therefore, four impact probability classes were obtained using the "natural breaks" function (considering a standard deviation equal to 6.56 and a mean equal to 1.5): low: 0-1.2; moderate: 1.2-3.5; high: 3.5-10; very high: 10-100. Shallow earth-debris rotational/translational slides developing from the gullies' slopes and head-cut apexes, following intense rainfall events can evolve into earth-debris flows that may channel within the hydrographic network. These latter phenomena occur more likely in the gullies of the central and southern hill sector like Mahamasina (Fig. 9a), Ankadilalana (Fig. 9b), and Tzimbazaza quarters (Fig. 9c). The earth-debris flow map obtained shows how these phenomena involve the lateritic soils as well as the anthropic waste deposits (since often the gullies are used as dumps), enhancing the gully head-cut erosion as well as the retreat of the scarps (Fig. 10a-f).

Rockmass instability and rock fall analysis
The 3D kinematic analysis was performed where the rock mass widely outcrops: the quarry front in the hill's south-eastern sector (Manakambahiny-Ankerakely quarter; "sector 1"), and the rock cliffs cropping along the western slope (from north to south: Mahamasina, Tsimialonjafy, Ambanin'Ampamarinana, Ankadilalana, and Tzimbazaza quarters; "sector 2") ( Fig. 11). In both areas, the 2D kinematic analysis shows that the most relevant instability   mechanisms detected in the field are topples (Fig. 11a, d), followed by wedge (Fig. 11b, f) and plane failures (Fig. 11c-f). The maps in Fig. 11g, h show the 3D spatial probability of occurrence (expressed in GKI %), for the combined above-mentioned instability mechanisms, with colors ranging from green which represent stable areas to yellow (low GKI) and red (maximum GKI = 26%). The quarry fronts in Fig. 11g show that the most unstable areas are in its central-southern sectors, while in correspondence to the western slope rock cliffs the most unstable areas are instead homogeneously distributed along the steepest sectors (Fig. 11h). These unstable areas (selected considering a value of GKI over the 95th percentile of its distribution curve) were used as sources for the rock fall trajectories, together with the source areas reported by BNGRC (2015BNGRC ( , 2019, and the slope areas exceeding a steepness of 70°. A new methodology was implemented in a CloudCompare software environment (DGM 2015), by interpolating on these source areas contour lines at a selected distance, based the operator's judgement. Along these contour lines, the code selected more potential block sources. The rock fall analysis was calibrated using the outcomes of the field surveys (block arrivals in Fig. 7), which allowed to recognize how the rock cliffs are frequently affected by rockfalls of various size (from decametric to plurimetric; Fig. 12d-g). The blocks' trajectories show that a large part of both sectors 1 and 2 are affected by the run out of rock blocks. Sector 2 is affected by rockfalls all along the foot of the rockcliff, from the Mahamasina quarter to Noth to Ankadilalana area in the area central sector, to the Tsimbazaza quarter to the south (Fig. 13a, c, e).

Preliminary landslide frequency assessment
A preliminary analysis on the landslide frequency, as reported by the surveys carried out among the local population and the rainfall data from 2015 to 2020, showed that the main landslide events occurred after major cyclones during periods of intense precipitations (Fig. 14). This analysis was followed by the definition of the most representative rainfall events connected to landslides, based on the cumulated rainfall associated to each landslide over different time intervals. As shown in Fig. 15a the cumulated rainfall associated with landslides grew constantly from 1 to 7 days, while the increase of rainfall from 15 to 30 days is lower than the one from 7 to 15 days, even if it is referred to a longer time interval; this could suggest that longer rain events have a lower influence on landslide triggering. The analysis of mean rainfall intensity (Fig. 15b) helped to highlight some differences between the rain events, in particular the intensity referred to 1-day duration shows rather scattered values, while the distribution of intensities for longer durations is more constant. This analysis also revealed that the rainfall intensity of events with 15 or 30 days duration is sensibly lower than that of shorter events.
On the basis of these analyses, it is possible to hypothesize that rainfall events ranging from 2 to 7 days are the ones mostly responsible for triggering landslides. The definition of I-D (intensity duration) rainfall thresholds is usually based on hourly rainfall data, but a recent paper (Abraham et al. 2021) explored the possibility of using daily data for the same purpose. To identify the threshold, rainfall data of events associated with landslides have been used; in particular, the intensity (mm/day) of events with a duration between 2 and 7 days has been considered, based on the outcomes of the previous analyses. For each landslide event, the maximum intensity has been extracted, along with its duration and the I-D couples have been plotted as shown in Fig. 16. In this figure, it is possible to notice that some rainfall events follow a common trend, while others, characterized by lower intensity, appear random. Due to the paucity of data, no robust statistical analysis was possible, so the threshold was defined as the lowest limit of the events characterized by a common trend, obtaining the equation:

Conservation problems of the High City and constraints of the analysis
Flooding and rainfall-triggered landslides are widespread hazards for the communities living on hillslopes of the High City. The rapid urban development and the lack of urban planning of Antananarivo has determined several environmental issues, such as intense slash and burn deforestation, illegal quarrying, slope cutting and terracing, garbage dumping, a lack of proper sewer and drainage systems, narrowing, culverting, and swamping of the creeks (Ciampalini et al. 2019). This urban pressure did not consider the fragile geological and geomorphological context of the Analamanga Hill, making the area (particularly the slopes and footslopes) particularly prone to landslide phenomena (e.g., shallow landslides, earth-debris flows and rock falls). This instability is particularly evident throughout the heavy cyclonic rainfall events often affecting the central highlands, as shown by the winter 2015 and the 2018-2019 landslide events (Frodella et al. 2020). In this framework, the expected future increase of extreme meteorological events can only exacerbate soil erosion and landslide events. The identification of the winter 2015 landslide processes in the field was not easy, given the rapid landform modifications caused by anthropic activity (the survey was carried out at the end of 2017). Nevertheless, the integration with high-resolution data (dating back to 2015 and 2016) and the helicopter-based pictures taken directly after the event (BNGRC 2015) allowed to detect the landslide processes with a satisfactory resolution. The rockmass 3D kinematic analysis surely suffered from a lack of resolution since the adopted 3D model was the 1-m Lidar DEM. The rockfall simulation unfortunately did not consider the volume of the involved material, this was due to the removal of the rock blocks and their use for construction materials. Future activities will have to involve the acquisition of a high-resolution 3D model of the rock cliffs, as well as a more accurate inventory of the rockfall events. The paucity of data influenced the temporal probability assessment as well, since it must be developed with a simple empirical approach, partially based on the expert judgment. Even with limitations, a validation of the rainfall threshold has been attempted, checking how many rainfall events (whether associated or not with landslides) overcame the thresholds in the analyzed period. This analysis showed that in 5 years the threshold was overcome 28 times, 10 of which are associated with landslides; this means that 18 rainfall events passed the threshold, but without triggering landslides. At first glance, this result might could look unreliable, but in the same period 30 landslides with unknown triggering dates were reported (some of them may have been triggered on the same day). Furthermore, in the same period of 298 days with at least 1 mm rainfall were reported. The definition of a rainfall threshold could be used to estimate the probability that a rainstorm is likely to cause landslides. If proper countermeasures are taken, a consequent reduction of social losses can be expected.

Identification of the areas characterized by the higher landslide susceptibility
The shallow landsides connected with the 2015 event affected the UNESCO core zone only slightly, except for its southern area, in correspondence of the pathway joining the Rova palace with Ambohipotsy Cathedral (Figs. 8 and 9). The cultural heritage buildings are not in landslide-prone areas, except for the Trano-Gasy Houses, which are situated in an unstable area on the topmost western slope, between the Rova and the Andafiavaratra Palace, as revealed by the cracks in the pavement and the surrounding signs of slope instability (Fig. 17). The most exposed linear structures are the tracks running on both mid-slopes in the central-southern sector of the Analamanga Hill. Within the UNESCO buffer zone, the most exposed areas (poor housings and pathways) are in the Middle City area, particularly the foot of the western hillslope (Tsimialonjafy-Ambanin'Ampamarinana-Ankadilalana-Tzimbazaza quarters), where shallow earth-debris rotational/translational slides take place (Fig. 17). In this area, linear creek erosion of the soft lateritic soil cover creates large gullies, which are rapidly expanding, as shown by the diffused damage of the road pavement and buildings (Fig. 18). The middle-southern eastern slope is characterized by a lower susceptibility, except for the creek basin east of the Rova. The widespread urban cover and the scarce slope steepness determines a general stability of the hill slope northern area, while the southern sector shows some prone areas at the top of the abandoned quarry niches. Rock falls severely affect sector 1 (the western footslope) and sector 2 (southeastern slope) (Fig. 18c). Landslide susceptibility assessment revealed that NDVI was the most important predictor to locate landslides; NDVI obviously cannot be considered as a predisposing factor, but its importance matches the outcomes of multitemporal analysis, where it was noticed that most of landslides were triggered along deforested areas, indicating that this practice can be considered the main predisposing factor to landslide. Slope orientation is considered an important parameter for landslide susceptibility, as well as the slope curvature, this may be since most of the known shallow landslides occur inside small creek catchments located in both the east and the west slopes. The lithologic parameter was not very important, likely because there is a low lithological variability in the study area. The north hill slope sector shows a general stability due to the widespread urban cover and the lower slope angles, except for the area around the roadway tunnel entrance, in the Ambohijatovo quarter. Shallow earth-debris rotational/translational slides can develop from the gullies' slopes and headcut apexes following intense rainfall events and can evolve into earth-debris flows that may channel within the hydrographic network. These phenomena involve the lateritic soils as well as the anthropic waste deposits (often the gullies are used as dumps), enhancing the gully headcut erosion as well as the scarps retreat. More likely they occur in the gullies of the central and southern hill sector (Ankadilalana-Tzimbazaza quarters), not affecting the High City's cultural heritage (only the hovels built along the scarps and the pathways connecting with the Low City; Figs. 10 and 18). Rockmass instability was assessed in correspondence to the fractured granite cliff on the western slope, while in the south-eastern hill sector in correspondence to the dismissed quarries. The carried-out simulation showed that these areas can be sources of hazardous rockfalls, which may impact the populated area of the Middle City, especially affecting poor housing.

Suggestions for mitigation strategies and landslide disaster risk reduction
The obtained landslide maps can be used not only for hazard zoning, but also for prioritizing landslide mitigation measures, e.g., slope reprofiling, bioengineering methods using local material (wood barriers and check-dams, bio-mats, and grids to enhance vegetation regrowth in correspondence to the bare soil sectors),   creek filtering dams, nets, and fences for the rockfall/slides. Even if some results are based on simple approaches and show a degree of uncertainity, this work can represent a starting point for creating structured activity aimed at landslide risk reduction. This activity should start with the creation of a more detailed and updated landslide inventory, where the triggering date, morphological-kynematic characteristics and magnitude of the landslides are defined. This landslide inventory could be used to produce a more accurate and detailed landslide hazard assessment, which, in turn, would result in more effective landslide risk mitigation. In the near future, the disaster risk management (DRM) activities in the High City should also focus on the installation of an automated weather station and sensor network, in order to collect accurate rainfall data and implement an early warning system (EWS). Risk assessment and spatial planning activities should be performed by the municipality of Antananarivo, BNGRDC and technical experts, in order to avoid the construction of new buildings in hazardous areas, manage waste materials, and avoid bad practices such as the waste burning, quarrying, and slope cutting. These activities should harmonize with existing on-going projects, while public awareness must be raised to understand the hydrogeological risks and reduce social vulnerability. All these activities in the future must focus on the sustainable management of the High City in the framework of its enrollment in the UNESCO World Heritage List. Original Paper

Conclusions
In this work, an integrated landslide hazard assessment was carried out in the High City of Antananarivo, one of the most important cultural heritage sites of Madagascar, part of UNESCO Tentative List. In detail the following activities were carried out: • As a first step a multitemporal slope scale landslide inventory spanning from 2015 to 2020 was collected by combining field survey campaigns, interviews with local people and analyzing remote sensing data (high-resolution DEM and satellite orthophotos). Various landslide phenomena were identified: shallow landslides in correspondence to loose slope deposits, debris flows involving both lateritic residual soil and anthropic debris, and rock falls. • Preliminary laboratory analyses were carried out on the soilrock samples collected during the field surveys to geotechni- cally characterize the involved material; this data was then used as input for the simulations carried out. • The spatial probability of occurrence and the impact of these slope instabilities phenomena was then assessed by means of several analyses carried out using a combination of codes and software in order to perform shallow landslide and rockmass susceptibility analysis, as well as debris flows and rockfall runout simulations. • The analysis carried out show that most of the studied area is affected by moderate to very high landslide susceptibility. The areas characterized by the higher susceptibility are located along the Middle City quarters, particularly at the foot of the central sector of the western hillslope (from north to south: Mahamasina, Ankadilalana, Tzimbazaza quarters), where shallow debris rotational/translational slides involve mainly the eluvial-colluvial and the lateritic deposits and rockfalls can occur from the granite cliffs. The latter phenomena only involved shacks and hovels of the Middle City. The soft lateritic soil cropping out in the western footslope are affected by intense linear creek erosion and large active gullies. The eastern slope is less prone to shallow landslides, except for two creek basins the middle-southern sector (one located east of the Rova in the Manjakamiadana quarter, the second is located east of the Ambohipotsy quarter). The southern sector shows ususceptible areas in correspondence to the abandoned quarry slope cuts (Anjahana-Ambohipotsy quarters). The High City cultural heritage buildings are not located in shallow landslide-prone areas, except for the Trano-Gasy Houses, situated in an unstable area on the top of the western slope, between the Rova and the Andafiavaratra Palace. The most exposed linear structures are the pathways running on both mid-slopes in the central-southern sector of the Analamanga Hill. Earth-debris flows do not involve the High City's cultural heritage, but only the hovels and shacks built inside the gullies, as well as the pathways connecting to the Low City. • A preliminary analysis on the rainfall-induced landslide frequency reported by the surveys and the rainfall data from 2015 to 2020 showed that the main landslide events occurred after major cyclones. The analysis of cumulated rainfall and intensity on a 30-day interval suggests that longer rainfall events have a lower influence on landslide triggering, and that rainfall events ranging from 2 to 7 days are more likely to trigger landslides. An intensity (I)-duration (D) rainfall threshold for landslide triggers was defined, despite the paucity of data. • The performed landslide hazard assessment represent a first step for land-use management and planning towards a geo-hydrological risk reduction strategy to be carried out by the scientists, practitioners, and decision makers involved in the High City protection and conservation. This study represents an important contribution for improving the knowledge on landslide processes in Madagascar (in particular in an area with limited landslide data such as Antananarivo) and can be reproduced in cultural heritage sites characterized by a similar complex geomorphological and urban scenarios.