1 Introduction

Over the last decades, the frequency of landslides has increased (Lee et al. 2021) due to increasing human use of the environment (Fell 2018) by communities with a lack of environmental knowledge or loss of cultural heritage (Harmsworth and Raynor 2005; Pradhan et al. 2011; Hadji et al. 2013; Pourghasemi et al. 2019; Anis et al. 2019). Urban communities disconnected from the surrounding environment pose new environmental risks, landsides being one of the starkest. Some researchers have argued that humans carry some responsibility for increasing the landslide risk by living in hazardous regions (Achour et al. 2017; Manchar et al. 2018; Ercanoglu et al. 2018). Land cover change dynamics such as deforestation and expansion of human settlements in exposed areas have defined a new risk aspect (Brander et al. 2018; Karim et al. 2019; Hadji et al. 2013). The strong relationship between landslide risk and land cover change can explain the increase of landslide events related to the loss of vegetated surfaces, being land cover one of the key factors in landscape and land management (Gonzalez-Ollauri and Mickovski 2017; Pisano et al. 2017), but at the same time climate change impacts have contributed to slope instability (Crozier 2010). Landslide risk expresses the relation between rapid economic–demographic changes and increase in environmental risk. Fast economic changes define a spread of urban and agricultural land use to support a growing population, resulting in more people being exposed to hazards (Smith 2003). The loss of ecosystem services related to land cover change which leads to vegetation denudation defines a loss of human welfare that is not immediately perceived and accounted for by society (Costanza 1997; Crespin and Simonetti 2016). Forest and landscape restoration are not often considered for supporting local communities to produce ecological benefits and avoid landside events (Paudyal et al. 2017), and fast economic-demographic growth generally leads to an increase in environmental hazards such as landslides due to a no planned urbanization (Cui et al. 2019). Above the natural triggers, anthropogenic influences on the environment is sometimes the primary cause of a landslide activation, dramatic land cover changes, alteration of morphological aspects, and hydrological changes are all disturbances that increase landslide risk (Jaboyedoff et al. 2018).

At the same time, landslides are a phenomenon in which Disaster Risk Reduction (DRR) can make a difference and avoid catastrophic impacts. DRR planning can take action in all phases of the landslide disaster cycle: Mitigation, Preparation, Response, and Recovery. Information on landslides can be crucial for policy to focus efforts on the most vulnerable regions. In this sense, it is important to highlight hotspot risk areas at global level to understand where to prioritize DRR practices. Susceptibility to landslides has been studied at global scale (Hong et al. 2006; Hong et al. 2007; Stanley and Kirschbaum 2017), also considering overall risk assessment (Yang et al. 2015; Nadim et al. 2006).

This paper aims to highlight where and when people could be affected by a landslide, considering monthly time granularity in order to facilitate a landslide DRR process. A multi-scale workflow is adopted, considering a global risk approach in order to define the most exposed areas and providing a downscaling methodology for the risk maps to sub-regional level. A Geographic Information System (GIS) approach is needed to cope with the strong spatial dimension of landslide events. The number of people living in hazard exposed areas is strictly related to the urban and human infrastructure distribution.

From a geophysical perspective, a landslide can be defined as material (rock or soil) movement due to gravity, activated by environmental trigger factors such as a whether event, rain, snow, or seismic energy release. In a landslide area, we can recognize a separation zone, a movement zone and an accumulation zone (Varnes, 1978). In all the above zones, topography must be considered in relation to its exposure to the hazard.

Our approach considers the landslide human exposure as a result of the hazard frequency with the presence of population considering rain incidence as spatial–temporal triggers. The risk is expressed as the number of affected people per month by the Global Administrative Unit Layer (GAUL) administrative boundaries (FAO 2015). The landslide exposure can be reduced based on a hazard factor, avoiding the landslide events through sustainable land management, and avoiding people’s falling within exposure hazardous areas. Figure 1 conceptualizes the landslide exposure triggered by rain on monthly basis.

Fig. 1
figure 1

Conceptualization of landslide human exposure. Pop (population density) characterizes the exposure and rainfall the landslide trigger event

Full size image

