Badland Erosion and Its Morphometric Features in the Tropical Monsoon Area

Climatically driven processes are important controls on the Earth’s surface and on interactions between the hydrological cycle and erosion in drainage basins. As a result, landscape forms such as hillslope topography can be used as an archive to reconstruct historical climatic conditions. Recent progress in the Structure-from-Motion (SfM) photogrammetric technique allows for the construction of high-resolution, low-cost topography data using remote-controlled unmanned aerial vehicle (UAV) surveys. Here, we present the climatic effects on the hillslope erosion rate that can be obtained from the drainage frequency of hillslopes. We quantify the centimeter-scale accuracy of surveys across 72 badland hillslopes in SE Taiwan, which is a tropical monsoon area with an annual precipitation of over 2 m. Our observations indicate that climatic erosion results in a higher drainage frequency and the number of furrows, instead of drainage density. Additionally, the morphometric slope index (MSI) has a strong positive correlation with erosion and its rate but shows a negative correlation with drainage length and a positive correlation with inclination. This suggests that the erosion pattern is due to gravitational mass wasting instead of hydrological erosion. MSI should always be calculated relying on the normalized slope length and is less applicable to landslide-dominated erosion. We, therefore, suggest that UAV-driven digital elevation models (DEMs) are integrated into erosion mapping to aid in identifying erosion patterns. We highlight the unique opportunity for cross-climate zone comparative studies offered by badland landscapes and differential rainfall patterns, with remote sensing techniques and the morphometric slope index.


Introduction
Landscapes evolve under the influence of external drivers and can be used as an archive for reconstructing historical climatic conditions. Badlands are a common landscape formation in softer sedimentary rocks and clay-rich sediment and are distributed across climate zones, e.g., Midwestern North America, as well as the Mediterranean and Asian tropical monsoon areas. Meanwhile, badlands are especially sensitive to environmental change due to dense drainage systems and sharp slopes in barren areas. Therefore, badlands lend a unique opportunity to extend the understanding of climatic signatures in landscapes.
Nevertheless, studies comparing high-resolution digital elevation model (DEM)-based morphometric analyses of badland hillslope development in different climate zones are still limited. Meanwhile, the hydrological erosion in badlands also causes high denudation rates and corresponding geohazards globally [1]. For example, gully erosion in Te Weraroa Stream, New Zealand affected ∼6% of the drainage basin and produced 28.7 Mt of sediment within 30 years [2]. In contrast, reforestation as a protection of the land surface can reduce the area of badlands and the expected erosion rate [3,4]. In tropical monsoon areas, badlands in Taiwan caused extensive landslides, mudflow [5], and denudation rates of around 9-13 cm per year [6], resulting in the Erren river, with its drainage area mostly in the mudstone badlands, having the highest mean sediment yield in the world of over 10 5 t/km 2 /yr [7]. Therefore, simple and systematic constraints on erosion help to improve the understanding of interactions of the lithosphere and the hydrosphere in topography.
The morphometric slope index (MSI) was developed to analyze Calanchi erosion in northern Sicily, Italy. Calanchi is a specific term for badlands in Italy which has dense drainage systems and steep clayey slopes [36]. The methodology builds a relationship among morphometric features, erosion type, process and volume, in slope scale. However, the Calanchi area is mainly distributed in the Mediterranean climate. In this study, we would like to apply this method to the mudstone badland area in SE Taiwan, which is a tropical monsoon area with an annual precipitation of over 2 m. Therefore, it provides an opportunity to verify the applicability of MSI to a specific alternate climate region with spring monsoons and typhoons in summer. Furthermore, comparing the two types of badland hillslope evolutions, in Italy and Taiwan, could also extend the understanding of climatic signatures in badland areas. In this contribution, the data from UAVs are used to create high-resolution aerial images and DEMs for the calculation of MSI and morphometric features. This study aims to evaluate the development of badland hillslopes in Taiwan using MSI, and establish the relationship between MSI and other topographic features. The quantitation of slope-scale erosion volumes during typhoon events would extend the understanding of erosion patterns induced by high-intensity rainfall.

