Bank gully extraction from DEMs utilizing the geomorphologic features of a loess hilly area in China

Abstract

As one of most active gully types in the Chinese Loess Plateau, bank gullies generally indicate soil loss and land degradation. This study addressed the lack of detailed, large scale monitoring of bank gullies and proposed a semi-automatic method for extracting bank gullies, given typical topographic features based on 5 m resolution DEMs. First, channel networks, including bank gullies, are extracted through an iterative channel burn-in algorithm. Second, gully heads are correctly positioned based on the spatial relationship between gully heads and their corresponding gully shoulder lines. Third, bank gullies are distinguished from other gullies using the newly proposed topographic measurement of “relative gully depth (RGD).” The experimental results from the loess hilly area of the Linjiajian watershed in the Chinese Loess Plateau show that the producer accuracy reaches 87.5%. The accuracy is affected by the DEM resolution and RGD parameters, as well as the accuracy of the gully shoulder line. The application in the Madigou watershed with a high DEM resolution validated the duplicability of this method in other areas. The overall performance shows that bank gullies can be extracted with acceptable accuracy over a large area, which provides essential information for research on soil erosion, geomorphology, and environmental ecology.

Keywords

bank gully / DEMs / topographic features / loess shoulder line / relative gully depth

Introduction

Gully erosion is a serious environmental problem and primary source of sediment loss in the Chinese Loess Plateau (^{Luo, 1956}; ^{Zhu, 1956}; ^{Liu, 1985}; ^{Valentin et al., 2005}; ^{Li et al., 2015}). Data collected in different parts of the world showed that soil loss rates by gully erosion represent 10% to 94% total sediment yield caused by water erosion (^{Poesen et al., 2003}) and contribute 60% to 90% total sediment production on agricultural land in the hilly areas of the Chinese Loess Plateau (^{Li et al., 2003}). In the Chinese Loess Plateau, bank gullies mainly develop below shoulder lines, that is, the boundary between positive (excavated or interfluves) and negative (valley) terrains (^{Wu and Cheng, 2005; Zhou et al., 2010}). As one of the most active types of gullies, they lead to the evolution of hillslope-channel coupling related to flow erosion, which is the most significant aspect of coupling in headwater fluvial systems (^{Harvey, 2002}; ^{Chen et al., 2007; Evans and Lindsay, 2010}).

Some bank gully-oriented studies have been carried out with a focus on erosion rate (^{Jing, 1986; Li et al., 2012}), sediment production (^{Li et al., 2014}), geomorphologic thresholds (^{Hu and Wu, 2005}), and developmental prediction (^{Sidorchuk, 1999}; ^{Hessel and van Asch, 2003}; ^{Li et al., 2012}, ^{2015; Zhang et al., 2015}) to understand the contribution of bank gullies to soil loss and landscape evolution. These studies have relied on the explicit presentation and efficient extraction of bank gullies. Mapping gullies helps in understanding their characteristics, such as location, length, longitudinal slope, morphology, density, and distribution, which are necessary information for modeling gully head retreat rate, calculating critical catchment area, exploring gully developmental stage, and predicting gully erosion and sediment yield (^{Harvey, 2001}; ^{Shruthi et al., 2014}; ^{Na et al., 2016}). Compared with sheet or rill erosion, gully erosion is more difficult to monitor accurately in the field. Therefore, this study proposed a semi-automatic DEM-based extraction method for the large-scale mapping of bank gullies.

Currently, remote sensing images and DEMs are the two main datasets for regional gully extraction. Except for visual interpolation (^{Hessel, 2002; Li et al., 2015}), pixel-based image analysis is a reliable technique (^{Metternicht and Zinck, 1998}) that applies spectral or textural heterogeneity to extract gullies. However, spectral heterogeneity is sensitive and easily affected by different variables such as local land cover, shadows, and the atmosphere, which strongly influence the performance of these methods. The advancement of technology and increased convenience of acquiring high-resolution images (HRI) facilitated a shift from traditional pixel-based methods to object-oriented analysis (OOA) methods. Some OOA-based studies have effectively monitored gullies over a large area (^{Knight et al., 2007}; ^{Shruthi et al., 2011}; ^{Shruthi et al., 2015}; Li et al., 2017). However, the results of these studies have failed to differentiate bank gullies from the gully system and could not be used for bank gully studies. Furthermore, although the detailed shape of gullies can be extracted from an HRI, OOA-based results merely present geometric information in two dimensions such as area and shape. Some important geomorphological properties (e.g., depth, longitudinal slope, curvature, and flow accumulation) and topographic structural features (e.g., gully head and gully channel line), which are necessary for understanding the development process of bank gullies, are difficult to extract directly.

Most DEM-based methods for mapping gullies focus on channel extraction, which could derive some 3D gully information but still fail to differentiate bank gullies from other gullies. The channel networks, that is, the gully networks in this study, can be generated by simulating surface runoff as part of the hydrological process (^{O’Callaghan and Mark, 1984}; ^{Fairfield and Leymarie, 1991}; ^{Tarboton et al., 1991}; ^{Lashermes et al., 2007}; ^{Afana and Del Barrio, 2009}). In this process, the gully network is defined by a threshold segmentation of flow accumulation. Small gullies include bank gullies that can be extracted with minimal flow accumulation thresholds. However, some problems are still unresolved. First, classical channel extraction produces a large amount of pseudo or parallel channels. Thus, numerous parallel channels still appear in flat areas when using a small flow accumulation threshold even when a multi-flow direction or priority-flood algorithm is employed (^{Hellweger, 1997}; ^{Saunders, 1999}; ^{Turcotte et al., 2001}; ^{Barnes et al., 2014}; ^{Zhou et al., 2016}). These parallel channels can be found in the main and branch channels. Second, the position of the gully head is usually inaccurate. For DEM-based hydrologic methods, the location of the gully head is controlled by the flow accumulation threshold, which is significantly influenced by subjectivity. In addition, fixed thresholds are unsuitable for bank gullies with diverse lengths. A bank gully head is generally located just below the loess shoulder line, which is the boundary between positive and negative terrains (P–N terrain) (^{Tang et al., 2007}; ^{Zhou et al., 2010}; ^{Song et al., 2013}; ^{Yan et al., 2014}; ^{Zhu et al., 2014}). If this spatial relationship is disregarded, then the extracted gully head usually failed to match its actual position. Third, bank gullies still cannot be readily differentiated from derived networks. Thus, this study addresses these problems relating to the lack of an automatic and robust method for a large-scale bank gully extraction.

