Evaluation of Rainfall-Triggered Debris Flows under the Impact of Extreme Events: A Chenyulan Watershed Case Study, Taiwan

Abstract

This study examined the conditions that lead to debris flows, and their association with the rainfall return period (T) and the probability of debris flow occurrence (P) in the Chenyulan watershed, central Taiwan. Several extreme events have occurred in the Chenyulan watershed in the past, including the Chi-Chi earthquake and extreme rainfall events. The T for three rainfall indexes (i.e., the maximum hourly rainfall depth (I_m), the maximum 24-h rainfall amount (R_d), and RI (RI = I_m $\times$ R_d)) were analyzed, and the T associated with the triggering of debris flows is presented. The P–T relationship can be determined using three indexes, I_m, R_d, and RI; how it is affected and unaffected by extreme events was developed. Models for evaluating P using the three rainfall indexes were proposed and used to evaluate P between 2009 and 2020 (i.e., after the extreme rainfall event of Typhoon Morakot in 2009). The results of this study showed that the P‒T relationship, using the RI or R_d index, was reasonable for predicting the probability of debris flow occurrence.

1. Introduction

In recent years, the frequency and magnitude of disasters associated with climate extremes have increased [1,2,3,4,5,6]. Extreme rainfall has given rise to massive sediment transport in mountainous areas; in particular, debris flows have caused significant damage in many countries [7,8,9,10,11,12,13,14].

Studies of rainfall that triggers debris flows, for the purpose of issuing hazard warnings, have focused on empirical thresholds based on rain gauge observations [15]. The properties of cumulative rainfall depth, rainfall intensity, and rainfall duration have been widely used to develop an empirical threshold for triggering of debris flows [16,17,18,19,20,21]. However, debris flow initiation depends not only on local rainfall properties, but also on sediment availability in the catchment area and local terrain conditions, such as topography, lithology, and soil cover [15,18,22]. The sediment supply can change significantly after extreme events, such as a major earthquake and extreme rainfall, and the triggering threshold or critical rainfall for debris flow may change accordingly [23,24,25,26,27,28]. This leads to several uncertainties when determining rainfall thresholds, and makes debris flow monitoring and prediction a challenge [29,30]. However, most of the existing empirical models were developed for real-time warning or monitoring of debris flows. In the face of climate change and increasing extreme events, there are few studies on how to reflect the long-term variation in rainfall characteristics of a region (composed of multiple debris flow gullies, townships, and villages) and the impact of extreme rainfall events on debris flow monitoring. After extreme rainfall events, the rainfall conditions that can trigger a debris flow will change, and this unstable period lasts for several years [26,28]. Thus, it is necessary to establish an empirical model for critical rainfall which can be combined with the effects of long-term rainfall changes and extreme events (on an annual scale) for debris flow monitoring and prediction. There are many rainfall indexes, such as cumulative rainfall depth and rainfall intensity, which researchers have used to develop empirical models of critical rainfall. This study aims to integrate these rainfall indicators and establish a model for evaluation of the probability of debris flow occurrence.

The Chenyulan watershed, located in central Taiwan, has experienced extreme events including the Chi-Chi earthquake (CCE) and extreme rainfall, and was therefore selected as the study area. The conditions that lead to debris flows, and their association with the rainfall return period (T) and the probability of debris flow occurrence (P), were examined. This study had two main purposes: (1) to develop an empirical model for how the P‒T relationship was affected by extreme rainfall and the CCE, using different rainfall indexes, and (2) to apply the P–T relationship to evaluate P after the recent extreme event of Typhoon Morakot.

2. Study Area and Extreme Events

The study area was the Chenyulan watershed (from 23°24′30″–23°53′30″ N latitude and 120°42′10″–121°02′30″ E longitude) in central Taiwan, as shown in Figure 1. This area is conducive to the study of debris flows caused by extreme events owing to its weak geological characteristics (many faults accompanied by fracture zones; enormous landslides and an abundant source of rock debris [31]), abundant rainfall (the average annual rainfall in the watershed is approximately 3500 mm), as well as the occurrence of the Chi-Chi earthquake (a moment magnitude of M_W 7.6, in 1999) [25,32] and many other extreme events [26,28]. Thus, many studies associated with landslides and debris flows have been conducted in this watershed [25,26,28,32].