To gain an initial overview of landslide exposure, the National Aeronautics and Space Administration (NASA) Global Landslide Catalogue (GLC) was considered (Kirschbaum et al. 2009, 2015). This catalogue has been collected based on landslide events from 2007 with updates by the Google alert system, considering online documents like news and reports that mention a recent event. The collection of a standardized landslide geospatial database is not common practice, but considering the NASA global landslide catalogue it is possible to recognize a landslide incidence in countries such as the USA, India, Philippines, Nepal, and Indonesia. There is not a clear relation between the number of landsides and the number of human fatalities (Fig. 2). Some countries such as the USA have a very high number of events (1661) and a low number of fatalities (61). However, in some other countries such as Afghanistan, 12 events with 2231 fatalities were registered. This can be explained in part by the intensity of the events and the population density in the hazard area, but also because some countries have developed a robust DRR framework. There is a clear influence of the national wealth of nations in relation to their DRR preparedness. Another aspect that must be considered is the under-reporting of events in emerging countries with respect to technologically advanced countries (Kirschbaum et al. 2010) or underestimation in remote or rural areas compared to more urban areas (Haque et al. 2016). The NASA GLC has been already implemented in testing landslide susceptibility analysis confirming that less than 25% of reported points are in low landslide susceptible areas (Emberson et al. 2020). The lack of reporting affects the report-based landslide inventories but at same time the NASA GLC remains a global reference data set being available, open, an already distributed in GIS format.

Fig. 2
figure 2

Landslide incidence and fatalities on global countries considering the NASA catalogue between 2007 and 2018

Full size image

The aims of this study were as follows:

  1. (i)

    Performing a monthly disaggregation of landslide exposure considering rainfall as a trigger factor;

  2. (ii)

    Defining a multi-scale approach to landslide exposure assessment and monitoring;

  3. (iii)

    Identifying explanatory variables for landslide events at local scale.

2 Material and methods

A probabilistic approach starting from Varnes’ landslide risk definition (Varnes 1984) was adopted. In the case of Varne’s equation the landslide risk (R) is defined as the probability that an extraordinary event such as rock or soil movement (hazard, H) affects some terrain elements (vulnerability, V) in a specific time window and spatial area. This probabilistic approach is summarized by Varnes’ equation:

$$R \, = \, V*H$$
(1)

According to the above formula the landslide risk is defined by the combination of morphological instability, trigger events, that define the Hazard, and human presence on territory that describes the Vulnerability. The vulnerability is a complex dimension depending on socio-politico-economical context (Peduzzi et al. 2009), difficult to quantify over large areas. In this study, the Varne’s concept was modified defining the human exposure to landslide as a combination of the hazard frequency pondered by rain and the population presence (Eq. 2).

2.1 Global approach

We applied the following formula to obtain a global risk analysis expressed as the ‘number of people exposed to the event’ considering the landslide event triggered by rainfall:

$${\text{Landslide}}\,{\text{Human}}\,{\text{Exposure}}\left( i \right) = \left( {F*{\text{Pi}}} \right)*\frac{{{\text{Pop}}}}{10}$$
(2)

where: Landslide human exposure of the month i coincides with the number of people living in hazard areas. \(F\) is the frequency of a landslide event expressed in percentage of occurrence. Pi is the average rainfall in month i.

POP/10 is the population divided by 10 in the spatial unit of 1 km2; based on empirical findings by Arnold et al. (2006) and Nadim et al. (2006) that about 10% of the population living in 1 km2 of landslide event areas is affected.

In order to calculate global landslide exposure, we have chosen a set of databases with global coverage and from open datasets:

  • Landslide frequency from the Geohazards Global Assessment Risk (GAR) representing the Hazard (UNEP 2018) with a spatial resolution between 1 and 5 km (UNISDR 2015);

  • Global landslide georeferenced catalogue by NASA, to integrate into the hazard for past landslide events.

  • Landscan 2018 to consider the population presence expressed in number of people/per grid cell, with spatial resolution of 1 km, representing population (Pop) (Bright et al. 2016; Dobson et al. 2000);

  • Average monthly rainfall (1970–2000) from ‘WorldClim’ at a spatial resolution of 1 km representing the trigger factor of the hazard (Fick and Hijmans 2017).

The priority coincided with the search and use of the landslide frequency data. We have chosen the global layer frequency elaborated by GAR. The GAR frequency has been elaborated considering the cumulative effect of a set of environmental features, susceptibility factors. The GAR has been elaborated based on the European continent, this being widely studied already and having one of the most complete and viable datasets on past landslide events (Jaedicke et al. 2014), serving as an approximation for a global analysis. The set of environmental elements considered in GAR are slope factor; lithological features; soil moisture; and precipitation. The result expresses the frequency of landslide events by year.