The accurate identification of bank gullies is impossible without expert knowledge on auxiliary information such as geometric properties (e.g., shape and orientation) and spatial relationships with surrounding features (^{Wu et al., 2008}). In this study, we propose an extraction method for bank gullies utilizing their geomorphologic features. First, the gully network containing bank gullies can be extracted by the proposed algorithm called iterative channel deepening based on a small threshold of flow accumulation, and pseudo or parallel gullies can consequently be eliminated. The extracted gully heads can then be positioned correctly considering the spatial relationship between bank gullies and the loess shoulder line. Finally, a new measurement called the relative gully depth (RGD) is proposed based on the geomorphologic features of bank gullies; this measurement can be employed to distinguish bank gullies from other gully types. This study aims to: (i) describe a DEM-based semi-automatic method of large-scale bank gully extraction considering geomorphologic features; (ii) propose a method of eliminating pseudo and parallel channels when using minimal threshold values of flow accumulation for regional gully network extraction; and (iii) establish the parameters of this method in the study areas and assess the accuracy of the results.

Study area

This study was carried out in the Linjiajian watershed at the middle reaches of the Wuding river in the Suide county of China. This watershed covers an area of 12.35 km² between 10°18′04″E–110°22′04″E and 37°32′47″N–37°34′47″N. The altitude ranges from 867 m to 1186 m, and the average surface slope is 30°. A training area of 1.8 km² was selected for parameter optimization, and six sample areas were used for field investigations (Fig. 1). High-density loess gullies are among the most distinguishing geomorphologic features of the selected area (^{Wu and Cheng, 2005}). The morphology of these gullies consistently has similar characteristics with those distributed in other areas such as an evident loess shoulder line, which separates the terrain surface into gentle upslope and steep downslope areas. This area is representative of the fragmented terrain in loess hilly landforms, where thousands of gullies develop, including rills as well as ephemeral, bank, and valley gullies, with intense soil erosion. The area has a semi-arid continental climate, with an average annual temperature of 8°C, average precipitation of ~450 mm/yr, and 1615 sunshine hours annually. The rainy season spans from July to September, and accounts for 64.4% of the total annual precipitation, often concentrated in one rainstorm. Field investigations show that the sediment yield from a single rainstorm is often greater than 60% of the annual yield.

Materials and methods

Materials

DEMs and digital orthophoto maps (DOMs) are utilized in this study. The former, generated from the contours of a 1:10,000 scale topographic map, is used for bank gully extraction, whereas the latter is employed to generate reference data for assessment. Both datasets are created by the National Administration of Surveying, Mapping, and Geoinformation of China with resolutions of 5 m and 1 m, respectively.

Methods

The bank gully extraction procedure consists of five steps: field survey, gully network extraction, gully head revision, bank gully identification, and accuracy assessment (Fig. 2). First, field surveys help establish the geomorphologic knowledge of bank gullies, such as critical catchment area, width, and depth. The second step is extracting a full gully network in which bank gullies are included, while eliminating parallel and pseudo-gullies through an interactive channel deepening algorithm. The third step is to correctly position gully heads, and in the process, the loess shoulder line acts as a critical reference. The fourth step is identifying bank gullies from the extracted gully network by RGD, a newly proposed indicator. Finally, the method is assessed for validation.

Field surveys of bank gullies

The basic characteristics of bank gully were first investigated through field surveys in the sample area to help establish criteria for extraction and validation. During the field work, the heads of bank gullies were mapped on 1:10,000 topographic map, and the topographic features of typical bank gullies were estimated by a laser rangefinder. The upslope catchment areas of the gully head measured were calculated from the DEMs. The selected watershed was mainly investigated relative to the following aspects:

1) Position of bank gully heads. A total station technique was applied to accurately measure the position of bank gully heads. The prism-free mode was also used when gullies were inaccessible. Those measured gully head positions were utilized with DEM to calculate the critical upslope catchment area of bank gullies. The critical area could set the threshold in the following iterative process of gully network extraction. Furthermore, gully head positions were used to validate the bank gullies extracted from DEM.

2) Width and depth. A laser rangefinder was employed for width and depth measurement. Each gully was measured at the top, middle, and bottom for three to five times for the averaged values to represent the whole gully. The width and depth of bank gullies could contribute to implementing the criterion of bank gully identification in the artificial bank gully delineation process.

All spatial data obtained from field surveys were digitized into a GIS database with corresponding DEM and DOM images. These results contribute to the identification and parameter determination of extractions.

Extraction of gully networks

A typical solution to the aforementioned parallel channel problem is the modification of DEMs using digitized rivers, lakes, or other hydrologic objects with techniques such as the Agree-DEM algorithm (^{Hellweger, 1997}), the Burn-in algorithm (^{Saunders, 1999}), or the Digital River and Lake Network (DRLN) method (^{Turcotte et al., 2001}). The basic principle of these algorithms is increasing the elevation difference near the river by deepening and decreasing the elevation of the river position. By comparing the multi-flow direction of the Dinf algorithm (^{Tarboton, 1997}) and the priority-flood algorithm, we found that a combination of the D8 (deterministic eight-node) algorithm (^{O’Callaghan and Mark, 1984}) and burn-in algorithms could effectively extract a clean channel network that conforms with the digital river data (^{Zhang et al., 2012}). Therefore, the burn-in algorithm is used in this study to eliminate parallel or pseudo gullies.

Digital river data are required to modify the DEMs in the burn-in algorithm. However, no perennial rivers exist in most upstream areas, which mean that no river data can be used to modify the DEMs. When using a greater threshold of flow accumulation in hilly areas, the extracted main channel networks emerged cleanly and naturally without any parallel channels. Thus, this result could be used in place of perennial rivers to modify DEMs. If we repeat this procedure with smaller threshold, then small gullies can be extracted meanwhile parallel gullies are eliminated gradually. Therefore, an iterative channel deepening algorithm was proposed (Fig. 3). In this process, gully heads continuously move upstream for the bank gully extraction.

The extraction of gully networks using ArcGIS Hydro analysis tools involves the following steps:

