Towards a Transferable Antecedent Rainfall—Susceptibility Threshold Approach for Landsliding

Abstract

Determining rainfall thresholds for landsliding is crucial in landslide hazard evaluation and early warning system development, yet challenging in data-scarce regions. Using freely available satellite rainfall data in a reproducible automated procedure, the bootstrap-based frequentist threshold approach, coupling antecedent rainfall (AR) and landslide susceptibility data as proposed by Monsieurs et al., has proved to provide a physically meaningful regional AR threshold equation in the western branch of the East African Rift. However, previous studies could only rely on global- and continental-scale rainfall and susceptibility data. Here, we use newly available regional-scale susceptibility data to test the robustness of the method to different data configurations. This leads us to improve the threshold method through using stratified data selection to better exploit the data distribution over the whole range of susceptibility. In addition, we discuss the effect of outliers in small data sets on the estimation of parameter uncertainties and the interest of not using the bootstrap technique in such cases. Thus improved, the method effectiveness shows strongly reduced sensitivity to the used susceptibility data and is satisfyingly validated by new landslide occurrences in the East African Rift, therefore successfully passing first transferability tests.

1. Introduction

Rainfall-triggered landslides pose a severe threat to societies on all continents [1,2]. Rainfall thresholds are therefore essential for characterizing landslide hazard and developing early warning systems [3,4,5]. Empirical approaches define thresholds on scales ranging from local [6,7] to regional and global [8,9], based on the observed relation between dated landslides and rainfall characteristics such as intensity, accumulation, duration, or antecedent rainfall (AR) conditions [10]. However, rainfall is only a proxy for what is regarded as the main trigger of landslides, i.e., the development of high pore-water pressure in the subsurface, constrained by water infiltration [11,12]. Interacting with retention and drainage processes [12], infiltration is a highly complex process affected by a myriad of factors such as soil physical properties (e.g., soil suction head, porosity, hydraulic conductivity) and their variations through the soil column [13,14,15], presence of cracks [16], hillslope morphology [17], vegetation [11,18,19], antecedent rainfall conditions [15,20,21,22], and rainfall intensity [23,24]. In contrast to the empirical threshold definitions, process-based approaches incorporate such hydrophysical parameters through a spatially extended infinite-slope stability model [25]. However, the large required data input for well-calibrated process-based thresholds explains their current limitation to mostly applications at the hillslope scale or through numerical simulations [4,21,25,26,27].

The estimation of empirical rainfall thresholds is also associated with additional sources of uncertainty. Firstly, landslide inventories are inherently biased towards high-impact landslide events and regions that are most accessible, while their accuracy is constrained by the scientific validity of the reporting sources, especially in data-scarce low-capacity environments [1,28,29,30,31]. Secondly, rainfall data comprise uncertainties related to the spatial representativeness of rain gauges or biases in satellite-derived estimates [32,33]. Thirdly, the definition of rainfall parameters, with intensity and duration forming the most frequently used parameter couple [3,5], varies strongly across studies [3]. Finally, the latter parameters’ interdependence is problematic, obscuring the physical processes associated with the calculated thresholds [34].

In order to account for and characterize threshold uncertainties, a growing number of reproducible statistical techniques have been developed [3]. A weakness of such methods is, however, that they are generally tailored to a specific area and available data sets, which often prevents straightforward transferability to other regions and data sets [35]. Nevertheless, transferability is not only essential for evaluating and comparing landslide hazard over different regions of the world [10,36], but also valuable in the context of the increasing availability of ever higher-resolution data relevant for threshold analysis, such as rainfall estimates from global-scale satellite data [32].

The most influential statistical threshold techniques include the probabilistic approach through Bayesian inference [10,37], the use of receiver operating characteristics (ROC) analysis with different optimization metrics [38,39], and the frequentist approach developed by [40]. The Bayesian and ROC approaches compare conditions that resulted or not in landsliding, the former fundamentally relying on prior and marginal probabilities [37] and the latter attempting to balance the true and false positive rates derived from a confusion matrix [39]. When rainfall data are only available for conditions that triggered landslides, the frequentist method provides a quantitative way to exploit it and calculate thresholds. This method, as developed by [40] for the (intensity, duration) parameter couple of rainfall, calculates the least-square fit of the log-transformed data and fits a Gaussian function to the probability density function of its residuals. Next, the Gaussian curve is used to adjust the intercept of the best fit equation to the desired threshold, expressed in terms of exceedance probability [40]. Practically this means that for a threshold at, e.g., 5% exceedance probability level, there is a 0.05 probability that any landslide be triggered by rainfall conditions below the threshold. The quality of the thresholds obtained by this method depends on the size of the data set and its good covering of the whole range of the parameters used [40]. An improvement of the frequentist method lies in the adoption of a bootstrapping statistical technique to assess the parameters’ uncertainty in the power-law threshold model [9]. Here, the bootstrap procedure involves many threshold calibrations (e.g., 5000 [9]), each of which based on n randomly sampled data (with replacement) from a data set of size n. The final threshold parameters and associated uncertainties are calculated as the mean and their standard deviations, respectively, of their many estimates. This approach has proved to be transferable over different regions where abundant information on landsliding and rainfall was available [9,41,42].