The GAR hazard layer was integrated with the georeferenced past landslide events from the NASA catalogue (2007–2018). A buffer of 3 km was generated around the NASA landslide point locations, adding to this buffer the highest frequency value (Table 1). The buffer was applied to consider not only the accumulation areas but also the movement area. The run-out distance has a direct relation with the height of the fall (Devoli et al. 2009), however considering landslide/rain events the run-out distance of 3 km seems to match with overview statistics (Brunetti 2014).

Table 1 Hazard susceptibility levels
Full size table

The GAR hazard, the frequency (F), was disaggregated at monthly timescale considering the monthly rain incidence (Ri) in order to obtain the hazard frequency per month. The monthly disaggregation was performed considering the monthly average rainfall by WorldClim in a time series of 1970–2000. This set of raster layers were reclassified into a susceptibility value according to Table 2. The GAR hazard layer was multiplied by the reclassified monthly rain raster obtaining 12 monthly hazard raster.

Table 2 Rain susceptibility levels.
Full size table

The 12 monthly hazard layers, reclassified as a susceptibility value according to Table 1, were multiplied by the number of people in each grid cell divided by 10. The population information was divided by 10 as an approximation of the number of people affected by a serious event in the spatial unit (1 km2) (Nadim, 2006).

The landslide physical exposure equation was calculated by GIS analysis based on raster analysis according to the workflow in Fig. 3.

Fig. 3
figure 3

Workflow to calculate the Exposed population according to the hazard frequency classes

Full size image

In Fig. 4 the global workflow has been schematized based on a cartographic representation considering Nepal as example. The hazard components, defined by GAR landslide susceptibility layer and the monthly rain layers, have been combined with the population layer (Pop). From this combination It has been obtained the population living in hazard landslide areas. The final raster output can be aggregated by administrative layer considering the zonal statistics.

Fig. 4
figure 4

Physical Exposure Equation in GIS workflow considering a specific month on Nepal

Full size image

2.2 National analysis by logistic regression

From the global landslide exposure described in the previous section, 4 countries were selected to carry out a Logistic Regression (LR) analysis. The aim of this step is to identify which environmental variables explain landslide probability best. Bangladesh, India, Indonesia, Nepal have been selected based on a long-standing landslide problem (Fig. 2) and based on the high number of landslides points that have been collected by the global landslide catalogue (Kirschbaum et al. 2015) for these countries.

The LR method was chosen to calculate the hazard probability based on the presence and absence of the landslide points as the dependent variable. The landslide probability was related to a series of independent variables: slope, rainfall, land cover, soil properties, and soil moisture.

The LR methodology expresses the relationships between dependent and independent variables on the basis of the Receiver Operating Characteristic (ROC) curve. This indicator expresses the explanatory power of a LR model ranging from 0 to 1, where 0.5 is commonly considered the threshold value for an explanatory model. The following equation describes the LR (Kleinbaum and Klein 2002):

$$f\left( z \right) = \frac{1}{{\left( {1 + e^{ - z} } \right)}}$$
(3)

where z:

$$z = \, \beta_{0 \, + } \beta_{{1}} {\text{X}}_{{{1} + \ldots + }} \beta_{{\text{p}}} {\text{X}}_{{\text{p}}}$$
(4)

and the exponent of e is a linear regression defined by the behavior of the independent variables.

Equation 3 relates a dependent variable to a set of environmental variables. The dependent binary variable is defined by a set of landslide points from the NASA global catalogue, and from an equal number of random points representing the absence of the phenomena (point with value 1 indicating presence, and 0 indicating absence of an event). In this stage the LR has been run considering the independent variables of the month with highest exposed population (Fig. 7).

Rainfall data were obtained from the WorldClim database (Fick and Hijmans 2017). For each country, the average monthly rainfall considered in the national LR model has been the one for the month with the most exposed population to landslide hazard, (Fig. 7). High rainfall is related to high hazard probability.

Soil moisture was obtained from the NASA-USDA Global Soil Moisture Data (Bolten et al. 2010) and aggregated to monthly average soil moisture on Google Earth Engine (Sazib et al. 2018) considering the single month with the most exposed population (Fig. 7). High soil moisture is related to high landslide probability.