1) Fill the sinks of the DEMs, calculate the flow direction based on the D8 algorithm, and generate the flow accumulation dataset.

2) Set a large threshold of flow accumulation to generate the initial gully network, which mainly includes downstream channels. If these downstream channels are confirmed without parallel channels through a comparison with terrain hillshade images, then they can be regarded as digital rivers and used to modify DEMs based on the burn-in algorithm.

3) After calculating flow direction and accumulation based on the modified DEMs, extract the new gully network using a small threshold of flow accumulation. The small channels in this network that can be extracted are increased.

4) If these small channels still cannot describe bank gullies, then remodify DEMs based on the burn-in algorithm using the new network derived in the above step as digital rivers.

5) Repeat Steps 3) and 4) until the small channels in the extracted network could describe bank gullies. In this iterative computation process, the flow accumulation threshold is set to a large number at first, and then gradually decreased. The parallel channels in the main gully are initially eliminated, followed by those in the branch gullies. Small bank gullies are retained.

In the above procedure, the series of threshold values is generally based on whether the extracted channel network has parallel channels. First, the initial threshold can be set to a large value (e.g., 50% of the watershed area). This value may increase if parallel channels still appear. Second, the values will continuously decrease. Finally, the termination threshold value is determined as the critical catchment area of the bank gullies measured by field surveys.

Revision of gully heads

The gully heads obtained in Section 3.2.2 cannot reach their virtual positions completely as they typically either overshoot or undershoot. Hence, a revision process for gully head position is necessary (Fig. 4). A gully head is usually located where a sharp gradient change exists from a gentle to a steep slope, related to the highly developed vertical joints in the loess stratum. Horizontally, gully heads usually tend to approach the loess shoulder line; this line is the most critical terrain structure line in loess landforms, and it divides the surface into a gentle inter-gully and a steep gully area in the positive terrain and negative terrain, respectively (Fig. 1).

Therefore, this spatial relationship is key for revising bank gully heads. Our comparison with existing DEM-based loess shoulder line extraction algorithms found that the P–N terrain segmentation algorithm (^{Zhou et al., 2010}) achieved satisfactory accuracy and continuity. Hence, gully heads could be modified by trimming the overshot part (red solid line) or extending the undershot gully line (yellow dashed line) by tracking an upstream cell with the maximum flow accumulation in a neighboring 3 × 3 cell window until the head reaches the loess shoulder line (^{Jiang et al., 2013}) (Fig. 4).

Identification of bank gullies

Field observation showed that the majority of bank gullies have a linear development despite few of them containing tributaries. This finding means that an overwhelming number of bank gullies could be labeled as first-order gullies according to Strahler’s classification (^{Strahler, 1963}), making it possible to distinguish bank gullies from gully network. Therefore, the first-order gullies from the Strahler’s classification are treated as candidates for bank gully identification. Note that bank gully is conceptually different from the first-order gully. The former is defined in a geomorphological sense, whereas the latter is from Strahler’s law of channel networks in basin hydrology. Therefore, not all first-order gullies are bank gullies, and some bank gullies may have branches. However, under the appropriate flow accumulation threshold in Section 3.2.2, numerous bank gullies are first-order gullies. Hence, an index of RGD is proposed to identify bank gullies.

The depth of valley gullies is generally much greater than that of bank gullies, showing a difference in their geomorphologic development stage (^{Liu et al., 1988}). The RGD is the average elevation difference between a gully skeleton line and its buffer boundary. The identification process of bank gullies from valley gullies uses an appropriate buffer size and segmentation threshold of RGD (Fig. 5).

The procedure details are as follows:

1) Classify the gully network using the Strahler classification method and extract the first-order gullies;

2) Generate a buffer zone of this first-order gully, and calculate the average elevation of the buffer boundary;

3) Calculate the average elevation of the first-order gully;

4) Calculate RGD, which is the difference between the results of Steps 2) and 3);

5) Set a segmented threshold of RGD. Gullies with RGD over the threshold are identified as bank gullies.

Accuracy assessment

The reference data for the assessment can be generated by manually identifying bank gullies from the gully network according to previous knowledge acquired through field surveys. Several standards including length, width, depth, and transverse profile are established for distinguishing bank gullies from the gully network (see Section 5.1.2). Bank gullies could be identified artificially one by one based on the above information. Furthermore, those missed in the extraction by DEM can be delineated by the DOM (aerial image) dataset. The gullies derived by artificial mapping from DOM may not be more accurate than that extracted using DEM because the image interpretation is inevitably influenced by individual subjectivity. In addition, gully lines derived from DEM rely on rigorous computation of well-designed algorithms with detailed description of the topographic feature, i.e., high resolution of DEMs (^{Martz and Garbrecht, 1992}; ^{Benaïchouche et al., 2016}; ^{Persendt and Gomez, 2016}). Therefore, accuracy can be assessed by comparing our result with i) field survey in a small sampling area, and ii) reference data by manual identification.

The accuracy of bank gully extraction could be acquired by measurements of producer’s accuracy (PA) and user’s accuracy (UA), which were designed by ^{Congalton (1991)} and widely applied in the accuracy assessment of imagery classification. We introduced the measurements and re-designed them for the assessment. PA is the positive precision for bank gully identification, whereas UA is the ability of the method to identify bank gullies from first-order gullies. The two measurements can be defined as follow:

(1)\mathrm{P}\mathrm{A}={\displaystyle \frac{\mathrm{I}\mathrm{C}}{\mathrm{I}\mathrm{C}+\mathrm{M}\mathrm{I}}},

(2)\mathrm{U}\mathrm{A}={\displaystyle \frac{\mathrm{I}\mathrm{C}}{\mathrm{I}\mathrm{C}+\mathrm{I}\mathrm{W}}},

where IC is the amount of correctly identified bank gullies; IW is the amount of wrongly identified bank gullies; MI is the amount of misidentified bank gullies. These two measurements are also used to select optimal parameters of gully buffer size and RGD threshold in Section 3.2.4 during the comparison of the identified results and reference data in the sampling areas.

Results

Gully system