Extreme rainfall events in Taiwan have not yet been clearly defined; however, they generally can be considered as rainfall events with accumulated rainfall (R) > 130 mm within 3 h, R > 200 mm within 6 h, or R > 350 mm within 24 h [33]. The rainfall index (RI) is defined as the product of the maximum hourly rainfall depth I_m and the maximum 24-h rainfall amount R_d (i.e., RI = R_d × I_m). Extreme rainfall events with an RI of >365 cm²/h have triggered multiple debris flows in the Chenyulan watershed [26]. The five most extreme rainfall events in the Chenyulan watershed were Typhoon Herb (TH) in 1996, Typhoon Toraji (TT) in 2001, Typhoon Mindulle (TMi) in 2004, a heavy rainstorm (HR) in 2006, and Typhoon Morakot (TM) in 2009.

3. Debris Flow Occurrence Associated with Rainfall Return Period

As shown in Figure 1, there are three meteorological stations in the study area, namely Sun Moon Lake, Yushan, and Alisan, which have recorded long-term hourly rainfall since 1963. The hourly rainfall recorded from the three meteorological stations was used to determine the representative rainfall for the Chenyulan watershed, which is the regional or average rainfall for the entire watershed. The regional hourly rainfall in the Chenyulan watershed can be determined [6] using the following equation:

I = 0.099I₁ + 0.387I₂ + 0.514I₃

(1)

where I₁, I₂, and I₃ represent the hourly rainfall records from the Sun Moon Lake, Yushan, and Alisan meteorological stations, respectively. Equation (1) calculates the areal average rainfall at the centroid of the Chenyulan watershed area, using the reciprocal-distance-squared method [34]; the weighting factors in Equation (1) represent a distance function, being inversely proportional to the square of the distance measured between the meteorological station and the centroid point. The reciprocal-distance-squared method may not actually reflect the rainfall characteristics at specific locations where local rainfall varied significantly due to elevation changes [6]. Although the different locations of meteorological stations may influence the average rainfall result, this is a simple method which can be used to directly compute the regional average rainfall characteristics of a watershed [6]. Moreover, Equation (1) can easily represent the long term variation trend in the Chenyulan watershed [6] and it was therefore used to determine the regional hourly rainfall in this study. However, it should be noted that the estimated average rainfall contains high uncertainty. Thus, three rainfall indexes, described in Section 3.1, and the probability concept were introduced in this study in the analysis of debris flow occurrence. Rainfall events in this study were identified by the following criteria: a rainfall event occurred when the hourly rainfall depth was >4 mm and ended when this value remained at <4 mm continuously for 6 h. These criteria have generally been used in official Taiwanese warning models to identify rainfall events that trigger debris flows [35].

3.1. Rainfall Return Period

Three rainfall indexes associated with debris flow occurrence, I_m, R_d, and RI, have been proposed and discussed by [26,28], and were used to analyze T in the Chenyulan watershed. The rainfall return period can be evaluated using many methods, such as the Weibull, Jenkinson, and Gringorten formulae, various computational methods, and the modified Gumbel method [34]. The Weibull formula was used in this study because it is the most efficient and commonly used for most sample data [36]; moreover, it can predict rainfall events with much shorter return periods than the other methods [37]. T can be estimated using the Weibull formula:

T = (n+1) / m

(2)

where n refers to the number of years in the record and m is the rank of the observed rainfall value in a list arranged in descending order. The rainfall index of the annual maximum series collected in the Chenyulan watershed between 1963–2017 (with the exception of 1964; no rainfall record) was used to determine T using Equation (2). A total of 53 years of rainfall index data (I_m, R_d, or RI) and T are shown in Figure 2. The relationship between the rainfall index and T can be expressed in the following form:

x_i = a (T − b) ^c

(3)

or T = (x_i / a)^1/c + b

(4)

where x_i is the variable of the rainfall index (i.e., I_m, R_d, or RI); and a, b, and c are empirical coefficients, which can be determined by regression analysis based on x_i and T data.

Table 1. Empirical coefficients a, b, and c for the rainfall indices: maximum hourly rainfall depth (I_m), maximum 24-h rainfall amount (R_d), and RI (I_m × R_d).

