A Landslide Probability Model Based on a Long-Term Landslide Inventory and Rainfall Factors

: The prediction and advanced warning of landslide hazards in large-scale areas must deal with a large amount of uncertainty, therefore a growing number of studies are using stochastic models to analyze the probability of landslide occurrences. In this study, we used a modiﬁed Thiessen’s polygon method to divide the research area into several rain gauge control areas, and divided the control areas into slope units reﬂecting the topographic characteristics to enhance the spatial resolution of a landslide probability model. We used a 2000–2015 long-term landslide inventory, daily rainfall, and e ﬀ ective accumulated rainfall to estimate the rainfall threshold that can trigger landslides. We then employed a Poisson probability model and historical rainfall data from 1987 to 2016 to calculate the exceedance probability that rainfall events will exceed the threshold value. We calculated the number of landslides occurring from the events when rainfall exceeds the threshold value in the slope units to estimate the probability that a landslide will occur in this situation. Lastly, we employed the concept of conditional probability by multiplying this probability with the exceedance probability of rainfall events exceeding the threshold value, which yielded the probability that a landslide will occur in each slope unit for one year. The results indicated the slope units with high probability that at least one rainfall event will exceed the threshold value at the same time that one landslide will occur within any one year are largely located in the southwestern part of the Taipei Water Source Domain, and the highest probability is 0.26. These slope units are located in parts of the study area with relatively weak lithology, high elevations, and steep slopes. Compared with probability models based solely on landslide inventories, our proposed landslide probability model, combined with a long-term landslide inventory and rainfall factors, can avoid problems resulting from an incomplete landslide inventory, and can also be used to estimate landslide occurrence probability based on future potential changes in rainfall.


Introduction
Taiwan is a relatively new island formed by plate movements. Due to its high mountains, steep slopes, and relatively unstable geological conditions, as well as frequent typhoons and torrential rains, slopeland disasters are common in mountainous areas. Thus, slopeland hazard prevention and mitigation projects are necessary. In slopeland hazard prevention work, landslides have a high level of unpredictability. In particular, estimating the likelihood of landslides in large watersheds using deterministic models is difficult when no detailed geomorphological and hydrological data have been collected for the whole area. Therefore, the use of a stochastic model to assess landslide probability is more feasible. According to the definition, landslide hazard involves both spatial and temporal probability [1]. The analysis of landslide spatial probability is generally seen as a landslide landslides [42,43], and antecedent rainfall can affect both of these factors. Accordingly, antecedent rainfall can be used to determine when landslides may occur [36]. In research on rainfall thresholds incorporating the antecedent rainfall conditions, large differences exist in the number of days of antecedent rainfall that were employed in each study. For example, daily rainfall was employed in conjunction with 15-day antecedent rainfall [39], both daily and 3-day cumulative rainfall were used [44], and 3-day and 30-day antecedent rainfall were employed [45]. Guzzetti et al. [36] suggested that the large variability in the antecedent rainfall may be attributed to three types of factors concerning the research area: diversity in lithological, morphological, vegetation, and soil conditions; differences in climatic regimes and meteorological circumstances; and the heterogeneity and incompleteness in the rainfall and landslide data used to establish the rainfall thresholds. As a consequence, the local conditions and availability of data in the research area must be assessed when choosing the number of days of antecedent rainfall.
Since the rainfall threshold determined using a single rain gauge for a large area constitutes one value for the entire area, as soon as rainfall reaches or exceeds the threshold, landslides may occur anywhere in that area, and knowing their precise locations is impossible. As a consequence, a denser array of rain gauges can be employed to acquire rainfall data with finer spatial resolution [21], and the research areas can be subdivided into analytical units with a smaller area, which can better account for the spatial variability of rainfall patterns in the analytical units and the spatial resolution of landslide prediction. However, 19.1% of recent studies on this subject failed to subdivide their research areas, and those studies that did subdivide their research areas had resulting analytical units with an average area of 302.0 km 2 [21]. For instance, a research area of 4660 km 2 was subdivided into 12 analytical units with an average area of 388.3 km 2 [39], but excessively large analytical units make it impossible to identify the precise possible locations of landslides. In addition, the subdivision approaches employed in some studies run into the problem of incomplete coverage. For instance, although a 25 km 2 research area was subdivided into eight analytical units, the landslide prediction results only represented the paths of roads in the subdivisions and not the entire subdivisions because most landslides (94%) in the study occurred on roadside slopes [26]. Althuwaynee et al. [28] divided the research area into six circular analytical units with their centers at rain gauges, but the analytical units did not cover the entire research area and also overlapped. Although these studies subdivided their research areas into different analytical units, the units could not provide a landslide probability distribution with a finer spatial resolution because they were excessively large, or experienced problems such as incomplete coverage and overlap. If the method of subdividing a research area into analytical units is improved so that the units are smaller in area, the spatial resolution of the landslide probability estimation results could be improved. There are seven types of analytical units subdivided in research areas: grid cell, terrain unit, unique condition unit, slope unit, geo-hydrological unit, topographical unit, and administrative unit [46,47]. The slope units are suitable for landslide probability analysis because they express topographic features and slope characteristics. In this study, we consequently selected slope units as our analytical unit.