As mentioned in Section 3.2.1, the potential of an iterative channel deepening algorithm for extracting gully networks which includes small gullies without pseudo or parallel channels was analyzed. A series of flow accumulation thresholds was adopted in this iterative process. These thresholds were tested according to the principle of an order from large to small. It is the fact that large threshold of flow accumulation will result in large gullies without parallel channels, which could be used for the following channel deepening algorithm. We found that a main channel network without any parallel channels in the watershed could be extracted with a threshold of 400,000 cells, that is, a catchment area of 10 km² used as the initial threshold. The second threshold would be smaller than the first to increase the number of small channels that could be extracted. The appropriate value for the second threshold is the critical one with no parallel channels generated at this level. The principle is the same for the next series of threshold determination. Finally, we obtain the threshold series for iterative deepening process using the Agree-DEM algorithm (see Section 3.2.2) as 300,000; 100,000; 50,000; 10,000; 1000; 500; 200; 100; 50; 40; and 20 cells, representing 7.5; 2.5; 1.25; 0.25; 0.025; 0.0125; 0.005; 0.0025; 0.00125; 0.001; and 0.0005 km², respectively. Among them, the last threshold depends on the survey of upslope catchment area from 60 bank gullies heads. Field surveys indicated that the average upslope catchment area of bank gully heads in this watershed is approximately 500 m², that is, 20 cells (Table 1).

Figure 6(a) illustrates the gully network extracted by the original DEMs with a threshold value of 20. Pseudo and parallel channels are evidently abundant in the main and branch gullies. Figures 6(b) show three enlarged windows of comparison between the result without channel deepening and those after 4 and 10 iterations, respectively. The parallel channels in the main and branch gullies were progressively eliminated. After 14 iterations, the bank gully network without pseudo and parallel channels was extracted.

The gully heads of the derived gully network were then revised using the spatial relationship restriction from the shoulder lines. As revealed in Fig. 7, the terrain surface was divided into gentle upslope and steep downslope surfaces by the loess shoulder line. The gully heads were revised to their actual positions at immediately below the loess shoulder line.

Bank gully identification by RGD

We examined the potential of RGD to obtain accurate classification results. In the bank gully identification procedure based on RGD, the buffer size of first-order gullies and the threshold of RGD are two key parameters. A total of 368 gullies in the training area were classified into 357 bank and 11 valley gullies based on image features and field surveys, which were used as reference data for parameter optimization. A series of buffer sizes and thresholds were employed as input parameters to systematically test the proposed model and to obtain the most appropriate results for determining the most effective parameter values. The radius of the buffer was set to 5, 10, 15, and 20 m, given that most bank gully widths are within 20 m, and the DEM resolution is 5 m; RGD was subsequently calculated. Meanwhile, different segmentation thresholds were set at each buffer size to check identification accuracy by comparing against the reference data of the training area.

Table 2 shows an increasing trend in PA and a consistent trend in UA with increasing buffer sizes using Eqs. (1) and (2). With a certain buffer size, increasing the RGD segmentation threshold leads to an upward trend in PA and a decline in UA. This finding can be attributed to more first-order gullies being selected as bank gullies with a high RGD threshold. Consequently, more valley gullies were mistakenly identified as bank gullies to restrain UA. Therefore, the most reasonable parameters should be a balance between PA and UA. After comparison with the reference data in the training area, a 15 m buffer radius and a 2.5 m/m RGD segmentation threshold were selected as optimal parameters for bank gully identification in the study area. The identification procedure was performed on the whole study area using this set of parameters, and 1340 bank gullies were identified automatically (Fig. 9).

Geomorphologic characteristics of bank gullies

Different types of gullies reflect various developmental stages and soil erosion intensities. Distinguishing bank gullies from other types of gullies will help understand their contribution to modern active soil erosion. Furthermore, exploring their characteristics will reveal their geomorphologic features, development mechanisms, and role in the complex integrated sequence of gully erosion. Therefore, some statistical indicators for bank gullies were calculated (Table 3).

Among the 1340 bank gullies identified, the longest was 236.7 m, the shortest was 17.5 m, and the mean was 65.7 m. As discussed above, the accuracy of bank gully extraction depends on the DEM resolution. The higher the DEM resolution is, the smaller the bank gullies that can be extracted. In this study with 5 m resolution DEMs, the smallest bank gully identified was 17.5 m long, which is 3.5 times of the DEM grid size. This experiment failed to identify any significantly shorter bank gullies due to the limited DEM resolution. Consequently, 5 m DEMs can be used to extract bank gullies over 20 m in length with reasonable accuracy.

Gully density is usually used to describe the fragmentation degree of the land surface and the intensity of soil erosion. The statistics of gully density generally consider permanent gullies but not bank gullies. Bank gullies make the terrain more fragmented, thereby causing considerable surface erosion; hence, bank gullies are indispensable for calculating gully density. Table 3 shows that the value of gully density including bank gullies is approximately twice as that without bank gullies. When the entire area is divided into sunny and shady slopes, bank gully distribution maintains a balance in both aspects, which reflects uniform environmental conditions. The longitudinal slope is the ratio of bank gully elevation difference against its horizontal length. The mean longitudinal slope is 62% with a range of 20% to 147%. This metric measurement indicates a relatively steep slope of loess bank gullies, which is consistent with their early stage of development.

Discussion

Validation

Validation with field survey

Evaluating the accuracy of extracted bank gullies is necessary although a visual assessment of the results already confirmed a highly accurate identification. Field surveys are the most fundamental means of validation. The field investigation confirmed that the skeleton line of bank gullies could be represented effectively based on DEM. Therefore, the bank gullies extracted from DEM are regarded as accurate if the measured gully heads from the field survey correspond to the ones extracted. We surveyed bank gully heads using total stations in six sampling areas. Our survey result showed a satisfactory matching result by comparing gully heads between our results and the surveyed ones (Fig. 8). Distance variations in the two results are all less than 4 m, which is acceptable for the DEM data in the 5 m resolution. It also showed that some small bank gullies could not be extracted owing to the missing simulation of small gully morphology in this DEM resolution.

Validation with image-based interpretation

Accuracy assessment is typically based on the comparison of analytical results and reference data, but it is insufficiently meaningful to verify each individual gully by field investigation (^{Shruthi et al., 2011}) because generating reference data for the whole study area through field surveys is challenging. Reference data were produced by manual interpretation based on DOMs and DEMs. The interpretation criteria came from a field investigation of 60 bank gullies in the sample area (Table 4), and are as follows: i)<200 m long,<20 m deep, and<30 m wide; ii) v-shaped transverse profiles with steep slopes; and iii) just below the shoulder line. A total of 1531 gullies were classified as bank gullies and 111 as valley gullies, which can be regarded as validation data for the accuracy assessment.