xi	a	b	c	R2
Im	40	0.94	0.23	0.91
Rd	410	0.96	0.26	0.95
RI	180	0.97	0.42	0.93

lists the empirical coefficients a, b, and c, and the coefficient of determination R² for rainfall indices I_m, R_d, and RI. All rainfall indexes have high R² to T.

Equations (3) and (4) represent the long-term rainfall characteristics of the Chenyulan watershed. The return period for the critical rainfall needed to trigger a debris flow can be determined using Equation (4). The x_i ‒ T relationship (Equations (3) or (4)) was developed based on 54 years of data, and the annual maximum x_i was selected. Thus, the optimal range for T was 1–54 years.

3.2. Debris Flow Events Related to the Rainfall Return Period

There are no field investigation data or scientific reports on debris flows in the Chenyulan watershed prior to 1985; therefore, debris flow events prior to 1985 were not considered in the analysis of P. Figure 3 shows rainfall events, using the rainfall indices of I_m, R_d, or RI, and rainfall-triggered debris flow events between 1985–2017. The blue lines in Figure 3 indicate various return periods (T), corresponding to three rainfall indies (I_m, R_d, or RI), which were determined by Equation (3). Data on debris flow events triggered by the five extreme rainfall events (TMi, TT, HR, TH, and TM) and the CCE are also shown in Figure 3.

The extreme rainfall events that triggered numerous debris flows had T values > 5 years. Within approximately 5 years of the CCE, debris flows could be triggered at lower I_m, R_d, or RI values, corresponding to rainfall events with lower T values. After the CCE, the critical rainfall for debris flows decreased significantly, and the critical rainfall returned to the pre-earthquake level within approximately 5 years of the earthquake. This result is consistent with those of previous studies [25,26,28]. In this paper, the period encompassing approximately 5 years after CCE is referred to as the CCE-affected period (CEAP) and is shown in Figure 3. After the CCE, debris flows could be triggered at low RI values with a T of <1 year, which was much smaller than after the other extreme rainfall events. Excluding debris flow events in the CEAP, the majority of debris flow events occurred with T > 1 year. Many rainfall-triggered debris flow events coincided with maximum annual I_m, R_d, or RI values.

4. Models for Evaluating Probability of Debris Flow Occurrence

The rainfall return period T is the average interval in years between events equaling or exceeding a certain magnitude and represents the long-term rainfall characteristics of an area. The probability of debris flow occurrence P considers rainfall events that do not trigger debris flow and represents the uncertainty of debris flow occurrence. The development of an empirical model for the relationship between P and T is presented in this section.

4.1. Probability of Debris Flow Occurrence (P)

Debris flow initiation depends on local rainfall properties, sediment availability, lithology, topography, and soil cover. These properties have a nonuniform spatial distribution and may lead to uncertainty in the evaluation of debris flow occurrence. The average hourly rainfall in the study area was obtained from the long-term rain records of the three meteorological stations. The rainfall was assumed to have a uniform distribution. The benefits of ignoring the spatial variability of hydrogeological properties are as follows: (1) more debris flow events can be identified and analyzed over the whole watershed, and (2) it assists in determining the long-term rainfall characteristics (or T) associated with the triggering of debris flows. However, ignoring spatial distribution characteristics can result in uncertainties. Thus, to reflect the resulting uncertainty, the probability concept was introduced in this study to analyze the probability of a rainfall event that triggered one or more debris flows in the entire watershed.

The P at a rainfall threshold of x_i is referred to as the probability that one or more debris flows are triggered when the rainfall amount is > x_i. P can be calculated using the following equation:

P = N_D/N_R

(5)

where N_R is the number of rainfall events for rainfall amounts > x_i. N_D is the number of those rainfall events that triggered debris flows. P is then determined at various rainfall thresholds x_i, from the maximum annual I_m, R_d, or RI or their corresponding T values. For example, if there were eight rainfall events and all of them triggered debris flows when RI > 285 cm²/h or T > 5 years (see Figure 3c), thus, N_D = N_R = 8 and P = 1.

4.2. Relationship between Probability of Debris Flow Occurrence (P) and Rainfall Return Period (T)

After obtaining datasets for the historical maximum annual I_m, R_d, or RI and their corresponding T, the empirical relationships between P and T can be determined for I_m, R_d, or RI at different periods.

