Extraction of thalweg networks from DTMs: application to badlands

. To study gully spatial patterns in the badlands us-ing a continuous thalweg vector network, this paper presents methods to extract the badlands’ thalweg network from a regular grid digital terrain model (DTM) by combining terrain morphology indices with a drainage algorithm. This method will delineate a thalweg only where the DTM denotes a signiﬁcant curvature with respect to DTM accuracy and relies on three major steps. First, discontinuous concave areas were detected from the DTM using morphological criteria, either the plan curvature or the convergence index. Second, the concave areas were connected using a drainage algorithm, which provides a continuous, thick, tree-structured scheme. We assumed that these areas were physically signiﬁcant and corresponded to a gully ﬂoor. Finally, the thick path was reduced to its main course and vectorised to obtain a thalweg network. The methods were applied to both virtual and actual DTM cases. The actual case was a LiDAR DTM of the Draix badlands in the French Alps


Introduction
The badlands are usually defined as an intensely dissected landscape, where vegetation is sparse or absent, and on which agriculture is useless (Bryan and Yair, 1982). A combination of various weathering and erosion processes affects the side slopes and channels at different spatial and temporal scales. This results in the modelling of a particular landscape. To understand the geomorphological and hydrological functioning over highly eroded areas, different scale cartographies are needed. The badlands morphology consists of numerous gullies of variable size (from 10 cm to 100 m wide) organised in hierarchical networks. The high drainage density and steep slope make fieldwork and detailed cartography difficult. Consequently, remote sensing is the most appropriate source of data for mapping the badlands. Badlands mapping provides essential information for studying the spatial pattern of this environment. The quantitative description of the pattern contributes to a better understanding of the functioning and scale relationships (Horton, 1945). In addition, it allows for accurate comparisons between different eroded areas.
Representing and characterising the gully network from remotely sensed data is not straightforward. One possibility is to represent the gullies by one of their constitutive elements: thalweg networks. A thalweg is defined as the line joining the lowest points of the gully. In badlands topography, these networks are characterised by tree-structured topologies. To study the spatial pattern, digital indices describing the network geometry (fractal dimensions) and topology (as Horton ratios and Tokunaga's indices) must be computed for different extents and resolutions. These latter Published by Copernicus Publications on behalf of the European Geosciences Union.
indices are very sensitive to noise, especially for small tree structures (Dodds and Rothman, 2001), so the network representations must fit the data information to provide stable results.
Airborne or high-resolution satellite sensors supply the source data for the representation of the thalweg. We need an automatic and easily repeatable method with the least subjective parameter possible to extract the thalweg network. Even if manual image interpretation is possible and provides consistent networks, the result depends on the operator, so it is not reproducible (Molloy and Stepinski, 2007). Moreover, this operational and fastidious technique is difficult to apply to extensive areas. On the other hand, landforms can be extracted directly from digital terrain models (DTMs). Indeed, DTMs are spatial, digital representations of the relief that allow one to stake out morphological features. Nevertheless, extracting features such as thalwegs from a DTM raises some fundamental issues, such as which methods permit extraction of complete and stable networks from spatial data.
The most-used methods are based on regular grid DTMs. Various drainage algorithms, including D8 (O'Callaghan and Mark, 1984), kinematic routing algorithm (Lea, 1992), FD8 (Quinn et al., 1991) or D∞ (Tarboton, 1997), offer the possibility of computing drainage networks all over the grid surface. However, the transition from one drainage flow path to a vector thalweg network is not self-evident (Martz and Garbrecht, 1995;Tarboton and Ames, 2001;Turcotte et al., 2001). Most of the time, the process uses a unique contributing area threshold beyond which thalwegs, valleys or hydrographical networks are depicted. Because this criterion has no direct physical significance, a relevant question to ask is where to start the thalwegs upstream? Some authors have proposed using topological or physical reasoning to establish the threshold. Montgomery and Dietrich (1994) have suggested distinguishing dissected areas, where fluvial transport influence is dominant, from undissected areas, where diffusional, hillslope transport processes are dominant, to locate channel heads. Tarboton et al. (1991) rely on searching the smallest scale, measured in terms of support area, where the constant drop property breaks. This property assumes that the average drop is approximately constant along the Strahler order. The authors argue that this point represents the physical transition from channel erosion to hillslope erosion. However, because we plan to analyse the network topology and relate it to morphology, we cannot use this reasoning. The question, then, is how to get thalweg networks to fit, given the amount of information about landforms contained in a DTM? In other words, how should we compute a thalweg network from a DTM grid to minimise algorithm artefacts and altimetrical noise effects?
The morphological index computed from the DTM permits one to locate concave areas represented in the data. A large number of parameters exist that identify the terrain convergences (Wilson and Gallant, 2000;Peuker and Douglas, 1975;Yokoyama et al., 2002), and some have been used for thalweg extraction (Tarboton and Ames, 2001;Molloy and Stepinski, 2007;Tarolli and Dalla Fontatana, 2009).
In this paper, we present a method that combines an existing drainage algorithm and a morphological index to delineate the thalweg network. Compared with the plan curvature, the convergence index seems to be a good choice due to its ability to underscore the morphological features, such as valleys. This method attempts to extract thalwegs objectively, considering the significant landforms included in a DTM. This report is divided into four main sections: Sect. 2 presents the study areas and the data; Sect. 3 describes the method proposed to extract and assess thalweg network from DTMs; Sect. 4 presents the results obtained with both virtual and actual DTM; and Sect. 5 discusses the method and the result and compares them with another extraction method.

Virtual catchment and simulated data
First, several virtual catchment DTMs were used to assess the benefit of the proposed method. The use of virtual catchment DTM allows for the representation and comparison of networks extracted using different methods, while the initial shape of the valley is known and controlled. Our virtual DTMs are simple shapes embedded into tilted planes. The chosen shapes are either virtual valleys or idealised mathematical shapes. The virtual valleys are presented as various confluences (mainly varying angles), width and network patterns. The most interesting mathematical shape used is the corrugated sheet.

Draix badlands catchments and data
The study area is located in the Southern Alps in France, at approximately 6 • 20 ′ E and 44 • 08 ′ N. The test site represents one of the upper Durance catchment's badland areas. These badlands were formed from the highly erodible, black, Jurassic marls. In addition, the Mediterranean climate and precipitation served to intensify the active weathering and erosion processes. This particular test site, known as the Draix experimental catchments, is a field laboratory for the study of mountain erosion processes (Richard and Mathys, 1999). Five small catchments, from 1000 m 2 to 1 km 2 , have been monitored since 1984 and are now labelled, "Observatory for Research on the Environment" (ORE). For this work, we focused on a particular catchment called the Moulin. The Moulin catchment is spread over 0.08 km 2 with a low vegetation cover rate (40%) and a very high drainage density. The steep slopes, an average of 33 • , make the fieldwork difficult and the use of remote sensing data, such as the digital terrain model (DTM), essential.
Because of the ORE status and the work of several research teams, ample data is available on the study area.
Three DTMs of the Moulin catchment were used: a 50-cmresolution, a one-meter-resolution and a 2-meter-resolution. The three DTMs come from the same LiDAR survey (detailed in Bretar et al., 2009). They were derived from the point cloud resulting from the last echo. The DTM RMSE is around 0.17 m. The xyz points acquired with a DGPS (Jacome et al., 2008) permitted an evaluation of the spatial distribution of the DTM's altimetrical error.

Method
The presented method is intended for the extraction of one morphological feature, the thalweg network, reflecting the information contained in the DTM. More specifically, the objective is to delineate thalwegs from where the DTM denoted a concave curvature in the terrain where flow converges.
Our method consisted of three steps, based on the combination of a drainage algorithm and a morphological index (Fig. 1). The morphological index permits to locate convergent features. Two indices were tested: the plan curvature and the convergence index (Köthe and Lehmeier, 1994;Kiss, 2004). First, areas with significant convergence were derived from the DTM. We assumed that these areas corresponded to gully floor elements. But, as these areas can be discontinuous, a second step was needed. The discontinuous convergent areas were connected using a drainage algorithm to obtain a directed tree structure, considered as the entire gully floor features. The result was a thick network that in the third step was thinned and filtered to get a refined vector thalweg network.

Detection of significant convergence areas
Several morphological indices are available to distinguish the flow convergence and divergence areas from a DTM (Wilson and Gallant, 2000). Two recent studies (mentioned below) have proposed to integrate the curvature parameters derived from the DTM to extract the thalweg network. For both studies, the core idea is to flag convergence cells; however, they each compute convergence differently. One uses a proxy of curvature stemming from the Peuker and Douglas (1975) algorithm (Tarboton and Ames, 2001), and the other uses the tangential curvature (Molly and Stepinski, 2007). The curvature parameter integration aims to improve the ground-truth of the delineated thalwegs and introduce variability in the drainage density of a catchment. In our work, we compared the propensity of the plan curvature and the convergence index (CI) to extract convergence cells. The plan curvature is a popular terrain parameter that seems relevant, as it should directly inquire about the presence of a thalweg. Since the plan curvature and the tangential curvature have the same spatial distribution of convex and concave areas (Molly and Stepinski, 2007), these two parameters have the same impact in this work. The CI showed great potential to reveal the gully floor features in a previous work (Bretar et al., 2009).

Using the plan curvature
The plan curvature (or horizontal convexity; Wood, 1996) is defined as the second derivative of the surface elevation (Evans, 1972;Rana, 2006). The plan curvature (PC) measures the contour curvature and distinguishes the flow convergent and divergent cells of the DTM (Wilson and Gallant, 2000). The local surface is approximated by a quadratic function, and PC is found by differentiating the surface equation in a direction transverse to the slope (Zevenbergen and Thorne, 1987). The PC is computed using a 3×3 moving window.

Using the convergence index
The convergence index (CI) calculation principle (Köthe and Lehmeier, 1994;Kiss, 2004) differs from the PC. The CI is based on the aspect, which is defined as the slope direction (Zevenbergen and Thorne, 1987) and computed from the local surface function.
Using a 3×3 moving window (Fig. 2) for each external cell i, θ i is the angle, in degrees, between the aspect of cell i and the direction to the centre (i.e., the direction of the vector joining the centre of cell i and the centre of the window). The convergence index is the average of θ i minus 90 • (Eq. 1).
CI can range from −90 • to 90 • (Köthe and Lehmeier, 1994; Kiss, 2004). In the resulting grid, positive values represent divergent areas, and negative values represent convergent areas. Null values relate to areas without curvature as planar slopes.

Significant values
Some CI or PC values close to zero may not be significant: the DTM noise impacts the results. To highlight the significant cells, which should effectively correspond to gully floors from the proposed morphological criteria, a threshold is applied on each grid (CI T for the CI grid and PC T for the PC grid). Lashermes et al. (2007) and Tarolli and Dalla Fontana (2009) have proposed methods to apply a threshold curvature grid to identify hillslope-valley transition. However, the method proposed in this paper considered an objective threshold, dependent on the DTM noise. To objectively determine the thresholds for CI and PC, we used the DTM altimetrical error spatial distribution for the main parameter.
To determine the CI T , a 50×50-cell, tilted plane DTM (with a slope of 33 • , which is the Draix mean slope) was simulated with a noise respecting the altimetrical error spatial structure of the actual DTM considered (Fig. 3). Since spatially correlated errors were found and modelled on the DTM, the noise simulation process was based on geostatistical Gaussian field simulation, with the spatial covariance function and parameters inferred in the same manner as those used in Bretar et al. (2009). The computation of CI from the simulated DTM provided a threshold CI T , beyond which CI values can be considered significantly different from a plane landform. The Gaussian distribution of the CI values allows setting the CI T as the mean value minus twice the standard deviation (σ ). We accepted that 2.5 percent of CI values due to hazard on the DTM noise can be classified as significant (Bretar et al, 2009). Following this, a binary image attributed a value of one to significant convergent cells (i.e., for CI>CI T and 0 to the non-convergent cells). This same process was applied to determine the threshold PC T .
At the end of this step, the disconnected convergence areas were identified and considered as gully floor features. This identification was objective and repeatable. Moreover, the threshold computation method is valid for different resolution DTMs.

Calculation of a tree-structured gully floor area
As the objective was to compute a continuous, tree-structured thalweg network, the disconnected convergence areas needed to be connected to create a continuous area. Drainage algorithms allowed for upstream-downstream reconnection, as in Tarboton and Ames (2001) and Molly and Stepinski (2007). Most of the drainage algorithms implemented are based on local analyses and follow the same steps: (i) determination of flow direction; (ii) computation of flow accumulation; (iii) and delineation of the network using a unique Hydrol. Earth Syst. Sci., 14, 1527-1536, 2010 www.hydrol-earth-syst-sci.net/14/1527/2010/ contributing area threshold (CA T ). The concept of flow accumulation allows for the generation of a continuous flow path. The choice of a unidirectional flow algorithm (each cell has a unique output) versus a multi-directional flow algorithm (Tarboton, 1997;Quinn et al., 1992) ensures a treestructure.
The main ideas are: (i) to start the drainage area upstream only where a convergence is denoted and (ii) to force the flow to pass through the convergence areas detected. To achieve the connection, we integrated the terrain parameter into the flow accumulation computation. Usually, the flow accumulation relies only on the flow direction grid, but to take into account the convergence, the flow accumulation must be weighted by the threshold convergence grid. As a consequence, the cells having an accumulation value other than zero become part of the network. Cells that have a zeroaccumulation value are cut out of the network, as this represents areas without concavity. This method differs from the usual accumulation, which needs the application of a contributing area threshold to delineate a network.
Several algorithms are relevant to achieve the weighted flow accumulation computation; we compared the widelyused D8 algorithm (O'Callaghan and Mark, 1984) and the kinematic routing algorithm (KRA) (Lea, 1992). The major difference between D8 and KRA is the flow direction determination: the D8 directions are multiples of 45 • , whereas the KRA has a 1 • accuracy (Wilson, 2007). Given that Lea's flow directions, based on the aspect, are more flexible, we mainly used this algorithm. Lea (1992) routes the flow as if it were a rolling ball on a plane released from the centre of each cell. A plane is fit to elevations of cell corners, wherein the corner elevations are estimated by averaging the elevations of adjoining cell centre elevations (Tarboton, 1997). In contrast, the D8 method routes flow from each cell to one of its eight neighbours that has the steepest slope.
Besides, the DTMs often contain sinks (cells without an output due to an elevation lower than its neighbours). Whichever algorithm is chosen, the DTM needs to be filled before the direction is calculated to obtain a continuous flow path (Martz and Garbrecht, 1999;Planchon and Darboux, 2001) or the flow path may be interrupted. Nevertheless, the morphological indices are computed on the unfilled DTM so the terrain parameters fit to the DTM morphology representation.
As a result, we obtained a binary raster map of gully floors that are represented as a continuous tree structure.

Vector thalweg network extraction
The last step was to extract a vector thalweg network to perform a topological analysis of the network extracted. The network was obtained from the raster gully floors area: we consider that the central line represents the thalweg network. As a consequence, the thalweg sinuosity will be erased unless the area also meanders. In the Draix badlands, this assump-tion seems reasonable; the thalweg path fit well to the gully morphology.
Extracting a vector thalweg network from the map of gully floors involves thinning the gully floor (Molloy and Stepinski, 2007) and vectorising and filtering the thinned network. The conversion from a raster area to a vector network is not trivial. The thinning of the area results in a one-pixel-wide network with a course led by the shape of the area. This raster network is then vectorised. The vectorisation corresponds to the classical operation that joins the pixel centres. Because we want a continuous, tree-structure vector network, the vectorisation process integrates both the main information -the raster network -and the flow direction grid. The incorporation of the flow direction in the calculation creates a continuous network with an upstream/downstream consistency and without closed loops. Lastly, a filtering process removed the shortest thalwegs. We assumed that removing segments shorter than twice the resolution was a reasonable option. Using the DTM planimetric error and the badlands morphology, these segments seemed insignificant and unstable.

Network delineation and errors quantification
We decided that a quantitative assessment of the extracted network was necessary to validate the extraction method. While there is no standard method to assess the quality of an extracted network (Molloy and Stepinski, 2007), Heipke et al. (1997) have proposed an assessment based on a comparison of the extracted with a reference network. This global assessment focuses on network lengths. The reference network is a field-mapped network. We consider ground-truth network as reference because it comes from a direct vision that allows a multi-scale representation. The assessment addresses the following questions: how much of the extracted network is over-detected (extracted segments that cannot be matched with the reference) and how much is under-detected (ground-truth segments missing on the extracted network)?
The assessment, which focuses on network lengths, allows for an estimation of the delineation error based on an overlap of the networks (Heipke et al., 1997;Molloy and Stepinski, 2007). The overlap permits a quantification of the match between the two networks. More specifically, this method measures the length of the extracted network included in the reference domain and inversely the length of the reference network included in the extracted domain. The domain is defined as a buffer around the network. The width of the buffer is determined by the operator, and the choice depends on the accuracy required considering the DTM resolution and the type of terrain. Two simple parameters are computed from the matching: the false negative (FN), which is the length of the reference not included in the extracted network domain, and the false positive (FP), which is the length of the extracted network not included in the reference domain. In other words, FN and FP should correspond to the length of under-detection and the length of over-detection, respectively. For more significant results, the parameters are normalised by the total length of the reference network. FN and FP values close to zero represent minimal under-or overdetection, whereas the highest values (around 1) represent a maximal under-or over-detection. Two error indicators stemming from these parameters provide a quantitative assessment of the network extraction, which makes an objective comparison between several extraction methods possible. The first indicator, referred to as error balance (Eq. 2), reflects the geometric balance of the extraction between the over -and under -detection values. It ranges between 0 and 1; the closer the error balance is to zero, the more compliant the network is. The second indicator, referred to as the error sum (Eq. 3), summarises the total of the errors; the value interpretation is the same as the error balance.
The overlap method provides valuable information on the network's completeness and geometric accuracy (Heipke et al., 1997); however, some anomalies can disrupt the meaning of the calculated length (Fig. 4). For example, the length measure cannot be justified if the streams do not have the same orientation. In this case, the measured thalweg does not necessarily correspond to the buffered thalweg. Thus, the results of the assessment should be discussed before the method is validated.

Results
The proposed method was applied to both virtual and actual terrain cases. Figure 5 compares the results obtained using different extraction methods for two virtual DTMs. The classical extraction method used was the kinematic routing algorithm (Lea, 1992) with a unique threshold approach (detailed in Sect. 3.2).

Virtual valleys DTM
For the tilted plane DTM (Fig. 5a), we observed that with the classical method, the main stream was divided into two branches, whereas the algorithms using a morphological index (PC or CI) returned a single thalweg segment. This unique downstream line reveals that the morphological index-based methods provide a more realistic result in wide valleys. Moreover, with the classical method, artefact lines were created on the sides of the DTM, even though it is a perfect plane; these lines do not appear when a curvature Hydrol. Earth Syst. Sci., 14, 1527-1536, 2010 www.hydrol-earth-syst-sci.net/14/1527/2010/ 1 Fig. 6. Comparison of thalweg networks obtained using the three methods evaluated in this paper with a reference network of the Moulin catchment test area.
index is taken into account. The side artefacts are the direct consequence of the unique contributing area threshold: a thalweg is delineated once the threshold is reached. Differences were also observed between the PC and CI based methods. The CI-based network effectively fit the virtual valley: all the thalweg segments of the valley were represented and no over-detected segments were observed, whereas the PC-based network included under-detected segments. The corrugated surface (Fig. 5b) shows that when a morphological index was used, the flow only passed in the concave part of the wave, whereas with the classical method, the flow went up on the side where no curvature was represented. In this case, no distinctions were observed between the use of PC and CI.
In summary, it seems essential to integrate a morphological index for the thalweg's delineation. Furthermore, the CIbased method provided a unique solution that fits the virtual valley.

Comparing the extraction methods (PC vs. CI)
As with virtual terrain case, we compared the networks obtained by applying the CI-based method, the PC-based method and the network extracted with unique contributing area criterion methods. The comparison was mainly geometric. We present here the numerical results of the network assessment.
For the CI-based method, the threshold CI T resulting from the geostatistical calculation, considering the error spatial distribution of the LiDAR DTM, was −9.54 • . For the plan curvature, the PC T was −0.66. For the classical extraction method, the kinematic routing algorithm was used with a contributing area threshold CA T of 50 m 2 .
To understand which method would provide more accurate results (more compliant with the ground-truth), we compared the extracted networks to a field reference (Fig. 6). The reference network came from a manual mapping in the field. The quantitative assessment of the three networks extracted with different methods is shown in Table 1. According to the quantitative assessment (Table 1), the CI-based method resulted in the most compliant network. Among all of the networks, it had the highest accuracy statistic: its indicators (EB and ES) were the closest to zero. The two other networks' indicators gave higher values of EB and ES, which means that the extractions were unbalanced and possessed many errors. Still, the PC-based extraction provided better results than the classical method.
In details, we observed that for both the classical and PCbased methods, the unbalanced extraction came from high normalised FNs (i.e., high under-detection). On the other hand, the classical and PC-based methods resulted in low normalised FPs (i.e., low over-detection). In these two cases, most of the extracted network corresponded to the field, so the geometry was correct, but the network admits a lot of lack. It seems that for the same DTM, the PC-based method extracted less information than the CI-based method. Nonetheless, CI-based normalised FP is slightly higher than for the other methods. This can be explained by the error location. Figure 7 shows that the main over-detection recorded was located in the same place, a forested area.
Each network domain was defined as a two-meter-wide buffer, considering that twice the resolution of the base DTM approaches the data's planimetric noise. The more the buffer is reduced, the higher FN and FP. This quality measure helps to assess the geometric accuracy of the extraction. The effect that the choice of buffer width and the length of shortest thalweg segments to be removed during network extraction have on this assessment is not known.
Besides, we must also take into account the possible scale differences between the assessed network and the reference. Defining a scale for fieldwork can be difficult, and it can be heterogeneous depending on the terrain's accessibility. Indeed, some areas are more detailed in the field network. The resulting interpretation must incorporate this scale consideration.

Testing the CI approach on other resolutions
Because the CI-based method appears to provide more compliant results than the other two methods, we tested the method on other DTM resolutions to explore the method's sensitivity to data resolution.
To test the sensitivity of the CI-based extraction, we applied the method to two other DTM resolutions (50 cm and 2 m). The results for the test area of the Moulin catchment are presented Fig. 8.
The pattern difference between the three networks extracted using the CI-based method is due to the change in resolution because the main part (the highest order streams) is similar between the networks, while the lower orders grow as the resolution gets higher (closer to the field reference). This shows that the networks extracted using the CI-based method accurately express the DTM landform information, verifying the method's robustness.
The quantitative assessment of the three resolution networks in the test area reveals error indicators relatively close to zero (Table 2). Theoretically, as the resolution gets higher, the extracted networks should possess a constant over-detection (FP normalised ) and lower under-detection because the network pattern is getting closer to the field ref-erence. Consequently, EB and ES should be reduced. However, the results are not so clear. The 50-cm-resolution network possesses slightly higher EB and ES than the 1-mresolution network. Because the main parts of the network patterns look similar, we assumed that this inconsistency comes from an inaccuracy of network positioning due to the resolution. As previously, we used a 2-m-wide buffer to determine the network domains.

Discussion and conclusions
In this paper, we proposed a method to extract tree-structured gully floors and their corresponding continuous thalweg networks. While classical methods of extraction are based on a unique threshold and a flow routing algorithm (e.g., D8 or kinematic routing approach), the presented method extracts thalwegs only where the DTM expresses a landform. Thus, this method allows for a more robust geometric positioning that can be disrupted otherwise, compared with the reference network. The advantages of this method are the objectivity and the accuracy of the thalwegs that are extracted.
To construct this method, we used terrain parameters reflecting the convergence. Two main parameters were compared: the plan curvature and the convergence index. The results showed the enhanced ability of the convergence index to reveal the landform. The CI-based method provided a more realistic and complete network, especially for bare soils. This method is suitable for the badlands area or any other bare landscape. In other areas, such as forested environments, the results are less accurate.
This method is data-dependent and data-driven because it relies on DTM noise parameters. These parameters give an objective threshold for the delineation of significant flow convergence areas. This dependence requires that the DTM quality be evaluated first. The main limit of the CI-based method comes from the DTM noise determination. The spatial distribution of altimetrical error is valid only for baresoil areas, which explains the key source of over-detection (Fig. 7). We observed that thalwegs tended to be delineated between sets of trees. In addition, we used only one noisy tilted plane, corresponding to the average slope of the Draix. The threshold would be more accurate if it were determined using several planes with different slopes, in accordance with the badlands' side slope distribution, or using differencing Hydrol. Earth Syst. Sci., 14, 1527-1536, 2010 www.hydrol-earth-syst-sci.net/14/1527/2010/ error spatial structure with regards to terrain cover (e.g., vegetated/non-vegetated areas). The CI-based method was applied to three different DTM resolutions (2 m, 1 m and 50 cm) to test the sensitivity of the approach to the resolution. However, for different resolutions (far from the tested values), local neighbourhood convergence indices can limit network extraction and be unsuitable for network extraction in wide or rough valleys. In that case, a pyramidal approach, based on similar principles, should be developed. At a given resolution, the use of pyramidal moving windows for CI computation should smooth the resulting network and provide a more accurate result. In this work, for both the plan curvature and convergence index, the kernel used was a 3×3 moving window. Pirotti and Tarolli (2010) showed that the use of other kernels can provide better results. According to the authors, the best combination for detecting thalwegs is a curvature window size of 15×15 pixels. Because we intend to extract different-sized morphological elements, a pyramidal analysis would assist in testing different kernels for the convergence index.
The CI-based method applied to DTMs with a resolution on the order of few meters provided an acceptable thalweg network. The application of this method to very high resolution DTMs will need to be investigated in the future. It is likely that a resolution threshold exists, beyond which the convergence index will not be the optimal parameter for revealing a thalweg network. Above this threshold, the convergence may only express the terrain roughness. More specifically, in the badlands, the convergence may report the micromorphology of the terrain rather than the continuous thalweg network. The link between the intended object's size and the resolution of the base DTM must be considered for a physically significant extraction.