Environmental Setting of Taipei Water Source Domain
Taipei Water Source Domain is located in the northeast part of Taiwan and supplies tap water for five million people in the greater Taipei area. The area is characterized by hilly and mountainous topography, as well as the Xueshan Range extending to the northeast and a subrange of Mt. Qilan extending to the northwest, both of which account for the area's high terrain in the south and low terrain in the north. Elevations in the area range from 12 to 2130 m ( Figure 1). Water 2020, 12, x FOR PEER REVIEW 4 of 17 Concerning the distribution of lithology, we followed the classification approach proposed by Lin et al. [48] by dividing the Taipei Water Source Domain into areas underlain by alluvium, loose sandstone and shale, hard sandstone and shale, and slate. Whereas alluvium found at the confluence of rivers and in downstream areas covers only a small part of the research area, hard sandstone and shale as well as slate underlay most of the research area. Of these types of lithology, areas underlain by slate had the highest number of landslides and the greatest landslide area. Wu et al. [49] indicated that the areas underlain by hard sandstone and shale as well as slate in the Kaoping River Watershed had the highest landslide ratios in 2008 and 2009. This indicates that the lithology condition of most areas is fragile. Typhoons and torrential rain events can readily wash away unconsolidated sand and gravel as well as trigger landslides, which deposit large loads of sediment in rivers and reservoirs.

Rainfall Data
The rain gauges employed in this study were located as shown in Figure 1, and rainfall data between 1987 to 2016 from these rain gauges were used. Average daily rainfall for the entire area during the same period as the 2000-2015 landslide inventory is shown in Figure 2. Figure 2 shows that apart from the eight typhoon events causing the corresponding landslide inventory, other events of high daily rainfall occurred without a significant increase in landslides. As a consequence, apart from calculating the exceedance probability that rainfall events will exceed the rainfall threshold, we also calculated the probability of landslides when rainfall events exceed the threshold. In addition, Figure 3 shows the average daily rainfall and standard deviation of the eight typhoon events during the 2000-2015 period in each control area of a rain gauge divided by a modified Thiessen polygon method, considering the morphology of the area, proposed by Salvaticic et al. [19]. Concerning the distribution of lithology, we followed the classification approach proposed by Lin et al. [48] by dividing the Taipei Water Source Domain into areas underlain by alluvium, loose sandstone and shale, hard sandstone and shale, and slate. Whereas alluvium found at the confluence of rivers and in downstream areas covers only a small part of the research area, hard sandstone and shale as well as slate underlay most of the research area. Of these types of lithology, areas underlain by slate had the highest number of landslides and the greatest landslide area. Wu et al. [49] indicated that the areas underlain by hard sandstone and shale as well as slate in the Kaoping River Watershed had the highest landslide ratios in 2008 and 2009. This indicates that the lithology condition of most areas is fragile. Typhoons and torrential rain events can readily wash away unconsolidated sand and gravel as well as trigger landslides, which deposit large loads of sediment in rivers and reservoirs.

