Direct Measurements of Bedrock Incision Rates on the Surface of a Large Dip-slope Landslide by Multi-Period Airborne Laser Scanning DEMs

This study uses three periods of airborne laser scanning (ALS) digital elevation model (DEM) data to analyze the short-term erosional features of the Tsaoling landslide triggered by the 1999 Chi-Chi earthquake in Taiwan. Two methods for calculating the bedrock incision rate, the equal-interval cross section selection method and the continuous swath profiles selection method, were used in the study after nearly ten years of gully incision following the earthquake-triggered dip-slope landslide. Multi-temporal gully incision rates were obtained using the continuous swath profiles selection method, which is considered a practical and convenient approach in terrain change studies. After error estimation and comparison of the multi-period ALS DEMs, the terrain change in different periods can be directly calculated, reducing time-consuming fieldwork such as installation of erosion pins and measurement of topographic cross sections on site. The gully bedrock incision rate calculated by the three periods of ALS DEMs on the surface of the Tsaoling landslide ranged from 0.23 m/year to 3.98 m/year. The local gully incision rate in the lower part of the landslide surface was found to be remarkably faster than that of the other regions, suggesting that the fast incision of the toe area possibly contributes to the occurrence of repeated landslides in the Tsaoling area.


Overview of Tsaoling Landslide
The Tsaoling area geologically belongs to a dip-slope terrain. Multiple historical landslide events occurred in this area, which were triggered by earthquakes, typhoons, and torrential rains in 1862, 1898, 1941, 1942, 1951, 1979, and 1999 [2,38-43,46]. Among these, the landslides in 1862, 1941, and 1999 were triggered by earthquakes. The magnitudes for the two earliest earthquakes were 6 to 7 and 7.1, respectively. On 21 September 1999, the Chi-Chi earthquake, with a ML magnitude of about 7.3 and a seismic moment magnitude Mw of 7.6, occurred in the area while additional slides were triggered by the continuous torrential rains that lasted for two to five days, with precipitation up to 776 mm [2]. The volume of the rock debris accumulated from the previous landslide was about 0.026-0.126 km 3 [2,42], which formed a landslide dam with a height of 70-170 m and blocked the Qingshui River. This brought disasters to the Tsaoling area adjacent to the landslide site. The downstream areas then suffered from sudden failures of the landslide dam. Several disasters in the past have also become important case studies for large landslides [41,[47][48][49][50][51]. Figure 3 shows the aerial view of the Tsaoling area taken in 2012. Before the incident of the Chi-Chi earthquake in 1999, there was already an obvious V-shaped scarplet in the landslide crown area, a precursor of sliding before the landslide event [43]. When the Chi-Chi earthquake occurred, a large-scale rock avalanche was triggered, forming a large-scale bare slope. Afterwards, there were incisions that initiated water gullies on the sliding surface after being flushed by the surface runoff. Up until 2012, there was obvious gully development caused by continual incisions on the slope. The Cholan Formation and the Chinshui Shale, as well as their stratigraphic boundary, could be seen from the terrain height difference in the Chunqiu Cliff. Figure 3 shows that the terrain of the Tsaoling landslide area is an obvious dip slope, with the upper part being the Cholan Formation and the lower part being the Chinshui Shale.
The Chunqiu Cliff line moves to the north over time because of the constant collapse of the sandstone in the Cholan Formation. Therefore, in this study, the gullies in the sliding surface were identified and interpreted according to the multi-period ALS DEMs and aerial images from 1999 to 2011 to gain an understanding of the ground surface erosion, especially the gully incision rate on the bedrock within the study area after the Chi-Chi earthquake triggered the landslide. of the sandstone in the Cholan Formation. Therefore, in this study, the gullies in the sliding surface were identified and interpreted according to the multi-period ALS DEMs and aerial images from 1999 to 2011 to gain an understanding of the ground surface erosion, especially the gully incision rate on the bedrock within the study area after the Chi-Chi earthquake triggered the landslide.

Materials and Methods
In this study, the ALS technique was used to produce three periods of DEMs by the Council of Agriculture and the Central Geological Survey together with five periods of orthoimages for interpretation. The orthoimages in 2011 and 2012 are the aerial images taken when the ALS survey was being conducted. The orthoimages in 1998, 1999, and 2003 are the historical images recorded in

Materials and Methods
In this study, the ALS technique was used to produce three periods of DEMs by the Council of Agriculture and the Central Geological Survey together with five periods of orthoimages for interpretation. The orthoimages in 2011 and 2012 are the aerial images taken when the ALS survey was being conducted. The orthoimages in 1998, 1999, and 2003 are the historical images recorded in

Materials and Methods
In this study, the ALS technique was used to produce three periods of DEMs by the Council of Agriculture and the Central Geological Survey together with five periods of orthoimages for interpretation. The orthoimages in 2011 and 2012 are the aerial images taken when the ALS survey was being conducted. The orthoimages in 1998, 1999, and 2003 are the historical images recorded in routine aerial photography campaigns for Taiwan forests by the Aerial Survey Office, Forest Bureau (ASO).

Airborne Laser Scanning (ALS)
Airborne laser scanning (ALS) obtains optimal resolution and accuracy in regional-scale terrain topographic surveying and mapping. The data used in this study were obtained in April 2002, August 2011, and September 2012. In April 2002, the Council of Agriculture used ALS to conduct the topographic survey of the Chi-Chi earthquake area to examine the feasibility of its use in surveying and the application of the technology [52]. In August 2011, the ALS survey was implemented to obtain an island-wide DEM with 1 m resolution in Taiwan by the Central Geological Survey, Ministry of Economic Affairs (MOEA), after serious damage was induced by Typhoon Morakot. A point cloud density of 1.5 pt/m 2 was targeted in the ALS survey for the hilly and mountainous areas above 800 m. In September 2012, the ALS survey for the Tsaoling area was conducted by this research team. The survey covered the same area as the DEM data obtained in 2011. All the relevant flight parameters for these three measurements are shown in Table 1. The ALS data for the three periods were obtained through the same processes. First, the data obtained after the flights were calculated with flight strip adjustment. Second, the original point cloud data were obtained for further classification of the point cloud data. The point clouds were divided into four categories: ground points, non-ground points, water body points, and outlier points. Subsequently, the digital surface model (DSM) and the bare-earth DEM can be generated. For the ALS surveys of the 2011 and 2012 periods, the aerial photos were taken and orthorectified into orthoimages of 25 cm resolution. As the ALS data in 2011 have a wide and complete coverage in all the areas, the terrain data of that period are employed as the benchmark for the accuracy assessment.
For comparison purposes, the TWD97 ground coordinate system was adopted for the ALS data in each period [53], and the Taiwan Vertical Datum 2001 (TWVD 2001) was adopted for the vertical coordinate system [54].

Aerial Photogrammetry
Aerial photogrammetry can be used to obtain the orthoimage and digital terrain model data, which are helpful in conducting topographic analysis and geomorphological interpretation. The Aerial Survey Office of the Forest Bureau (ASO) regularly takes aerial photographs across Taiwan. The images from previous years can be used for terrain change analysis. The historical orthoimages in the areas adjacent the Tsaoling landslide were selected in this work, including images taken on 21 July 1998, 24 September 1999, 21 June 2002, and 1 July 2012 as shown in Figure 4. In addition, the two periods of aerial photographs obtained during the ALS surveys were also used in this study. Therefore, orthoimages from six periods with a pixel resolution of at least 25 cm were used in the topographic feature interpretation in the Tsaoling landslide area. Additionally, in order to understand the condition of the slope terrain immediately after the Chi-Chi earthquake, the aerial images of 1999 were used for generating the DEM data. The six aerial photographs together with the ALS DEMs were used to conduct the topographic feature interpretation in the landslide area from the perspective of shallow landslides, landslide scars, erosion gullies, and landslide dams. topographic feature interpretation in the Tsaoling landslide area. Additionally, in order to understand the condition of the slope terrain immediately after the Chi-Chi earthquake, the aerial images of 1999 were used for generating the DEM data. The six aerial photographs together with the ALS DEMs were used to conduct the topographic feature interpretation in the landslide area from the perspective of shallow landslides, landslide scars, erosion gullies, and landslide dams.

Equal-Interval Cross Section Selection Method
The terrain features in the landslide area can be characterized by slope analysis, hillshade analysis, and other visualization calculations using the ALS DEMs. Subsequently, the terrain features of the landslide surface, such as gullies and landslide cliffs, can be interpreted. The

Equal-Interval Cross Section Selection Method
The terrain features in the landslide area can be characterized by slope analysis, hillshade analysis, and other visualization calculations using the ALS DEMs. Subsequently, the terrain features of the landslide surface, such as gullies and landslide cliffs, can be interpreted. The characterization of erosion gullies shown in Figure 3 is based on the ALS DEM data in 2012. The five gullies (A-E), which are located in the Cholan Formation, are interpreted according to the different lithological characteristics on the landslide surface. Gullies A, B, and C are continuously eroding the Chinshui Shale to the Qingshui River channel, and gullies D and E merge into Gully C. Using GIS spatial analysis and other tools, the gully cross sections are set at 20 m intervals. In each section, the high point and low point of the gully are determined to obtain the depth of the gully in each period. Moreover, the depth change of the gully incision is also obtained through the different periods of the ALS DEMs. In Figure 5, PH is the highest point of the gully section, while PL is the lowest point of the gully section. The depth of the gully, D, can be obtained by D = PH − PL. The gully incision rate, R, from t1 to t2 is  Figure 3 is based on the ALS DEM data in 2012. The five gullies (A-E), which are located in the Cholan Formation, are interpreted according to the different lithological characteristics on the landslide surface. Gullies A, B, and C are continuously eroding the Chinshui Shale to the Qingshui River channel, and gullies D and E merge into Gully C. Using GIS spatial analysis and other tools, the gully cross sections are set at 20 m intervals. In each section, the high point and low point of the gully are determined to obtain the depth of the gully in each period. Moreover, the depth change of the gully incision is also obtained through the different periods of the ALS DEMs. In Figure 5, PH is the highest point of the gully section, while PL is the lowest point of the gully section. The depth of the gully, D, can be obtained by D = PH − PL. The gully incision rate, R, from t1 to t2 is R = (Dt2 − Dt1)/(t2 − t1).

Continuous Swath Profiles Selection Method
The traditional single cross section method uses a single straight line to extract the topographic elevations. The swath profiles method [55,56], in contrast, implements continuous calculations for the gully terrain data. In this way, the high and low points of a continuous series of cross sections for the whole gully can be obtained from the source to the point where it flows into the main stream. The gully depth and incision rate can be calculated accordingly.
In general, the swath profiles method uses a rectangle to select the area to be drawn. The DEM data within the rectangle is rotated so that one side of the cross section to be drawn is parallel to the south-north direction. After calculation, the maximum value, minimum value, mean, and standard deviation of all the height values can be obtained in each line. The values obtained in each line are used to draw the profile, so the profile lines of the maximum value, minimum value, mean, and standard deviation within the rectangle [55,56] are obtained. This study selected irregular areas along the gullies to draw the swath profiles [11]. As described above, the DEM data within the irregular areas were rotated to make one side of the profiles be parallel to the south-north direction. Then, the highest and lowest values of the profile were determined. Subsequently, the elevation difference should represent the gully depth for each period and the difference in gully depth between two periods represents the gully incision depth. The gully incision rate can be obtained by dividing the gully incision depth by the time interval between the selected periods.

Image Interpretation Results
Orthoimages with at least 25 cm resolution from six different periods were employed to interpret the terrain changes, as shown in Figure 4. Figure 4a shows the image of the Tsaoling area before the Chi-Chi earthquake, in which a small settlement and several roads can be seen in the forefront of the landslide area. After the earthquake, the settlement slid to the opposite bank of the

Continuous Swath Profiles Selection Method
The traditional single cross section method uses a single straight line to extract the topographic elevations. The swath profiles method [55,56], in contrast, implements continuous calculations for the gully terrain data. In this way, the high and low points of a continuous series of cross sections for the whole gully can be obtained from the source to the point where it flows into the main stream. The gully depth and incision rate can be calculated accordingly.
In general, the swath profiles method uses a rectangle to select the area to be drawn. The DEM data within the rectangle is rotated so that one side of the cross section to be drawn is parallel to the south-north direction. After calculation, the maximum value, minimum value, mean, and standard deviation of all the height values can be obtained in each line. The values obtained in each line are used to draw the profile, so the profile lines of the maximum value, minimum value, mean, and standard deviation within the rectangle [55,56] are obtained. This study selected irregular areas along the gullies to draw the swath profiles [11]. As described above, the DEM data within the irregular areas were rotated to make one side of the profiles be parallel to the south-north direction. Then, the highest and lowest values of the profile were determined. Subsequently, the elevation difference should represent the gully depth for each period and the difference in gully depth between two periods represents the gully incision depth. The gully incision rate can be obtained by dividing the gully incision depth by the time interval between the selected periods.

Image Interpretation Results
Orthoimages with at least 25 cm resolution from six different periods were employed to interpret the terrain changes, as shown in Figure 4. Figure 4a shows the image of the Tsaoling area before the Chi-Chi earthquake, in which a small settlement and several roads can be seen in the forefront of the landslide area. After the earthquake, the settlement slid to the opposite bank of the Qingshui River at the foot of the slope and was covered by soil and rocks, causing the death of 29 residents. Figure 4b shows the image several days after the earthquake. After the rock slide, a flat sliding surface was revealed and a dip-slope terrain was formed. In the middle of the slope, the Cholan Formation can be seen in the top half and the Chinshui Shale can be seen in the lower half, with the Chunqiu Cliff as their boundary. The Qingshui River was nearly filled by the large amount of soil and rocks sliding to the river course after the landslide, initiating the formation of the landslide dam. Figure 4c shows the image from 2002. The landslide dam, formed by soil and rocks, still existed two years after the earthquake. Most areas on the slope were still bare ground, without full coverage of plants. Many gullies were developing along the slope. On the northeast side of the slope, there were several intercepting ditches parallel to the slope after a manual restoration project, as well as several artificial gullies going downward along the slope. Figure 4d shows the image from 2011 when the landslide dam had disappeared. Compared with Figure 4c taken in 2002, the slope was covered by plants. Several larger gullies could be identified from the aerial photo, while the remaining small gullies were mostly covered by plants. Terrain cliffs caused by the large lithological difference between the sandstone and shale, such as the Chunqiu Cliff, could be identified from the image. In Figure 4e, only the area circumscribed by the yellow dashed line was taken in July 2012. At that time, clouds covered a wide area, so only the cloud-free area was presented here to observe the location change of the Chunqiu Cliff. Figure 4f shows the image from September 2012. Compared with that of 2011, shown in Figure 4b, a further collapse of the Chunqiu Cliff, a position change of the slope, and sand and rocks after a slide can be observed. However, there was no significant difference in the other geomorphological characteristics.

DEM Topographic Interpretation
The ALS technique was used to produce three periods of DEM data for interpretation. In order to understand the slope terrain change shortly after the earthquake, aerial photogrammetry was also used to produce the DEM data shortly after the Chi-Chi earthquake in 1999. The four periods of DEM data are shown in Figure 6 in the form of a color-shaded relief map. This kind of terrain shadow map can be used to show the geomorphological characteristics of slopes. Figure 6a shows the DEM image shortly after the Chi-Chi earthquake in 1999. In the figure, the zone of depletion triggered by the earthquake and the zone of accumulation can be clearly seen. The sliding surface of the landslide area was an obvious planar slope, with no gully development. Sand and rocks from the landslide filled the channel of the Qingshui River, resulting in a landslide dam. Figure 6b shows the DEM image from 2002 after the landslide. A large amount of soil and rocks blocked the river course. In some parts of the zone of accumulation, new river channels were formed by the effects of erosion. The landslide dam existed for two years after the earthquake. Gullies gradually developed on the slope during this period. The grey areas represent non-measured areas or areas with no effective data. Figure 6c shows the DEM from 2011. Affected by erosion and incision, the rock blockage was eroded and cleared from the main course of the Qingshui River, and the landslide dam no longer existed. Significant gully development can be observed in this period. Figure 6d is the DEM measured in 2012, from which the topographical conditions of the Tsaoling landslide area 10 years after the earthquake can be observed, including the gully development of the collapsing slope, the Chunqiu Cliff, and the Qingshui River channel. In Figure 6a-d, the geomorphological evolution of the Tsaoling landslide area since the Chi-Chi earthquake in 1999 can be observed. This verifies the feasibility of using multiple-period DEM data to identify and interpret gullies, river course incisions, and the development and degradation of landslide dams.

Calculations of Incision Rates
To calculate the gully incision rate, five gullies on the slope of the Tsaoling landslide area, which can be interpreted by the orthoimages and the DEM data, were selected. After the interpretation, the three periods of ALS DEMs were used to calculate the incision rate using the 20 m equal-interval cross section method and the continuous swath profiles selection method.

The Equal-Interval Cross Section Selection Method
For each selected gully, the cross sections were evaluated every 20 m, which resulted in 42-78 sets of cross section data, depending on the length of each gully ( Table 2). Calculations were performed on each cross section, and the results are shown in Table 2

Calculations of Incision Rates
To calculate the gully incision rate, five gullies on the slope of the Tsaoling landslide area, which can be interpreted by the orthoimages and the DEM data, were selected. After the interpretation, the three periods of ALS DEMs were used to calculate the incision rate using the 20 m equal-interval cross section method and the continuous swath profiles selection method.

The Equal-Interval Cross Section Selection Method
For each selected gully, the cross sections were evaluated every 20 m, which resulted in 42-78 sets of cross section data, depending on the length of each gully ( Table 2). Calculations were performed on each cross section, and the results are shown in Table 2

The Continuous Swath Profiles Selection Method
The continuous calculation was also conducted using the swath profiles method for the selected gullies. Figure 7 shows the swath profile results for the five gullies in each DEM period, as well as the elevation difference between the highest and lowest values or the gully depth in each period. The grey dashed line in the figure shows the boundary between the Cholan Formation featuring the thick sandstone layer and the Chinshui Shale. Due to significant differences in lithological properties between the two areas, the Chunqiu Cliff was formed at the boundary with a large elevation difference. As shown in Figure 7, the sliding of the Chunqiu Cliff sandstone caused the cliff line to move gradually to the north. The sliding area was also included in the swath profiles calculation. However, since the elevation difference within the sliding area was not the result of gully erosion, the elevation difference values of the collapsed area of the Chunqiu Cliff were excluded in the calculation for Gullies A, B, and C. These results are shown in Table 3

Three-Period ALS DEM Data Comparison and Error Estimation
When comparing DEM data from different periods, the difference in data accuracy from the different periods affects the results of the comparison. Therefore, in this work, error estimation is implemented based on the statistics between the two periods of data [11]. Forty-one points located outside the Tsaoling landslide affected areas were adopted as the ground control points (GCPs), as they were easy to identify, flat, not covered by vegetation, and characteristic corners of buildings. The GCPs were selected based on the multi-period ALS data and orthoimages, where the selected points show no obvious signs of movement during the periods. These points were used to carry out the error estimation for the elevations from the three-period DEM data.  Table 5.
The comparative analysis of multiple-period DEM data has been widely used in studies on terrain change, such as geomorphological evolution, geological hazards, and river channel sediment transport [57][58][59]. The principle of this approach is to calculate the difference between the later elevation and the corresponding earlier elevation. Three periods of ALS DEMs constructed for the Tsaoling area were used in this study. Figure 8 shows the terrain change results of the Tsaoling area in two different periods, i.e., the nine-year long-term DEM data during 2002-2011 and the one-year short-term DEM data during 2011-2012. Figure 8a,b shows the elevation change diagrams for the two periods, respectively, where the color map represents the topographic change. The cool colors represent the decrease in terrain elevation, indicative of the effects of surface erosion. The warm colors represent the increase in terrain elevation, indicative of the effects of debris accumulation. The increasing and decreasing values were then enhanced by a continuous numerical method, the results of which are shown in Figure 8c,e to demonstrate changes during 2002-2011. It can be seen that, in the long term, the erosion effect was mainly exerted on the adjacent area of the Qingshui River, the Chunqiu Cliff, the shallow collapsing area of the upper edge of the slope, and the slope gullies. Accumulation occurred under the Chunqiu Cliff, as well as the location under the shallow sliding site on the upper edge of the slope. Figure 8d,f shows the changes during 2011-2012. Both erosion and accumulation were densely concentrated in the areas adjacent to the Qingshui River and the Chunqiu Cliff. It can be concluded from Figure 8 that the volume and range of erosion and accumulation are greater considering long-term geomorphological effects than when considering short-term geomorphological effects.  Figure 8d,f shows the changes during 2011-2012. Both erosion and accumulation were densely concentrated in the areas adjacent to the Qingshui River and the Chunqiu Cliff. It can be concluded from Figure 8 that the volume and range of erosion and accumulation are greater considering long-term geomorphological effects than when considering short-term geomorphological effects. Table 5. Elevation error estimations between every two periods of ALS DEMs using the same ground control points.

Three-Period Data Comparison and Error Evaluation
The high-resolution ALS DEMs generated from three different periods were employed in this study to identify the gully distribution on the sliding surface of the Tsaoling landslide area. With the cross section locations determined based on the lithology on the sliding surface, incision depths were calculated. The three periods of DEM data produce different errors related to different production times and flight measurement planning. When comparing the DEMs, one must have sufficient information about data quality; especially the flight mission plan for the particular measurement area, including flight height, flight speed, heading, and route overlap rate, which will all affect the DEM resolution and error. In general, the elevation error of ALS is about 0.1-0.3 m. In planar areas with no vegetation, this can be below 0.1 m [11,52]. Under general treatment procedures, the DEM error will be affected by terrain slope, plant type, and vegetation density. The DEM error for bare planar ground should be better than that for bare slope areas. Similarly, areas with sparse vegetation should have superior results in terms of error than areas with dense vegetation. The target of this study is mainly to investigate gully bedrock incisions, which can be considered as areas with sparse vegetation. In Table 5, the relative errors of the data in each period calculated with the same ground control points are compared. Using the 2011 DEM as the benchmark, the root mean square error of 2002 DEM and 2011 DEM is 0.31 m, and the root mean square error of 2011 DEM and 2012 DEM is 0.17 m. The ALS DEMs can be used to directly measure the terrain changes by comparing the incision depths obtained during different periods for each gully. Due to flight and equipment parameters, the resolution and accuracy of the 2002 DEM are not as good as in the two later periods. However, the longer observation time and more extensive terrain change during this period allowed the 2002 DEM to be appropriate for gully evolution calculations. Previous studies [30,35] have showed that the ALS DEMs can be applied to gully measurements as long as the point cloud density is sufficiently high. The present study also showed that the ALS DEMs could be applied to the direct observation and measurement of the erosion gullies, greatly reducing the amount of time-consuming fieldwork, such as configuring erosion pins or cross section measurements. The comparison of multiple periods of DEMs can directly illustrate the terrain evolution in different stages within the whole area. As shown in Figure 8, the extent of erosion or accumulation can be clearly interpreted by a certain time.

Gully Incision Rate
In this study, the gully bedrock incision rate was calculated using two different methods, i.e., the equal-interval cross section selection method and continuous swath profiles selection method ( Table 4). The results show that the incision rate of the bedrock within the Chinshui Shale was significantly higher than that of the Cholan Formation. For the five selected gullies in the study area, their water accumulation areas, slopes, and climatic conditions are similar, which suggests that the erosion resistance of the Chinshui Shale is weaker than that of the Cholan Formation. The resistance to erosion is closely related to the lithological strength. The overall rock strength of the Cholan Formation, which is composed of interlayers of thick sandstone and shale, is greater than that of the Chinshui Shale that consists mainly of shale. The finding that the incision rate of gully bedrock is affected by lithological composition is also in agreement with the results of river channel incision rate results from previous studies [16,17,60].
A field investigation was also performed to examine the lithological composition of the sidewall and bedrock of the gullies, as shown in Figure 9. The gully sidewalls within the Cholan Formation mainly consist of thin interlayers of sandy shale, while the bedrock primarily features thick sandstone. The thick sandstone of the bottom has significantly larger lithological strength than the sidewall. A wide and shallow U-shaped gully is currently formed within the Cholan Formation (Figure 9b,c), possibly due to the fact that erosion occurs faster on the sidewall than the gully bottom. Within the Chinshui Shale, the gullies are exposed as shale, forming a wide and deep V-shape due to fast incision ( Figure 9d). The differences in morphology and incision rate among the studied gullies should be associated with various mechanisms of riverbed incision and erosion, including plucking, abrasion, cavitation, and solution erosion [61]. In addition, these differences might also relate to the location and interval distance of the discontinuous cross sections, as well as the regional geological structure, which should be examined in future studies. Nevertheless, the ALS DEM is a useful analytical tool in erosion and incision research.  [61]. In addition, these differences might also relate to the location and interval distance of the discontinuous cross sections, as well as the regional geological structure, which should be examined in future studies. Nevertheless, the ALS DEM is a useful analytical tool in erosion and incision research.

Continuous Sampling and Intermittent Sampling
Most studies of the erosion rate of river channel incisions in the past employed direct cross section measurements in the river channel [15,16]. The advancements in remote sensing technology with improved DEM resolution offer efficient alternative methods for river channel cross section selection [7,30,34,35]. For the equal-interval method, the interval distance needs to be carefully determined based on the scale of the slope gully. If the distance is too large, the volume of the cross section data would be too small, and it would be doubtful whether the gully depth could be accurately reflected. If the distance is too small, the volume of data should be sufficient, but the extraction will be time-consuming, and the effectiveness might still be questionable in some cases. In this study, based on the slope scale of the study area, 20 m was set as the interval distance for this method so that the number of sections is statistically significant.
The application of the swath profiles method has mostly appeared in discussions associated with tectonic structures or landslides on a regional scale [11,62]. By using this method, the maximum value and minimum value of each cross section can be directly calculated in a continuous manner from the source of the gully based on the DEM resolution. With the obtained minimum values, the

Continuous Sampling and Intermittent Sampling
Most studies of the erosion rate of river channel incisions in the past employed direct cross section measurements in the river channel [15,16]. The advancements in remote sensing technology with improved DEM resolution offer efficient alternative methods for river channel cross section selection [7,30,34,35]. For the equal-interval method, the interval distance needs to be carefully determined based on the scale of the slope gully. If the distance is too large, the volume of the cross section data would be too small, and it would be doubtful whether the gully depth could be accurately reflected. If the distance is too small, the volume of data should be sufficient, but the extraction will be time-consuming, and the effectiveness might still be questionable in some cases. In this study, based on the slope scale of the study area, 20 m was set as the interval distance for this method so that the number of sections is statistically significant.
The application of the swath profiles method has mostly appeared in discussions associated with tectonic structures or landslides on a regional scale [11,62]. By using this method, the maximum value and minimum value of each cross section can be directly calculated in a continuous manner from the source of the gully based on the DEM resolution. With the obtained minimum values, the longitudinal section of the gully considered consistent with the actual topography can be drawn. The elevation difference between the highest and lowest values can be used to indicate the incision depth of the gully, as shown in Figure 7. This method can quickly obtain the data of the longitudinal cross sections of the gully, as well as the corresponding elevation change. In addition, it also has the advantage of automatically excluding a small number of special topographic variations affecting the overall calculation of incision depth. For example, between the Cholan Formation and the Chinshui Shale in the investigated slope, the Chunqiu Cliff is formed because of the remarkable lithological difference. The shale erosion below the Chunqiu Cliff results in a continuous collapse of the sandstone cliff, and the elevation difference observed should not be considered in the gully erosion depth calculation. If the erosion rates calculated by the two methods (Table 4) are compared, the difference between the two methods within the Cholan Formation is small, but large discrepancies can be seen for the Chinshui Shale. The reason for this probably lies in the ability of the continuous swath method to reflect accurately the changes in both the accumulation and erosion of the gully. In the area of the Chinshui Shale, Gullies A, B, and C are all affected by debris accumulation from the Chunqiu Cliff (Figure 8c,d).
The equal-interval method may have partially ignored the influence of debris accumulation at certain locations, leading to much smaller incision rates in the Chinshui Shale than those calculated using the continuous swath method. Therefore, it is reasonable and convenient to obtain complete gully evolution data using the continuous swath method before conducting targeted data analysis according to the actual geomorphological features.

The Effects of Typhoon Soula
Three periods of ALS DEMs were used in this study to discuss gully erosion evolution on the landslide slope. The interval of the first two periods of data is nine years, and the interval of the last two periods is about one year. Compared to previous studies, the present study should be considered as a short-term topographic erosion investigation. Figure 10 shows the data of the last two periods from August 2011 and September 2012, where the effects of erosion can be observed from the retreat of the Chunqiu Cliff outline. According to the rainfall data in Figure 2, during this three-day period, only Typhoon Soula caused torrential rain with more than 700 mm of rainfall. Figure 10a,c shows aerial photos taken in the last two ALS DEM measurements. Additionally, Figure 10b shows an orthoimage (selected from the ASO historical aerial photograph library) taken in July 2012, just prior to Typhoon Soula. According to the DEM-based terrain shadow map illustrating elevation changes during 2011-2012 (Figure 10d), along with the superimposed Chunqiu Cliff position change from Figure 10a-c at three time points, no major change in the cliff position can be observed before Typhoon Soula. Therefore, it can be inferred that the topographic elevation change during 2011-2012 was mainly due to the torrential rain brought by Typhoon Soula on 30 July 2012. This also suggests that, when using multiple-period ALS DEMs, the determination of the time cycle for data acquisition could be based on topographic changes induced by large-scale external forces, such as typhoons and earthquakes.

Gully Bedrock Incision Rate
The gully bedrock incision rate reflects terrain features, such as lithology, slope, discontinuous surface distribution, tectonic characteristics, river width, and basin size, as well as the impacts of regional tectonic structures and climatic conditions [63][64][65]. The long-term incision rate is often obtained using the sediment dating or river terrace dating methods. These data, covering thousands of years or more than ten thousand years, usually result in incision rates on the order of mm/year. The short-term incision rate of a few years or decades can be obtained using erosion pins or cross section measurements, with a typical incision rate ranging from a few to tens of mm/year. The short-term rate is much larger than the long-term rate [15][16][17]. In this study, the gully bedrock incision rates range from as small as 0.23 to as large as 3.98 m/year. With an incision rate several or dozens of times larger than that of a typical river, these gullies develop so rapidly that the pace can be analogous to direct cutting on the slope rock. In addition, with shale as the main lithological component, the toe area of the slope has been eroded by the mainstream of the Qingshui River for decades, and the slope has become a geologically unstable area. As a result, when earthquakes, typhoons, or torrential rains occur, landslides would be triggered again in this area.

Conclusions
Understanding the fundamental processes of landslides will help manage and mitigate landslide-related natural hazards. Recent advancements in ALS technology have provided a great opportunity for characterizing landslides through developing improved methods in multi-period

Gully Bedrock Incision Rate
The gully bedrock incision rate reflects terrain features, such as lithology, slope, discontinuous surface distribution, tectonic characteristics, river width, and basin size, as well as the impacts of regional tectonic structures and climatic conditions [63][64][65]. The long-term incision rate is often obtained using the sediment dating or river terrace dating methods. These data, covering thousands of years or more than ten thousand years, usually result in incision rates on the order of mm/year. The short-term incision rate of a few years or decades can be obtained using erosion pins or cross section measurements, with a typical incision rate ranging from a few to tens of mm/year. The short-term rate is much larger than the long-term rate [15][16][17]. In this study, the gully bedrock incision rates range from as small as 0.23 to as large as 3.98 m/year. With an incision rate several or dozens of times larger than that of a typical river, these gullies develop so rapidly that the pace can be analogous to direct cutting on the slope rock. In addition, with shale as the main lithological component, the toe area of the slope has been eroded by the mainstream of the Qingshui River for decades, and the slope has become a geologically unstable area. As a result, when earthquakes, typhoons, or torrential rains occur, landslides would be triggered again in this area.

Conclusions
Understanding the fundamental processes of landslides will help manage and mitigate landslide-related natural hazards. Recent advancements in ALS technology have provided a great opportunity for characterizing landslides through developing improved methods in multi-period DEM differencing and geomorphometry. This study used the ALS DEMs from three different periods to interpret the gullies formed on the slope after the Tsaoling landslide. The elevation change was used to calculate the gully bedrock incision rate. By estimating errors associated with the gully incision depth calculation, it was verified that the gully erosion related terrain change could be efficiently determined using ALS DEMs, greatly reducing the amount of time-consuming fieldwork, such as installing erosion pins and making direct cross section measurements. In this study, the incision rate was calculated using two different methods for cross section extraction, the 20 m equal-interval cross section selection method and the continuous swath profiles selection method. The two methods gave similar results for the Cholan Formation, but different results for the Chinshui Shale. Based on the findings, we suggest that the continuous swath method should be used as it can reflect the elevation change of the gullies on a continuous basis. The gully bedrock incision rates in the Chinshui Shale were found to be significantly greater than in the Cholan Formation, indicating that the erosion resistance of the former is weaker than that of the latter. In this study, the gully bedrock incision rates ranged between 0.23 and 3.98 m/year, remarkably higher than the typical results from the previous studies, suggesting that the high erosion rate is also one of the primary factors leading to repeated landslides in this region. By comparing the DEM data, aerial photos, and precipitation records of this area, the effects of erosion could be observed from the retreat of the Chunqiu Cliff outline during August 2011 to September 2012. It was inferred that the change in the topographic elevation during 2011-2012 was mainly due to the torrential rain brought by Typhoon Soula, which occurred on 30 July 2012. Thus, it can be said that if the sampled period is too short or there are no external events during this period, the degree of terrain change is likely to be in the range of the ALS DEM error, and it may not reflect the actual changes in the terrain.