Study Area
Our study area comprises about 1 km 2 of badlands landscape in Tainan County, SW Taiwan ( Figure 1a). This area hosts a military facility with high levels of security and, therefore, has undergone little intensive development in the past. Consequently, the area exhibits a natural landscape without extensive human disturbance and is, therefore highly suitable for a landscape-development study. It has a hot and wet tropical climate with major rainfalls from the seasonal monsoon and typhoons. The rainy season is between May and October and accounts on average for 89% of the total annual rainfall (2118 mm). The peak is in August with >400 mm of monthly rainfall on average. The average monthly rainfall is typically <40 mm from November to April.
The stratigraphy of the study area is dominated by clay and sandstone, which are prone to erosion caused by rainfall. As rainfall infiltrates the mudstone, the ionic bond link is destroyed and then releases strain energy, causing slaking. Illite (30.54%) and chlorite (28.70%) are the main minerals in the mudstone [5]. There is also a large proportion of soluble salt, which is easily dissolved by rainfall to form sodium ions. The sodium ions make soil particles disperse easily and result in slaking of the mudstone. After a long dry period, the mudstone surface shrinks and then forms cracks and crusts, 8-10 cm thick on the weathered surface layer and the cracks are preferential routes for infiltration. Plus, the alternate climate makes the annual erosional cycle in the the badland hillslope.
The typical badlands landscape is composed of steep slopes, dense rills, and gullies, as well as pipes or tunnels. Several studies have focused on the slope evolution of the mudstone areas in Taiwan. Yen [40] mentioned that the homogeneous lithology controlling the gradient tends to form hillslopes of 40-50 • during the erosion process, followed by a parallel slope retreat. However, the longitudinal sections of the mudstone badlands slope are not completely linear, consisting of a steep slope in the upper part and a gentle slope in the lower part. The two steps in the slope indicate that the badlands slope evolution is simultaneously driven by gravity and water erosion. This implies that it is the external force, rather than tectonic movement, that mainly controls the form of the hillslope [40]. Yen and Chen [41] developed the landscape evolution of the mudstone badlands slope-drainage density, roughness, and local slope gradient would increase in the young state and decline the relief by climatic erosion and then eventually reach stabilization.

DEMs Created from UAV Surveys
UAV surveys were conducted between July and August of 2017. The dates of the acquisition were chosen according to the weather conditions. The two consequent typhoons (July 29-August 01), Nesat and Haitang, packed 579 mm of rainfall with the maximum rainfall intensity of 74 mm/h likely causing severe erosion in the mudstone areas. We conducted UAV surveys to collect surface near-

DEMs Created from UAV Surveys
UAV surveys were conducted between July and August of 2017. The dates of the acquisition were chosen according to the weather conditions. The two consequent typhoons (July 29-August 01), Nesat and Haitang, packed 579 mm of rainfall with the maximum rainfall intensity of 74 mm/h likely causing severe erosion in the mudstone areas. We conducted UAV surveys to collect surface near-infrared (NIR) images for DEM generation and ortho-rectification with a Canon s110 camera and an eBee Classic drone [43]. Detailed information on this equipment is provided in Table 1. In order to achieve a high survey accuracy, we installed 11 ground control points (GCPs) across the 1 km 2 surveying area ( Figure 1b) and conducted high precision GCP surveys using a Leica RX1250XC differential GPS to obtain their coordinates. The GCPs were made of iron nails with a 5 mm diameter surrounded by aerial photogrammetric targets. All GCPs were positioned on solid sites, such as paved roads or bridges, to ensure stable anchoring for the survey. The coordinates of each GCP were measured over a period of more than 15 min using the real-time kinematic (RTK)-GPS system. The average flight altitude was 325 m and the average flight time to cover the entire study area was about 30 min. The resulting ortho-rectified images have a ground sample distance (GSD) of 11 cm. The data were transformed into 3D-models of the study area using the Acute 3D software (Bentley Systems).