Rainfall Data
The rain gauges employed in this study were located as shown in Figure 1, and rainfall data between 1987 to 2016 from these rain gauges were used. Average daily rainfall for the entire area during the same period as the 2000-2015 landslide inventory is shown in Figure 2. Figure 2 shows that apart from the eight typhoon events causing the corresponding landslide inventory, other events of high daily rainfall occurred without a significant increase in landslides. As a consequence, apart from calculating the exceedance probability that rainfall events will exceed the rainfall threshold, we also calculated the probability of landslides when rainfall events exceed the threshold. In addition, Figure 3 shows the average daily rainfall and standard deviation of the eight typhoon events during the 2000-2015 period in each control area of a rain gauge divided by a modified Thiessen polygon method, considering the morphology of the area, proposed by Salvaticic et al. [19]. from the satellite images, we often found that human mapping errors affected interpretation quality. We followed the recommended procedures proposed by Liu et al. [50] to map landslides in the research area. Table 1 shows the dates of the eight landslide events and landslide statistical data. The size of landslides ranged from 16 to 118,108 m 2 and the average area was 2474 m 2 . The resulting distribution of landslides caused by the eight typhoon events was shown in Figure 4, which revealed that landslide sites were concentrated in the southwestern portion of the research area.      area of landslides [47], which minimized the chance that any specific landslide site would be a part of different slope units, and thereby confuse the analysis results. We also divided the research area into rain gauge control areas ( Figure 3) based on rain gauge locations and using the modified Thiessen polygon method. The rainfall measured by each rain gauge was taken as representative of the control area in which that gauge was located, and we expected this approach to reflect the different rainfall distribution characteristics within the research area.

Analytical Units and Rain Gauge Control Areas
We employed slope units as analytical units due to their relatively well-defined topographic boundaries, as well as topographic and geological meaning. We employed the subdivision method used by Xie et al. [51] to divide the watershed into slope units. The original topography could be divided into sub-watersheds, and the combination of sub-watershed units before and after reversal yielded the slope units. We ensured that the smallest area of slope units was larger than the average area of landslides [47], which minimized the chance that any specific landslide site would be a part of different slope units, and thereby confuse the analysis results. We also divided the research area into rain gauge control areas ( Figure 3) based on rain gauge locations and using the modified Thiessen polygon method. The rainfall measured by each rain gauge was taken as representative of the control area in which that gauge was located, and we expected this approach to reflect the different rainfall distribution characteristics within the research area.

Analysis of Discrete Rainfall Groups
The two rainfall parameters considered in this study consisted of daily rainfall (I) and effective accumulated rainfall (R t ). After selecting rain gauges near the research area with rainfall data for recent years, we obtained daily rainfall data for the 1987-2016 period from the Water Resources Agency and Central Weather Bureau. This study calculated the effective accumulated rainfall based on rainfall for that day and rainfall during the previous 7 days using the method proposed by Jan [52]; this calculation was performed using Equation (1): where R 0 is the rainfall amount on that day, R 1 is the rainfall amount on the day before that day, and so on, and the weighting coefficient α = 0.7 proposed by Jan [52]. Adopting the concept proposed by Tsai [53], after using daily rainfall data to calculate effective accumulated rainfall (R t ), we obtained a group of daily rainfall and effective accumulated rainfall (I, R t ) for each day. The daily rainfall and effective accumulated rainfall were continuous variables and would not facilitate subsequent calculation of a joint cumulative distribution function, therefore we rounded off the daily rainfall and effective accumulated rainfall values to the 10th place and made them discrete variables. The group of daily rainfall and effective accumulated rainfall (I, R t ) for each day was termed as "discrete rainfall group" in this study.
We defined different rainfall events by the continuity of daily rainfall. Consecutive days of non-zero daily rainfall were considered to be the same rainfall event, and the number of the consecutive days varied from event to event. We then calculated the distance (d) from each discrete rainfall group to the origin (0, 0), and assumed that the greater the value of d, the greater the likelihood of landslides. The discrete rainfall group with the greatest d in each rainfall event was used to represent that rainfall event in subsequent analysis.