Land cover was obtained from the ESA Climate Change Initiative (CCI) database for 2015 (Defourny et al. 2012). After a land cover generalization to derive vegetation scarcity the following classes were selected: Wetland-0, Forest-1, Shrubland-2, Grassland-3, Cropland-4, Barren-5, Urban-6. The landslide probability is related to vegetation scarcity.

A proxy for soil properties was chosen by the K-factor elaborated with Wischmaier’s equation (Wischmaier 1958) based on the global soil map from the Food and Agriculture Organization-FAO (FAO et al. 2012). The K factor is an index used in erosion assessment to consider soil properties. The highest value of K has been related to high landslide probability values.

The slope was obtained from Google Earth Engine from the Shuttle Radar Topography Mission SRTM-NASA 500 m. This variable introduces gravity to the model and is directly related to the landslide probability.

In addition to the Rainfall monthly amount further three independents variable have been elaborated to introduce in the LR models the rain intensity aspect. Indeed for the estimation of the rainfall characteristics that could be regarded as triggering factors of landslides, the NOAA CDR Climate Prediction Center MORPHing technique (CMORPH) product was used (Xie et al. 2021, 2017). CMORPH is bias corrected and re-processed global precipitation product that covers area between 60°N and 60°S parallels with a spatial resolution of 8 km. The time resolution of this product is 30-min and it covers period from 1998 onwards. The method takes use of precipitation estimates derived based on the low Earth orbit satellite-based passive microwave observations (Kim et al. 2020). This products has been used for numerous applications around the globe and provided reasonable estimates for the precipitation patterns (e.g., Dis et al. 2018; Kim et al. 2020; Palharini et al. 2020). Based on the location of landslides, next variables were calculated at monthly time scale as average values in the period from 1998 until 2019:

  • Maximum 30-min precipitation rainfall sums (P30m);

  • 10 days moving precipitation rainfall sums (P10d);

  • Monthly rainfall erosivity values as calculated by Bezak et al. (2021) (R).

Monthly values of these three rain intensity variables were used as input to the LR model together with the other independent variables.

2.3 Subnational analysis

The last step of our methodology is the analysis at the sub-regional level, in the case study area in Bangladesh’s Chittagong district. This area was selected because it is easily accessible, and it was possible to collect a local landslide catalogue and it is a district bordering Cox’s Bazar. From a humanitarian perspective, Cox’s Bazar and the surrounding regions have crucial importance in hosting a high number of refugees and being strongly affected by the monsoon season (Ahmed et al. 2018).

2.3.1 Subnational study area

Chittagong district is located between 20′35° and 22′59° N latitude, and 91′27° to 92′22° E longitude. The landslide-prone hilly landscape has been formed by sedimentary rocks of the Tertiary (65–1.8 Ma) (Ali et al. 2014). The monsoon is responsible of the high monthly average maximum rainfall in July (74.70 mm), while the minimum precipitation is in January (0.66 mm) (BMD 2020). Almost 90% of the total yearly precipitation takes place between the months of June and October (BBS 2011). In the last four decades the demography experienced exponential growth, where the urban population increased by 583%, with a current 6992 people/km2 population density (Ahmed and Dewan 2017) (BBS 2011).

Between January 2001 and March 2017, 730 landslides occurred in Chittagong (Rabby and Li 2020).

Data suggest that Chittagong has a record number of deaths associated with landslides since 1999. From 1999 to 2018, landslides caused over 600 deaths in Chittagong (Alam 2020). For example, on 11 June 2007 landslide events alone caused 128 deaths and 100 injuries adjacent to hilly areas because landslides were triggered by heavy rainfall (610 mm) for eight consecutive days. Five years later, on 26 June 2012, another eight days of continuous rainfall (889 mm) triggered landslides that claimed 90 casualties (Ahmed and Rubel 2013). These landslide events occurred in hilly clear-cut areas with steep slopes. Slope failure in these fragile hilly areas occurred during the rainy season between June and September. It is to note that the population of Chittagong has increased by a factor of four since 1974 with a significant number of people living in highly vulnerable hill slope areas.

2.3.2 Dependent and independent local variables

To analyze the landslide probability in Chittagong district a specific LR model was created considering most of the independent variables of the national analysis.