Several extreme events occurred in the Chenyulan watershed during the study period, such as the CCE and five extreme rainfall events. The critical rainfall for debris flows drops significantly after extreme events [25,26,28]; however, it does recover within a certain period thereafter. Thus, periods affected by the CCE and extreme rainfall events were separated in this analysis. The CEAP in the Chenyulan watershed was defined as the 5 years after the CCE [23,32]. The maximum recovery period of critical rainfall is approximately 3 years after extreme rainfall [28]. Rainfall events occurring within the 3 years after extreme rainfall events, were selected to determine the relationship between P and T during the extreme rainfall-affected period (ERAP). Such events included the rainfall events for 3 years after Typhoon Herb and Typhoon Morakot, 2 years after the heavy rainstorm in 2006, and 1 year after Typhoon Mindulle.

Furthermore, the relationship between P and T for the whole period (WP) between 1985–2017 and the unaffected by extreme events period (UEEP) (i.e., the WP excluding CEAP and ERAP) can also be developed. The four periods (i.e., WP, CEAP, ERAP, and UEAP) and their analysis range of rainfall events are summarized and listed in

Table 2. Empirical coefficients $\mathsf{\alpha }$, $\mathsf{\beta },$ and, $\mathsf{\gamma }$ of Equation (6) during the four periods.

Period	Range of Analyzed Rainfall Events	$\mathit{\alpha }$	$\mathit{\beta }$	γ	R2
I. Whole period (WP)	Between 1985–2017	0.47	0.78	0.85	0.89
II. Chi-Chi earthquake-affected period (CEAP)	Between CCE and TMi (the period of approximately 5 years after the CCE)	0.23	0.87	8.65	0.56
III. Extreme rainfall-affected period (ERAP)	3 years after TH and TM, 2 years after HR, and 1 year after TMi	0.32	0.79	2.89	0.84
IV. Unaffected by extreme events period (UEEP)	WP excluding CEAP and ERAP	0.20	0.46	0.99	0.82

Notes: R² = coefficient of determination.

. Figure 4 shows the relationships between P and T during WP, CEAP, ERAP, and UEEP. Whether using I_m, R_d, or RI indices, the data describing the relationship between P and T showed a similar distribution in all periods. The P‒T relationship can be expressed in the form of a Weibull distribution, i.e.;

$$\mathrm{P}=1-{\mathrm{e}}^{-\mathsf{\alpha }{\left(\frac{\mathrm{T}}{\mathsf{\beta }}\right)}^{\mathsf{\gamma }}}$$

(6)

The empirical coefficients $\mathsf{\alpha }$, $\mathsf{\beta },$ and $\mathsf{\gamma }$ for the four periods were determined by fitting the given data and are listed in Table 2. The P–T relationship models of the four periods are compared in Figure 5. The WP model used long-term data between 1985–2017, and also included data that were affected by extreme events, that is, five extreme rainfall events and the CCE. The UEEP model excluded the influence of extreme events, and the P value predicted by the UEEP model was slightly smaller than that predicted by the WP model at the same T. Compared to the UEEP model, P increased significantly after an extreme rainfall event or the CCE, under the same rainfall conditions or at the same T. In particular, the P value affected by the CCE was markedly higher than that affected by extreme rainfall.

The benefits of determining the P–T relationship (Figure 5) are as follows: (1) all rainfall indexes can be integrated by the P–T relationship, and (2) the P values during periods affected by the CCE or extreme rainfall can be evaluated at various T values (or different rainfall conditions). For example, P is 48% at T = 1.5, when the area of the Chenyulan watershed is unaffected by extreme events (see light green curve in Figure 5), whereas P increases to 87% after an extreme rainfall event (see blue curve in Figure 5), and P = 100% after the CCE (see red curve in Figure 5). When affected by the CCE, the P value can reach up to approximately twice that of a period unaffected by extreme events, and P during a period affected by the extreme rainfall events can be approximately 80% higher than during a period unaffected by extreme events at the same T.

5. Discussion

The heavy rainfall caused by Typhoon Morakot in August 2009 was one of the extreme rainfall events in the Chenyulan watershed and had an I_m of 85.5 mm, R_d of 1192.6 mm, and RI = 1019.7 cm²/h. It caused numerous debris flows that buried more than 20 houses in Shenmu, Tongfu, and Xinyi villages. After Typhoon Morakot, there were seven debris flow events triggered by rainstorms or typhoons between 2009–2020, as shown in

