Mapping and Analyzing the Evolution of the Butangbunasi Landslide Using Landsat Time Series with Respect to Heavy Rainfall Events during Typhoons

: Large rainfall-induced landslides are among the most dangerous natural hazards in Taiwan, posing a risk for people and infrastructure. Thus, better knowledge about the evolution of landslides and their impact on the downstream area is of high importance for disaster mitigation. The aim of this study is twofold: (1) to semi-automatically map the evolution of the Butangbunasi landslide in south-central Taiwan using satellite remote sensing data, and (2) to investigate the potential correlation between changes in landslide area and heavy rainfall during typhoon events. Landslide area, as well as temporary landslide-dammed lakes, were semi-automatically identiﬁed using object-based image analysis (OBIA), based on 20 Landsat images from 1984 to 2018. Hourly rainfall data from the Taiwan Central Weather Bureau (CWB) was complemented with rainfall data from Climate Hazards Group Infrared Precipitation with Station data (CHIRPS) to examine the potential relationship between landslide area changes and rainfall as a triggering factor. The OBIA mapping results revealed that the most signiﬁcant landslide extension happened after typhoon Morakot in 2009. We found a moderate positive relationship between the landslide area change and the duration of the heavy rainfall event, whereas daily precipitation, cumulative rainfall and mean intensity did not present strong signiﬁcant correlations.