DEM of Difference (DoD)
DoD is the difference between two DEMs of the same area but acquired at two different times [44]. It reveals the geomorphic change and provides insight into the relationship between process and form. It can also be used to assess the simulation of numerical morphodynamic models. In this study, we constructed a DEM of each UAV survey and then obtained the DoD accordingly.
Wheaton et al. [44] and Williams [45] proposed that the DEM elevation value is likely to contain a vertical error component (or uncertainty). In DoDs, the combined errors are as follows: where δZ 1 and δZ 2 are the errors associated with Z 1 and Z 2 , respectively. We followed this method to calculate the error in our DoDs.

Parameter Extraction and Calculation of MSI
The geomorphological features of individual hillslopes were mapped manually using ArcGIS 10.5. The ortho-rectified UAV images were used to identify the hillslope divides and drainages. From the hillslope divides and drainages, we can calculate the plane tributary area (A 2D ), slope length (L), circularity ratio (R C ), drainage frequency (F), drainage density (D), and drainage length (l). The slope inclination (P) is derived from the constructed DEM and we then calculate the reconstructed surface area (Ar) accordingly. The eroded volume (V) can reflect the amplitude of the transient topography in response to climate forces. The former is defined as the maximum vertical distance between drainage divide and valley bottom; the latter can be derived by creating a TIN surface capping the hillslope and subtracting the hillslope DEM from this surface, normalized by drainage area for comparison among hillslopes of different size (V/A 2D ). The MSI is then calculated for the reconstructed slope of each hillslope [36]: Spearman's correlation test was conducted to obtain the relationship between topographic indices calculated from the DEM and other data mentioned previously, and the MSI. The significance values of Spearman's rank correlation coefficient (p-value) were set at 0.01 and 0.05, respectively.

DoD Uncertainty Estimation
The DoD of the slopes studied between July 1st and August 4th, 2017 show a mean of 0.07 m and a standard deviation of~0.01 m. The distribution of survey error in each grid is shown in Figure 2. As calculated, the uncertainty of DoD is ± 0.1 m. We, therefore, mark the change between ±0.1 m as uncertainty (error). The DoD of the slopes studied between July 1st and August 4th, 2017 show a mean of 0.07 m and a standard deviation of ~0.01 m. The distribution of survey error in each grid is shown in Figure  2. As calculated, the uncertainty of DoD is ± 0.1 m. We, therefore, mark the change between ±0.1 m as uncertainty (error).

Morphometric Slope Index (MSI)
A total of 72 hillslopes were identified and mapped from the DEM obtained from the UAV data in the study area, as shown in Figure 1c. Overall, the elevation of the slopes range from 70 m to 130 m, and 42% of the 71 hillslopes are about 90-100 m in elevation (Figure 1d). In addition, 75% of the 71 hillslopes are S-facing (including SE, S, and SW) (Figure 1e). This implies that the effect of orographic precipitation and the difference of tectonic activity is negligible. The MSI of the 72 hillslopes is shown in Table 2.

Morphometric Slope Index (MSI)
A total of 72 hillslopes were identified and mapped from the DEM obtained from the UAV data in the study area, as shown in Figure 1c. Overall, the elevation of the slopes range from 70 m to 130 m, and 42% of the 71 hillslopes are about 90-100 m in elevation (Figure 1d). In addition, 75% of the 71 hillslopes are S-facing (including SE, S, and SW) (Figure 1e). This implies that the effect of orographic precipitation and the difference of tectonic activity is negligible. The MSI of the 72 hillslopes is shown in Table 2.  The relationship between MSI and morphometric features is shown in Table 3. MSI has a strong positive correlation with total drainage length (R 2 = 0.77), total number of furrows (R 2 = 0.79), plane length (R 2 = 0.86), eroded volume (R 2 = 0.87), and mean eroded depth (R 2 = 0.70). In addition, it also has no correlation with inclination (R 2 = 0.32) and circularity ratio (R 2 = 0.41). In contrast, it has no correlation with drainage density (R 2 = −0.12) and mean drainage length (R 2 = −0.27).