Table 3. Debris flow event and related rainfall characteristics in the Chenyulan watershed between 2009–2020, after the extreme rainfall events of Typhoon Morakot (TM).

No.	Year	Date	Name of Event	Im (mm/h)	Rd (mm)	RI (cm2/h)	Remark
0	2009	6–11 August	Typhoon Morakot (TM)	85.5	1192.6	1019.7	Extreme rainfall event
1	2010	23–24 May	Rainstorm	35.8	227.2	81.3	Events within the period affected by extreme rainfall TM
2	2011	17–20 July	Rainstorm	33.6	256.2	86.1
3	2012	9–12 June	Rainstorm	33.6	384.6	129.2
4	2012	18–21 June	Typhoon Talim	22.2	243.4	54.0
5	2012	1–3 August	Typhoon Saola	36.4	502.2	182.8
6	2013	12–13 July	Typhoon Soulik	52.4	661.7	346.7	Events unaffected by TM
7	2017	1–4 June	Rainstorm	46.1	515.7	237.7

Notes: I_m (mm/h), maximum hourly rainfall in each rainfall event; R_d (mm), maximum 24-h rainfall in each rainfall event; RI, rainfall index for each rainfall event, which can be determined using the equation RI, R_d × I_m.

. The relationship between P and T was used to evaluate P after the rainfall events of Typhoon Morakot.

5.1. Procedures

The empirical model for evaluating P was applied through the following steps:

• The hourly rainfall data from three metrological stations was input and the regional hourly rainfall was evaluated using Equation (1).
• I_m, R_d, and RI were determined from regional rainfall data.
• T was determined using Equation (4) with I_m, R_d, or RI.
• P was determined using the P–T relationship. There are three P–T relationships with different empirical coefficients based on various periods (Equation (6) and Table 2).

That is, the P–T relationship of the WP is:

$$\mathrm{P}=1-e^{-0.47{\left(\frac{\mathrm{T}}{0.78}\right)}^{0.85}};$$

(7)

the P–T relationship of the ERAP is:

$$\mathrm{P}=1-e^{-0.32{\left(\frac{\mathrm{T}}{0.79}\right)}^{2.89}};$$

(8)

and the P–T relationship of the UEEP is:

$$\mathrm{P}=1-e^{-0.20{\left(\frac{\mathrm{T}}{0.46}\right)}^{0.99}},$$

(9)

where T can be determined using Equation (4). The coefficients in Equation (4) depend on the parameters I_m, R_d, or RI, i.e., T = (I_m/40)^4.35 + 0.94, T = (R_d/410)^3.85 + 0.96, and T = (RI/180)^2.38 + 0.97. P can be determined when I_m, R_d, or RI is given according to Equations (7), (8) or (9).

However, the three equations were developed based on different periods and datasets, and the valid conditions for the three equations may not be identical.

Equation (7) predominantly reflects the long-term characteristics of debris flow occurrence and may not reflect the short-term characteristics caused by extreme events. In contrast, Equation (8) focuses on the influence of extreme rainfall events. Equation (9) reflects P during periods that are unaffected by extreme events, such as significant earthquakes and extreme rainfall. Hence, field data of debris flow occurrence and rainfall events between 2009‒2017 were used to assess the suitability of the proposed equations and their adopted rainfall parameters (I_m, R_d, or RI).

5.2. Results

The model for the WP based on the return period of different rainfall indices was used to evaluate P. Figure 6 shows the variation in the predicted P derived from Equation (7) for rainfall events with return periods greater than 1 year (T > 1 year) from 2012‒2020, and the associated debris flow events are labelled (note: red dots and those labeled No. 1‒7 correspond to the events shown in Table 3).

In Figure 6a, using the T of I_m to evaluate P, 10 rainfall events that did not trigger debris flows were predicted with P ≥ 0.5, although most of the debris flow events (i.e., six of seven) were reasonably predicted by P values > 0.5. In Figure 6b,c, using the return periods of R_d and RI to evaluate P, respectively, only four of the seven debris flow events were predicted with occurrence probabilities of P ≥ 0.5. The occurrence probability P for the four rainfall events that did not trigger debris flows was overestimated with values of P ≥ 0.5.