Compared with the manual identification, 1340 bank gullies among a total of 1354 candidate gullies (orange line in Fig. 9) were correctly identified, whereas the 14 remaining gullies were classified as valley gullies. A total of 191 misidentified gullies (green line in Fig. 9) were generally located in the high-gully-density slopes. The gully buffer boundaries in these areas are influenced by nearby neighboring gullies, that is, crossing other gullies, which misestimated the gully depth in calculating GRD. Overall, the producer accuracy of identification was 87.5% with a user accuracy 99.0%, which can be regarded as reasonable and acceptable.

Comparison with other methods

Comparison with other flow direction algorithms

As mentioned in Section 3.2.2, the proposed iterative channel deepening algorithm focuses on removing the parallel flow patterns in flat areas of DEMs based on the stream burn-in method. However, note that many algorithms have been proposed to define flow direction over flat areas (^{Garbrecht and Martz, 1997}; ^{Liang and MacKay, 2000}; ^{Jones, 2002}; ^{Wang and Liu, 2006}; ^{Coppola et al., 2007}; ^{Barnes et al., 2014}). Some of them have reduced parallel channels significantly in flat areas. Therefore, we carried out a comparison to investigate the feasibility of the flat area-oriented flow direction definition algorithm for bank gully extractions. Here, we adopted the method of Barnes (^{Barnes et al., 2014}) to realize flow direction definition in flat areas for two reasons; 1) its effectiveness for flow direction definition in depressed or flat areas; 2) its relatively simple and easy implementation.

Figure 10 shows a comparison of channel networks generated by D8 and Barnes’ algorithms with the flow accumulation threshold of 20, which is the critical threshold for bank gully extraction in our study area. No difference was observed in the upstream, whereas great difference appeared in the downstream flat area. The majority of parallel channels in the downstream (flat bottom of main gullies in this study) based on the D8 algorithm were removed using Barnes’ algorithm. However, some parallel channels were still generated in main stream channel.

Although the flow direction in flat area could be defined and improved by flat area-oriented algorithms, the appearance of parallel channels could be attributed to the small threshold of flow accumulation in a relatively large flat area.

As described in Fig. 11(a) which shows the flow direction definition result in the flat area from the Barnes algorithm (^{Barnes et al., 2014}), parallel channels are still generated (Figs. 11(b) and 11(c)) when using a small flow accumulation threshold. Suppose that the position ○ in Fig. 11 is the outlet of the drainage area, nine parallel channels are derived using the threshold of 3 cells (Fig. 11(b)), which is much more than that of 3 parallel channels by using the threshold of 5 cells (Fig. 11(c)). The greater the threshold is, the fewer parallel channels. Therefore, parallel channels would not be eliminated with such a small flow accumulation threshold value even when using a flat area-oriented flow direction algorithm. In this study, a small flow accumulation threshold is required for bank gully extraction owing to small scale bank gullies in the upslope area. Hence, the proposed iterative channel deepening algorithm is necessary for removing parallel channels in large flat areas.

Interestingly, no parallel channel was generated when the flow accumulation threshold is greater than 500 cells based on Barnes’ algorithm (Fig. 12). The iterations will significantly decrease by Barnes’ algorithm as the iterative channel deepening algorithm starts at a flow accumulation value when a clean channel network is extracted without any parallel channels. According to the experiment result, the iteration of channel deepening algorithm could start with the flow accumulation of 500 cells, followed by 200, 100, 50, 40, and 20, similar to before. Therefore, the iterations of the channel deepening process could be reduced from 14 to 6. Furthermore, no channel network difference in the upslope area from the D8 and Barnes’ algorithms could be found except for the main channel network in the downstream area (Fig. 10). The series of flow accumulation thresholds for the channel deepening process in our method could be applied in the channel network result from Barnes’ algorithm, given that our aim is to identify the bank gullies located in the upslope area.

The iterative channel deepening algorithm is necessary to obtain a clean channel network by removing parallel channels for the bank gully extraction in this study. Our method could be easily duplicated without any programing, and its iterations can be reduced by the improved algorithm for flow direction definition in flat areas.

Comparison with an automatic adaptive gully extraction method

The upslope catchment area is one of the most important features to automatically extract bank gully accurately, and is used as the termination threshold of iterative channel deepening in Section 3.2.2. In this study, the catchment area can be settled through field survey. This finding means that our method depends on prior field surveys. Afana and del Barrio (2009) proposed an adaptive approach for channel network delineation based on the assumption that DEMs are self-contained structures for detecting drainage networks, and channel complexity is best reflected by its corresponding intrinsic properties. The method could detect the optimal channel network that best describes landscape dissection especially in a large area of heterogeneous landforms. Therefore, we incorporated the method into this study area to inspect whether it remains effective for bank gully delineation.

A comparison of the extracted channel network between Afana’s and our method showed a clear difference (Fig. 13(a)). According to his method, when using 170 cells as the critical threshold, which is the first maximum change point in the relationship curve between exterior and interior link lengths ratio (Fig. 13(b)), namely, R_A (^{Schumm, 1956}) and critical catchment area (A_S), respectively, the channel network could be delineated (yellow line in Fig. 13(a)). Although this method extracted the valley gullies successfully, it failed to extract the bank gullies (red line in Fig.13(a) in our result). Note that his original work was applied on DEM with 1 cm and 6 cm resolution, and thus, we suppose that the performance of Afana’s method is limited by the 5 m DEM resolution.

Application on high resolution data

To validate the applicability of our method in other areas and on high resolution data, an additional experiment was performed on a small watershed of a hilly loess area: Madigou, located in Jingbian County with an area of 0.23 km². A 0.5 m DEM of this area was generated through point clouds by the terrestrial laser scanner, Riegl VZ400, and aerial photos were acquired by an unmanned aerial vehicle in August 2014 for validation. Based on the method proposed in this study, the threshold flow accumulation (upslope catchment area) parameters during the iterative channel deepening process were the same as that of the Linjiajian watershed. Owing to the small area, the latter six parameters such as 0.0125, 0.005, 0.0025, 0.00125, 0.001, and 0.0005 km² were used. The loess shoulder line was derived accordingly. The buffer size of RGD identification was set to 5 m and the segmentation threshold to 2.5 m/m. Compared with the field investigation, the PA of the extraction results increased to 88.9% owing to the increased small bank gullies extracted with high resolution DEM.