Recently, this frequentist approach with bootstrapping [32,41,43] has been modified by [35] through coupling a dynamic rainfall variable (AR) with a static indicator of the spatially varying predisposing ground conditions (landslide susceptibility, S) (further referred to as the AR-S approach). The first step in AR-S threshold estimation is similar to the frequentist method developed by [40] and [9], calculating the residuals of the least-square fit on the log-transformed data. Then, it proceeds to select 2x% of the data with the largest negative residuals, on which a new least-square regression is applied, providing a threshold at the x% exceedance probability level. In this way, not only the intercept α but also the slope β of the threshold equation are based on the smallest AR data able to cause landsliding. In parallel, the following novel AR index, covering a period of 42 days (n) of antecedent rainfall [35], was proposed to account for the non-linear decay of the effect of rainfall on soil wetness.

$$A{R}_i={\displaystyle \sum }_{k=i}^{i-n}{e}^{\frac{-a\times \left(t_i-t_k\right)}{r_k^b}}\times r_k,$$

(1)

with t referring to time (here expressed in days), and the characteristic time $\tau =r_k^b/a$ varying non-linearly with daily rainfall r_k [35].

Identifying thresholds for rainfall-triggered landsliding in data-scarce environments is challenging with respect to information on landslide occurrence and hydrophysical parameters, resulting in the quasi-absence of research on this topic in regions such as Central Africa [3] despite high hazard potential [29,44,45,46,47]. The AR-S approach allowed defining the first regional threshold for landsliding in the western branch of the East African Rift (WEAR) [35]. To the authors’ knowledge, it has so far not been used in other regions. Moreover, the cited study relied on limited data available on landslide occurrence, global satellite-based rainfall estimates [48], and continental susceptibility data [45]. There is hence a strong need for testing the method’s robustness with other data sets. A regional S model is now available for the WEAR [49], which outperforms the global and continental models with regard to prediction accuracy and geomorphological plausibility [49]. Moreover, the landslide event database used in [35] has now grown by about 27%. In this paper, our aim is thus to use these new data and test the transferability of the AR-S threshold method as designed by [35] to these new data.

2. Study Area and Data

2.1. Underreported Landslide Events in the WEAR

The WEAR covers an area of ~350,000 km² in tropical Africa (Figure 1). This highly populated region is characterized by high rainfall intensities, recent seismicity, deeply weathered substrates, and a complex rift topography [33,47,49,50,51,52]. These factors render the area highly susceptible to landsliding [45,49]. Indeed, recent studies incorporating observations from satellite images and fieldwork reveal high landslide activity, with hundreds of recent landslide events over the last 10 years accounting for more than 5000 individual landslides mapped in the area [44,49,50,53]. An event is defined as a single landslide or a group of landslides with a common trigger over the same area [29]. Landslide types vary greatly from dominant shallow slides and earth and debris flows to less frequent rock avalanches and deep-seated rotational slides [44,47,49,50]. Rainfall has been identified as the prevailing trigger of these recent landslides [29,35,47,49,54]. Despite this large number of landslide events, usable information about the day of their occurrence is rare [29]. Actually, for a period of observation similar to the one covered by the satellite rainfall estimates, [35] compiled an unprecedented landslide event inventory that comprises 145 events for which the location and day of occurrence is known. This clearly highlights the fact that landslides are severely underreported in this poorly accessible region [29].

In this research, the calibration data set (‘CAL’ in Figure 1) used for the threshold approach is the same as in [35]. The 145 dated landslide events are located with an individual accuracy better than 25 km and a mean accuracy of 7.2 km. The inventory covers the period from 2001 to 2018, with most landslides occurring from March to May, after the second rainy season [35] (Figure 2).

The event inventory used for validation (‘VAL’ in Figure 1) consists of newly acquired information on 39 additional events located with similar individual accuracy and a mean accuracy of 2.3 km. Information about these events was obtained through field observations and newly identified online media sources, explaining their 2002–2019 temporal coverage overlapping with that of the calibration inventory. Their seasonal distribution confirms the temporal pattern previously observed in [35] (Figure 2).

Both the calibration and validation data sets are strongly biased in space and time, due to the severe constraints inherent to the political, economic, and environmental context in the WEAR [29,35,55]. Considering the restricted accessibility of many parts of the study area, fast vegetation regrowth or land rehabilitation, and areas with poor temporal (cloud-free) satellite coverage, uncertainties are too high to claim absence of landsliding in places and periods where no event has been reported. A distinction between landslide types cannot be asserted, owing to inadequate information in media reports which constitute the principal source of both landslide event inventories. Information on the sub-daily timing of the landslide occurrence is rare. Therefore, when the day of occurrence of the event is known, there is still an uncertainty mainly related to the ambiguous reporting of landslides that occurred during the night when the date sequence is not mentioned. For this reason, we consider that a reported landslide may have occurred randomly at any time over a 36-h period centered on the reported day [35].