Distribution of Hillslope Erosion
Over the surveyed periods, about 68% of the study area was affected by erosion. We measured 31 cm of average erosion depth during the typhoon. The two selected badlands hillslopes represent the erosion pattern ( Figure 3). Overall, Hillslope 33 exhibited the largest observed erosion rate, located on the steepest hillslopes and at sites with drainage areas under 100 m 2 . In addition, erosion and deposition occurred simultaneously in drainage areas under 100 m 2 s. Drainage areas between 200 m 2 and 300 m 2 predominantly exhibited a deposition of 0.4 m, and these are the areas where the main gully is distributed (Figure 3c). Furthermore, erosion and deposition occurred in drainage areas over 300 m 2 but resulted in a change of only ± 0.1 m. On Hillslope 59, the observed erosion was also located on the steepest hillslopes with gradients greater than 40 • , and at sites with drainage areas under 50 m 2 (Figure 3d). In contrast, deposition occurred in relatively gentle areas with slope gradients under 50 • , and the deposition rate increased with decreases in slope gradient. Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 15

The Effect of Rapid Incision on Morphology of Hillslopes in Badlands
In order to isolate the topographic effects of different climatic forcing it is essential to understand how landscapes work on badlands evolution. The calculated MSI in this study was compared to the results of Buccolini et al. [36] (Figure 4), who documented 65 morphological features of hillslopes and MSI in the central Sicily. Badlands in Taiwan have lower drainage length and drainage density than those in Italy, which may imply that the relatively small hillslope area controls the drainage network development. The badlands in the study area are caused by incision into weakly resistant mudrocks of the river terrace [46], and the formation of mudstone badlands is mainly controlled by climatic erosion, i.e., high rainfall events brought by either typhoons or summer monsoons. The strong climatic erosion also steepens the hillslope to an average hillslope inclination of 47.1°, which is 8.3° higher than that of the badlands in Sicily, Italy. Additionally, the origin of the river terrace is ascribed to strong uplift combined with a storm-induced incision in unconsolidated mudstone since 2000 YBP [46]. Such a young badlands landscape may explain why the morphology of badlands landscapes in

The Effect of Rapid Incision on Morphology of Hillslopes in Badlands
In order to isolate the topographic effects of different climatic forcing it is essential to understand how landscapes work on badlands evolution. The calculated MSI in this study was compared to the results of Buccolini et al. [36] (Figure 4), who documented 65 morphological features of hillslopes and MSI in the central Sicily. Badlands in Taiwan have lower drainage length and drainage density than those in Italy, which may imply that the relatively small hillslope area controls the drainage network development. The badlands in the study area are caused by incision into weakly resistant mudrocks of the river terrace [46], and the formation of mudstone badlands is mainly controlled by climatic erosion, i.e., high rainfall events brought by either typhoons or summer monsoons. The strong climatic erosion also steepens the hillslope to an average hillslope inclination of 47.1 • , which is 8.3 • higher than that of the badlands in Sicily, Italy. Additionally, the origin of the river terrace is ascribed to strong uplift combined with a storm-induced incision in unconsolidated mudstone since 2000 YBP [46]. Such a young badlands landscape may explain why the morphology of badlands landscapes in Taiwan is relatively short, small, and steep. However, the annual precipitation in the Taiwan badlands is three times higher than in those of Italy, which may still cause strong erosional feedback, i.e., rills and gullies. Thus, badlands in Taiwan have a higher drainage frequency and a greater number of furrows. Taiwan is relatively short, small, and steep. However, the annual precipitation in the Taiwan badlands is three times higher than in those of Italy, which may still cause strong erosional feedback, i.e., rills and gullies. Thus, badlands in Taiwan have a higher drainage frequency and a greater number of furrows.

Figure 4.
Boxplots showing that significant differences (P-value < 0.01) exist between Taiwan and Italy in terms of drainage frequency, drainage density, total drainage network length, total number of furrows in the drainage unit, mean drainage network length, plane length, inclination, mean eroded depth, and MSI.