Joint Cumulative Distribution Function
The joint cumulative distribution function was obtained from the joint probability mass function of the foregoing discrete rainfall groups. The probability (P I , Rt (I i , R tj )) of each discrete rainfall group (I i , R tj ) was defined [54] as shown in Equation (2): where i = 0, 10, 20, 30, . . . ; j = 0, 10, 20, 30, . . . ; the joint probability mass function has a probability value only when I and R t are multiples of 10 and the probability values in other places are 0. The foregoing joint probability mass function yielded a joint cumulative distribution function using: The joint cumulative distribution function was a monotonic increasing function with a range between 0 and 1, and had the form of a three-dimensional curved surface when plotted on coordinate axes. The farther the point (I i , R tj ) from the origin, the greater its probability value. The probability of a discrete rainfall group on the curved surface expressed the cumulative probability of all discrete rainfall groups, which were nearer to the origin than this discrete rainfall group (I i , R tj ).

Selection of a Rainfall Probability Threshold
After establishing a joint cumulative distribution function, taking each 0.05 as an interval, we set 20 rainfall probability thresholds ranging from 0.05 to 1.00, and employed the error matrix concept to calculate the true positive rate (TPR), true negative rate (TNR), and positive predictive value (PPV) for each rainfall probability threshold at each rain gauge control area. The rainfall probability threshold was treated as the threshold of cumulative probability of the discrete rainfall groups which was used to predict whether rainfall events could trigger landslides. Here, TPR expresses the ratio of discrete rainfall groups that correctly predicted landslide occurrence to discrete rainfall groups triggering landslides actually, TNR expresses the ratio of discrete rainfall groups that correctly predicted no landslide occurrence to discrete rainfall groups triggering no landslides actually, and PPV expresses the ratio of discrete rainfall groups that correctly predicted landslide occurrence to discrete rainfall groups predicting landslides. To capture the performance of each threshold, PPV and Youden's index were used for comprehensive consideration. The higher the PPV and Youden's index values, the more accurate the rainfall probability threshold at classifying landslide occurrence and landslide occurrence for discrete rainfall groups. The TPR, TNR, PPV, and Youden's index calculations were performed employing Equations (4) Youden s index = TPR + TNR − 1 (7)

Poisson Probability Model
A Poisson probability model relies on the past frequency of events to predict their occurrence probability in the future. The basic assumption underlying this type of model is that future events will occur with the same frequency as past events. In this model, the probability of at least one event occurring in the time interval (t) is given by Equation (8): where P(N(t) ≥ 1) indicates the probability of at least one event occurring within a period of t years; this probability is known as the exceedance probability.
We calculated the number of discrete rainfall groups exceeding the threshold at each rain gauge in the past using the optimal rainfall probability thresholds and then divided by the years of the rainfall data to obtain the occurrence frequency (λ), which was used to calculate the exceedance probability. The exceedance probability indicated the probability of at least one rainfall event exceeding the threshold of discrete rainfall groups within any one year.

Conditional Probability
We employed the concept of conditional probability in the analysis. We first used the Poisson probability model to calculate the exceedance probability of at least one rainfall event exceeding the threshold of discrete rainfall groups within any one year at each rain gauge control area. We divided the number of landslides occurring in each slope unit by the number of rainfall events exceeding the threshold of discrete rainfall groups to estimate the probability that a landslide would occur in that slope unit when the rainfall exceeded the threshold. Lastly, we multiplied the two probabilities together to obtain the probability that a rainfall event would exceed the threshold of discrete rainfall groups and at least one landslide would also occur in each slope unit within any one year, as shown in Equation (9): where R ≥ RT indicates rainfall events exceed the threshold of the discrete rainfall group and L indicates the occurrence of a landslide.