2.2. Satellite-Based Rainfall

Due to the scarcity and poor spatial representativeness of rainfall data from ground observations [33,35], we rely on satellite-based rainfall estimates (SRE) from TRMM Multisatellite Precipitation Analysis 3B42 Real-Time, version 7 (hereafter spelled TMPA-RT). TMPA-RT data are freely available with a latency of 8 h over 50° N–50° S, at 0.25° × 0.25° and three-hourly spatiotemporal resolution, covering the period 2000 to present without gaps in space and time [48]. A recent paper on SRE over Central Africa [56] shows that TMPA has overall good skills in detecting and estimating daily rainfall as compared to ARC, CHIRPS, CMORPH, PERSIANN, TAPEER, TARCAT (see [56] for these acronyms’ meanings). In their study, the Research Version of TMPA was used, a product that is outperformed by TMPA-RT over the WEAR with regard to rainfall detection skills and absolute errors [33]. Moreover, the short latency of TMPA-RT compared to the two-month latency of the Research Version, is of crucial importance in the context of early warning systems. Despite the relatively good performance of TMPA-RT in Central Africa, the WEAR is a challenging environment for SRE due to its complex topography, high rainfall variability, and presence of large lakes, with a resulting mean bias in daily rainfall estimates in the order of ~40% [33]. TMPA-RT three-hourly rainfall data have been downloaded from NASA Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/, last access: 14 April 2019) for the period 2000 to 2019 and accumulated to daily rainfall to maintain consistency with the temporal resolution of the landslide inventory.

2.3. Susceptibility Models

Two S models are used in this study. The continental-scale S model of [45] is calibrated for all landslides regardless of type at a spatial resolution of 0.0033°. This model is produced through logistic regression using a ~4:1 landslide to no-landslide (L/NL) ratio and is based on four predictor variables: maximum slope (~90 m SRTM [57]), mean local relief (~90 m SRTM [57]), peak ground acceleration [58], and lithology [59]. The landslide inventory used for the model contains more than 18,000 landslides, of which 765 are located in the WEAR. The second S model is the regional-scale model of [49] which was calibrated for a representative part of the WEAR and extrapolated within this study for the entire WEAR. This model includes all landslide types and is trained at a 0.0003° resolution using logistic regression with a 1:1 L/NL ratio based on a local inventory and 11 global/continental predictor variables [49]: slope (~30 m SRTM [57]), peak ground acceleration [58], distance to active faults and inactive faults [52,60], lithology [59], land cover [61], distance to drainage network (~30 m SRTM [57]), planar curvature (~30 m SRTM [57]), profile curvature (~30 m SRTM [57]), aspect (~30 m SRTM [57]), and two-day 15 mm rainfall accumulation threshold exceedance [62]. Note that the rainfall predictor was of minor importance in the model and had no significant impact on the susceptibility pattern in the study area [49]. The inventory contained more than 6000 landslides and the regional model shows predictive power and geomorphological plausibility that strongly outperform the continental model [49].

In order to exploit AR and S data at the same spatial resolution, both S models are resampled to the coarser 0.25° resolution of TMPA-RT data while assigning the 95th percentile of the original values to the coarser pixels (Figure 1). The S range of the continental-scale model for pixels containing calibration (validation) landslides is 0.38–0.97 (0.31–0.97) with mean and standard deviation equal to 0.80 ± 0.15 (0.79 ± 0.16). The regional-scale S data range is 0.10–0.72 (0.12–0.72) with mean and standard deviation equal to 0.57 ± 0.14 (0.49 ± 0.15). The difference in the data range between the two S models mainly results from their different sampling strategies (L/NL). Furthermore, S values are scaled for different geographical extents, with the continental-scale S model comprising areas that are not representative for the WEAR.

3. Problem Statement

We applied the AR-S threshold method according to [35] at the 5% and 10% exceedance probability levels, using the same calibration landslide data set, the same TMPA-RT-based AR data, but the new regional-scale S data of [49]. We obtained the following general AR-S relation and threshold equations:

$$AR=(\alpha \pm \Delta \alpha )\times S^{\left(\beta \pm \Delta \beta \right)}=\left(38.8\pm 1.6\right)\times S^{-0.06\pm 0.06} \left(R^2=0.00\right)$$

(2)

$$AR \left(5\%\right)=\left(13.1\pm 1.7\right)\times S^{0.24\pm 0.16} \left(R^2=0.05\right)$$

(3)

$$AR \left(10\%\right)=\left(17.2\pm 1.7\right)\times S^{0.22\pm 0.16} \left(R^2=0.03\right).$$

(4)

