Multitemporal relative landslide exposure and risk analysis for the sustainable development of rapidly growing cities

In recent decades, developing countries have experienced an increase in the impact of natural disasters due to ongoing climate change and the sustained expansion of urban areas. The intrinsic vulnerability of settlements, due to poverty and poor governance, as well as the lack of tools for urban occupation planning and mitigation protocols, has made such impacts particularly severe. Cuenca (Ecuador) is a significant example of a city that in recent decades has experienced considerable population growth (i.e. exposure) and an associated increase in loss due to landslide occurrence. Despite such effects, updated urban planning tools are absent, so an evaluation of multitemporal exposure to landslides and related risks is required. In this perspective, a potential urban planning tool is presented based on updated data depicting the spatial distribution of landslides and their predisposing factors, as well as population change between 2010 and 2020. In addition, a multitemporal analysis accounting for changes in exposure between 2010 and 2020 and an estimation of relative landside risk was carried out. Due to the absence of spatially distributed population data, energy supply contract data have been used as a proxy of the population. The results show that the current higher exposure and related relative risk are estimated for parishes (parroquias) located in the southern sector of the study area (i.e. Turi, Santa Ana, Tarqui, Nulti, Baños and Paccha). Moreover, the exposure multitemporal analysis indicates that most parishes located in the hilly areas bounding the city centre (i.e. Sayausi, San Joaquin, Tarqui, Sidcay, Baños, Ricaurte, Paccha and Chiquintad) are experiencing sustained population growth and will be potentially exposed to an increased risk with a consistently growing trend. The obtained relative risk map can be considered a valuable tool for guiding land planning, land management, occupation restriction and early warning strategy adoption in the area. The methodological approach used, which accounts for landslide susceptibility and population variation through proxy data analysis, has the potential to be applied in a similar context of growing population cities in low- to mid-income countries, where data usually needed for a comprehensive landslide risk analysis are non-existing or only partially available.