Joint Cumulative Distribution Functions of the Rain Gauges
In this study, we collected multi-year daily rainfall data from each rain gauge and calculated the effective accumulated rainfall (R t ) by employing Equation (1), which yielded rainfall and effective accumulated rainfall for each day. We then rounded off the daily rainfall and effective accumulated rainfall values to the 10th place, which yielded discrete rainfall groups including both daily rainfall and effective accumulated rainfall. The next step was establishing frequency tables for different discrete rainfall groups, which we used to show the frequency of the discrete rainfall groups. Figure 5 shows the frequency of discrete rainfall groups at the Bihu rain gauge with daily rainfall and effective accumulated rainfall (R t ) ranging from 0 to 100 mm. The depth axis represents daily rainfall, the horizontal axis represents the effective accumulated rainfall (R t ), and the vertical axis represents the frequency of a discrete rainfall group. We then calculated the cumulative frequency of each discrete rainfall group on this basis, and this represented the frequency of all discrete rainfall groups with values lower than that of any designated discrete rainfall group. The cumulative frequency was then divided by the total frequency of all discrete rainfall groups, which yielded the cumulative probability of each discrete rainfall group. rainfall values to the 10th place, which yielded discrete rainfall groups including both daily rainfall and effective accumulated rainfall. The next step was establishing frequency tables for different discrete rainfall groups, which we used to show the frequency of the discrete rainfall groups. Figure  5 shows the frequency of discrete rainfall groups at the Bihu rain gauge with daily rainfall and effective accumulated rainfall (Rt) ranging from 0 to 100 mm. The depth axis represents daily rainfall, the horizontal axis represents the effective accumulated rainfall (Rt), and the vertical axis represents the frequency of a discrete rainfall group. We then calculated the cumulative frequency of each discrete rainfall group on this basis, and this represented the frequency of all discrete rainfall groups with values lower than that of any designated discrete rainfall group. The cumulative frequency was then divided by the total frequency of all discrete rainfall groups, which yielded the cumulative probability of each discrete rainfall group. The joint cumulative distribution function of each rain gauge was then obtained from the cumulative probability of the discrete rainfall groups, and this function was used to plot a joint cumulative distribution chart. Figures 6 and 7 are joint cumulative distribution functions for the Bihu and Fushan (3) rain gauges, and daily rainfall and effective accumulated rainfall (Rt) are shown within a 0-300 mm range. The joint cumulative distribution functions have areas with gentler slopes indicating fewer and more dispersed discrete rainfall groups within a certain interval, and areas with steeper slopes indicating more and more concentrated discrete rainfall groups within a certain interval.  The joint cumulative distribution function of each rain gauge was then obtained from the cumulative probability of the discrete rainfall groups, and this function was used to plot a joint cumulative distribution chart. Figures 6 and 7 are joint cumulative distribution functions for the Bihu and Fushan (3) rain gauges, and daily rainfall and effective accumulated rainfall (R t ) are shown within a 0-300 mm range. The joint cumulative distribution functions have areas with gentler slopes indicating fewer and more dispersed discrete rainfall groups within a certain interval, and areas with steeper slopes indicating more and more concentrated discrete rainfall groups within a certain interval. and effective accumulated rainfall. The next step was establishing frequency tables for different discrete rainfall groups, which we used to show the frequency of the discrete rainfall groups. Figure  5 shows the frequency of discrete rainfall groups at the Bihu rain gauge with daily rainfall and effective accumulated rainfall (Rt) ranging from 0 to 100 mm. The depth axis represents daily rainfall, the horizontal axis represents the effective accumulated rainfall (Rt), and the vertical axis represents the frequency of a discrete rainfall group. We then calculated the cumulative frequency of each discrete rainfall group on this basis, and this represented the frequency of all discrete rainfall groups with values lower than that of any designated discrete rainfall group. The cumulative frequency was then divided by the total frequency of all discrete rainfall groups, which yielded the cumulative probability of each discrete rainfall group. The joint cumulative distribution function of each rain gauge was then obtained from the cumulative probability of the discrete rainfall groups, and this function was used to plot a joint cumulative distribution chart. Figures 6 and 7 are joint cumulative distribution functions for the Bihu and Fushan (3) rain gauges, and daily rainfall and effective accumulated rainfall (Rt) are shown within a 0-300 mm range. The joint cumulative distribution functions have areas with gentler slopes indicating fewer and more dispersed discrete rainfall groups within a certain interval, and areas with steeper slopes indicating more and more concentrated discrete rainfall groups within a certain interval.