Regarding the dependent variable, the landslide point was collected locally.

To develop Chittagong’s landslide catalogue, local histories, archives of institutional and administrative records, newspapers, reports, digital archives, and published peer-reviewed journal papers were consulted. Following the development of the landslide catalogue, field visits and investigation were conducted to validate the information collected from the secondary sources and to identify accurate locations of landslide occurrence. During the field survey, information on location name, coordinates (latitude, longitude), area of displacement mass, landslide mechanism (type of movement), causes of movement, and consequences (casualties, injuries, damages, impacts) were collected and validated. Coordinate values of landslide locations were determined using a handheld GPS. This data represented the dependent variable for the LR model at sub-national level.

An additional independent variable was introduced at sub-regional level, the distance from the topographic wetness areas. Topographic wetness area has been generated inside the GIS software considering the NASA digital elevation model of 30 m pixel spacing derived from the SRTM mission. This variable was implemented to introduce the local morphology and soil moisture conditions, indeed the previous soil moisture layer, with a resolution of 25 km, was considered not useful for a sub-regional scale.

3 Results

3.1 Global landslide exposure

From the landslide human exposure equation, equation number 2, it was possible to obtain the people living in hazardous areas globally. The landslide human exposure equation used average annual rainfall to obtain a first representation of the exposure layer without monthly disaggregation. Figure 5 shows the national clustering and Fig. 6 the regional clustering (map ‘a’ absolute value; maps ‘b’ normalized). From these two global maps, it is possible to identify some macro-scale hotspot areas like India, Nepal, China, Philippines, Indonesia, Ethiopia, Colombia, and Central America among others.

Fig. 5
figure 5

a People living in landslide hazard areas at national level according to landslide human exposure equation; b People living in landslide hazard areas normalized by to total population at national level according to landslide human exposure equation owner elaboration from equation number 2

Full size image
Fig. 6
figure 6

a People living in landslide hazard areas at regional level according to landslide human exposure equation; b People living in landslide hazard areas normalized by total population at regional level according to landslide human exposure equation, owner elaboration from equation number 2

Full size image

A second global analysis measured the landslide exposure monthly disaggregation by the monthly rainfall. From that step twelve landslide exposure layers were obtained according from monthly rainfall data. Figure 7 shows the mean monthly rainfall and the number of people exposed from landslides for Nepal, Philippines, India, Bangladesh, Afghanistan, which are systematically affected by landslides (Figs. 2, 5, 6).

Fig. 7
figure 7

Exposed population at national level by equation number 2, per Nepal, India, Afghanistan, Philippines, Bangladesh

Full size image

The monthly rainfall aggregation is useful to understand not only where but also when a landslide event is more likely to occur. Figure 7 shows the estimated number of people living in hazard areas related with the monthly rainfall amount.

3.2 National landslide exposure analysis by logistic regression

After the global analysis, a Logistic Regression model was applied to particular countries: Bangladesh, India, Indonesia, Nepal. Table 3 shows the explanatory variables for the country scale. Afghanistan was not considered because of the number of past landslide events in the Nasa catalogue was low (17 events reported since 2007). The Nasa catalogue is updated with Google alerts based on mass-media and news reports. For Afghanistan, the global landslide catalogue highlighted a limitation in reporting events in rural and sparsely populated areas where international news often fail to report on landslides. In the region of Badakhshan in the northeast of Afghanistan, an analysis of high-resolution satellite images recognized more than 600 past landslide events (Zhang et al. 2015).

Table 3 ROC values per countries per independent variables, plus final model variable (*R10 = cumulative 10 days rainfall; **R30 = maximum cumulative 30 min rainfall)
Full size table

Table 3 shows the ROC values per country and for each independent variable from the LR analysis. The LR model highlights the most important variable per country. Variables under the threshold of 0.5 were excluded from the model.

Rainfall and slope are significant in all four countries. Land cover has been considered according to the vegetation scarcity of its land cover types. The hypothesis that high vegetation scarcity leads to high landslide probability holds true except for India. Also, soil moisture has a direct correlation with landslide probability, except in Bangladesh.

A possible explication for the lack of explanatory power of land cover for landslide events in India is the large fragmentation that land cover has experienced in this country (Reddy et al. 2013). The global land cover layer in our analysis is the ESA land cover with 250 m pixel resolution which is unable to describe highly fragmented landscapes adequately.