Contrary to [35], the close to zero determination coefficients R² (averaged from 5000 bootstrap iterations) associated with the two calculated thresholds show no dependence of threshold AR values on S (Equations (3) and (4)). The meaningless character of these threshold estimates is further confirmed by the positive slope of the regression lines suggesting counter evidence that higher rainfall would be needed to trigger landslides in more susceptible areas (Figure 3). Analysis of the individual bootstrap iterations likewise uncovers a major issue lying in the estimation of parameter β, which is significant in only ~1 of 2 iterations, with relative uncertainties ∆β/β of 0.7 on average, much larger than the generally accepted 10% level [9].

Such poor thresholding cannot be ascribed to low-quality S data, the regional-scale data of [49] having been shown more accurate than the continental-scale S data of [45]. The reason for very weak and unrealistic positive correlation between AR and S has thus to be found elsewhere, most certainly in some hidden deficiency of the AR-S threshold method of [35]. We suggest and test hereafter that the problem arises from the way the data subset used in the threshold calibration is defined in the frequentist-based approach, based on the selection of the most negative residuals of the general fit. Indeed, in the case of the relatively small data set available in the WEAR and the unequal spread of the data across the S range, the frequentist method’s assumption that the data set is large and well-spread [40] is not satisfied. In particular, using the regional S data, the distribution of the data points within the AR-S space is such that the 10% and 20% subsets sample (i.e., 2x%) comprise almost no data in the domain of low S, due also to the quasi horizontality of the general fit that forces the location of the most negative residuals in the high-S region (Figure S1). This means that a large number of the bootstrap iterations are based on data belonging exclusively to a narrow range of high S values, biasing the threshold AR-S relation and degrading the method’s robustness. In any case, this failed test of the method highlights the need for improving it in order to overcome limitations imposed by heterogeneously distributed and relatively small-sized data sets. It also points to the possible role of the bootstrap procedure and calls for a critical evaluation of its use in such contexts. We thus propose two major methodological modifications of the AR-S approach in the next sections.

4. Improving the Data Distribution over the S Range of the Data Used for Threshold Calculation

4.1. Rationale

While sticking to the choice of [35] of defining x% thresholds from 2x% subsets of data, we first propose a major modification of the AR-S method aimed at optimizing the use of the information available over the entire S range. Fundamentally, the data that are now considered for inclusion in the calibration subsets rely no longer on residuals of an often non-significant fit over the whole data set but rather on minimum AR values. The best possible distribution of the latter is obtained by stratified sampling, dividing the actual S range of the data set in a number of slices from which (as much as possible) equal numbers of minimum AR data are selected. The slices are taken of equal size in log(S) and their optimal number was fixed at 10 on a trial-and-error basis. As an example, suppose you want to estimate a 10% threshold based on a data set containing 150 landslide events, i.e., 450 event dates (see [35] and A.1 below). Homogeneously distributing over 10 S classes the 90 data of the 20% subset required for this threshold calculation implies to select the nine events with lowest AR in each class. Obviously, some classes may contain less than nine data, thus contributing less to the composition of the data subset, whose final size will often be slightly smaller than expected. In addition, when an S class does not contain enough data to fully contribute to the subset, all its data will be selected, however far their AR values are from minimum. However, tested through down weighting of the data proportionately to the deficit in contribution of their provenance class, this possible bias appeared to insignificantly affect the threshold estimates. The modified method is described in detail hereafter (see also Figure 4). The source code is provided in the Supplementary Material (Code S1).

• Data preparation.

  A.1.

  AR values associated with each day of a reported landslide plus the days prior and after these dates are extracted from the AR time series of the corresponding pixels calculated according to Equation (1) and the parameterization adopted in [35], i.e., a = b = 1.2, n = 42 days, for which the index is relevant for landslide types ranging from shallow to deep-seated landslides [35,63]. Data with AR < 5 mm are discarded from the data set as unlikely to have been triggered by rainfall [35]. The size of the provisional data set Q is then q ≤ 3p, where p is the number of landslide events in the raw calibration set.

  A.2.

  The data are weighted to account for the event date uncertainty: w = 24/36 for the day a landslide was reported, w = 6/36 for the days prior and after the landslide was reported. This weighting is implemented by expanding the data set as described in [35]. The expanded set is noted R.

• Threshold calibration.

  B.1.

  The number t_C of data to be selected per S class is determined as

  $$t_C=\frac{2\times \mathrm{TPE}\times r}{10},$$

  (5)

  where TPE refers to the desired threshold probability of exceedance, r is the number of data in R, and 10 is the number of log(S) classes.

  B.2.

  The data of R are grouped by S class. For each S class, data with the lowest AR values are selected until they amount to t_C. The set of selected data points over all S classes is referred to as T and contains a number of data t ≤ (2 × TPE × r).

  B.3.

  Thresholds are then calculated through linear least-square regressions using the log-transformed AR and S data from T and the bootstrap technique as in [35] to obtain threshold relations in the form of Equation (2).