Selection of Rainfall Probability Thresholds of the Rain Gauges
Following the analysis results of the joint cumulative distribution functions of the rain gauges, we used rainfall data from the rain gauges during the eight rainfall events triggering landslides to select rainfall probability thresholds. The rainfall probability threshold was treated as the threshold

Selection of Rainfall Probability Thresholds of the Rain Gauges
Following the analysis results of the joint cumulative distribution functions of the rain gauges, we used rainfall data from the rain gauges during the eight rainfall events triggering landslides to select rainfall probability thresholds. The rainfall probability threshold was treated as the threshold of cumulative probability of the discrete rainfall groups which was used to predict whether rainfall events could trigger landslides. Starting with a rainfall probability threshold value of 0.05, we set a rainfall probability threshold at each interval of 0.05 until a value of 1.00 was reached, and then calculated the TPR, TNR, PPV, and Youden's index of each rainfall probability threshold. Here, the number of landslide events predicted correctly divided by the number of rainfall events triggering landslides actually equaled TPR, the number of no landslide events predicted correctly divided by the number of rainfall events triggering no landslides actually equaled TNR, and the number of landslide events predicted correctly divided by the number of rainfall events predicting landslides equaled PPV. Table 2 shows the results of these calculations for the Bihu rain gauge. In the analysis results for the individual rain gauges, the rainfall probability thresholds with the highest Youden's index were within the probability interval of 0.85-0.95, and the rainfall probability thresholds with the highest PPV were at the probability of 0.95 in all cases. We consequently opted to use a rainfall probability threshold value of 0.95 for the whole area. The TPR, TNR, PPV, and Youden's index for all rain gauges when the rainfall probability threshold was 0.95 are shown in Table 3.  Table 3. TPR, TNR, PPV, and Youden's index for all rain gauges at a rainfall probability threshold of 0.95.