Introduction
In recent decades, ongoing climate change, associated with global population growth and the related expansion of urbanized areas, has been responsible for an increase in the frequency and impact of natural disasters due to floods, landslides and wildfires (Knox 1993;Xu et al. 2013;Altan et al. 2015;Arnell and Gosling 2016;Gariano and Guzzetti 2016;Di Napoli et al. 2020a). Developing countries have experienced more severe impacts because people are often concentrated in high-hazard urban areas where housing is highly vulnerable due to poor construction, and early warning systems are commonly absent (Zorn 2018;Aguirre-Ayerbe et al. 2020). The significance of such an impact can be easily understood considering that between 1996 and 2015, approximately 90% of disaster-related deaths occurred in middle-to low-income countries (http:// relie fweb. int/ report/ world/ pover ty-death-disas ter-and-morta lity-1996-2015).
Vulnerability factors, such as poverty, poor governance and the lack of experience in facing natural disasters, are responsible for this disproportionate effect (i.e. the death concentration in developing countries) (Petley 2012). The lack of tools for urban occupation planning, recommendations and mitigation protocols is an additional element of vulnerability that particularly applies to rapidly growing urban areas prone to floods and landslides. For this reason, evaluating and preventing exposure to geohazards, in terms of susceptibility and hazards, are a fundamental step for adequate environmental planning and management, as shown by several scientific contributions in this field (Goetz et al. 2015;Guerriero et al. 2018Guerriero et al. , 2020aLombardo et al. 2020;Segoni et al. 2020;Di Napoli et al. 2021;Novellino et al. 2021;Allocca et al. 2021).
In Latin American countries, between 2004 and 2013, 611 landslides triggered by rainfall and earthquakes were responsible for approximately 12,000 deaths (Sepúlveda and Petley 2015). A relevant example is an event occurred in 2017 in the city of Mocoa in southern Colombia, which killed more than 300 people and destroyed 130 houses (García-Delgado et al. 2019). While the spatial distribution of such events is consistently related to a combination of slope morphometry, rainfall distribution and population density, poverty is a factor controlling the impact on people and is particularly relevant in urban areas of developing countries. Indeed, the presence of informal settlements and their localization greatly impact on the number of fatalities. In such conditions, landslide susceptibility and risk maps are useful tools to develop land planning strategies for preventing such impacts and supporting the sustainable development of cities (Musakwa and van Niekerk 2015;AlQahtany and Abubakar 2020).
Landslide susceptibility indicates the probability of a slope failure occurring in an area depending on its geological, geomorphological, structural and vegetation cover peculiarities (van Westen et al. 2003;Guzzetti et al. 2006;Reichenbach et al. 2018). It Landslides 20 • (2023) Original Paper differs from hazard since it does not directly consider any evaluation of the expected magnitude of an event and its recurrence time (Fell et al. 2008). A landslide susceptibility map spatially reproduces the landslide occurrence likelihood, providing an overview of areas that need recommendations from a perspective of settlement development (Chen et al. 2020;Di Napoli et al. 2020b;Zhang et al. 2020;Arabameri et al. 2021). Landslide risk depends on the characteristics of elements at risk, their vulnerability and the temporal-spatial probability of the occurrence of a damaging landslide event. Risk maps are powerful tools because they also consider the characteristics of exposed elements and provide potential damage scenarios (Bignami et al. 2018;Santangelo et al. 2021;Novellino et al. 2021). In general, landslide risk is evaluated through a multistep analysis: (i) hazard identification, (ii) hazard assessment, (iii) inventory of elements at risk and exposure, (iv) vulnerability assessment and (v) risk estimation (Dai et al. 2002;Glade et al. 2006;van Westen et al. 2008;Corominas et al. 2014). Due to the frequent lack of landslide occurrence timing data, the risk is often evaluated by adopting a simplified approach based on susceptibility scenarios rather than hazard (Ercanoglu 2008;Fell et al. 2008;Arabameri et al. 2017). As an alternative, hazard can be estimated based on the return period of landslide-triggering events (Grelle et al. 2014). A further element of simplification is generally related to vulnerability estimation: it could be very challenging to be correctly estimated and might cause diffuse under-or over-estimation of landslide risk, especially for heterogenous settlements (Glade 2003;Li et al. 2010;Mavrouli et al. 2014). The concept of relative risk applies well to regions where a reduced complexity estimation is preferred due to limits in data availability, such as in rapidly growing urban areas of developing countries (Andrejev et al. 2017).
The city of Cuenca (Ecuador) is an example of a rapidly growing city that in recent decades has experienced a sustained increase in population, having reached ~637,000 people. The magnitude of the change can be estimated by considering that 30 years ago, the total population was approximately 200,000 people. Due to this growth, the urban area has consistently expanded from its original position within the Tomebamba River floodplain and now occupies surrounding hilly slopes. Such slopes are prone to landslides, so damage to settlements and infrastructures is very frequent, preventing the sustainable and safe development of the city (Miele et al. 2021). On this basis, an analysis of the multitemporal variation of exposure, as a function of population growth, and the relative landslide risk in the city of Cuenca was completed to provide (i) an updated overview of the exposure and relative landslide risk to the population and (ii) a general tool for supporting land planning in rapidly growing cities suffering the effects of natural hazards due to the lack of planning instruments. Such a tool can be used to identify areas where the relative landslide risk is higher and, at the same time, the most suitable areas for urbanization. Indeed, for Cuenca city, the first attempt at natural hazard assessment dates back to the early 1990s with a pilot project called PRevention ECuador CUenca PAute (PRECUPA). From the perspective of the analysis, an integrated method consisting of (i) machine learning-based susceptibility assessment by using a maximum entropy algorithm, (ii) energy supply contract analysis for exposure quantification and (iii) relative risk estimation was used.

Study area
The study area comprises Cuenca city and the surrounding hilly area that has been increasingly occupied by settlements (Fig. 1). Cuenca is the capital of the Azuay province of Ecuador and extends over an area of approximately 124 km 2 . In 1999, the city's historic centre was included on the UNESCO World Heritage Site list due to its importance as a cultural and governmental centre of the Cañari and Inca civilizations,and as an example of Renaissance urban planning in the Americas during the Spanish colonial period.
The city lies within an inter-Andean valley, which formed following compressional deformation controlled by major NE-trending faults (Noblet et al. 1988;Hungerbühler et al. 2002). The geology of this area is represented by Mesozoic marine and subaerial sedimentary deposits that cover the Palaeozoic metamorphic basement (Noblet et al. 1988). The sedimentary series is more than 2400 m thick and is composed of two main sequences separated by a regional unconformity. The lower sequence consists of fluvial and brackish delta plain deposits containing ubiquitous metamorphic pebbles from the Cordillera Real. From the bottom, this series consists the Biblián, Loyola, Azogues and Mangan Formations, which include sandy clays, laminated shales with gypsum, tuffaceous sandstones, siltstones and conglomeratic sandstones. The upper sedimentary sequence is composed of volcanic clast-bearing rocks of the Turi Formation, divided into the Turi and Santa Rosa Members; it consists of tuffaceous coarse sandstones, volcanic clast-supported conglomerates, matrix-supported volcanic breccias and minor tuff layers. Furthermore, the late Miocene volcanic series of the Tarqui Formation crops out in the area that unconformably covers a wide range of volcanic and Tertiary sedimentary formations. The Tarqui Formation is composed of two members: (i) the Tarqui Member, which is formed by poorly consolidated and intensely weathered red volcanic airfall deposits, and (ii) the Llacao Member, represented by debris flow deposits of volcanoclastic materials (Steinmann 1997).
This area is known for frequent landslides involving settlements and infrastructures (Fig. 2). For instance, the buildings of the Faculty of Philosophy of the University of Azuay are consistently affected by slow-moving landslides and periodically damaged by slope deformation (Sellers et al. 2021a). On March 29, 1993, a large landslide (20 million m 3 ) occurred northeast of Cuenca city, damming the Paute river and causing the formation of an artificial lake that flooded fertile land and destroyed houses, roads, railways and a regional thermoelectric plant (Plaza et al. 2011). Moreover, the segment of the Pan-American Highway crossing the city is continually affected by landslides that induce damage and, in some cases, accidents (Miele et al. 2021). The high frequency of landslides is related to the considerable yearly rainfall amount (approximately 900 mm), the remarkably high frequency of earthquakes of relevant magnitudes (ranging from 4.0 to 4.9 Mw according to Ecuador's Geophysical Institute, https:// www. igepn. edu. ec/) and the predisposing effects of slope morphologies and geological characteristics or rocks forming slopes.

Data and methods
To evaluate the relative landslide risk to the population of the city of Cuenca and the multitemporal variation of exposure in terms of population growth, a procedure consisting of machine learningbased susceptibility analysis, population growth estimation through energy supply contract analysis and relative risk evaluation was executed (Fig. 3). Below, the detailed data and methods used for the analyses are provided.

Landslide inventory map
From the perspective of evaluating the relative landslide risk to people, the assessment of susceptibility to landslides was carried out using the available landslide inventory map (LIM) prepared by Miele et al. (2021). The LIM is composed of 710 landslides obtained by (i) field surveys, following the approach presented in Sellers et al. (2021b;(ii) photo interpretations; and (iii) interferometric data obtained from the processing and interpretation (Di Martire et al. 2016;Guerriero et al. 2019;Ammirati et al. 2020) of Sentinel-1A and Sentinel-1B images (ascending and descending passes) acquired between October 2016 and May 2019 using the coherent pixel technique-temporal phase coherence (CPT-TPC) approach (Mora et al. 2003;Iglesias et al. 2015). The landslide database contains valuable information regarding the type of movement according to Cruden and Varnes (1996), state of activity, location, triggering factor (i.e. precipitation or anthropic), geology, land use, velocity and further information (such as the field survey date, induced damages and fatalities) (Sellers et al. 2021b). In the database, it is possible to recognize rockfalls (72-10.1%), topples (3-0.4%), flows (8-1.1%), spreads (5-0.7%), rotational slides (550-77.1%) and translational slides (72-10.1%) (Fig. 6). The mass movements represent the principal hazard of the area since they affect the urban area damaging road networks and buildings. Landslide initiation is mainly related to intense or prolonged rainfall.  Table 1 for the code descriptions of lithologies). Red lines represent the four principal rivers that crossing Cuenca town

Original Paper
Geostructural discontinuities notably influence the occurrence of gravitational phenomena involving rock masses, such as falls, topples and planar slides. Rockfalls and topples affect steep artificial slopes around the main infrastructures (Miele et al. 2021). These phenomena mainly occur where anthropogenic actions have caused cuts linked to the construction of infrastructures. The high slope angles represent an essential element in favouring translational slides (Raso et al. 2020), whose action is often enhanced by diffuse and/or concentrated erosion.
The abovementioned covariates were tested for multicollinearity, representing the occurrence of high intercorrelations among two or more independent variables. When collinearity issues arise, the concerned variables should be removed because they imply a redundancy in a model (Gareth et al. 2013). The variance inflation factor (VIF), by using the R "usdm" package (Naimi et al. 2014), was employed to test multicollinearity. As a rule of thumb, a VIF value that exceeds 5 or 10 indicates severe issues (Gareth et al. 2013).
Landslide susceptibility was assessed using the MaxEnt modelling algorithm (Elith et al. 2011). MaxEnt is a presence-only (PO) spatial distribution method that addresses only locations where landslides are present (Phillips et al. 2006;Phillips and Dudík 2008). It uses occurrence data and a large (typically 10,000) number   (2023) of points throughout the study area, which are referred to as background points. To reconstruct the potential distribution of an event, MaxEnt calculates two probability densities: for all presence and background points (Guillera-Arroita et al. 2014). The resultant probability occurrence map ranges from 0 (no landslide probability) to 1 (highest landslide probability). Since MaxEnt predictions are sensitive to initial modelling settings (Merow et al. 2013), different MaxEnt implementations were evaluated through the ENMeval R package (R Core Team 2021) to detect the settings that optimize the trade-off between goodnessof-fit and overfitting (Muscarella et al. 2014). In fact, in MaxEnt, it is possible to set up two main parameters: (1) feature classes and (2) regularization multipliers (for further details, consult Elith et al. (2010) and Merow et al. (2013). For the analysis, regularization values between 0.5 and 10, with 0.5 steps, were investigated, and the following feature classes were considered: linear, linear + quadratic, hinge, linear + quadratic + hinge, linear + quadratic + hinge + product and linear + quadratic + hinge + product + threshold (Muscarella et al. 2014).
Model evaluation was executed using a spatial block crossvalidation scheme (Muscarella et al. 2014). This method converts part of the occurrence and background points into evaluation bins and uses them to reduce the spatial -autocorrelation between training and validation points (Hijmans 2012;Wenger and Olden 2012). The block cross-validation scheme could assess model transferability, i.e. the ability to extrapolate predictions into new areas (Roberts et al. 2017) and penalize models based on meaningless predictors (Fourcade et al. 2018).
Because no consensus currently exists regarding the most appropriate metric or approach to evaluate the performance of models (Fielding and Bell 1997;Warren and Seifert 2011;Peterson et al. 2011), different statistical approaches were adopted for the models' predictive performance with presence-background data (Muscarella et al. 2014). In this case, the best model reliability combination was chosen following three criteria: (i) the lowest delta Akaike information criteria (ΔAICc) (Burnham and Anderson 2002), (ii) the area under the curve plot based on the training data (AUC train ) (Hanley and McNeil 1982) and (iii) the difference between training and testing AUC (AUC diff ) (Warren and Seifert 2011).
Moreover, the landslide ratio of each predicted landslide susceptibility class (LR class ) was employed as a further performance evaluation of the landslide model. The LR class formula is as follows (Eq. 1): (1) LR class = % of contained sites in each susceptibility class % of predicted landslide areas in each susceptibility class This index was developed specifically to deal with situations when landslide boundaries are not available but their locations are known. The advantage of the LR class index is to consider both the predicted stable and unstable areas and thus significantly decrease over-prediction. A higher value of LR class corresponds to a lower over-prediction by the model (Park et al. 2013).
Predictive performance estimation is only a partial metric of model goodness. Predictor contributions represent a further key step that should be assessed to comprehensively estimate the validity of a model. In this contribution, the investigation has been carried out considering (1) predictor importance and (2) percentage contribution (also called permutation importance; Oke and Thompson 2015). Predictor importance represents the degree to which single environmental variables contribute to the final model so that the percentage contributions for all predictors in a model sum to 100% (Phillips and Dudik 2008). A large permutation importance decrease indicates that the model depends heavily on that variable (Phillips 2017). A very useful and detailed explanation was given by Bradie and Leung (2017).

Energy supply contract analysis
The exposure and the relative landslide risk evaluation were completed by considering the population as the element at risk. To account for population growth, a multitemporal assessment was executed considering data representative of the population distribution in the study area in 2010, 2018, 2019 and 2020 (the only available years). Although population distribution and projections are available for the entire city of Cuenca, this is not available in a discrete way for each parish. Such data are notoriously essential in any landslide risk evaluation. To provide proxy data on inhabitant allocation, the distribution of energy supply contract (in the following specified as ESC) data derived from the Institute of Studies on the Sectional Regime of Ecuador (IERSE, https:// gis. uazuay. edu. ec/ visor es/ info-z6/) were used. Although it does not correspond to the number of people living in each area of the city, the utilities' availability (i.e. ESC) was exploited as a proxy of households and can be considered as alternative data from a relative risk assessment perspective. ESC data consist of georeferenced points corresponding to different addresses and therefore to a different family unit (household). Each point corresponds to a single family unit, and the potential of such data is related to their public availability. For each parish, such information was aggregated and associated with the specific area and used for relative risk evaluation. Figure 4A

Exposure and relative risk analysis
Considering the risk definition introduced by Varnes (1984), landslide risk can be assessed qualitatively (Wang et al. 2013) or quantitatively (Chang et al. 2021). Commonly, for a wide area, where the quality and quantity of available data are inadequate for quantitative analysis, a qualitative risk evaluation may be more appropriate (Andrejev et al. 2017). In the presence of data that might support a simplified quantitative exposure estimation, the relative risk could be assessed. In the context of the proposed analysis, data from susceptibility analysis and ESC, and their variations between 2010 and 2020 (i.e. 2010, 2018, 2019 and 2020), were used as a basis for multitemporal landslide exposure evaluation and risk estimation over the study area. The evaluation was completed by considering only the risk for people, so that exposure analysis accounted for household distribution. For the city of Cuenca, the use of energy supply contracts as a proxy of people distribution fit the choice of evaluating the multitemporal exposure and landslide risk relatively. Indeed, while exposure analysis is complete considering the distribution of elements at risk for each landslide susceptibility class, relative risk is estimated using the matrix of Fig. 5 which relates susceptibility to landslide and exposure (Andrejev et al. 2017). The outcomes of exposure analysis were expressed in the form of histograms, representing the number of ESC located in each susceptibility class. The histograms were obtained by using the Sturges method (Sturges 1926), which highlights how different Cuenca sectors underwent an increase in contracts over time and, therefore, also in an exposure. For the relative risk estimation, a matrix-based approach, allowing the determination of a specific risk class for each association of susceptibility class and exposure class, was adopted. Specifically, a 5 by 5 classification matrix (Fig. 5) was developed accounting for classified exposure data, in terms of concentration in each cell of the susceptibility map (see Supplementary Material 3), through the natural break criterion. Exposure classes considered for the analysis are E5, 0 to 1; E4, 1 to 2; E3, 2 to 3; E2, 3 to 6; and E1, 6 to 11. In particular, this analysis employed the last available information referred to in 2020. As a final product, a relative risk raster map (10 m pixel resolution) was produced. In the context of this analysis, the vulnerability was implicitly considered equal to the unit. Table 2 shows the results of multicollinearity analysis carried out through VIF estimation and its comparison with a predefined threshold value. Concerning that point, it must be noticed that there is no definite and approved threshold in the scientific literature. However, it is generally accepted that VIF values higher than 10 indicate severe collinearity (Hair et al. 2010), even though this rule of thumb lacks a theoretical basis (Gómez et al. 2016).

Multicollinearity examination
The employment of many environmental covariates might lead to overfitting problems, but, in this work, the individual predisposing factor values differ considerably from the aforementioned threshold. Based on Table 2, the highest VIF value is 2.02, corresponding to the topographic wetness index, while the smallest VIF values are 1.01 and 1.02, which are associated with northness and eastness, respectively. Therefore, no environmental variables exceed the critical   value, and thus, these results satisfy the criterion (VIF < 10) indicating that there is no multicollinearity among the landslide PFs. Figure 6 shows the landslide susceptibility analysis outcome subdivided into five classes through the natural breaks classification (Jenks 1967). This is accomplished trying to reduce the standard deviation within each class and maximizing that between the classes themselves (Basofi et al. 2018;Novellino et al. 2021). The percentage of susceptibility classes is summarized in Table 3.

Landslide susceptibility
As shown in the map, the most susceptible areas in Cuenca are placed on mountainsides that border the city. In these portions, slopes are steep and concave, and roads create local discontinuities. The central part of the map is characterized by very low and low susceptibility zones representing the Cuenca urban area located in the plain characterized by the presence of alluvial deposits and different terrace orders. Along the Tomebamba riverbanks, as well as the other rivers of the area, due to their erosive action that affects the foot of the slopes, there are sectors predisposed to landslides with a medium to very high susceptibility. Peri-urban and rural areas, instead, are in medium to very high susceptibility zones where there are steeper mountainsides.

Susceptibility model validation
According to Swets (1988), the obtained models achieved fair-togood predictive performance, with AUC values ranging from 0.763 to 0.866 (Fig. 7). The lower value is associated with models with linear or linear + quadratic features and high regularization values (i.e. 9.5 and 10). The higher value is associated with a model that contains all features and low regularization values (i.e. 0.5 and 1). AUC diff values scored from 0.06 to 0.14. Among the resulting 120 combinations, the one reporting the lowest ΔAIC c was chosen. The selected model is characterized by the following: linear + quadratic + hinge + product + threshold features, AUC value of 0.82, average AUC difference value of 0.08 and ΔAIC c value equal to 0 (for further information, refer to Supplementary Material 2).
The LIM availability allows an assessment of the model performance that considers field data. By intersecting the landslide detachment points and the final susceptibility map, it is possible to achieve information about the number of landslides and their distribution into each different landslide susceptibility class. Subsequently, the areal extent of the susceptibility classes was calculated. Areas characterized by high and very high susceptibility involve 15.5% and 10.6% of the total study area, respectively (Table 3). Very low and low susceptibility classes cover approximately 55% of the study area, particularly in the central sectors, namely, Cuenca city. The remaining portions are assigned to the moderate (19.4%) class. Furthermore, the highest concentration of landslides can be found within the highest susceptibility class values (high, 25.2%, and very high, 48.9%). In addition, approximately 3% are in the very low susceptibility class, and the remaining 22.9% are distributed in the low class (i.e. 9.4%) and moderate class (i.e. 13.5%). The last column of Table 3 highlights the LR class percentage for each susceptibility class. The highest LR class value corresponds to the very high susceptibility class: more than 80% are in the high and very high susceptibility classes. This evidence roughly implies that, if a landslide occurs, then the predicted susceptible area has approximately an 80% chance of including the landslide itself.
The spatial aggregation of the susceptibility map confirmed that the largest part of the study region has a low susceptibility to landslide events. Approximately 75% of actual landslides were localized in the high and very high susceptibility classes. Moreover, the higher susceptibility classes showed higher values of LR class percentages. These outcomes show significant quantitative agreement between the simulated scenario and the landslide inventory map. All the produced analyses permit zoning of the complex territory of the Cuenca area to identify the spatial probability of landslide initiation in areas characterized by specific conditions that materialize by the considered environmental variables. Moreover, landslide density is characterized by an increasing trend when passing from the lowest to the highest susceptibility classes. These observations underline that, despite the limited area extension of the very high susceptibility class, most of the surveyed landslides are within the latter.

Factors predisposing to slope instability
As the last outcome, the variables' contribution was evaluated. This result allows us to understand which variables have greater importance in the final models' implementation. Table 4 presents the impact of each variable. In particular, the results reveal that the main conditioning variables (i.e. the variables that assume a fundamental role in the final landslide susceptibility) are slope steepness, distance to roads and planform curvature. Therefore, a model with a higher fit is achieved through the aforementioned variables. In general, observing the per cent contribution, a good spread of values is desirable. Conversely, if the contribution of a variable is high in the model (i.e. higher than 70%), something is not correct and that variable is not encompassing many variations or is correlated to other variables. Other variables, such as distance to streams, land   (2023) use, geology and relative slope position show noteworthy values in the model. Last, low values close to zero are assigned to eastness, northness, profile curvature and the topographic wetness index. The primary roles in slope stability are related to slope steepness, which influences water infiltration and modulates gravity force (Huat et al. 2006). Particular conditions are related to the presence of road cuts where steep slopes and slope unloading related to the excavations for the road construction influence slope stability. An additional fundamental covariate contribution to the slope stability of this area is represented by the planar curvature that adjusts the convergence or divergence of water in the direction of landslide movement and landslide material (Ohlmacher 2007). Furthermore, in the Cuenca territory, hillsides with low planar curvature values are the most susceptible to earth and debris flows and earth and debris slides. Indeed, this factor is used to identify gullies (Wieczorek et al. 1997), and debris flow initiation areas can be recognized where the curvature values are negative (Park et al. 2016). Finally, other considerable predisposing factors to slope stability are geology and land use (e.g. Glade 2003;Cevasco et al. 2014). Unconsolidated material, such as alluvial and colluvial deposits or pyroclastic lithologies derived from the near volcanic activity, covers the surrounding mountainous landscape, resulting in a very high susceptibility to sliding. Figure 8 provides an overview of the exposure in each parish updated to 2020 (i.e. ESC number for each susceptibility class). In this way, the parishes with higher exposure to risk were identified. As observed from the inset graphs, the exposure is higher in boroughs surrounding the southern portions of Cuenca, namely, Turi, Santa Ana, Tarqui, Nulti, Baños and Paccha. Such boroughs show Original Paper a higher number of ESC (i.e. an element at risk) in high and very high landslide susceptibility classes. In contrast, the northern sectors of the study area report a lower number of ESC included in the higher landslide susceptibility classes. Last, the central portion of the study area, in which the city of Cuenca is located, is characterized chiefly by very low and low exposure to landslides. However, recently, as shown in Figs. 4 and 8, the city of Cuenca and its neighbouring areas have experienced a substantial exposure increase. Similar problems can be easily found in other rapidly expanding cities where the construction of new boroughs takes place in hilly and mountainous areas more prone to instability (Di Martire et al. 2012). This condition is marked by the histograms shown in Fig. 8: it is possible to note that 11 out of 15 parishes have several contracts distributed almost uniformly, therefore highlighting a large number of ESC located in the high and very high susceptibility classes.

Exposure and relative risk assessment
Since the city has experienced consistent growth of population between 2010 and 2020 (Fig. 4A, B), with a sustained occupation of peri-urban hilly areas characterized by higher susceptibility in comparison to the centre of the city, a detailed exploration of the spatial and temporal rates of change for both energy supply contracts and related exposure for each parish is reported in Fig. 9. In Fig. 9, graphs "a" highlight ESC for each landslide susceptibility class, yielding the percentage of the change in the exposure to landslides from the reference year 2010 to 2020 (i.e. 2018, 2019 and 2020), while the same rates in terms of the total exposure (i.e. normalized to the total 100%) are given in graphs "b". In addition, the relative areal extensions of different susceptibility classes are also reported on graph "a" (violet bar). In the reference period, the parishes of Sayausi, San Joaquin, Tarqui, Sidcay, Baños, Ricaurte, Paccha and Chiquintad experienced a higher increase in ESC located in the very high susceptibility class corresponding to an increase in the exposure to landslides (Fig. 9). Such an increase, in comparison with 2010, ranged between 33% of Paccha and 300% of Sayausi. Conversely, the parishes of Turi, Nulti and Cuenca experienced a higher increase in ESC located in the very low to medium susceptibility class. Such an increase, in comparison with 2010, ranged between 20% of Cuenca and 100% of Nulti.
However, in absolute terms, the districts of Tarqui, Turi, Nulti and Santa Ana show the highest exposure, with the highest number of energy supply contracts located in high susceptibility areas between 2018 and 2020. Conversely, the boroughs of Sayausi, Sidcay, Ricaurte, Cuenca and Chiquintad show the lowest exposure, with the lowest number of energy supply contracts located in high susceptibility areas in the same period. In Supplementary Materials 3, it is possible to consult the tables that quantitatively represent the graphs shown in Fig. 9. Figure 10 depicts the outcome of the relative risk estimation in the form of a risk map. It shows five classes of landslide risk for the studied area. Specially, the centre of Cuenca city is located within a very low-risk zone because of a significant number of buildings and infrastructures built up in a flat area characterized by very low landslide susceptibility. However, some areas with higher landslide risk correspond to the banks of the rivers that pass through the historic centre of the city. Overall, the low-risk class (green pixels) is observable in correspondence to low inclined slopes of the hills surrounding Cuenca where the landslide susceptibility is generally low. On the other hand, steep slopes are generally characterized by medium to high degrees of landslide susceptibility. In these areas, the analysis pointed out a large number of ESC and, consequently, a final estimate of elements at risk. These slopes are generally Fig. 9 Multitemporal evolutionary perspective of exposure to landslides from 2010 to 2020. Graph "a" highlight the ESC variation, in percentage terms, for each landslide susceptibility class for the considered years (i.e. 2018, 2019 and 2020); graph "b" represents the same rate in terms of total exposure, namely, normalized to the total, i.e. 100%. The coloured bars of the histograms indicate the different susceptibility classes ranging from very low to very high, and the violet bar, included in histogram "a", specifies the landslide susceptibility class areal extent. In the middle, maps of exposure variation for each parish of the study area are given during the time span considered coloured with orange and red pixels indicating the high to very high risk. In particular, the southern sector of the study area is characterized by the highest relative landslide risk.
In the last decade, landslide susceptibility analyses have also been conducted in emerging and developing countries (O'Hare and Rivas 2005; Klimeš and Rios Escobar 2010;de Listo and Carvalho Vieira 2012;Jamalullail et al. 2021). Several attempts have been made to estimate landslide risk in various contexts where demographic growth is very pronounced (Rahman 2012;de Listo and Carvalho Vieira 2012;Rojas et al. 2013; Alcántara-Ayala and Moreno 2016). Considering the high growth rate recorded in the last few years, it is appropriate to assess the exposure of landslide risk over time. This type of analysis makes it possible to identify the most critical areas and sectors of the city of Cuenca in terms of risk evolution due to population growth. The outcome of the analysis represents a significant land planning tool for the definition of urban occupation plans, land use recommendations and mitigation protocols that should be applied to reduce the impact of a landslide occurring in urban areas (Klimeš et al. 2020;Sultana and Tan 2021). The problem of landslides involving settlements and claiming human lives is particularly relevant in low-to mid-income countries because people are often concentrated in high-hazard urban areas and vulnerability factors, such as the poor building of housing, poor governance, the lack of experience in facing natural disasters and the absence of early warning systems, consistently exacerbate landslide impacts (Petley 2012;Zorn 2018;Aguirre-Ayerbe et al. 2020).
This qualitative procedure for evaluating the landslide exposure in Cuenca attempts to provide information for risk assessment, which is useful in a preliminary stage of regional planning or more detailed studies on high-exposure areas. Therefore, the procedure proposed in this study can be implemented when not all the useful information for risk assessment is attainable. Exposure quantification, which is a basic input in spatial and risk reduction planning, is the main objective of this study. It is important to mention that, as not all the information is available, landslide risk values are not expressed in absolute terms, but relative landslide risk may be a good proxy for districts that have encountered, in the last decade, a population growing into very high-risk territories. With limitations, this study has allowed us to estimate which areas are more prone to instability, which areas have a high relative landslide risk and to establish future risk by consulting the trends in the different parishes.

Conclusions
In this paper, an analysis of the exposure to landslides, its multitemporal variation between 2010 and 2020 and related relative risk for the city of Cuenca in Ecuador (Latin America) has been presented. This study provides insights into important issues, such as (i) the effect of the sustained expansion of urban areas due to population growth on landslide exposure variation and (ii) the reduced complexity method of risk assessment in the presence of partial data only (i.e. landslide susceptibility rather than hazard and energy supply contracts rather than population distribution or exposure). The results indicate that the current higher exposure and related relative risk are estimated for districts located in the southern sector of the study area (i.e. Turi, Baños, Santa Ana, Nulti, Tarqui and

Original Paper
Paccha). In addition, by carefully observing the study area, it is possible to understand how most boroughs of the city are located in the hilly areas bounding the centre (i.e. Sayausi, San Joaquin, Tarqui, Sidcay, Baños, Ricaurte, Paccha and Chiquintad); these boroughs, which are experiencing sustained population growth, could be exposed to an increased risk with a sustained growth trend. This is also related to the presence of modest buildings and the absence of early warning systems.
The obtained results can be considered a relevant tool for future land management in the town of Cuenca. The proposed method, using potentially available data for mid-and low-income countries (i.e. landslide inventory and a proxy of population distribution), has the potential to be applied in many contexts where a minimum dataset is available or can be developed based on either field or remotely sensed data. The results of the risk analysis are useful for ranking the parishes based on increasing risk and for supporting decisionmakers in prioritizing funding for risk mitigation measures.