Regarding the lack of correlation with the July Bangladesh soil moisture, this country has experienced a general increase in rainfall events with strong intensity in the last decades (Shahid 2011), which is not represented by monthly soil moisture averages. It is not the aim of this paper to perform a specific analysis on each independent variable but without any doubt, global change could have a dramatic domino effect where only biophysical factors are not sufficient to explain the probability of hazardous phenomena in particular landslide-prone areas.

The soil properties, represented by the K factor, were a less important variable and only had a high explanatory power for Indonesia with an ROC of 0.59.

In Fig. 8, the ROC values represented by the graph show the sensitivity and the specificity, which are defined as the true positive rate of landslide-prone areas, and the false positive rate of areas without landslides, respectively. The goodness-of-fit of the model is defined by how well the LR classifies the presence/absence of landslides and that goodness-of-fit is expressed by the ROC value (Cantarino et al. 2019). In Fig. 8, we can observe the variable with the largest Area Under the Curve (AUC). The amplitude of the area between the reference line (threshold) and the curve drawn by each variable represents the explanatory power of the model (Zweig and Campbel 1993). For example, we can easily recognize in India how important the slope is, as the variable with a ROC value of 0.82 is represented by a well-defined curve. At the same time, we can note for Indonesia and Nepal that there is not one strong explanatory variable. With the sum of the different independent variables, it is possible to obtain an explanatory model with ROC over 0.75 for Indonesia, and 0.7 for Nepal.

Fig. 8
figure 8

Area Under the Curve representing ROC values: a India; b Bangladesh; c Indonesia; d Nepal

Full size image

In Fig. 9, the behaviors of the most significant independent variables versus the landslide probability of the corresponding LR model are shown. The curves show the probability trends by the independent variables and the analytic thresholds. Figure 9 shows that according to the LR models, values with >50% probability of a landslide occurring present the following independent variables with the strongest influence:

  • India–Slope > 10 degrees;

  • Bangladesh–mean monthly rainfall > 600 mm;

  • Indonesia–land cover classes of grassland, cropland, barren, and urban land;

  • Nepal–mean monthly rainfall > 400 mm.

Fig. 9
figure 9

Probability trend considering the most significant variable behavior in the LR models per country

Full size image

3.3 Sub-national results of logistic regression

Table 4 shows the LR model results at local level for Chittagong district. The most significant variable is the erosivity index, and the variable with the weakest relationship is soil moisture. The distance from the topographic wetness areas has been introduced to compensate the coarse spatial resolution of soil moisture at sub-national scale because of the inadequacy of soil moisture data with 25 km spatial resolution. At district level, a spatial resolution of 25 km does not have enough capacity for an adequate spatial discrimination of patterns. The final model is defined from all the independent variables with the exclusion of soil moisture (ROC < 0.5). The model with all significant variables has a ROC value of 0.98.

Table 4 ROC values per Chittagong District per independent variables, plus final model with significative variables (*R10 = cumulative 10 days rainfall; **R30 = maximum cumulative 30 min rainfall)
Full size table

The rainfall intensity variables, at district level, reflect the national LR models results showing their importance, as well as rainfall erosivity (ROC 0.93), cumulative 10 days rainfall (ROC 0.81), and maximum 30 min rainfall (ROC 0.89) calculated for the rainiest month. It is also interesting the ROC (0.69) of the land cover variable, reflecting the local landscape reality.

Figure 10 shows the probability trend for erosivity index (a), and the cumulative 10 days rainfall (b), the significant variables for Chittagong district.

Fig. 10
figure 10

Probability trend considering a Erosivity index; b 10 days accumulated rainfall

Full size image

The erosivity index is an expression of rainfall intensity and energy. The part of the monthly maximum 30 min rainfall meanings records can also be considered as a proxy of different rainfall intensities across the year. Accordingly, both variables are related to the short period term rainfall behaviors, and in that sense, it was interesting across to observe the effects of the 10 days the cumulative rainfall threshold (560 mm).

For this subnational scale, a landslide probability layer (Fig. 11) was defined based on the most significant model from Table 4 and related independent variable layers (Fig. 12). A landslide probability map was created using kriging interpolation of the probability values and the data points where landslides occurred / did not occur. The landslide probability map discriminates in the space the probability expressed by ROC value according to the prediction capacity of the independent variable. The map shows the area around Chittagong City has high landslide probability, while it unfortunately is the most populated area, too.