Landslide Probability Analysis Employing a Rainfall Probability Threshold and a Long-Term Landslide Inventory
After determining a rainfall probability threshold for the rain gauges, we calculated the number of discrete rainfall groups exceeding this threshold at each rain gauge during the 1987-2016 period. We then divided these values by the years of statistics at each gauge, which yielded the λ values in Equation (8). Substituting t = 1 year into Equation (8) allowed us to calculate the probability of at least one rainfall event exceeding the threshold of discrete rainfall group within any one year (i.e., P(R ≥ RT) in Equation (9)) under the assumption that future rainfall conditions will be the same as past conditions. Figure 8 shows the exceedance probability value calculated for each rain gauge overlaid on each rain gauge control area. The Quchi rain gauge control area had the highest probability of 0.76502 that at least one rainfall event will exceed the threshold of the discrete rainfall group within any one year, whereas the Bihu rain gauge control area had the lowest probability of 0.43886.
Water 2020, 12, x FOR PEER REVIEW 12 of 17 overlaid on each rain gauge control area. The Quchi rain gauge control area had the highest probability of 0.76502 that at least one rainfall event will exceed the threshold of the discrete rainfall group within any one year, whereas the Bihu rain gauge control area had the lowest probability of 0.43886. Figure 8. The exceedance probability that at least one rainfall event will exceed the threshold of discrete rainfall group within any one year in each rain gauge control area.
In this study, we also divided the number of landslides occurring in each slope unit during the 2000-2015 period by the number of rainfall events exceeding the threshold of discrete rainfall group at the rain gauges to which the slope units were assigned during the same period to estimate the landslide probability in the slope units when the rainfall exceeded the threshold, which is P(L│R ≥ RT) in Equation (9). The resulting probability distribution is shown in Figure 9. Figure 9 shows that the different slope units within a single rain gauge control area have different landslide probabilities, and these differences should be attributed to different geomorphological conditions in the slope units.
Lastly, employing Equation (9), we multiplied the probability P(R ≥ RT) that at least one rainfall event will exceed the threshold of discrete rainfall group within any one year in each rain gauge control area by the landslide probability P(L│R ≥ RT) in each slope unit when rainfall exceeds the threshold, which yielded the probability that at least one rainfall event exceeds the threshold of discrete rainfall group at the same time that one landslide will occur in each slope unit during the future one-year period ( Figure 10). The two probability maps shown in Figures 9 and 10 were validated by the landslide inventory data respectively. The landslides were mainly distributed in the slope units where the landslide probability values were greater than 0.01. The top 2% of slope units ranked with landslide probabilities included 50.40% of slope units where landslides occurred while the top 6% of slope units ranked with landslide probabilities included 100.00% of slope units where landslides occurred in Figure 10. The results indicated these maps had reasonable landslide probability distributions. Figure 10 reveals that the Fushan (3) rain gauge control area, which is located in the southwest part of the research area, contained relatively many slope units with high landslide probability, and the highest probability value was 0.26. Apart from having fragile lithology consisting of hard sandstone and shale as well as slate, this area has a higher elevation and steeper slopes than other control areas, which suggests that elevation and slope have a definite correlation with landslide occurrence. Figure 8. The exceedance probability that at least one rainfall event will exceed the threshold of discrete rainfall group within any one year in each rain gauge control area.
In this study, we also divided the number of landslides occurring in each slope unit during the 2000-2015 period by the number of rainfall events exceeding the threshold of discrete rainfall group at the rain gauges to which the slope units were assigned during the same period to estimate the landslide probability in the slope units when the rainfall exceeded the threshold, which is P(L|R ≥ RT) in Equation (9). The resulting probability distribution is shown in Figure 9. Figure 9 shows that the different slope units within a single rain gauge control area have different landslide probabilities, and these differences should be attributed to different geomorphological conditions in the slope units.
Lastly, employing Equation (9), we multiplied the probability P(R ≥ RT) that at least one rainfall event will exceed the threshold of discrete rainfall group within any one year in each rain gauge control area by the landslide probability P(L|R ≥ RT) in each slope unit when rainfall exceeds the threshold, which yielded the probability that at least one rainfall event exceeds the threshold of discrete rainfall group at the same time that one landslide will occur in each slope unit during the future one-year period ( Figure 10). The two probability maps shown in Figures 9 and 10 were validated by the landslide inventory data respectively. The landslides were mainly distributed in the slope units where the landslide probability values were greater than 0.01. The top 2% of slope units ranked with landslide probabilities included 50.40% of slope units where landslides occurred while the top 6% of slope units ranked with landslide probabilities included 100.00% of slope units where landslides occurred in Figure 10. The results indicated these maps had reasonable landslide probability distributions. Figure 10 reveals that the Fushan (3) rain gauge control area, which is located in the southwest part of the research area, contained relatively many slope units with high landslide probability, and the highest probability value was 0.26. Apart from having fragile lithology consisting of hard sandstone and shale as well as slate, this area has a higher elevation and steeper slopes than other control areas, which suggests that elevation and slope have a definite correlation with landslide occurrence. Water 2020, 12, x FOR PEER REVIEW 13 of 17  The probability that at least one rainfall event exceeds the threshold of discrete rainfall group and one landslide will also occur during the future one-year period within the Taipei Water Source Domain.

Discussion
In comparison with a landslide probability model based solely on the use of landslide inventories, our landslide probability model based on the use of landslide inventories and rainfall factors reflect different basic assumed conditions. The assumption of the landslide probability model incorporating rainfall factors is that the frequency of future rainfall events exceeding the threshold and the frequency of landslides occurring when the threshold has been exceeded are the same as in the past. In contrast, the assumption of a landslide probability model based solely on the  The probability that at least one rainfall event exceeds the threshold of discrete rainfall group and one landslide will also occur during the future one-year period within the Taipei Water Source Domain.