• Threshold evaluation
  • Threshold quality is evaluated through the correspondence between the obtained false negative rate (FNR, actual ratio of data in R below the calculated threshold) and the nominal TPE. Differences may result from t significantly smaller than (2 × TPE × r), large outliers in T, and possibly also from bootstrap issues (see Section 5).

4.2. Increased Efficiency of the Method

We test the modified AR-S method using the regional-scale S data [49] and the calibration landslide data set from [35] to calculate thresholds with 0.05 and 0.10 exceedance probability. The 145 landslide events constituting the calibration set yield 435 weighted event dates, of which eight are discarded from the analysis because they do not meet the AR ≥ 5 mm requirement, thus 427 data instances remain in the threshold analysis (constituting Q, Figure 4). The 5% probability level is most frequently used in landslide hazard and early warning studies [9,41,64]. We also include the 10% level because threshold estimation relies then on a larger data subset T (at 5% probability level, t = 85; at 10% level, t = 171). Significance measures mentioned throughout the paper are associated to the significance level p = 0.05. The R open-source software, release 3.4.3 (http://www.r-project.org, last access: 14 April 2019) was used for all analyses. AR thresholds at the 5% and 10% exceedance probability levels were estimated as (Figure 5)

$$AR \left(5\%\right)=\left(4.8\pm 0.6\right)\times S^{\left(-1.16\pm 0.08\right)} \left(R^2=0.69\right)$$

(6)

$$AR \left(10\%\right)=\left(6.4\pm 0.6\right)\times S^{\left(-1.08\pm 0.07\right)} \left(R^2=0.62\right)$$

(7)

Contrary to the unrealistic results obtained from the original AR-S approach (Equations (3) and (4)), we get here plausible marked inverse relations between S and AR [65,66]. Moreover, the threshold equations are now associated with meaningful average R² coefficients of 0.69 and 0.62. All bootstrap iterations provide significant α and β parameters for both thresholds. We remind that here the bootstrap procedure consists in repeating the threshold calibration phase 5000 times, each iteration being based on a random sampling (with replacement) out of the R data set until the number of sampled data equals that of the r points of the data set. The subset of lowest-AR data is then selected from the random sample before threshold estimation. The mean and standard deviations of the 5000 estimates of α and β define the parameter values and uncertainties (∆α and Δβ). The results indicate an excellent performance of the modified AR-S threshold approach where the spread of the data subset for threshold calibration is forced over the entire S range. Obviously, strongly negative slopes result in decreased values of intercept α in Equations (6) and (7) as compared to Equations (3) and (4), respectively.

Though performing satisfyingly well, the modified AR-S threshold method leaves two minor issues open. The first one is related to the very close parameterization of the 5% and 10% thresholds and finds its cause in the similar actual FNRs of 0.05 and 0.07 obtained from 5% and 10% thresholding, respectively. In particular, the too low actual FNR associated with the 10% threshold equation betrays the real nature of the problem, which lies in the insufficient number of data in the low-S classes preventing the constitution of a complete data subset T to estimate the desired threshold. This issue is independent of the size of the original data set because, however large the number of recorded events might be, their distribution across the S range will remain similarly unequal, with low-S classes relatively deficient in data, especially for thresholds with higher exceedance probability demanding larger calibration subsets. Owing to the specific distribution of the data in the AR-S space, the AR-S approach inevitably implies to make a trade-off between high exceedance probability levels and degraded distribution of the data from which the threshold is estimated. Fortunately, more conservative low-exceedance probability thresholds (typically 5%) are the least affected by this issue.

High relative uncertainties (in the order of 10%) on parameter α might be another source of concern. However, beyond being subjective, the criterion chosen by [9] to qualify the threshold quality, namely a > 10% relative uncertainty, is barely usable in the AR-S space, where the many outliers in data distribution alter the efficiency of the bootstrap technique of uncertainty estimation (see Section 5). Moreover, in addition to the fit uncertainty, the bootstrap-based errors on the parameters obtained here from our weighted approach include the event date uncertainty and are also affected by the effect of the partly erratic character of the data distribution, inherent to the combination of ground (S) and meteorological (AR) variables on which the method relies. We thus conclude that the benefits of a method yielding thresholds directly modulated by the environmental conditions greatly outweigh the shortcomings of slightly higher uncertainty mainly on the threshold line intercept.

5. Bootstrapping Called into Question

The non-parametric bootstrap statistical techniques, including that introduced by [9] in the frequentist approach of threshold estimation, were designed to estimate the sampling distribution of a variable based on an empirical data set and assign measures of accuracy to statistical estimates [67]. While [9] acknowledge that, owing to the use of the same data for calculating the regression and estimating its parameters’ uncertainties, the bootstrap may yield optimistic estimates of the latter, other possible drawbacks are not discussed in studies having incorporated the bootstrap technique in threshold estimations [35,41,64]. However, the bootstrap may fail when the data set is incomplete, resulting in overestimation of the uncertainty, or when there are outliers in the data set, to which least-square regression estimates are highly sensitive [67]. Therefore, in the light of the observed uncertainty level and hints of variability in the bootstrap results, we decided to evaluate the pros and cons of applying this technique by performing a run of threshold estimation without using it.

Performing a single threshold calculation (no bootstrap), we obtained the following AR thresholds (Figure 6a):

$$AR \left(5\%\right)=4.6\times S^{-1.18} \left({\mathrm{R}}^2=0.70\right)$$

(8)

$$AR \left(10\%\right)=6.2\times S^{-1.10} \left(R^2=0.65\right).$$

(9)

Parameters α and β are significant for both threshold levels, with α barely smaller and β barely larger compared to the thresholds obtained using the bootstrap method (Equations (6) and (7)), thus well within the bootstrap-defined uncertainty boundaries. Opposed changes in α and β might be anticipated from the inverse correlation that links coefficient and exponent of power law fits to a given data set. Therefore, the two parameter changes damp each other, thus inducing almost no difference in thresholds calculated with or without bootstrap (

Table 1. AR threshold values (in mm) at 5% and 10% exceedance probability with (Equations (6) and (7)) and without (Equations (8) and (9)) bootstrap for the extreme susceptibility values S observed in the data set.

Threshold	with Bootstrap	without Bootstrap
5% threshold, S = 0.10	69.4	69.6
5% threshold, S = 0.72	7.0	6.8
10% threshold, S = 0.10	76.9	78.1
10% threshold, S = 0.72	9.1	8.9

). Using no bootstrap, only the information about fit uncertainty is lost because date uncertainty is still accounted for through data weighting. Moreover, in the case of the AR-S approach, the inherent poor S-spread of the data and the presence of large outliers in the data subset used for threshold estimation imply that the bootstrap procedure, which, sampling with replacement n data from a set of size n, is nothing more than a kind of random data weighting, includes a number of iterations with oversampled outliers. These iterations yield erratic results and may alter the final mean threshold estimate and exaggerate the fit uncertainty to an unknown extent. This is highlighted here by the better coefficients of determination of the AR-S thresholds obtained from the no-bootstrap approach. Furthermore, with or without bootstrap, the AR-S method does not account for crucial uncertainties affecting AR and S data themselves, so that providing bootstrap-derived uncertainties is actually misleading. We thus conclude that the AR-S threshold procedure is more meaningful when no bootstrap is applied. The corresponding source code of the AR-S threshold method is provided in the Supplementary Material (Code S2). As for the other issue affecting the modified AR-S approach mentioned in the previous section, namely the bias in higher exceedance probability threshold estimates (FNR < TPE), it is essentially linked to the lack of data in the low-S classes. It is thus independent of the use of a bootstrap technique and cannot be solved by discarding the latter.

6. Robustness of the Modified AR-S Threshold Method

6.1. First Test: Sensitivity to New Data on Landslide Occurrence

The modified AR-S method with no bootstrap is tested firstly by using the recent addition to the WEAR data set of dated landslide events. Taking into account their date uncertainty, the 39 landslide events constituting this validation set yield 117 new weighted event dates, of which four are discarded from the analysis because they do not meet the AR ≥ 5 mm requirement. The 113 remaining data instances (constituting Q, Figure 4) are distributed in the log(AR)–log(S) space in such a manner that 6% of them are located below the 5% threshold line derived from the calibration (Equation (8)) and 8% below the 10% threshold line (Equation (9)), indicating a good performance of the calculated thresholds (Figure S2) considering the small sample size.

Another test using the validation set, which in the same time should improve the accuracy of the calibrated thresholds, has consisted in combining the data of the calibration and validation sets into a larger data set of 540 event dates in order to recalculate the thresholds. The new thresholds read as

$$AR \left(5\%\right)=4.5\times S^{-1.14} \left(R^2=0.66\right)$$

(10)

$$AR \left(10\%\right)=6.1\times S^{-1.08} \left(R^2=0.59\right)$$

(11)

and do not much differ from those derived from the calibration set only (Equations (8) and (9)) (Figure 6b), confirming the relevance of the modified AR-S method. Though slightly decreased by additional noise brought in the middle- to low-S classes by the new data (Figure S2), their coefficients of determination remain highly significant. Likewise, their FNRs (0.04 and 0.06 for the 5% and 10% thresholds, respectively) are slightly degraded mainly as a result of an increased deficit in data in these S classes. Owing to the larger size of the data set, we nevertheless consider these thresholds (Equations (10) and (11)) more reliable than those based only on the calibration set, especially the 5% threshold, for which FNR ≈ TPE.

6.2. Second Test: Robustness to Different S Data Sets

We test the modified AR-S approach for the adoption of a different data set for S, using the continental-scale S data [45] and the complete (calibration+validation) data set of landslide events, obtaining the following AR thresholds (Figure 7)

$$AR \left(5\%\right)=5.7\times S^{-2.10} \left({\mathrm{R}}^2=0.73\right)$$

(12)

$$AR \left(10\%\right)=7.6\times S^{-2.08} \left({\mathrm{R}}^2=0.61\right)$$

(13)

with significant and meaningful values for R² and threshold parameters α, and β. Moreover, these thresholds show a stronger relation between threshold AR values and S with increased values for parameters (α,) β and R², explained by the increased dispersion of the data over the S range (Figure 7) relative to when the regional-scale S data was applied (Figure 6b). Where the AR-S approach developed by [35] posed problems for adopting a different S model than that used for its development (Figure 3), these results show that the modified AR-S approach proved to solve this matter. The threshold at the higher exceedance probability remains affected by a bias similar to that in Equation (11) with the actual FNR lower than the TPE (FNRs equal 0.05 and 0.07 for the 5% and 10% thresholds, respectively).

Another indicator for the robustness of the modified AR-S method is the generally satisfying correspondence between AR threshold results for the continental- and regional-scale S models (

Table 2. AR threshold values (in mm) calculated using continental- (Equations (12) and (13)) vs. regional-scale (Equations (10) and (11)) S data, provided at 5% and 10% exceedance probability for the extreme susceptibility values S observed in the data sets. The estimations are based on the complete calibration+validation data set of landslide events.

Threshold	Continental S Data	Regional S Data
5% threshold, min S	66.7	62.1
5% threshold, max S	6.1	6.5
10% threshold, min S	86.9	73.3
10% threshold, max S	8.1	8.7

). The single main discrepancy is observed for the 10% threshold in the low S data range, related to the actual FNRs of these thresholds being significantly smaller than their TPE. The latter is explained by the sensitivity of the threshold slope to the deficient number and exact location of data in the low-S classes, causing the largest threshold difference to appear for the 10% threshold (implying a greater lack of data in low-S classes) at the low end of the S range. By contrast, the intercept of the threshold equations, being located in the AR-S space with the highest density of data, remains quasi stable for different S models (Figure 6b, Figure 7).

Because of the enhanced relation between threshold AR values and S in Equations (12) and (13), it is tempting to suggest that thresholds based on the continental-scale S data would be more efficient when adopted in a landslide early warning system. However, the spatial pattern of the AR thresholds based on the regional-scale S data are closer to the reality, given that this regional-scale S model has a higher predictive power and geomorphological plausibility as compared to the continental-scale model [49]. The respective AR threshold maps are presented at the 5% probability of exceedance level in Figure 8. In general, we observe lower AR thresholds within the Rift. Nevertheless, there are some major differences between the two threshold maps, caused by differences in the quality of the susceptibility models. First, the threshold model using the continental susceptibility map of [45] assigns low AR thresholds to the rainforest in DR Congo south of the equator, despite the fact that the area is characterized by high amounts of rain and few landslides [29,35]. Second, the regional susceptibility data of [49] overall shows a much lower threshold in Uganda. In conclusion, we confirm the earlier observation that the AR thresholds based on the regional-scale S data and the currently most extensive landslide event inventory are to date the most accurate available thresholds for landsliding in the WEAR (Equations (10) and (11)).

7. Relevance to Landslide Hazard and Early Warning Studies

The modified AR-S approach is relevant for the increased accuracy of the resulting AR thresholds, which is partly also on account of the regional-scale S data [49] and the enlarged landslide inventory. In particular, the conservative low-exceedance probability thresholds are most reliable, being least affected by a degraded distribution of data used for the threshold calibration (FNR ≈ TPE). Depending on the local susceptibility, thresholds at the 5% exceedance probability range from AR = 62 mm in the least susceptible areas to AR = 7 mm in the highest susceptibility pixels, respectively, for which landslide have been reported (Equation (10)). These triggering AR conditions might seem low at first sight when compared to values obtained in other studies that look into antecedent rainfall conditions based on gauge measurements, e.g., a required minimum of 139 mm cumulated over 20 days to trigger landslides in the NE Himalaya [68]; a mean triggering rainfall accumulation of 376 mm for periods ranging between 15 and 40 days in NW Spain [69]; a critical rainfall amount of 450 mm over a two-week period in the greater Durban region in South Africa [70]. However, the triggering values obtained in our study are conceivable given the following main factors contributing to their relative lower values: (1) the exponential decay function applied in our AR calculation (Equation (1)) in contrast to the values obtained in the above cited studies through mere accumulation; (2) the high weathering conditions in the tropical context of the WEAR that may increase the sensitivity to landsliding [50]; and (3) the underestimation of the area-averaged SRE [33,56] used in the calculation of AR (Equation (1)). The latter is not necessarily an issue when thresholds are evaluated with the same SRE used for their calibration [32]. To date, SRE-based studies form only a small fraction in landslide threshold research [3,32]. A TMPA-RT-based threshold was established for Italy at the 20% exceedance probability, obtaining a critical rainfall accumulation of 189 mm for an extrapolated duration of 42 days [32]. However, the extrapolation is doubtful, given the limited threshold calibration range of ~11 days [32]. In addition, no decay function is applied to this extrapolated value for accumulated antecedent rainfall and the higher exceedance probability level obviously renders an elevated threshold compared to the 5% and 10% levels deployed in our study. On the other hand, the 5% thresholds calibrated for central Italy by [71] based on TMPA (Research Version) data, estimate critical accumulated rainfall in the order of 30 mm over an extrapolated duration of 42 days, which falls in our estimated triggering range at the 5% exceedance probability.

Furthermore, the modified AR-S threshold method is relevant in the context of landslide hazard analysis when rainfall data are only available for conditions that triggered landslides, for it proved to be a robust alternative for frequentist-based threshold approaches [9,40] when the method’s assumptions (i.e., large and well-spread data set [40]) are not met. To the authors’ knowledge, it is the first time that a stratified data selection technique is adopted in the threshold calibration approach, which effectively showed to enhance the data distribution over the whole range of the causative threshold variable (S). This allows the method to be transferable not only to other data sets for S (and hence other study areas), but to any parameter that might be considered as a possible cause for landsliding [34] without the requirement of the data to be homogeneously distributed, to which further research should be carried out. This is significant in the context of the ‘trigger-cause’ conceptual framework of threshold definition as proposed by [34], in response to the shortcomings of the rainfall-only thresholds with regard to their limited physical meaning [34]. The framework was designed to introduce hydrological information on different timescales, with the choice of the parameters and timescales depending on their expected significance for slope failure given the physiographic context and considered landslide types [34]. In the AR-S approach, the causative hydrological status of the slope is substituted by information on spatially varying predisposing ground conditions, while AR presents the progressive build-up of the landslide trigger [35]. We could anticipate, however, that improvements in satellite-based soil moisture data, with regard to their spatial resolution and performance over dense vegetated areas or complex topography [72,73,74], would allow in the future to replace the static S variable by a dynamic causative hydrological factor over data-scarce regions.

In addition, the enhanced relation between AR threshold values and S renders a potential landslide early warning system more efficient. In this respect, the integrated spatial component of the AR-S approach (Figure 8) has a clear advantage over traditional thresholds, such as rainfall intensity-duration thresholds, the latter only informing ‘when’ the probability of a landslide occurrence increases but not ‘where’ [65]. The spatial component furthermore avoids data to be partitioned according to homogeneous physiographic units to enhance the accuracy of thresholds [9,75], which is of particular relevance in data-scarce contexts. Finally, the method was designed to use area-averaged SRE, allowing its adoption in regions where a dense rain gauge network is absent and evaluate hazard in near-real time.

However, the efficient use of the AR-S approach in hazard or early warning studies is hampered by a number of factors that might contribute to an obscured relation between AR and S, including: (1) the level of accuracy and completeness of the landslide inventory [28,42]; (2) a lack of differentiation in landslide processes whose triggering conditions are different [63]; (3) the accuracy of S and SRE data [33,49]; (4) the resampling of S data to the coarser SRE data resulting in inaccurate S-classifications of landslide data; (5) the anthropogenic influence on the environment, such as in the context of exponential demographic pressure in the WEAR [76,77]; and (6) the empirically defined parameters in the AR equation (Equation (1)) by [35]. The highlighted obstacles and limitations serve as pathways for further investigation and improvements in the AR-S threshold approach.

8. Conclusions

We propose a modified antecedent rainfall–susceptibility (AR-S) threshold approach that improves on the initial AR-S method of [35], being transferable to other data sets for landsliding and S. For its development and evaluation we exploit the most current and extensive landslide inventory for the western branch of the East African Rift comprising 184 dated landslide events from 2001 to 2019, satellite-based rainfall estimates from TMPA 3B42 RT, and two S models, i.e., the continental-scale model of [45] and the regional-scale S model of [49]. The main novelty in the modified AR-S approach is the stratified selection of data associated with the lowest AR values able to cause landsliding, allowing to deploy data sets that are not necessarily homogenously distributed over the S range. Furthermore, we highlight that the threshold procedure is more meaningful when no bootstrapping statistical technique is applied, as the uncertainties in the parameters that define the power-law threshold model are mainly introduced by the bootstrapping related random sampling in combination with the presence of outliers in the data set. We obtain improved AR thresholds with an increased susceptibility-dependent gradient, and AR threshold maps with a higher accuracy through the use of the regional-scale S model in the modified AR-S approach. The improved AR threshold values at the 5% exceedance probability range from 7 mm in the most susceptible areas (S = 0.72) to 62 mm in the lowest susceptible areas (S = 0.10) where landslides have been recorded (uncorrected for underestimation by TMPA). Our approach is foremost relevant in data-scarce regions, where the lack of abundant data from rain gauges and in particular on landslide occurrence hampers the use of homogenously distributed data sets. Moreover, we suggest that this modified method is transferable not only to other data sets for S, but to any parameter that might be considered as a possible cause for landsliding.