The model for the WP mainly presents the long-term characteristics for P and is unable to respond to the occurrence of a debris flow that is affected or unaffected by extreme events. Thus, the predicted P for the period of 3 years after Typhoon Morakot was determined using the ERAP model (Equation (8)), and the predicted P after the ERAP was determined using the UEEP model (Equation (9)). As shown in Figure 7, three rainfall indexes (I_m, R_d, and RI) were used to evaluate the T values that were associated with the ERAP and UEEP models.

Using the return period T of I_m, R_d, and RI (Figure 7a–c), six of the seven debris flow events were predicted with P ≥ 0.5. One debris flow event occurred at P ≥ 0.5. Most rainfall events that did not trigger debris flows were predicted with P ≥ 0.5. There were six debris flows within a total of ten rainfall events for T of I_m, six debris flows within eight rainfall events for T of R_d, and six debris flows within eight rainfall events for T of RI when P ≥ 0.5. Figure 7 shows that 30% (Figure 7a), 60% (Figure 7b), and 75% (Figure 7c) of rainfall events triggered debris flows for T of I_m, R_d, and RI, respectively, when predicted with P ≥ 0.5. That is, the predicted P using the T of RI and R_d would be more accurate than that using T of I_m. Furthermore, [25] examined the period affected by extreme rainfall that triggers debris flow and the modification of critical rainfall for debris flows after extreme events in the Chenyulan watershed. Three rainfall indices, I_m, R_d, and RI, were used. The modifications of critical rainfall and recovery period have higher correlations with RI driven by extreme rainfall than with I_m and R_d. Thus, the RI index was suggested for convenience. When applying the model, one must initially check whether extreme rainfall (i.e., with RI > 365 cm²/h) has occurred in the Chenyulan watershed within 3 years. If so, the model that responded to the critical reduction RI (Equation (8)) is suggested; if not, Equation (9) should be used.

6. Conclusions

This study investigated the occurrence of debris flows after extreme events such as the CCE and extreme rainfall in a 449 km² region of the Chenyulan watershed. The probability of debris flow occurrence was analyzed to determine the uncertainty of the hydrogeological conditions in the study area. The innovation of this study was the combining of the concepts of T and P to establish a probabilistic model of debris flow occurrence that can be adjusted following extreme events. The findings of this study are as follows:

• The rainfall indices of I_m, R_d, and RI associated with T were analyzed. The advantages of using T are as follows: (1) it can reflect the long-term rainfall variation characteristics of a region and adjust with rainfall variability, and (2) the concept of T can be easily combined with hydrological analysis to facilitate the subsequent engineering design. Most debris flow events occurred with T > 1 year, excluding debris flow events affected by the CCE, while the T of extreme events triggering numerous debris flows could exceed approximately 5 years. T affected by the CCE for debris-flow-triggering rainfall was < 1 year, (i.e., much smaller than that affected by other extreme rainfall events).
• The rainfall indices of I_m, R_d, and RI can be used to determine the relationship between P and T. Four empirical models of the P‒T relationship were developed (i.e., CEAP, ERAP, WP, and UEEP models). The P values, affected or unaffected by extreme events (such as the CCE or extreme rainfall), can be evaluated at various T values (or different rainfall conditions) using the P‒T relationship.
• P increased significantly after extreme rainfall events or the CCE, with the same T. In particular, the P value of periods influenced by CCE was higher than that during periods influenced by extreme rainfall events. The P during a period affected by the CCE can reach up to approximately twice that of a period that is unaffected by extreme events, while P of periods affected by extreme rainfall events can be approximately 80% higher than the P of a period that is unaffected by extreme events at the same T.
• A model relating P and T was applied to estimate P during recent rainfall events (2009–2020) after the extreme rainfall of Typhoon Morakot, which showed that a model using the R_d or RI index was reasonably accurate at predicting debris flow occurrence.
• An empirical model for evaluating P was developed based on the regional characteristics of the Chenyulan watershed, and may not be applicable to areas with different hydrogeological properties. The suitability of this model must be assessed, and empirical coefficients will most likely be required for calibration if the model is applied to other areas. In addition, the RI index was derived empirically, and further studies on the association of the RI index with physical mechanisms are needed.