Discussion
In comparison with a landslide probability model based solely on the use of landslide inventories, our landslide probability model based on the use of landslide inventories and rainfall factors reflect different basic assumed conditions. The assumption of the landslide probability model incorporating rainfall factors is that the frequency of future rainfall events exceeding the threshold and the frequency of landslides occurring when the threshold has been exceeded are the same as in the past. In contrast, the assumption of a landslide probability model based solely on the Figure 10. The probability that at least one rainfall event exceeds the threshold of discrete rainfall group and one landslide will also occur during the future one-year period within the Taipei Water Source Domain.

Discussion
In comparison with a landslide probability model based solely on the use of landslide inventories, our landslide probability model based on the use of landslide inventories and rainfall factors reflect different basic assumed conditions. The assumption of the landslide probability model incorporating rainfall factors is that the frequency of future rainfall events exceeding the threshold and the frequency of landslides occurring when the threshold has been exceeded are the same as in the past. In contrast, the assumption of a landslide probability model based solely on the use of landslide inventories is that the frequency of future landslides occurring is the same as in the past. As a consequence, landslide probability models incorporating rainfall factors possess the following advantages: (1) This model can reflect the differences in landslide probability between the rain gauge control areas that have different rainfall conditions. (2) When rainfall data were added in the analysis, the probability model we obtained yielded more reliable results because the rainfall data were collected from a longer period (29-45 years) than the landslide inventory (16 years). (3) If we know how the probability of at least one rainfall event exceeding the threshold will change in the future, the incorporation of rainfall factors in the landslide probability model will allow the effect of possible rainfall changes on the landslide probability to be assessed.
However, several aspects connected to the application of this landslide probability model still require further investigation: (1) The method of analyzing landslide probability proposed in this study requires a long-term landslide inventory and rainfall data, therefore attention must be paid to the completeness of rainfall data for the research area and handling methods when data are incomplete. (2) Whereas the rainfall factors used in this study reflect daily rainfall and effective accumulated rainfall, the use of different rainfall factors will yield different analysis results, which may be explored further in future research. (3) Apart from the modified Thiessen polygon method, the division of rain gauge control areas can be performed using other methods, such as the height-balance polygon method. A finer division method should yield more precise results of a landslide probability distribution, therefore future research can also compare the applicability of different methods of division into rain gauge control areas. (4) We obtained a long-term landslide inventory consisting of only eight events, therefore all events collected were used in the process of building the model. The landslide inventory covering the period of other events may be collected to verify the predictive ability of this landslide probability model.

Conclusions
In this study, we employed joint cumulative distribution functions to calculate the TPR, TNR, PPV, and Youden's index for different rainfall probability thresholds, selected a threshold of 0.95 as suitable for the research area, and used this rainfall probability threshold to calculate the Poisson probability of at least one rainfall event exceeding the threshold of discrete rainfall groups at each rain gauge within the future one-year period. We then combined this probability with the landslide probability in individual slope units when rainfall exceeded the threshold value, which allowed us to estimate the probability that a landslide will occur in individual slope units during the future one-year period. Many of the slope units with a high landslide probability are located in the Fushan (3) rain gauge control area, and the highest probability is 0.26. Apart from fragile lithology, this area is characterized by high elevations and steep slopes, which indicates that the elevation and slope have a significant influence on the occurrence of landslides. This finding suggests that this area should be a focal area for landslide prevention and mitigation efforts. The landslide probability model established based on the use of a long-term landslide inventory and rainfall factor had a finer spatial resolution and data for a longer period, which yielded more reliable results and enabled the effect of possible rainfall changes on the landslide probability to be assessed. The effects of the completeness of rainfall data for the research area, the use of different rainfall factors, as well as the different methods of division into rain gauge control areas on the landslide probability analysis results can be significant and still require further investigation.

Conflicts of Interest:
The authors declare no conflict of interest.