Fig. 11
figure 11

Chittagong district landslide probability map from LR model

Full size image
Fig. 12
figure 12

Independent variables for the district level RoL model, Table 4

Full size image

In Table 5, the landslide points and points for which no landslides were recorded are classified according to the landslide probability layer ROC values. The table shows the odds ratio from the raster interpolation. From the table, it is possible to affirm that the model explains all past landslide events, 34 landslide event matching with ROC value 0.5 > , with just one misleading false positive.

Table 5 Probability map values for Chittagong district coincident to points with and without landslide occurrences
Full size table

4 Discussion

According to Chae et al. (2017), the landslide studies can be classified in (i) landslide susceptibility analysis, expressing the measurement of the probability considering a set of explanatory variables; (ii) landslide runout analysis, expressing the travel distance to outline the potentiality hazardous areas; (iii) landslide monitoring, analysis carried out to recognize the variation in slope stability measuring changes in parameters; and iv) landslide early warning, analysis strictly related to the monitoring focusing on the overpass of a certain threshold for early recognition of the hazardous event. According to this classification, the present study can be classified within the landslide susceptibility group, aiming at performing an analysis of the probability of landslide events across time and space. Concerning the LR models, they indicate the threshold values according to the overpass of ROC 0.5, an information useful to reinforce the early warning system at country and regional scale.

The workflow adopts a multi-scalar approach to spatially explicit landslide human exposure analysis. In the first step (the global overview) the output is an exposure map based on a combination of hazard (landslide frequency activated by high rainfall) and presence of people.In the second and third steps, the analysis considers landslide hazards at national and sub-national levels without crossing the result with the population.

The processing scheme builds up on previous global landslide analysis (Nadim et al. 2006), providing a global overview on landslide exposure disaggregated across 12-months. This considering the rain temporal distribution as a pondering factor able to highlight both where and when the landslide exposure is higher. A temporal disaggregation, considering the rain as trigger factor, has been considered already in literature (Froude and Petley 2018), putting in correlation the number of landslide events. In our study, the landslide hazard monthly disaggregation has been crossed with the population layer in order to estimate the exposed population on monthly basis. This approach aims to be innovative in the disaster risk reduction management suggesting not only where but also when a community is exposed to landslide events.

The monthly aggregation in the global approach (step 1) proved to be useful for understanding when the population is most exposed to landslides over the year. However, the national approach is a generalization that omits regional meteorological dynamics. From the histograms in Fig. 7, ‘exposed population at the national level in India’, the rainfall peak does not match perfectly with the exposed population peak. This is due to the different climatic areas that characterize a large country like India. Instead, in Nepal (Fig. 7) the coincidence between the exposed population with the mean monthly rain values is stronger.

Regarding the LR models, they have been defined considering the past landslide events. In the early warning landslide analysis to compare the previous conditions to understand future hazard situation is a recognized practice (Glade et al. 2000; Aleotti 2004).

The LR results defined the rainfall and slope as explanatory variables in all countries, while soil properties had the weakest influence on landslide hazard. This can be explained by the strong spatial characterization of the soil properties and because the global FAO soil data is too generalized for use at local scales (Batjes 1997).

The LR models identify the most important explanatory variables of landslide probability and their threshold values after which they become significant. For Bangladesh, mean monthly rainfall is the most significant variable and it indicates that if rainfall exceeds 600 mm per month the probability of landslides passes the threshold. Also, for Nepal rainfall is significant and the threshold in this case is 400 mm per month. These values are not a perfect match with the histograms derived at national scale (Fig. 7). This is because the histograms are remarking the values coming from the zonal statistic of the average monthly rainfall over the entire national area, while LR uses rainfall values centered on the points with and without a landslide. The national analysis uses the average amount of rain to estimate the number of people living in hazard areas and the related exposed month considering the entire country. Instead, the LR defines the statistical rainfall thresholds that indicate the threshold value at landslide point scale.

For Indonesia, the most significant variable is land cover. In this case, the threshold suggests that the land cover classes such as grassland, cropland, barren, urban are most prone to landslides. For India, the slope is the most significant independent variable, with a threshold of slope > 10 degrees.