Introduction
The mountains of Taiwan, with the highest peaks rising to almost 4000 m a.s.l., are characterized by fractured rock formations, high relief and steep stream gradients. These mountains, particularly the central mountain range (CMR), influence the tracks and intensity of typhoon events [1]. During summer and autumn, Taiwan is regularly affected by typhoons (tropical cyclones), three to four per year on average, which bring heavy rainfall [2]. Some studies indicate that the number of typhoons hitting Taiwan increased after the year 2000, also resulting in heavier precipitation in recent years [3,4].
Heavy rainfall is the main landslide triggering factor in Taiwan [5]. The rainfall events are usually associated with typhoons, which account for almost 50% of the island's total rainfall [6,7]. Studies have found relationships between long-duration and moderate intensity rainfall events and large and deep-seated landslides [8], sediment yield and debris flows [5,9].
Large, rapid landslides and debris flows frequently lead to fatalities. Major endeavors are necessary to protect people and settlements in areas at risk, and to implement prevention and early warning measures [10][11][12][13]. In general, as well as in Taiwan, such catastrophic landslides also cause severe damages to infrastructure, and efforts both in time and money are needed to recover and maintain the transportation infrastructure such as roads and bridges [14][15][16][17]. Beyond the direct landslide hazard, large landslides can initiate natural hazard cascades by damming rivers and inducing catastrophic flash floods and debris flows [18][19][20]. The large amount of mobilized debris that originates from landslides significantly affects the drainage system, for example, resulting in an increase in erosion and sediment discharge in rivers, and changes in channel size and shape [5,21,22]. According to Chen et al. [9] about 384 Mt y −1 of sediment is transported into the ocean in Taiwan, whereby the high proportion of large landslides significantly contributes to this high annual sediment yield.
In August 2009, typhoon Morakot caused a record-breaking cumulative rainfall (more than 2000 mm in three days), which led to debris flows and mudflows, flooding in coastal areas, and massive landslides [23,24]. The rainfall also triggered one of the most famous and fatal landslides in Taiwan, the Xiaolin landslide [25][26][27][28]. The Butangbunasi landslide [29][30][31] is another example of a large rainfall-triggered landslide in Taiwan. Reactivation and extension of this landslide have been resulting in repeated sediment delivery to the Laonong River, especially during torrential rainfall. The river course has been frequently affected, leading to the formation of a landslide-dammed lake several times during the past three decades [29,[32][33][34][35][36]. The magnitude of the Butangbunasi landslide is even significantly larger than of the disastrous Xiaolin landslide [29]. A deeper knowledge of the evolution of landslides and their triggering factors is crucial for hazard mitigation [37,38]. Therefore, mapping and analyzing the evolution of such large landslides over time helps to better understand their reactivation rates and their impact on downstream areas.
Remote sensing plays a key role in studying landslides and provides an adequate and cost-effective source to derive information about landslide distribution and types [39][40][41][42]. The use of remote sensing data also helps to investigate the potential impacts of landslides such as the damming of rivers, particularly in difficult to access and remote mountain regions [43]. The value of remote sensing for landslide studies becomes more and more evident with the increasing amount of freely available Earth observation (EO) data, which provides remarkable opportunities to map and monitor landslides over time [44][45][46].
Object-based image analysis (OBIA) provides a suitable methodological framework for efficient landslide mapping, as well as landslide change analysis [47,48]. By working on the object-level instead of the pixel-level, OBIA allows considering spectral, spatial, textural, morphometric and hierarchical properties for the classification of landslides [49][50][51]. Moreover, it is argued that using OBIA yields better classification accuracies than pixel-based classifications [51][52][53][54][55]. Several studies employed OBIA for landslide mapping and landslide change detection in Taiwan [32,49,53,[56][57][58][59], but none of them used time series of images for investigating the evolution and reactivation of an active large landslide.
The aim of this study is to analyze the evolution of the Butangbunasi landslide in south-central Taiwan using OBIA and time series of freely available Landsat images and to investigate the potential correlation between changes in landslide area and heavy rainfall events during selected typhoon (including tropical storm) events.

Study Area
For this study, we selected the Butangbunasi landslide, which is located in the Taoyuan District of Kaohsiung City, south-central Taiwan (R.O.C.; Figure 1). The study area comprised of two geological strata, the Changchihkeng formation, which is composed of sandstone and shale, and the Tangenshan formation, which consists of massive sandstones and is characterized by joints that can form precipitous scarps and deep ravines [29]. The area experiences between 2000 and 4000 mm of annual rainfall within a subtropical monsoon climate [34]. Reactivation and extension of this large landslide area have been taking place since the 1980s and have been resulting in repeated pulses of sediment delivered to the Laonong River, especially during torrential rainfall brought by typhoons.

Study Area
For this study, we selected the Butangbunasi landslide, which is located in the Taoyuan District of Kaohsiung City, south-central Taiwan (R.O.C.; Figure 1). The study area comprised of two geological strata, the Changchihkeng formation, which is composed of sandstone and shale, and the Tangenshan formation, which consists of massive sandstones and is characterized by joints that can form precipitous scarps and deep ravines [29]. The area experiences between 2000 and 4000 mm of annual rainfall within a subtropical monsoon climate [34]. Reactivation and extension of this large landslide area have been taking place since the 1980s and have been resulting in repeated pulses of sediment delivered to the Laonong River, especially during torrential rainfall brought by typhoons. The accumulation of debris and sediments affects the river course and can lead to the damming of the river as it happened several times in the past. The main landslide area is difficult to access and hardly visible from the Laonong River valley, where major construction efforts are needed to maintain and rebuild the transportation infrastructure such as roads and bridges. Figure 2 gives an impression of the Butangbunasi landslide.

Optical Satellite Data
We used time series of optical satellite data, i.e., Landsat 5, Landsat 7 and Landsat 8 imagery (20 images overall) with 30 m spatial resolution, from 1984 to 2018 to semi-automatically map the evolution of the Butangbunasi landslide (Table 1). In particular, we selected the first cloud-free image of the area of interest available after a typhoon (including tropical storm) event that shows a noticeable change in the landslide area compared to pre-event images. This was done based on a visual inspection of the Landsat database and considering the historical hurricane/typhoon tracks of the National Oceanic and Atmospheric Administration (NOAA) Office for Coastal Management [60]. In addition, all major hurricanes crossing Taiwan within a radius of 100 km around the Butangbunasi The accumulation of debris and sediments affects the river course and can lead to the damming of the river as it happened several times in the past. The main landslide area is difficult to access and hardly visible from the Laonong River valley, where major construction efforts are needed to maintain and rebuild the transportation infrastructure such as roads and bridges. Figure 2 gives an impression of the Butangbunasi landslide.

Optical Satellite Data
We used time series of optical satellite data, i.e., Landsat 5, Landsat 7 and Landsat 8 imagery (20 images overall) with 30 m spatial resolution, from 1984 to 2018 to semi-automatically map the evolution of the Butangbunasi landslide (Table 1). In particular, we selected the first cloud-free image of the area of interest available after a typhoon (including tropical storm) event that shows a noticeable change in the landslide area compared to pre-event images. This was done based on a visual inspection of the Landsat database and considering the historical hurricane/typhoon tracks of the National Oceanic and Atmospheric Administration (NOAA) Office for Coastal Management [60]. In addition, all major hurricanes crossing Taiwan within a radius of 100 km around the Butangbunasi landslide, which reached category H3 or above according to the Saffir-Simpson hurricane wind scale (SSHWS), were identified and the respective post-event Landsat imagery added to the database. The Landsat 5 scene from 1984 is the first available image for the study area and served as a starting point for the analysis. The Landsat 8 scene from 2018 does not follow a specific typhoon event but was used as the final image for the analysis to see if any changes in landslide area are also identified following a period without a heavy rainfall event. Moreover, this image temporally coincides with the field visit in November 2018 (cf. Figure 2c). landslide, which reached category H3 or above according to the Saffir-Simpson hurricane wind scale (SSHWS), were identified and the respective post-event Landsat imagery added to the database. The Landsat 5 scene from 1984 is the first available image for the study area and served as a starting point for the analysis. The Landsat 8 scene from 2018 does not follow a specific typhoon event but was used as the final image for the analysis to see if any changes in landslide area are also identified following a period without a heavy rainfall event. Moreover, this image temporally coincides with the field visit in November 2018 (cf. Figure 2c). All Landsat scenes were downloaded as Level-1 products from the EarthExplorer user interface of the United States Geological Survey (USGS). In order to increase the comparability across the Landsat imagery captured by different Landsat sensors and/or at different times, a top-of-atmosphere (TOA) calibration was performed on all satellite images [61]. To do so, the apparent reflectance function implemented in the ArcGIS Desktop version 10.7 was applied. The algorithm uses sun elevation, satellite position, acquisition date and sensor properties such as gain and bias settings for each spectral band [62]. In this way, the spectral differences between images from different dates and sensors were decreased. In addition to the Landsat images, the derived slope layer from the Advanced Land Observing Satellite (ALOS) Palsar digital elevation model (DEM; 12.5 m spatial resolution), acquired in 2008, was used as ancillary data to support the landslide mapping. All Landsat scenes were downloaded as Level-1 products from the EarthExplorer user interface of the United States Geological Survey (USGS). In order to increase the comparability across the Landsat imagery captured by different Landsat sensors and/or at different times, a top-of-atmosphere (TOA) calibration was performed on all satellite images [61]. To do so, the apparent reflectance function implemented in the ArcGIS Desktop version 10.7 was applied. The algorithm uses sun elevation, satellite position, acquisition date and sensor properties such as gain and bias settings for each spectral band [62]. In this way, the spectral differences between images from different dates and sensors were decreased. In addition to the Landsat images, the derived slope layer from the Advanced Land Observing Satellite (ALOS) Palsar digital elevation model (DEM; 12.5 m spatial resolution), acquired in 2008, was used as ancillary data to support the landslide mapping.

Typhoon and Rainfall Data
The International Best Track Archive for Climate Stewardship (IBTrACS) data repository [63] compiles the hurricane/typhoon best-track position and intensity from different data sources and makes them available in one consolidated archive. We obtained the fourth version of this data in shapefile format to extract the spatial information corresponding to 19 typhoon events that had an impact on the landslide area or had a SSHWS category three (H3) or higher within a 100 km radius from the Butangbunasi landslide (Table 2). As for the rainfall data, Taiwan's Central Weather Bureau (CWB) provides the open Typhoon Database [64], which compiles various sources of information relevant to the study and the monitoring of typhoon events in the western North Pacific. Among their data sources, CWB has made station weather data available, which includes hourly precipitation for each typhoon event for all the automatic rain gauge stations in Taiwan. For this study, we identified the three closest stations to the Butangbunasi landslide, all located within the Laonong River basin, which are summarized in Table 3. The locations of the three selected CWB stations are shown in Figure 3. Since station C1V210 was upgraded to a new weather station coded C0V210, which continues the data registration, we will refer to it as a single station called Fuxing. As for the rainfall data, Taiwan's Central Weather Bureau (CWB) provides the open Typhoon Database [64], which compiles various sources of information relevant to the study and the monitoring of typhoon events in the western North Pacific. Among their data sources, CWB has made station weather data available, which includes hourly precipitation for each typhoon event for all the automatic rain gauge stations in Taiwan. For this study, we identified the three closest stations to the Butangbunasi landslide, all located within the Laonong River basin, which are summarized in Table 3. The locations of the three selected CWB stations are shown in Figure 3. Since station C1V210 was upgraded to a new weather station coded C0V210, which continues the data registration, we will refer to it as a single station called Fuxing. Since the entire study period was not covered by the CWB rain gauge station data available on the Typhoon Database, we also collected data from the CHIRPS (Climate Hazards Group Infrared Precipitation with Station data) dataset. CHIRPS combines information from the Tropical Rainfall Measuring Mission (TRMM) 3B42 product, the NOAA Climate Forecast System (CFS) and other precipitation data sources to provide daily precipitation in grid format at a resolution of 0.05°, for latitudes between 50° S and 50° N [65,66]. We accessed the data through the Google Earth Engine, where we obtained the mean rainfall in the area of interest for the 19 typhoon events selected for Since the entire study period was not covered by the CWB rain gauge station data available on the Typhoon Database, we also collected data from the CHIRPS (Climate Hazards Group Infrared Precipitation with Station data) dataset. CHIRPS combines information from the Tropical Rainfall Measuring Mission (TRMM) 3B42 product, the NOAA Climate Forecast System (CFS) and other precipitation data sources to provide daily precipitation in grid format at a resolution of 0.05 • , for latitudes between 50 • S and 50 • N [65,66]. We accessed the data through the Google Earth Engine, where we obtained the mean rainfall in the area of interest for the 19 typhoon events selected for analysis. For each event, we obtained the daily precipitation on the day the typhoon passed closest to the Butangbunasi landslide area.

Semi-Automated Landslide Mapping
For semi-automatically mapping the Butangbunasi landslide area on each Landsat image, OBIA was used. The analyses were conducted using the eCognition (Trimble) software. We defined a set of knowledge-based classification rules for mapping the landslide area per image. Since images from 20 different points in time were used, efforts were made to design a transferable classification routine that could be applied to all images without or with only minor adaptations of the classification thresholds ( Table 4). The classification rules were developed based on the Landsat 5 image from 1984 and then transferred to the other images. Table 4. Segmentation and classification parameters used for object-based landslide mapping.

Bands for Segmentation Classification Parameters
Landsat 5 (1984,1989,1990,1992,1994,1996,1998 As a first step, the normalized difference vegetation index (NDVI) and the modified soil-adjusted vegetation index (MSAVI), as well as a brightness layer were calculated. Next, we applied the multiresolution segmentation [67] for creating image objects as a basis for the classification. The multiresolution segmentation in eCognition is a bottom-up region merging technique. It starts with single pixel objects and merges them stepwise to larger objects based on local homogeneity criteria that describe the similarity of adjacent image objects [67]. The scale parameter and homogeneity criteria (shape vs. color weighting, compactness vs. smoothness weighting) determine the maximum allowed heterogeneity and control the average size and shape of the resulting image objects. The segmentation parameters (Table 4) were selected based on an expert-guided trial and error approach and a visual assessment of resulting image objects.
The knowledge-based classification primarily relied on the usage of the calculated spectral indices. The main indication for mapping the landslide area was the absence of vegetation, which leads to a distinctive spectral contrast between the landslide-affected area and its surroundings, especially in densely vegetated regions [49,58,68]. This change in land cover related to landslide occurrence can be well represented with spectral indices such as the NDVI [51]. The same classification parameters were applied to all Landsat 5 images, except the March 2008 image, where the thresholds for the NDVI and MSAVI were slightly adapted. The reason for this was that this image is the only one acquired in spring, i.e., before the rainy season, and the seasonal variation leads to a slightly different spectral reflection. Minor adaptations were also needed for the Landsat 7 and Landsat 8 images, while the same classification rules and thresholds were used per sensor.
Since the DEM used in this study only represents one point in time, the topographic signature of the landslide at different dates cannot be accurately represented by the DEM data [49]. Thus, the DEM and the slope were mainly used as ancillary data to avoid the classification of obvious false positives, for example, debris accumulation areas with low slopes at low elevations in the Laonong river bed.
Finally, a few true but unwanted positives such as landslides on slopes not belonging to the Butangbunasi landslide area were manually removed, and small landslide objects in shadow areas that were obviously missed were manually added.
Some of the images revealed a temporary landslide-dammed lake where the debris from Butangbunasi reaches and blocks the Laonong River. The lake area was additionally mapped with OBIA by using relatively low NDVI and near-infrared values.
Finally, changes in landslide and lake area were calculated as the area change in the respective class between two successive images within the time series.

Accuracy Assessment of OBIA Results
To assess the classification accuracy, we compared the OBIA results to results from visual image interpretation. Reference data were exemplarily created based on three selected images, whereby one image from each Landsat sensor was used (Landsat 8 from 4 December 2016; Landsat 7 from 24 August 2008; Landsat 5 from 27 September 2000). Another criterion for the selection of these images was the presence of a lake close to the landslide area so that the accuracy of both classes could be assessed. Manual digitizing from the satellite imagery was carried out in ArcGIS 10.7 at a scale of 1:25,000. We assessed the classification accuracy based on the overlapping area between the manual and the OBIA mappings. Producer's accuracies were computed by dividing the correctly classified area (overlap area) by the total area of the reference data (i.e., the manual mapping results), and user's accuracies were obtained by dividing the correctly classified area by the total area mapped by OBIA [37,45,69]. The producer's accuracy is the map accuracy from the point of view of the creator of the map (here: the producer of the OBIA classification). It is a measure of omission error (error of exclusion) [69]. The user's accuracy is the accuracy from the point of view of a map user, not the map creator. It is a measure of commission error (error of inclusion) [69].

Analysis of Rainfall Data during Typhoon Events
Rainfall data registered during the selected typhoon and tropical storm events were analyzed to find potential relationships between the reactivation and evolution of the landslide area and heavy rainfall. From the CHIRPS data, we extracted daily precipitation for the date when the typhoon was closest to the landslide. For the CWB rain gauge station data, hourly precipitation data allowed us to identify rainfall events during the selected typhoons. We followed the definition of a rainfall event as explained in Chang et al. [5] and Chen et al. [8], where an event starts when there is more than 4 mm of rain registered in the gauge during an hour, and ends when the rain registered is lower than 4 mm for six consecutive hours. We calculated the duration, cumulative rainfall and intensity for each rainfall event. Rainfall duration is defined as the hours the event lasted; cumulative rainfall is the amount of precipitation registered during the entire rainfall event; rainfall intensity is the average amount of precipitation registered per hour during a rainfall event. The rainfall event with the highest duration and intensity within each typhoon event was selected for further analysis.

Correlation between the Change in Landslide Area and Rainfall
An exploratory analysis of the variables derived from the rainfall data, as well as of the landslide area change, showed a non-parametric distribution of the data with the presence of outliers. Hence, we tested for rank correlation after Spearman between the landslide area change and the derived rainfall parameters for each CWB rain gauge station and for the daily precipitation extracted from the CHIRPS data. The Spearman rank-order correlation coefficient (ρ) determines whether there is a monotonic relationship between variables, that is, two variables increasing or decreasing together, or as one increases, the other decreases. Unlike the Pearson correlation coefficient (p), which is a measure of the linear correlation between two variables, Spearman does not test for a linear relationship with a constant rate. The Spearman coefficient is more robust to outliers and appropriate for skewed distributions [70,71], as observed for our data. Finally, the correlation results were evaluated at a 90% confidence level [72]. Figure 4 shows the OBIA landslide mapping results for each Landsat image. In addition to the Butangbunasi landslide area, we detected a temporary landslide-dammed lake at the confluence of the Butangbunasi tributary with the Laonong River on six images.

Semi-Automated Mapping Results
Each selected typhoon event was linked to the respective Landsat image date when the OBIA mapping of the landslide area was performed. Figure 5 shows the time series where landslide area evolution was tracked along with the preceding typhoon events. After typhoon Morakot in 2009, an abrupt increase in landslide area was identified, keeping a steady or even slightly decreasing trend for the following years. Section 3.3 provides additional information about the evolution of the Butangbunasi landslide and the temporary landslide-dammed lake in the Laonong river bed with respect to the heavy rainfall events.

Comparison of OBIA Results with Visual Interpretation
The accuracy of the semi-automated OBIA results was assessed by comparison to results from visual interpretation. Table 5 summarizes the results of the accuracy assessment. Similar accuracy values were achieved for the OBIA classification among the selected images, whereby higher accuracies were reached for the landslide class compared to the lake class. Mixed pixels-or in our case mixed objects as a result of undersegmentation [73]-are an important issue in the identification and proportion estimation of classes in Landsat satellite scenes, since they cover more than one ground cover type and thus decrease the separability of classes [74]. We faced this problem especially in areas where small patches of vegetation within the landslide area exist or where revegetation leads to a sparse vegetation cover that influences the spectral reflectance on post-event images. This results in small differences in the OBIA landslide mapping compared to the visual interpretation. For the detected landslide-dammed lakes, classification uncertainties are mainly associated with mixed objects, partly shallow water areas with high sediment load or wet areas in the river bed.
When interpreting the accuracy values it has to be considered that results from visual expert interpretation, especially for natural phenomena such as landslides, cannot constitute an entirely 'true' reference, as their creation depends on various factors such as the data used or the expertise of the interpreter [49].

Comparison of OBIA Results with Visual Interpretation
The accuracy of the semi-automated OBIA results was assessed by comparison to results from visual interpretation. Table 5 summarizes the results of the accuracy assessment. Similar accuracy values were achieved for the OBIA classification among the selected images, whereby higher accuracies were reached for the landslide class compared to the lake class. Mixed pixels-or in our case mixed objects as a result of undersegmentation [73]-are an important issue in the identification and proportion estimation of classes in Landsat satellite scenes, since they cover more than one ground cover type and thus decrease the separability of classes [74]. We faced this problem especially in areas where small patches of vegetation within the landslide area exist or where revegetation leads to a sparse vegetation cover that influences the spectral reflectance on post-event images. This results in small differences in the OBIA landslide mapping compared to the visual interpretation. For the detected landslide-dammed lakes, classification uncertainties are mainly associated with mixed objects, partly shallow water areas with high sediment load or wet areas in the river bed.
When interpreting the accuracy values it has to be considered that results from visual expert interpretation, especially for natural phenomena such as landslides, cannot constitute an entirely 'true' reference, as their creation depends on various factors such as the data used or the expertise of the interpreter [49].

Rainfall Data Analysis for Each Typhoon Event
Two data sources of rainfall data were individually analyzed for this study. CHIRPS daily precipitation was obtained for the 19 typhoon events selected, where the highest daily precipitation registered on the day of the typhoon passing over Taiwan was 208 mm for typhoon Fung-Wong, while the lowest daily precipitation registered was 43.7 mm for typhoon Sarah. The mean daily precipitation for the selected events was 120 ± 48.2 mm. Nevertheless, a comparison of this satellite-derived data with rain gauge data from the CWB showed that the CHIRPS dataset under/overestimates the rainfall amount for different events (Figure 6). This observation is supported by the findings of Chen et al. [24], who found that remote sensing rainfall products underestimate rainfall compared to rain gauges for typhoon Morakot over Taiwan.

Rainfall Data Analysis for Each Typhoon Event
Two data sources of rainfall data were individually analyzed for this study. CHIRPS daily precipitation was obtained for the 19 typhoon events selected, where the highest daily precipitation registered on the day of the typhoon passing over Taiwan was 208 mm for typhoon Fung-Wong, while the lowest daily precipitation registered was 43.7 mm for typhoon Sarah. The mean daily precipitation for the selected events was 120 ± 48.2 mm. Nevertheless, a comparison of this satellitederived data with rain gauge data from the CWB showed that the CHIRPS dataset under/overestimates the rainfall amount for different events (Figure 6). This observation is supported by the findings of Chen et al. [24], who found that remote sensing rainfall products underestimate rainfall compared to rain gauges for typhoon Morakot over Taiwan. Given that the CHIRPS data are freely available and span the entire study period, we decided to include it in the study. However, the correlation results of the mapped landslide area with CHIRPS daily precipitation data should be interpreted carefully.
For the CWB station data, we analyzed rainfall event parameters derived from hourly precipitation records for three separate rain gauges. Several studies used spatial interpolation for the analysis of rainfall parameters in comparison to landslide area [8,75,76]. These studies investigate several landslide locations but do not focus on a detailed analysis of the progressive evolution of one major landslide. Given that spatial interpolation can introduce uncertainty to the resulting rain measurements [77] and that the rain gauges were located within a range of 2.8-8 km from the landslide area, we decided to analyze the three closest rain gauges separately.
The rain gauge data are available since 1992, and thus cover only 17 of the selected 19 typhoon events. In addition, rainfall data for the Longwang and Haitang typhoons were not available for station Xiaoguanshan. Cumulative rainfall, duration and mean intensity were derived for each rainfall event within the typhoon period for each station. Summary statistics of the derived parameters are presented in Table 6. In general, station Xiaoguanshan recorded the highest cumulative rainfall, duration and intensity during rainfall events, which could be explained by the altitude of the station, as well as its location downstream of the valley. Given that the CHIRPS data are freely available and span the entire study period, we decided to include it in the study. However, the correlation results of the mapped landslide area with CHIRPS daily precipitation data should be interpreted carefully.
For the CWB station data, we analyzed rainfall event parameters derived from hourly precipitation records for three separate rain gauges. Several studies used spatial interpolation for the analysis of rainfall parameters in comparison to landslide area [8,75,76]. These studies investigate several landslide locations but do not focus on a detailed analysis of the progressive evolution of one major landslide. Given that spatial interpolation can introduce uncertainty to the resulting rain measurements [77] and that the rain gauges were located within a range of 2.8-8 km from the landslide area, we decided to analyze the three closest rain gauges separately.
The rain gauge data are available since 1992, and thus cover only 17 of the selected 19 typhoon events. In addition, rainfall data for the Longwang and Haitang typhoons were not available for station Xiaoguanshan. Cumulative rainfall, duration and mean intensity were derived for each rainfall event within the typhoon period for each station. Summary statistics of the derived parameters are presented in Table 6. In general, station Xiaoguanshan recorded the highest cumulative rainfall, duration and intensity during rainfall events, which could be explained by the altitude of the station, as well as its location downstream of the valley. Table 7 presents the individual parameters per typhoon event next to the semi-automated mapping results. For every station, the maximum cumulative rainfall and duration correspond to typhoon Morakot, whereas the maximum intensity corresponds to typhoon Toraji, and the minimum duration corresponds to typhoon Nepartak. Minimum cumulative rainfall and intensity varied between typhoon events and stations.

Relation between Landslide Area Change and Rainfall-Derived Parameters during Typhoon Events
No cloud-free satellite image was available between typhoons Megi and Nepartak (cf. Figure 5). Hence, the landslide area mapped for December 2016 could represent either of those events. To avoid any uncertainty, both typhoon events were removed from the correlation analysis.
CHIRPS daily precipitation data for the exact date when the typhoon passed over Taiwan were correlated with the landslide area change for 17 typhoon events. The results show that there was not enough evidence of a significant correlation between the variables (ρ = 0.186, p = 0.47).
For the CWB rain gauges, we tested for correlations between the landslide area change and the cumulative rainfall, duration and mean intensity of the rainfall events during 15 typhoon events for the stations Meishan and Fuxing, and for 14 typhoon events for the Xiaoguanshan station. Figure 7 shows the rainfall parameters per CWB station plotted against the landslide area change and the resulting Spearman's rank correlation coefficients. Moderate positive correlations were found at a 90% confidence level between the landslide area change and the duration of the heavy rainfall event for stations Fuxing and Xiaoguanshan. This means, there is less than 10% chance that the found relationship was due to chance. Cumulative rainfall and mean intensity did not present strong significant correlations with the landslide area change. For the Meishan station, no significant correlation was found, which may be explained by the rain gauge location within a tributary valley.

Discussion
Semi-automated techniques can limit the subjectivity in landslide mapping and can contribute to improving the reproducibility of landslide maps [78]. OBIA is such a technique and provides a set

Discussion
Semi-automated techniques can limit the subjectivity in landslide mapping and can contribute to improving the reproducibility of landslide maps [78]. OBIA is such a technique and provides a set of suitable tools to semi-automatically map the evolution of landslides with time series of satellite images. The developed OBIA workflow was designed to be transferable across images, whereby only minor modifications for each Landsat sensor were necessary. This reduces the analysis time and increases the transferability of the approach. We used spectral indices for the landslide classification supported by ancillary DEM data to avoid the classification of specific false positives. While a DEM was available for only one point in time for our analysis, using DEMs acquired after each triggering event would increase the classification accuracy and would also allow a volume change estimation. In practice, however, multi-temporal DEM data are rarely available. The semi-automated mapping led to reasonable results. However, the determined accuracy values need to be considered with care, since any reference data created by manual mapping includes a certain degree of uncertainty and subjectivity [45].
The Landsat archive offers multispectral imagery since the 1980s suitable to identify recurrent changes in the area of large landslides such as the Butangbunasi landslide. Even if the spatial resolution of 30 m does not allow to identify very small changes, trends over time can be well depicted. However, the exact timing of landslide extension/reactivation following typhoons or tropical storms remains difficult. We employed the first cloud-free Landsat image acquired after such an event, but in some cases, the time span between a rainfall event and image acquisition date was up to several months or even longer. The accuracy of the OBIA mapping probably also depends on the time elapsed between a landslide triggering event and the acquisition of the next satellite image [37]. A short delay would allow deriving more detailed information about the landslide reactivating, the revegetation time and the formation of landslide-dammed lakes. The time span of over two years between the last image used and the last identified typhoon is probably the reason for the slightly decreasing trend in landslide area since initial revegetation happens quite fast. Several studies investigated the vegetation recovery after the occurrence of landslides. For example, Lin et al. [79] estimated a vegetation recovery rate of approximately 60% two years after landslides and Chou et al. [80] found a vegetation recovery rate of approximately 90% six years after landslides in central Taiwan. The new generations of satellites, for example, freely available data such as Sentinel-2 or EO data from commercial data providers, already provide a higher temporal and spatial resolution for recent years. This offers great opportunities for improving studies similar to the presented one in the future when longer time series of very high resolution (VHR) images will be available.
Rainfall-triggered landslides are particularly frequent in areas heavily affected by typhoon events, such as Taiwan. Freely available rainfall data are essential to improve early-warning models, as new technologies and research opportunities emerge. In this study, we focused on two free and open rainfall data sources. The CHIRPS data are globally available since 1981, however, its coarse spatial and temporal resolution limits its usage for local studies. In addition, the CHIRPS data under/overestimate daily precipitation. In our case, rain gauge data available from the CWB typhoon database were a more suitable alternative. Our results indicate that the duration of the heavy rainfall event is the main parameter linked to the landslide area change, while cumulative rainfall and mean intensity did not show significant correlations with the extension of the Butangbunasi landslide. Further analyses should be performed to find a direct causation between the rainfall event duration and the landslide area change. However, the freely available CWB data are limited to those hours when the typhoon passed Taiwan, and hence a complete historical rainfall database cannot be analyzed together with these extreme events. Several parameters could be computed from such a database, as indicated by Guzzetti et al. [81], which could lead to a more robust evaluation of the relationship between rain events and landslide evolution.
Better knowledge about the reactivation of large landslides and the recurring impact on downstream areas is of high importance for disaster mitigation. Even if our results did not indicate a direct relationship between the extension of the Butangbunasi landslide and the strength of the typhoon event, it became evident that also comparatively small typhoons or tropical storms could cause landslide reactivation. This is, for example, of high relevance for implementing early warning measures. The repeated sediment delivery after rainfall events frequently impacts the rivers system, which can result in the formation of landslide-dammed lakes and debris flows, and eventually poses a risk for people, settlements and infrastructure downstream. Major efforts are taken to maintain the transportation infrastructure in the study area and to avoid the repeated damming of the river (Figure 8). In particular, the Southern Cross-Island Highway, which is a popular tourist route that crosses the CMR and that provides a connection to the remote Yushan National Park, was severely affected by typhoon Morakot and following events, and the associated debris and sediment from the Butangbunasi landslide [29,30]. Information about the evolution of the Butangbunasi landslide is thus also important for planning and implementing maintenance and reconstruction activities of roads and bridges. In this study, we focused on one large landslide. By studying larger areas and relating spatiotemporal landslide hotspots to rainfall events brought by typhoons or tropical storms, more robust correlations between landslide extensions and triggering events might be found. Further research should also emphasize on the combination of OBIA and machine learning approaches for automated landslide time series analysis.

Conclusions
Large rainfall-induced landslides are among the most dangerous natural hazards in Taiwan, putting people and infrastructure at risk. Thus, better knowledge about the evolution of large landslides, their triggering factors and their potential to initiate cascading hazards is important in several respects.
Often insufficient information exists on landslide occurrence and reactivation intervals. Findings from the analysis of time series of satellite imagery, as provided for example by the Landsat missions, can provide useful information for supporting hazard mitigation and spatial planning. At the same time, the constantly increasing amount of satellite imagery at higher spatial and temporal resolutions In this study, we focused on one large landslide. By studying larger areas and relating spatio-temporal landslide hotspots to rainfall events brought by typhoons or tropical storms, more robust correlations between landslide extensions and triggering events might be found. Further research should also emphasize on the combination of OBIA and machine learning approaches for automated landslide time series analysis.

Conclusions
Large rainfall-induced landslides are among the most dangerous natural hazards in Taiwan, putting people and infrastructure at risk. Thus, better knowledge about the evolution of large landslides, their triggering factors and their potential to initiate cascading hazards is important in several respects.
Often insufficient information exists on landslide occurrence and reactivation intervals. Findings from the analysis of time series of satellite imagery, as provided for example by the Landsat missions, can provide useful information for supporting hazard mitigation and spatial planning. At the same time, the constantly increasing amount of satellite imagery at higher spatial and temporal resolutions implies the need for efficient and robust landslide (change) mapping methods. OBIA provides a suitable methodological framework for addressing these challenges. In this study, we semi-automatically mapped the evolution of the Butangbunasi landslide using Landsat time series data. The OBIA mapping results showed that the most significant landslide extension happened after typhoon Morakot in 2009. Freely available rainfall data were analyzed to find potential relationships between the reactivation and evolution of the landslide area and heavy rainfall during typhoon events. Our results indicate that the duration of the heavy rainfall event is the main parameter linked to the landslide area change, whereas daily precipitation, cumulative rainfall and mean intensity did not present strong significant correlations.
While landslides and associated hazards are a significant problem under present-day climate regimes, it is likely that climate change will lead to more frequent and extreme landslide-triggering events such as typhoons and tropical storms in Taiwan. Consequently, even more landslides may occur in the future. With this in mind, the relevance of studies that investigate and analyze the evolution of large landslides with respect to triggering events becomes even more important. Respective results can serve as input for hazard and risk analysis and the implementation of prevention and mitigation measures. Funding: This research has been supported by the Austrian Academy of Sciences (ÖAW) through the project RiCoLa (Detection and analysis of landslide-induced river course changes and lake formation).