The Jingbian test area results showed that the process of iterative channel deepening could eliminate parallel channels in the main stream (Fig. 14(a)) and that gully heads could be revised correctly by the shoulder line (Fig. 14(b)). All bank gullies could be extracted with this high-resolution DEM (Fig. 14(b)).

Interestingly, the upslope catchment area parameters used in iterative channel deepening with the 5 m DEM can also be applied to the 0.5 m DEM, which proves its applicability in hilly loess areas. The minimal upslope catchment area of 500 m² is significant and can be utilized in bank gully extraction in hilly loess areas. We are optimistic that this parameter of the interactive channel deepening procedure can also be used in other hilly loess areas.

Note that the parameters of RGD identification differ from those for 5 m DEMs because the distance between adjacent bank gullies is closer than that in Linjiajian; thus, a large buffer size of RGD may cover two or more gullies, thereby decreasing accuracy. We believe that undertaking a preliminary field survey of a small training area can facilitate the successful extraction of large-scale bank gullies with this method.

Implications of bank gully extraction

With the continuously increasing interest in bank gullies, extracting and mapping them over a large area is crucial for studies on soil erosion and geomorphologic development.

In the gully head retreat research, the accurate position of a gully head is essential to model the erosion process (Oostwoud ^{Wijdenes et al., 2000; Nyssen et al., 2002}; ^{Martınez-Casasnovas, 2003}; ^{Castillo and Gómez, 2016}; ^{Rengers et al., 2016}). Previous works have been usually conducted in a limited area like a certain hillslope or a runoff plot, and erosion rates have been estimated using this information (^{Castillo and Gómez, 2016}). However, whether the gully head retreat rate in a large area is homogeneous remains unclear. In this study, the head retreat rate of each gully in the whole watershed could be calculated by our extraction, which could provide detailed information and help to understand the development and heterogeneity of gullies in the spatial distribution from a large scale.

From the perspective of gully development, information such as amount, length, position, slope, and spatial distribution which could be acquired through the method proposed in this paper, are of great help for gully developmental predictions, such as the statistic model based on gully density (^{Zhao et al., 2016}), and for simulations based on cellular automatic models (^{Dunkerley, 1997; Cao et al., 2013}). When separating the bank gully from the gully system, it will help to understand the role and contribution of the bank gully to the drainage system development.

Landform, terrain complexity, and fragments are mainly determined by gully development within a watershed (^{Strahler, 1957}; ^{Wolman and Gerson, 1978}; ^{Wondzell et al., 1996}). The gully density or area could present general information on watershed characteristics without any interior structure content. A gully system usually develops within a drainage area, in which various types of gullies represent different developmental stages and evolve into an organic system (^{Seginer, 1966}; ^{Betts et al., 2003}; ^{Stolte et al., 2003}). The compositional proportions and spatial distributions of these gullies depict geomorphologic features and gully erosion processes directly. For instance, various types of gullies with the same gully density may cause different soil loss and form different landforms, which would require further research. The classification and extraction of different gullies could provide basic information which would contribute to understanding their roles in the loess deposition process and geomorphologic evolution mechanism.

Conclusions

Accurate and comprehensive information on erosion features are of critical importance to farmers, land managers, and scientists. Extracting and distinguishing bank gullies from the gully system provide basic information to the bank gully head retreat research and gully development prediction and simulation, which will help understand the role and significance of this kind of active gully in gully erosion, soil loss, and geomorphologic evolution. To achieve these, we first need to determine gully location and extent, and our study introduces a method to provide this essential information. In this study, a semi-automatic method of bank gully extraction was proposed, which considers the topographic features and spatial relationships of bank gullies. The experiment in the loess hilly areas showed that the overwhelming majority of bank gullies could be extracted in a large area with this method, although a field investigation in the sampling area should be conducted first. This condition shows an evident advantage of improving efficiency. The application in the Jingbian test area presented its reproducible in other areas and high resolution DEMs.

Digital river datasets are usually involved in the burn-in algorithm to eliminate the parallel channels in the DEM-based extraction of channel network. The “iterative channel deepening” procedure proposed in this paper could solve this problem and achieve a clean network especially in the area with no applicable digital river. The application of different resolution DEMs proved the practicality of this method. A comparison with another advanced flow direction algorithm over flat areas proved the necessity and validation of the iterative channel deepening algorithm for bank gully extraction in this study.

The accuracy of gully head location can be improved by a revision using the shoulder line, which corresponds to the relationship between shoulder line and bank gully head. Only under the constraint of a shoulder line could each gully head be repositioned correctly.

The RGD can be used to identify bank gullies from valley gullies according to a relatively large difference in gully depth between bank and valley gullies. The segmentation threshold of RGD varies in different landforms, and its optimization could be achieved by an experiment in the sample area.

The 5 m resolution DEM may be coarse for identifying small bank gullies, but its effectiveness was confirmed for extracting the majority of bank gullies. The additional experiment in the Jingbian test area indicated that the methodology and parameter were valid on other loess hilly areas and high-resolution DEM data. Thus, we are optimistic about the potential application of these techniques to other areas.

References

[1] Afana A, del Barrio G (2009). An adaptive approach for channel network delineation from digital elevation models. In: Proceedings of Geomorphometry. Zurich, 224–230

[2] Barnes R, Lehman C, Mulla D (2014). Priority-flood: an optimal depression-filling and watershed-labeling algorithm for digital elevation models. Comput Geosci, 62(1): 117–127

[3] Benaïchouche A, Stab O, Tessier B, Cojan I (2016). Evaluation of a landscape evolution model to simulate stream piracies: insights from multivariable numerical tests using the example of the Meuse basin, France. Geomorphology, 253: 168–180

[4] Betts H D, Trustrum N A, Rose R C D (2003). Geomorphic changes in a complex gully system measured from sequential digital elevation models, and implications for management. Earth Surf Process Landf, 28(10): 1043–1058