The last step of the workflow in our methodology is the landslide probability analysis of Chittagong district, representing the analysis at the sub-regional level. In this administrative area, we introduced a new geomorphological independent variable to replace the coarse spatial resolution soil moisture dataset, which is the distance from the topographic wetness areas. This independent variables achieved a correlation of 0.61 respect the 0.25 from the moisture. The topographic wetness represents the local geomorphological conditions well enough to be useful in landslide probability assessment, and with a high concentration of soil moisture it identifies the area that captures a lot of water after a rain event (Dahal and Hasegawa 2008). Topographic wetness is a GIS, DEM based, derived product. Chittagong’s LR model results are in line with other authors highlighting the importance of topography variables at local level (Chen et al. 2019; Sevgen et al. 2019; Goyes-Peñafiel and Hernandez-Rojas 2021).

The intensity rain variables are the real asset of the Chittagong’s model confirming the importance of local climatic condition, these variables were not so important at national level, indeed.

The contribution of this work to the knowledge in the field can be so summarized as follows: (i) Multiple scale dimensions; (ii) Workflow replicable in a GIS environment considering open data set and a parsimonious approach; (iii) Temporal disaggregation of the global human exposure on monthly basis.

The landslide probability has been considered by multiple scale approach where the relative outputs can support an integrated approach to the disaster risk reduction interesting a wide range of disciplines remarking the importance of reducing the distance between scientific dimension and policy making territorial subject (Alcántara-Ayala 2020).

5 Conclusions

Landslide hazards have a strong spatial–temporal dimension. Disaster risk reduction and disaster management practices can benefit significantly from a spatio-temporal landslide risk assessment and monitoring system such as the one presented here. The multi-scale approach of the method presented here supports the spatial and temporal quantification of the landslide probability. This study highlights are as follows:

  1. (1)

    Rainfall is a trigger of landslides at the national level. Where the country has more than one bioclimatic zone, different analyses can be conducted.

  2. (2)

    The multi-scale approach is useful to get an overview of the hazard and exposition and subsequently to identify the local landslide probability and regional climatic realities

  3. (3)

    The LR method represents a methodology to identify the most important variables related to the landslide probability. However, the LR models can only be as good as the quality and the spatial resolution of the input data allow.

  4. (4)

    The output from the landslide physical exposure assessment system in the form of statistical information and multi-temporal exposure maps are useful for disaster risk reduction programs to decrease the hazard. A simplified summary of these results could be sent to the local community in order to raise awareness of landslide risks.

  5. (5)

    At national level biophysical variables like rainfall, and soil moisture are relevant to landslide exposure, but at sub-regional level geomorphological variables such as slope, land cover, topographic wetness, have greater explanatory power. At the same time specific analysis considering the rain intensity variables such as rainfall erosivity, 30 minutes maximum rainfall and 10 days accumulated precipitation rainfall, reveal the importance in the preparedness phase highlighting empiric threshold values and crucial geostatitic correlation, in light of climate change scenario too.

  6. (6)

    The workflow of landslide exposure mapping and monitoring introduced in this paper can be used to produce thematic maps. From overlaying the territorial elements (e.g., urban settlements) and susceptibility factors such as topographic wetness areas and past landslide events it is possible to delineate local exposed areas.

The global approach of this work uses landslide frequency tables (Tables 1, 2) based on Nadim (2006) which provide vital information for the global landslide analysis and correlating the global landslide frequency map (UNEP 2018) with the mean monthly rainfall (Fick and Hijmans 2017). The global landslide catalogue from NASA (Kirschbaum 2010) has been crucial for this work and this data represents the state of the art. The methodology to identify landslide points represents, however, a generalization that could affect the quality of the LR model to a certain degree. The landslide physical exposure analysis performed in the first and second steps at global and national level could be recalculated considering the regional or district level to disaggregate the exposed population to a different administrative unit.

The multi-scale approach presented here can facilitate humanitarian support to disaster risk reduction. The outputs from different scales of analysis could be used to support different actors involved in the disaster management cycle. These outputs could suggest to policy-makers and donor agencies where to focus the disaster risk reduction interventions. The landslide probability analysis can pinpoint where to decrease the hazard by improving the territory with geomorphological assets and infrastructures, and where to decrease the final risk avoiding the proximity between hazards, exposure, and vulnerable areas.