Distinct Erosion Patterns in the Badlands Landscapes of Taiwan
In Taiwan, the higher mean eroded depth (Figure 4h) suggests that more sediment can be produced from the tiny hillslopes, which has low drainage density. To our knowledge, gully erosion in Taiwan is a dominant contributor to sediment production during the rainy season. However, no correlation between drainage density and mean eroded depth (R 2 = −0.18) are consistent with the results from Italy, and this may suggest that the drainages in Taiwan could play the specific role of shredding the landscape, instead of transporting eroded material, such as gully erosion (Figure 5a) and, therefore, high drainage frequency likely cause more free-surface, making hillslopes more prone to erosion. Moreover, no correlation between drainage density and mean eroded depth (R 2 = −0.26) can explain why the relatively short drainage length in Taiwan is not the main contribution to soil loss. Besides, the observed difference of mean eroded depth, which reflects slope inclination ( Figure  5c), may suggest that the erosional pattern is mainly due to gravitational mass wasting, which preferentially occurred during high-intensity rainfall events [47][48][49][50][51]. In addition, the mean circularity ratio is 0.52 in Taiwan, which suggests that the shape of the badlands is not affected by the evolution of a single gully but expands through many short drainage systems. Thus, the mean erosion depth is not influenced by the circularity ratio, which may also suggest that hydrological erosion is not the main erosional process.
During the study period, the total rainfall was 748.1 mm over 35 days (typhoons brought 579 mm of rainfall, with a maximum intensity of 74 mm/h). The high rainfall intensity causes severe erosion Boxplots showing that significant differences (p-value < 0.01) exist between Taiwan and Italy in terms of drainage frequency, drainage density, total drainage network length, total number of furrows in the drainage unit, mean drainage network length, plane length, inclination, mean eroded depth, and MSI.

Distinct Erosion Patterns in the Badlands Landscapes of Taiwan
In Taiwan, the higher mean eroded depth (Figure 4h) suggests that more sediment can be produced from the tiny hillslopes, which has low drainage density. To our knowledge, gully erosion in Taiwan is a dominant contributor to sediment production during the rainy season. However, no correlation between drainage density and mean eroded depth (R 2 = −0.18) are consistent with the results from Italy, and this may suggest that the drainages in Taiwan could play the specific role of shredding the landscape, instead of transporting eroded material, such as gully erosion (Figure 5a) and, therefore, high drainage frequency likely cause more free-surface, making hillslopes more prone to erosion. Moreover, no correlation between drainage density and mean eroded depth (R 2 = −0.26) can explain why the relatively short drainage length in Taiwan is not the main contribution to soil loss. Besides, the observed difference of mean eroded depth, which reflects slope inclination (Figure 5c), may suggest that the erosional pattern is mainly due to gravitational mass wasting, which preferentially occurred during high-intensity rainfall events [47][48][49][50][51]. In addition, the mean circularity ratio is 0.52 in Taiwan, which suggests that the shape of the badlands is not affected by the evolution of a single gully but expands through many short drainage systems. Thus, the mean erosion depth is not influenced by the circularity ratio, which may also suggest that hydrological erosion is not the main erosional process. of the steepest areas in the badlands slopes. The eroded material was then transported downwards and deposited in the lower part of the gully or downstream area (Figure 3). The observed erosion pattern is consistent with the argument that gravitational mass wasting dominates the steep areas of the slope during storm events.