[5] Cao M, Tang G A, Zhang F, Yang J (2013). A cellular automata model for simulating the evolution of positive–negative terrains in a small loess watershed. Int J Geogr Inf Sci, 27(7): 1349–1363

[6] Castillo C, Gómez J (2016). A century of gully erosion research: urgency, complexity and study approaches. Earth Sci Rev, 160: 300–319

[7] Chen L, Wei W, Fu B, Lü Y (2007). Soil and water conservation on the Loess Plateau in China: review and perspective. Prog Phys Geogr, 31(4): 389–403

[8] Congalton R G (1991). A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens Environ, 37(1): 35–46

[9] Coppola E, Tomassetti B, Mariotti L, Verdecchia M, Visconti G (2007). Cellular automata algorithms for drainage network extraction and rainfall data assimilation. Hydrol Sci J, 52(3): 579–592

[10] Dunkerley D (1997). Banded vegetation: development under uniform rainfall from a simple cellular automaton model. Plant Ecol, 129(2): 103–111

[11] Evans M, Lindsay J (2010). High resolution quantification of gully erosion in upland peatlands at the landscape scale. Earth Surf Process Landf, 35(8): 876–886

[12] Fairfield J, Leymarie P (1991). Drainage networks from grid digital elevation models. Water Resour Res, 27(5): 709–717

[13] Gan Z M (1990). Geomorphology and Soil Erosion Research in Loess Plateau of China. Xi’an: Shaanxi People’s Publishing House (in Chinese)

[14] Garbrecht J, Martz L (1997). The assignment of drainage direction over flat surfaces in raster digital elevation models. J Hydrol (Amst), 193(1–4): 204–213

[15] Harvey A (2001). Coupling between hillslopes and channels in upland fluvial systems: implications for landscape sensitivity, illustrated from the Howgill Fells, northwest England. Catena, 42(2): 225–250

[16] Harvey A M (2002). Effective timescales of coupling within fluvial systems. Geomorphology, 44(3–4): 175–201

[17] Hellweger F (1997). AGREE-DEM surface reconditioning system. Dissertation for PhD degree. University of Texas at Austin

[18] Hessel H (2002). Modelling Soil Erosion in a Small Catchment on the Chinese Loess Plateau. Dissertation for PhD degree. University Utrecht, 61

[19] Hessel R, van Asch T (2003). Modelling gully erosion for a small catchment on the Chinese Loess Plateau. Catena, 54(1): 131–146

[20] Hu G, Wu Y Q (2005). Progress in the study of geomorphic threshold theory in channel (gully) erosion. Mountain Research, 23(5): 565–570 (in Chinese)

[21] Jiang L, Tang G A, Zhao M W, Song X D (2013). Extraction and analysis of loess gully heads considering geomorphological structures. Geogr Res, 32(11): 2153–2162 (in Chinese)

[22] Jing K (1986). A study on gully erosion on the Loess Plateau. Scientia Geographica Sinica, 6: 340–347 (in Chinese)

[23] Jones R (2002). Algorithms for using a DEM for mapping catchment areas of stream sediment samples. Comput Geosci, 28(9): 1051–1060

[24] Knight J, Spencer J, Brooks A, Phinn S (2007). Large-area, high-resolution remote sensing-based mapping of alluvial gully erosion in Australia’s tropical rivers. In: Proceedings of the 5th Australian Stream Management Conference. New South Wales:199–204

[25] Lashermes B, Foufoula‒Georgiou E, Dietrich W E (2007). Channel network extraction from high resolution topography using wavelets. Geophys Res Lett, 34: L23S04

[26] Li B, Huang L, Feng L, Ma D (2014). Uncertainty of gully sediment budgets based on laser point cloud data. Transactions of the Chinese Society of Agricultural Engineering, 30(17): 183–191 (in Chinese)

[27] Li Y, Poesen J, Yang J C, Fu B, Zhang J H (2003). Evaluating gully erosion using ¹³⁷Cs and ²¹⁰Pb/¹³⁷Cs ratio in a reservoir catchment. Soil Tillage Res, 69(1): 107–115

[28] Li Z, Zhang Y, Yao W J, Zhu Q K (2011). Estimating gully development rates in hilly loess region of western Shanxi Province based on Quickbird images. Transactions of the Chinese Society of Agricultural Engineering, 28(22): 141–148 (in Chinese)

[29] Li Z, Zhang Y, Zhu Q K, He Y M, Yao W J (2015). Assessment of bank gully development and vegetation coverage on the Chinese Loess Plateau. Geomorphology, 228: 462–469

[30] Li Z, Zhang Y, Zhu Q K, Yang S, Li H J, Ma H (2017). A gully erosion assessment model for the Chinese Loess Plateau based on changes in gully length and area. Catena, https://doi.org/10.1016/j.catena.2016.04.018

[31] Liang C, MacKay D S (2000). A general model of watershed extraction and representation using globally optimal flow paths and up-slope contributing areas. Int J Geogr Inf Sci, 14(4): 337–358

[32] Liu D S (1985). Loess and Environment. Beijing: Science Press

[33] Liu Y B, Zhu X M, Zhou P H, Tang K L (1988). The law of hillslope channel erosion occurrence and development on loess plateau. Research of Soil and Water Conservation, 7: 9–18 (in Chinese)

[34] Luo L X (1956). A tentative classification of landforms in the loess plateau. Acta Geogr Sin, 22(3): 201–222 (in Chinese)

[35] Martínez-Casasnovas J (2003). A spatial information technology approach for the mapping and quantification of gully erosion. Catena, 50(2): 293–308

[36] Martz L W, Garbrecht J (1992). Numerical definition of drainage network and subcatchment areas from digital elevation models. Comput Geosci, 18(6): 747–761

[37] Metternicht G, Zinck J (1998). Evaluating the information content of JERS-1 SAR and Landsat TM data for discrimination of soil erosion features. ISPRS J Photogramm Remote Sens, 53(3): 143–153

[38] Na J M, Yang X, Li M, Ding H, Tang G A (2016). Progress in geomorphology research on young gully in loess plateau. Geography and Geo-Information Science, 32(4): 68–75 (in Chinese)