Benefitting from the UAV-Driven High-Resolution Topography in Application of MSI
Traditionally, the investigation of hillslope erosion, i.e., landslides or gully erosion, has been conducted directly through on-site human surveys or through satellite imagery technology (e.g., [52][53][54]). The results of this study demonstrate the practicality and feasibility of using UAV data for topographical analysis of landscape change. A combination of UAV and RTK-GPS surveys can facilitate the efficient acquisition of high-precision and accurate information about slope erosion. This particularly benefits the construction of high-resolution DEMs and DoDs for studying erosion patterns and volumes. The detailed topographic data also further extends the understanding of erosional processes, which is essential to understanding the response of a landscape to climatic forces. Our findings suggest that high-resolution topographic UAV surveys may guide further investigations in other badland areas or active-hillslope erosion regions across climatic zones. This may provide opportunities for clarifying the reaction of morphometric features to different environments. For further application of MSI, two related issues need to be considered. Firstly, while slope length is a function of MSI, the scale of hillslopes can range from tens to hundreds of meters. Therefore, the normalized slope length should always be used. Secondly, MSI is representative of hydrological erosion but is not suitable for landslide-dominated areas. Thus, we suggest that using UAV-driven DEMs to construct the DoD is beneficial for identifying erosion distribution, especially in areas of high rainfall intensity. During the study period, the total rainfall was 748.1 mm over 35 days (typhoons brought 579 mm of rainfall, with a maximum intensity of 74 mm/h). The high rainfall intensity causes severe erosion of the steepest areas in the badlands slopes. The eroded material was then transported downwards and deposited in the lower part of the gully or downstream area (Figure 3). The observed erosion pattern is consistent with the argument that gravitational mass wasting dominates the steep areas of the slope during storm events.

Benefitting from the UAV-Driven High-Resolution Topography in Application of MSI
Traditionally, the investigation of hillslope erosion, i.e., landslides or gully erosion, has been conducted directly through on-site human surveys or through satellite imagery technology (e.g., [52][53][54]). The results of this study demonstrate the practicality and feasibility of using UAV data for topographical analysis of landscape change. A combination of UAV and RTK-GPS surveys can facilitate the efficient acquisition of high-precision and accurate information about slope erosion. This particularly benefits the construction of high-resolution DEMs and DoDs for studying erosion patterns and volumes. The detailed topographic data also further extends the understanding of erosional processes, which is essential to understanding the response of a landscape to climatic forces. Our findings suggest that high-resolution topographic UAV surveys may guide further investigations in other badland areas or active-hillslope erosion regions across climatic zones. This may provide opportunities for clarifying the reaction of morphometric features to different environments. For further application of MSI, two related issues need to be considered. Firstly, while slope length is a function of MSI, the scale of hillslopes can range from tens to hundreds of meters. Therefore, the normalized slope length should always be used. Secondly, MSI is representative of hydrological erosion but is not suitable for landslide-dominated areas. Thus, we suggest that using UAV-driven DEMs to construct the DoD is beneficial for identifying erosion distribution, especially in areas of high rainfall intensity.

Conclusions
Through a morphometric analysis of UAV-acquired DEMs and DoDs, we have reached the following conclusions. First, by conducting UAV surveys, we were able to quantify the landform change to centimeter-scale accuracy across 72 badlands hillslopes and construct high-resolution DEMs and DoDs. These provide reliable data constraint on the relationship between hillslope-morphology and erosion and will benefit the monitoring of land surface change as a result of exogenic forces.
The topographic effects of the climatic difference are important to reveal the mechanism of landscape evolution. Our observations indicate that badlands hillslopes in Taiwan are relatively short, small, and steep compared with those of Italy. Because the annual precipitation of over 2 m causes strong incision, the mudstone badlands in Taiwan are characterized by higher drainage frequency and more furrows. Thus, we suggest that MSI can be used in the badlands areas of Taiwan. Our findings suggest that MSI has a strong positive correlation with erosion, which is consistent with results from Italy. However, the mean eroded depth shows no correlation with mean drainage length and a positive correlation with inclination, suggesting that the erosion pattern is due to gravitational mass wasting, rather than hydrological erosion.
MSI combined with UAV-driven DEMs, provides a relatively convenient and rapid method for assessing the severity of soil loss and evaluating landscape maturity. We highlight the unique cross-climate zone comparative study offered by badlands landscapes, with differential rainfall patterns, using simple survey techniques and morphometric slope index. The results of this study also suggest that UAV techniques are able to provide quantitative evidence for hillslope morphological analysis, and through simple morphological indicators, we can achieve comparative studies in different climate zones.