[39] Nyssen J, Poesen J, Moeyersons J, Luyten E, Veyret-Picot M, Deckers J, Haile M, Govers G (2002). Impact of road building on gully erosion risk: a case study from the northern Ethiopian highlands. Earth Surf Process Landf, 27(12): 1267–1283

[40] O’Callaghan J F, Mark D M (1984). The extraction of drainage networks from digital elevation data. Comput Vis Graph Image Process, 28(3): 323–344

[41] Oostwoud Wijdenes D J, Poesen J, Vandekerckhove L, Ghesquiere M (2000). Spatial distribution of gully head activity and sediment supply along an ephemeral channel in a Mediterranean environment. Catena, 39(3): 147–167

[42] Persendt F, Gomez C (2016). Assessment of drainage network extractions in a low-relief area of the Cuvelai Basin (Namibia) from multiple sources: LiDAR, topographic maps, and digital aerial orthophotographs. Geomorphology, 260(s 1–2): 32–50

[43] Poesen J, Nachtergaele J, Verstraeten G, Valentin C (2003). Gully erosion and environmental change: importance and research needs. Catena, 50(2–4): 91–133

[44] Rengers F, Tucker G, Mahan S (2016). Episodic bedrock erosion by gully-head migration, Colorado High Plains, USA. Earth Surf Process Landf, 41(11): 1574–1582

[45] Saunders W (1999). Preparation of DEMs for use in environmental modeling analysis. In: ESRI User Conference. Beijing: China Cartographic Publishing House, 24–30

[46] Schumm S A (1956). Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey. Geol Soc Am Bull, 67(5): 597–646

[47] Seginer I (1966). Gully development and sediment yield. J Hydrol (Amst), 4(4): 236–253

[48] Shruthi R B, Kerle N, Jetten V (2011). Object-based gully feature extraction using high spatial resolution imagery. Geomorphology, 134(3): 260–268

[49] Shruthi R B, Kerle N, Jetten V, Abdellah L, Machmach I (2015). Quantifying temporal changes in gully erosion areas with object-oriented analysis. Catena, 128: 262–277

[50] Shruthi R B, Kerle N, Jetten V, Stein A (2014). Object-based gully system prediction from medium resolution imagery using Random Forests. Geomorphology, 216: 283–294

[51] Sidorchuk A (1999). Dynamic and static models of gully erosion. Catena, 37(3–4): 401–414

[52] Song X D, Tang G A, Li F Y, Jiang L, Zhou Y, Qian K J (2013). Extraction of loess shoulder-line based on the parallel GVF snake model in the loess hilly area of China. Comput Geosci, 52(1): 11–20

[53] Stolte J, Liu B, Ritsema C, Van den Elsen H, Hessel R (2003). Modelling water flow and sediment processes in a small gully system on the Loess Plateau in China. Catena, 54(1): 117–130

[54] Strahler A (1963). The Earth. NY: Sciences, Harper and Row, 557

[55] Strahler A N (1957). Quantitative analysis of watershed geomorphology. Eos (Wash DC), 38(6): 913–920

[56] Tang G A, Xiao C C, Jia D X, Yang X (2007). DEM based investigation of loess shoulder-line. In: Geoinformatics 2007: Geospatial Information Science. Nanjing: SPIE, 67532E

[57] Tarboton D G (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resour Res, 33(2): 309–319

[58] Tarboton D G, Bras R L, Rodriguez-Iturbe I (1991). On the extraction of channel networks from digital elevation data. Hydrol Processes, 5(1): 81–100

[59] Turcotte R, Fortin J P, Rousseau A, Massicotte S, Villeneuve J P (2001). Determination of the drainage structure of a watershed using a digital elevation model and a digital river and lake network. J Hydrol (Amst), 240(3–4): 225–242

[60] Valentin C, Poesen J, Li Y (2005). Gully erosion: impacts, factors and control. Catena, 63(2): 132–153

[61] Wang L, Liu H (2006). An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int J Geogr Inf Sci, 20(2): 193–213

[62] Wolman M G, Gerson R (1978). Relative scales of time and effectiveness of climate in watershed geomorphology. Earth Surf Processes, 3(2): 189–208

[63] Wondzell S M, Cunningham G L, Bachelet D (1996). Relationships between landforms, geomorphic processes, and plant communities on a watershed in the northern Chihuahuan Desert. Landsc Ecol, 11(6): 351–362

[64] Wu Y, Cheng H (2005). Monitoring of gully erosion on the Loess Plateau of China using a global positioning system. Catena, 63(2–3): 154–166

[65] Wu Y, Zheng Q, Zhang Y, Liu B, Cheng H, Wang Y (2008). Development of gullies and sediment production in the black soil region of northeastern China. Geomorphology, 101(4): 683–691

[66] Yan S J, Tang G A, Li F Y, Zhang L (2014). Snake model for the extraction of loess shoulder-line from DEMs. J Mt Sci, 11(6): 1552–1559

[67] Zhang W, Yang X, Tang G A, Zhu S J, Li C L (2012). DEM-based flow direction algorithms study of stream extraction and watershed delineation in the low relief areas. Science of Surveying and Mapping, 37(2): 94–96 (in Chinese)

[68] Zhang Y, Yang S, Li Z, Li H, He Y (2015). Effect of narrow terrace on gully erosion in Northern Shaanxi loess area. Transactions of the Chinese Society of Agricultural Engineering, 31(7): 125–130 (in Chinese)

[69] Zhao J, Amerce M, Chen L, Govers G (2016). Vegetation cover and topography rather than human disturbance control gully density and sediment production on the Chinese Loess Plateau. Geomorphology, 274: 92–105

[70] Zhou G, Sun Z, Fu S (2016). An efficient variant of the priority-flood algorithm for filling depressions in raster digital elevation models. Comput Geosci, 90: 87–96

[71] Zhou Y, Tang G A, Yang X, Xiao C C, Zhang Y, Luo M L (2010). Positive and negative terrains on northern Shaanxi Loess Plateau. J Geogr Sci, 20(1): 64–76

[72] Zhu H C, Tang G A, Qian K J, Liu H Y (2014). Extraction and analysis of gully head of Loess Plateau in China based on Digital Elevation Model. Chin Geogr Sci, 24(3): 328–338

[73] Zhu X M (1956). Classification on the soil erosion in the loess region. Acta Pedologica Sinica, 4: 99–115 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature