Intermittent Channel Systems of a Low‐Relief, Low‐Gradient Floodplain: Comparison of Automatic Extraction Methods

Intermittent floodplain channels are low‐relief conduits etched into the floodplain surface and remain dry much of the year. These channels comprise expansive systems and are important because during low‐level inundation they facilitate lateral hydraulic connectivity throughout the floodplain. Nevertheless, few studies have focused on these floodplain channels due to uncertainty in how to identify and characterize these systems in digital elevation models (DEMs). In particular, their automatic extraction from widely available DEMs is challenging due to the characteristically low‐relief and low‐gradient topography of floodplains. We applied three channel extraction approaches to the Congaree River floodplain DEM and compared the results to a channel reference map created through numerous field excursions over the past 30 years. The methods that we tested are based on flow accumulation area, topographic curvature, and mathematical morphology, or the D8, Laplacian, and bottom‐hat transform (BHT), respectively. Of the 198 km of reference channels the BHT, Laplacian, and D8 extracted 83%, 71%, and 23%, respectively, and the BHT consistently had the highest agreement with the reference network at the local (5 m) and regional (10 km) scales. The extraction results also include commission “error”, augmenting the reference map with about 100 km of channel length. Overall, the BHT method provided the best results for channel extraction, giving over 298 km in 69 km2 with a detrended regional relief of 1.9 m. Further, these analyses allow us to shed light on the meaning and use of the term “low‐relief landscapes”.


Introduction
Floodplains are low-relief and low-gradient landscapes (LRLGs) at valley bottoms and they vary in width from tens of meters to hundreds of kilometers (Dunne & Aalto, 2013;Latrubesse, 2015). Floodplains may develop expansive, along-valley surface flow networks (Fagan & Nanson, 2004) that include flow conduits ranging from topographically connected depressions to channels with well-defined banks (David et al., 2017;Kupfer et al., 2015;Rowland et al., 2009). These surface water conveyance systems are important because they facilitate material exchange between the main river and the floodplain surface (Day et al., 2008;Dunne et al., 1998;Park & Latrubesse, 2017). Further, an increased understanding of these systems provides deeper insight into biotic and abiotic processes that maintain high ecosystem productivity (e.g., Amoros & Bornette, 2002). Analyses of floodplain channel networks have been mostly limited to larger systems with perennial channels (e.g., Dunne et al., 1998;Park & Latrubesse, 2017;Trigg et al., 2012), although smaller intermittent systems have been receiving greater attention (e.g., David et al., 2017).
Floodplains may experience frequent low-magnitude inundation under high, but below bankfull river stages (Kupfer et al., 2015;Rowland et al., 2009). Floodplain inundation under these conditions occurs via "through-bank" channels and their networks (Park & Latrubesse, 2017;Rowland et al., 2009). A through-bank channel can be observed along the main river at low flow; it occurs as a channel that cuts through the main bank or levee and appears as a "hanging valley". However, at higher but subbankfull flow conditions through-bank channels are the hydraulic openings to expansive albeit low-relief surface flow networks that give rise to low-level floodplain inundation. Hence, floodplain material cycling results from ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2020WR027603
Key Points: • Comparison of D8, Laplacian, and bottom-hat transform for automatic channel extraction in a low-relief and low-gradient floodplain • Bottom-hat transform requires minimal preprocessing of DEM and yielded the highest accuracy of extracting channels with complex morphology • Congaree River floodplain is highly channeled with a drainage density of 4.3 km/km 2 Supporting Information: • Figure S1 • Figure S2 • Figure S3 • Figure S4 Correspondence to: H. Xu, xu79@email.sc.edu
With increased availability of light detection and ranging (lidar) data, automatic identification of subtle flow systems or other surface features becomes tractable. However, the characteristically low-relief attributes of floodplain landscapes complicate the development of automatic and unbiased channel extraction procedures, and the extraction of low-relief landscape features remains challenging (e.g., Passalacqua et al., 2012;Schwanghart et al., 2013). Given their important role in floodplain surface processes, there is an urgent need for a reliable and robust method of extracting subtle surface water conveyance features from floodplains and low-relief landscapes in general.
The extraction of channels and channel networks from DEMs has been studied for decades (e.g., Tarboton et al., 1991). A commonly used method is the "D8" and its modifications (Jenson & Dominque, 1988;O'Callaghan & Mark, 1984;Quinn et al., 1991;Seibert & Mcglynn, 2007;Tarboton, 1997). D8 is based on topographically determined flow direction and flow accumulation area. It has been successfully applied to channel network extraction with dendritic channel patterns, or in moderate-to high-relief landscapes (e.g., James et al., 2007;Persendt & Gomez, 2016;Turcotte et al., 2001). However, in LRLGs typical D8 assumptions may not apply because areas contributing to channelized runoff may not be well defined, and channel flows have been observed to experience abrupt changes in direction and velocity, and in some cases complete reversal (e.g., Rowland et al., 2009;Rudorff et al., 2014). Hence, the LRLG topography may not always dictate flow conditions (e.g., Alsdorf et al., 2007).
Alternatively, elevation-based methods may be suitable for channel analyses. For example, if channels and all surface water conveyance conduits have lower elevations relative to their immediate surroundings, then elevation thresholds can be applied for their extraction (Chirol et al., 2018;Fagherazzi et al., 1999). Inherent in this approach is the assumption that relief contributed by regional gradient is negligible. In a similar fashion, channel identification using relief, the difference between a DEM and a focally averaged surface of the DEM, may be suitable. In this manner, channels, sloughs, or connected depressions can be identified as contiguous negative relief features (Cazorzi et al., 2013;Liu et al., 2015). Liu et al. (2015) enhanced this approach through convolution of a Gaussian function for channel cross sections with local relief, in combination with adaptive thresholding. However, the assumption of Gaussian-like channel cross sections might limit such an application to landscapes with complex channel morphology (e.g., Chirol et al., 2018).
Channel extraction based on local topographic structure is another viable alternative. Local slope can augment channel identification and extraction, especially in the presence of well-developed banks. For example, Mason et al. (2006) extracted channels using multilevel and knowledge-based image processing where channel banks were identified as steeply sloped areas. Another topography-based approach relies on curvature where curvature in a DEM is defined as the second derivative of topography, and one can calculate the curvature along planform contours or along an elevation profile (Moore et al., 1991). The Laplacian is taken as the difference between planform curvature and profile curvature, and it accounts for changes in slope along both horizontal directions (Moore et al., 1991;Passalacqua et al., 2012). DEM curvature has been applied to channel extraction in high-relief (Lashermes et al., 2007;Pirotti & Tarolli, 2010) and low-relief (Fagherazzi et al., 1999;Schwanghart et al., 2013;Sofia et al., 2014) landscapes. In applying curvature analyses to DEMs with channels, an important assumption is that the bottom of the cross-channel profile has the largest curvature (Passalacqua et al., 2012). Thus, applying thresholds to DEM curvature analyses can aid in the identification and extraction of channels, especially in LRLGs.
Another approach is based on "mathematical morphology" (a technique of processing images or DEMs), such as the bottom-hat transform (BHT, Gonzalez & Woods, 2018). BHT is a fundamental method in that it only assumes channels are lower than the banks (e.g., Cho, Slatton, Cheung, & Huang, 2011;Cho, Slatton, Krekeler, & Cheung, 2011;Schwanghart et al., 2013). In fact, BHT requires only two user-defined parameters, the structuring element (SE) which describes the shape of the target feature, and the threshold value of the transformed DEM. The BHT of a DEM can highlight depressions narrower than the SE (e.g., Cho, Slatton, Cheung, & Huang, 2011;Cho, Slatton, Krekeler, & Cheung, 2011), and it is distinct from the elevation thresholding approach described above because it is less influenced by the regional gradient. The size of target features can directly translate to the size of SE. For example, in floodplains with both narrow channels and wide oxbows or lakes, a smaller SE can be chosen so that only narrow channels are highlighted while precluding oxbows.
Above we summarized several well-known techniques for extracting subtle or low-relief topographic features from a DEM. The objective of this study is to evaluate the performances of three widely accepted channel extraction methods, the D8, Laplacian, and BHT within a topographically complex floodplain terrain. Products from these efforts will help fill a knowledge gap related to the alignment and structural organization of landscapes with subtle relief, and they will shed light on potential pathways of floodplain inundation and connectivity during subbankfull conditions. Hereafter, we use the term "channel" to refer to any conduits or structural depressions on the surface. These features include conduits with well-defined banks, much wider slough-like features or contiguous depressions with poorly defined banks, or any confined pathway that can provide a focused transfer of water over the surface. Hence, we use the term in a general sense, and we do not mean to imply that channels reported here arise from the typical channel forming processes (e.g., Leopold & Maddock, 1953). Lastly, the term "low relief" is widely used throughout the literature but a clear definition is lacking. In other words, "low relief" means different things to different investigators. In section 5 we attempt to provide some clarity on the issue.

Study Area and Data
The study site is in the North American southeastern coastal plain and it is a 93 km 2 floodplain within the Congaree National Park (CNP), South Carolina, USA ( Figure 1a). CNP is bound by the Congaree River to the south/southwest, by bluffs to the north, and by the Wateree River to the east. The site is densely forested and optical remote sensing, regardless of spatial resolution, fails to detect or represent the fine-scale topographic features (Figure 1a), even during leaf-off months. However, lidar technology provides viable through-canopy observation of the subtle ground features. Analyses of recent lidar data establish an opportunity to greatly improve the understanding of the subtle floodplain surface features, and we contend that the methods presented here are applicable to other LRLGs, in floodplains and beyond.
We Elevation of the floodplain gradually decreases downstream from 35.0 to 23.5 m over 46 km river reach ( Figure 1b). Topographic flow pathways generally include narrow (intermittent channels) to wide conduits (infilled oxbows, see Appendix A). Most oxbows have widths comparable to the modern Congaree River, although a relict oxbow of~400 m occurs further downstream (Suther et al., 2018). Most oxbows are forested (Figures 1a) and remain as surface depressions that may intermittently retain water after pluvial or fluvial events. A total of seven lakes represent the unvegetated areas, and their depths range from 0.6 to 2.7 m ( Figure 1b); of the latter, only two are perennial lakes.
Local systems with perennial and intermittent channels have been reported by several authors (Kaase & Kupfer, 2016;Kupfer et al., 2015;Shelley, 2007;Wohl et al., 2011). Channels and oxbows of CNP are highlighted in the unpublished Congaree National Park surface feature map produced by John E. Cely (Appendix B). The Cely map was created by interpretation of color infrared aerial imagery (1989 to 1991) and detailed field observations over the past 30 years, and it continues to be updated. Cely's mapping was particularly arduous given the dense forest canopy that inhibits orientation (Figure 1a), a lack of affordable GPS at the start of the effort, the very subtle relief and frequently inundated terrain. Although the mapping effort was monumental, the map provides only a few opportunities to establish geodetic control. This notwithstanding, the map does provide detailed insight on landscape surface variability, surface feature shapes, and positioning. More importantly, the map is a valuable resource for the purpose of validating our automatic extraction methods.
The Cely map was scanned to the digital format and georeferenced to the DEM with 10 ground control points at artificial features. We digitized the Cely channel centerlines and scroll bars directly connected to the channel network with two criteria. (1) Channel length greater than 1 km (e.g., the width of the Congaree River levees), and (2) for Tom's Creek ( Figure 1b) which does not clearly emerge on the DEM, we used 10 cm × 10 cm Pictometry true color images from 2009 to help identify the position. While the Cely map does not include all channels apparent in the DEM, it nonetheless serves as a viable reference that reflects the complex channel systems. To validate our channel extractions that yield highly variable width, we also created channel polygons based on the Cely map in three subareas (large rectangles in Figure 1b).
The first subarea is adjacent to the Congaree River and has a typical near-bank topography where lateral migration of the main river has formed levees, oxbows, crevasse channels, and batture channels (Rowland et al., 2009;Saucier, 1994). Also, levees are typically associated with through-bank channels ( Figure 1b). The second site highlights Cedar Creek, the longest tributary channel as it enters the floodplain, and one might expect D8 to perform well here. It also shows the transition from a single to a multithreaded system. Further downstream, the third subarea contains a reticulating reach of a larger channel surrounded by a variety of narrower and shorter channels or depressions. This site shows a typical landscape of the interior floodplain where the elevation is lower than the bank/levee and experiences more frequent inundation.

Methods
Three independent methods were applied to the extraction of channels ( Figure 2) from the lidar DEM. The system of extracted channels was reported as polylines (D8) or as polygons (Laplacian and BHT), which account for the highly variable channel width. Section 3.1 introduces the D8 method. Sections 3.2 and 3.3 introduce the DEM transform by Laplacian and BHT, respectively, and the thresholding which highlights potential channel areas from the background. Section 3.4 introduces the morphological filtering and supervised classification used by both the Laplacian and BHT. Section 3.5 introduces the accuracy assessment of three methods.

The D8
Application of the D8 approach requires the DEM of the targeted watershed (HUC 03050110). The first step was to remove topographic sinks; this step is necessary because sinks have no outlet and create errors in the flow direction. We chose to fill the sinks with the priority-flood algorithm (Barnes et al., 2014), and the topographically driven flow directions were simultaneously obtained. Flow accumulation was then calculated based on flow directions. A 1 km threshold of channel length was chosen because the Congaree River levees, where D8 might work well, are up to 1 km wide. The corresponding catchment area should be 6.69 × 10 5 m 2 , following Hack's Law (Hack, 1957). Hence, pixels with flow accumulation area greater than 6.69 × 10 5 m 2 were extracted as channels and converted to polylines.

The Laplacian
The Laplacian operator (e.g., Gonzalez & Woods, 2018) of a function f is a linear and isotropic operator defined as Its discrete form on a raster image f (x,y) is In this study, f (x, y) is the elevation at (x, y) and a is the grid size equal to 1 m. The Laplacian of the DEM highlights areas with sharp changes in local gradient (e.g., the bottom of V-shaped channel cross sections). It only accounts for the local five pixels (center plus the two ordinate directions) and it is sensitive to DEM noise that can produce spurious results (Lashermes et al., 2007;Passalacqua et al., 2012). Therefore, the first step of this method was to smooth the DEM using a 2-D Gaussian filter with standard deviation (σ, typically referred as kernel size) as a user defined variable ( Figure 2). To extract channels of variable width, we applied σ ¼ 1, 5, 10, and 15 m. The smallest σ ¼ 1 m was chosen to retain the narrowest channels, and the remaining σ values have an increment of 5 m. Most wide channels in the study area can be extracted with σ ¼ 15 m. A value of σ larger than 15 m (e.g., 20 m) does not extract additional channels, and the extracted channels are substantially wider than the real channels. After smoothing the DEM, the Laplacian was calculated with Equation 2 at each pixel ( Figure 2) and thresholds were applied to the Laplacian result as determined by quantile-quantile plots (Lashermes et al., 2007). The threshold is the value where the quantile-quantile curve deviates from the straight fitting line near the mean value. Pixels with curvature value higher than the threshold were classified as 1 and the remaining pixels were assigned 0. The resulting binary images emphasize features with high curvature (presumably channels) across the various scales defined by σ.

The BHT
Mathematical morphology is a well-established branch of image processing, and it has been used for extracting image components such as boundaries and the basic structures, or skeletons (Gonzalez & Woods, 2018). All morphological operations require a structuring element (SE), a user-defined image shape similar to spatial convolution kernels. The two basic operations are dilation (⊕) and erosion (⊖, see Appendix C for detail). The higher-level operation closing is defined as dilation followed by erosion. The BHT is defined as the difference between the closing of an image and the original image. In this study, it highlights depressions narrower than the SE: A circular SE was used because channels can occur in any direction, and hereafter we use R to represent the SE radius (Cho, Slatton, Cheung, & Huang, 2011). The diameter of the SE represents the width of the widest channels to be extracted  and users should know the distribution of channel width before choosing R. A preliminary survey of channels in the DEM and in the field showed most CNP channels are narrower than 40 m, hence, we chose R ¼ 19 m. However, this radius was not wide enough to extract some wide reaches of Cedar Creek (Figure 1b), and thus we used a second value of R ¼ 49 m. Finally, we also used R ¼ 5 m to include a few very narrow (~3 m) channels.
We first attempted to process the DEM with a global threshold (e.g., Cho, Slatton, Cheung, & Huang, 2011;Otsu, 1979;Rishikeshan & Ramesh, 2018), but found this approach was not suitable because of the high variability in channel depth from less than 20 cm to over 1.5 m. Consequently, we invoked an adaptive thresholding (e.g., Liu et al., 2015). Therefore, adaptive thresholding is defined as where μ(BHT, 2R) is the focal mean of the BHT with the radius equal to 2R of the circular SE. Ideally, the constant c should be zero if the spacing between channels is over 2R. However, the proximity of channels and noise in the DEM together led to a positive optimal c. The optimal c values are 0.05, 0.1, and 0.2 m for R ¼ 5, 19, and 49 m, respectively (see Appendix C for detail).

Morphological Filtering, Supervised Classification, and Union of Multiple Scales
To remove noise and to fill in small holes in the binary images of sections 3.2 and 3.3, we applied morphological filtering (Rishikeshan & Ramesh, 2018; Appendix C). The filtered images were converted to polygons, and then a supervised two-class classifier was used to separate channel polygons from nonchannel polygons distinguished as occurring outside of the channels. Many supervised multiclass classifiers can serve this

10.1029/2020WR027603
Water Resources Research purpose and we used the support vector machine (SVM) approach, a classic machine learning classifier (Vapnik, 2013) available as open source software LibSVM v3.23 (Chang & Lin, 2011). Polygons classified as channels were retained and the remainings were discarded (Appendices D and E). For both Laplacian and BHT, the above extraction was repeated on different scales (Figure 2; Appendices D and E) similar to other studies (Isikdogan et al., 2015;Liu et al., 2015). For each method, classified channel polygons were combined as the union of all scales.
where FC i is the classified result of channel polygons on scale i, and FC is the combined polygons. The Laplacian method includes four scales and the BHT includes three scales. The boundary of each overlapping and connecting polygon cluster in FC was then updated by the dissolve tool in GIS and the lengths of the dissolved channels were recalculated. In order to compare these results with D8, we removed from the FC the isolated channels shorter than 1 km.

Accuracy Assessment
Extraction was applied to the whole 93 km 2 study area, but accuracy was only assessed in the upstream area (69 km 2 ), the area with higher quality lidar data (Figure 1b), hereafter referred to as "the floodplain". The downstream section (24 km 2 ), mostly blue in Figure 1b, had only 36% area with lidar ground returns due to the presence of standing water at the time of the lidar mission, and an even lower percentage from channels. We define the channels that match the reference map, miss the reference, or added to the reference map as "match", "omission", and "commission", respectively. Extraction accuracy is presented as the percentage of match channel length (or area) divided by reference channel length (or area). It is certain that any omission should be omission error; however, commission includes both commission error and unmapped channels by Cely.
Given the variable pattern of the D8 polylines due to the high-resolution DEM there are discrepancies between the result and the reference channels. Discrepancies in channel position of up to 10 m were allowed, and this accounts for the close but discordant channel positioning. Reference channels within 10 m of the D8 channels were considered as match, and the remaining reference channels were omission. The D8 channels outside 10 m buffer of the reference channels belong to the commission channels. Therefore, length-based accuracy can be calculated for the D8 result.
For the Laplacian and BHT the output format is polygon. We define the reference channels inside the polygons as match channels, and outside the polygons as omission channels. The commission length was calculated as the polygon-based length (see Appendix E) minus the match channel length. Therefore, length-based accuracy, omission, and commission can be calculated for the Laplacian and BHT. In the three subareas (Figure 1b), the overlapping area between the reference channel polygons and the extracted polygons is defined as match channels, and reference area not extracted is the omission. Therefore, area-based extraction accuracy and omission error can also be calculated. Additionally, we also reported the percentage of correctly classified training polygons referred to as training accuracy (Appendix E).

Channel Properties
Field observations and the DEM reveal that channels or structured conduits that may confine surface water flows vary substantially in width and depth (Figures 3a-3c). Also, images from the field show that typical channels are unvegetated (Figures 3d and 3e), although larger systems may develop small tree islands (Figure 3f). Cross-section profiles help illustrate the very subtle nature of channel expression and the underlying challenges for automatic, unbiased, and representative extraction (Figures 3g-3i). Channels show a range of bank development ranging from broad U-shaped swales (Figure 3g) to quasi-trapezoidal cross sections (Figures 3h and 3i). In Figures 3g-3i, channel banktop widths are about 7, 30, and 48 m, and maximum depths are 0.3, 1.2, and 1.6 m, respectively. These metrics give a relatively uniform width-to-depth ratio of 25:1 and they highlight the subtleties of channel relief.  Figures 4h and 4i). These results show that the Laplacian with a single σ value does not provide a complete description of the channel network, and a more effective approach requires a combination of smoothing kernels (e.g., Liu et al., 2015).
A major advantage of the BHT approach is that it works on the original data; it does not require DEM smoothing. For example, with R ¼ 5 m, the BHT extracted A-A′ (Figure 4j) but missed the wider channels (Figures 4k and 4l). However, with R ¼ 19 m, all three channels were successfully detected (Figures 4m-4o). This result arises because the lack of smoothing preserved many of the subtle channel features of this low-relief terrain, especially with the narrower conduits. Hence, with a judiciously selected R the BHT approach can successfully extract channels of various sizes. The BHT products include a representation of channels as inverted cross sections to aid in their visualization (e.g., Figures 4m-4o).

Extraction Results: The Floodplain
The 198 km of channels from Cely includes a range of conduit types, and these channels serve as the reference or "ground truth" for the automatic extractions ( Figure 5a; Table 1; Appendix B). These channels also topographically connect oxbows (Appendix A), giving a more expansive system of conduits.

Water Resources Research
The D8 approach extracted only 23% of the reference channels ( Figure 5b; Table 1), and overall D8 failed to represent the channel network. It performed better in the upstream areas along some through-bank channels which cut across the wide natural levee of Congaree River, most likely because the levees serve as a type of "upland" source area. Fewer through-bank channels were extracted in the downstream where levees are

10.1029/2020WR027603
Water Resources Research less developed. Omission errors occur throughout the floodplain and commission channels mostly occur to the north. Tom's Creek was extracted (Figures 1b and 5b) due to its 131 km 2 upland watershed area but shows as a commission channel due to the large spatial discrepancy.
The Laplacian effectively extracted channel networks which mimic the reference channels with a success of 71% (Figure 5c; Table 1). Omission channels include long segments of Cedar Creek and Tom's Creek, and channels occurring in the inner floodplain with a total of 58 km (Figure 5c; Table 1). Most omission channels occur in areas with shallow oxbows (<0.5 m depth) and overall lower local relief, which can be further influenced by the requisite DEM smoothing (Figure 4). Most commission channels occur in the upstream half of Figure 5c, where local relief is higher

10.1029/2020WR027603
Water Resources Research and oxbows generally have over 1 m depth. Some commission channels seem to be curved extensions of the match channels; later, we will see that those mainly include spurious channels along some inner boundaries of oxbows.
BHT extracted 165 km or 83% of the reference channels, achieving the highest accuracy of three methods ( Figure 5d; Table 1). Extracted channel networks also have similar patterns to the reference system ( Figure 5a), but with substantially fewer omission channels. The main omissions include the upstream and middle reaches of Cedar Creek and Tom's Creek. The 100 km commission channels occur throughout the floodplain. Conversely, the BHT did not extract spurious channels along oxbow inner boundaries, giving a cleaner channel representation of the channel system.
In summary, the D8 method did not perform well because the assumption that channels occur in areas with high topographical flow accumulation does not hold in this LRLG. Both Laplacian and BHT performed substantially better because they rely on local topography instead of flow accumulation. In the next three subsections we will thoroughly compare the challenges, advantages, and disadvantages of three methods in different topographic settings.

Extraction Results: Main River Edge
The three test subareas for extraction methods are depicted in Figure 1b and in greater detail in Figures 6-8, and accuracies are reported in Tables 2-4. The river edge test site (Figure 6a) slopes northward to the floodplain interior, from the well-developed modern levee to the lower central floodplain. A variety of channels exist, including through-bank, crevasse, and batture (tie) channels (Figure 6b). Through-bank channels usually have distinctive levees and, thus, have relatively higher local relief of up to 1.7 m (e.g., Dead River

Water Resources Research
Gut), whereas field surveys found the depth to be~3 m. Meanwhile, the reference channel widths range from 2 to 50 m (slough).
The D8 approach achieved 16% accuracy (Table 2) and extracted only five segments of channels ( Figure 6c). Three "first-order" channels merge at Big Lake and then enter Dead River Oxbow, and finally leave this subarea via a relict crevasse channel (Figure 6c). Omission channels are ubiquitous, and this result seems to be unrelated to channel depth because the deepest channel, Dead River Gut, is not extracted. Commission channels in the oxbows are also spurious based on field observations.
The Laplacian extracted 86% channels by length (Table 2). Note the accuracies by area are up to 8% different for both Laplacian and BHT methods (Tables 2-4) because they are more influenced by wider channels. Omission channels are mainly limited to a couple of channels and sloughs in the north (Figure 6d). Commission channels mainly include the inner boundaries of oxbows where curvatures are higher. Such behavior can be explained by the cross section of Figure 4f and is also confirmed in Appendix D, where inner boundaries have been extracted at all four scales. However, those commission channels are false positives because no such channel was confirmed by field observations. Additionally, Laplacian channels are typically  Water Resources Research wider than the reference channels due to smoothing and they appear to match best with channels of 15 to 30 m width.
BHT again performed best of three and achieved 91% accuracy by length ( Table 2). The extracted channels preserve the fine details (e.g., the small island inside Horsepen Gut, Figure 6e) and usually match the reference width for both narrow and wide channels; omission is limited to two channels. Commission channels include very narrow channels which were not extracted by Laplacian due to smoothing. For example, see the horizontal channel connecting Dead River Gut and Horsepen Gut, and northern part of A-A′ (Figures 6b  and 6e). Of greater significance however is that the BHT did not extract oxbow inner boundaries (Figure 6e), a distinct advantage for floodplains with lakes or relict oxbows that are not modern channels.
In summary, the D8 failed along the river edge with a relatively poor accuracy. Both Laplacian and BHT more effectively extracted the channels and in particular, performed perfectly on channels with high levees. BHT also identified more narrow channels not included in the reference map.

Extraction Results: Upland Edge
Cedar Creek is the largest tributary channel entering CNP and takes a meandering course for 26.7 km to Congaree River (Figure 1b). Upon entering the floodplain, the single-channel transitions to an anastomosing system (Figures 7a and 7b), and further downstream it enters an oxbow, Wise Lake. Cedar Creek has a variable width from 3 m near the bluff to 55 m. The narrow reach of Cedar Creek has very low local relief of 0.25 to 0.4 m.
Here, the D8 achieved a maximum accuracy of 45% relative to its respective test areas (Table 3). Most of the main trunk of Cedar Creek has been extracted with some alternating between match and commission channels due to the high resolution ( Figure 7c). Interestingly, the extraction did not follow Cedar Creek near the bluff but joined another channel that originated from the main river edge (Figures 5b and 7c); therefore, D8 missed a short reach of Cedar Creek. D8 is also not capable of extracting multithread channels. Nevertheless, D8 did not produce as many spurious channels as in Figure 6c, giving the overall better performance.
The Laplacian achieved 85% accuracy by length (Table 3) and effectively extracted the multithread channel system (Figure 7d). Omission mainly includes the upstream Cedar Creek due to the low local relief. The most obvious commission channel is the periphery of Wise Lake, and the extracted narrower channels appear considerably wider than the reference channels.
The BHT approach also achieved 85% accuracy by length (Table 3) with an overall extracted network similar to the Laplacian (Figure 7e). Omission is mainly limited to upstream Cedar Creek, same as the Laplacian. Commission channels include slightly more narrow conduits, and the whole of an oxbow lake because R ¼ 49 m covers the 85 m wide oxbow. Again, extracted narrow channels are closer to the reference widths, causing less commission area.
In summary, the notable improvement with the D8 over the main river edge likely resides with the upland source area where flow accumulation assumptions are more applicable. Both the Laplacian and BHT highlighted the complex anabranching channel network with an 85% extraction accuracy, although some omission occurs along Cedar Creek due to the very low local relief, for example, poorly developed banks along this reach.

Extraction Results: Inner Floodplain
Cedar Creek with a "loop channel" up to 40 m wide is the major channel in the third test site (Figure 8a).
Lidar data indicate that only 21% area in the reference channels had ground returns (Figure 8b), therefore producing a DEM with locally diminished quality. For example, Cedar Creek has an overall patchy appearance likely due to lidar voids (Figure 8a). We contend that this low laying region of Figure 8 and the low-laying southeast part of CNP ( Figure 1b) are likely inundated more frequently than other locations.
The D8 achieved 8% accuracy and performed worst of all test sites and among all methods (Table 4). Only one segment in Cedar Creek and several short segments were extracted and substantial omission error occurs (Figure 8c). Commission channels include a few subtle channels and sloughs not included in the reference map. However, the poor performance is not linked to data quality, because the downstream part of Cedar Creek with higher ground return density and better DEM quality was not extracted (Figure 8a). Instead, the lack of wide natural levee near this site does not satisfy the flow accumulation assumption; therefore, D8 performs worse than on the upstream river edge (Figures 5b and 6c) and on the upland edge ( Figure 7c).
The Laplacian achieved 44% accuracy by length (Table 4). Only three segments of Cedar Creek and a set of narrow channels have been extracted (Figure 8d). The remaining trunk of Cedar Creek, the loop channel, and some narrow channels are omitted. The patchiness of the channels negatively influenced the Gaussian smoothing and therefore, channels became fragmented after thresholding. Those short channel segments were then discarded by the classification. Therefore, the DEM quality decreased the extraction accuracy of Laplacian. Moreover, an effect of DEM quality is apparent with the overall successful extraction of the downstream reach of Cedar Creek (Figure 8d), an area with a more representative DEM.
The BHT achieved 69% extraction accuracy by length and outperformed D8 and Laplacian in this challenging site (Table 4). It extracted most of the trunk of Cedar Creek and half of the loop channel ( Figure 8e). The most obvious omission error occurs along the western part of the loop, where local relief is only 0.4 m. The BHT result also suffered from the patchiness of the DEM; however, it does not require smoothing and the channel segments are still interconnected after thresholding and filtering (Figure 8e). The BHT also extracted the most commission channels of all method (Figure 8e; Table 4).
Overall, this test site demonstrated a typically challenging extraction case. The patchiness in channels caused by inundation can be expected in similar floodplain interior. Such prevailing conditions are particularly challenging for most extraction methods and they created the largest difference in our test results. Further, the products from all sites highlight the effectiveness of Laplacian and BHT in channel extraction. Both methods perform well on channels with well-developed banks and higher local relief. On the other hand, the representativeness of both declines where the local relief is typically less than 0.4 m. This notwithstanding, the BHT still outperforms the Laplacian because it can consistently highlight the fine detail of channels and extract more channels from degraded DEM quality. Overall, the BHT is suited for extracting intermittent channels from LRLGs. The D8 method is not suitable for channel extraction in floodplains because the concept of flow accumulation area is generally not applicable (Jones et al., 2008). Given the overall weak performance, the D8 is not further discussed.
Besides the extracted narrow oxbows (e.g., Figure 7e), four important subtypes of commission channels have been observed. The first subtype includes channels with steep banks that are absent in the reference map; these are not commission errors but true channels (e.g., see ellipses in Figures 6-8). The second subtype includes connected scroll bars (irregular polygons in Figures 6-8) with overall higher relief that may disrupt or connect channel networks and thereby influence floodplain connectivity (e.g., Trigg et al., 2012). The third subtype is connected sloughs (rectangles in Figures 6-8) where sloughs may function as channels if both ends are connected to channels (e.g., the slough in northern Big Lake in Figure 6). The last subtype includes pond-like depressions connected to channels (circles in Figures 6-8). These features likely serve as water storage sites instead of conduits and should be considered as commission errors. Overall, the topographic depressions associated with the channel networks complicate the commission channels; nevertheless, a substantial number of those channels are true channels. These subtle channels, either commission or reference channels, represent an important facet of surface water hydrology, especially at this study site (e.g., Kupfer et al., 2015).

Discussion
Previous studies of channel extraction often merge the results of two or more distinct methods. For example, Fagherazzi et al. (1999) extracted intertidal creeks as the union of results from the application of thresholds in elevation, and in curvature. Later, Lashermes et al. (2007) extracted channels by applying a threshold of curvature and in the derivative of slope direction. The intersection of the two extractions was taken as the final product, and it significantly reduced the noise within the individual results. Based on our work above we propose that it is not necessary to merge results of different methods to produce a robust and representative network of subtle channel features. This is important because a single method reduces the overall level of requisite DEM data conditioning and the number of arbitrary decisions, for example, as with establishing a threshold or filtering of extraction results.
On the other hand, merging results of different scales is necessary to achieve the best extraction result, especially for the Laplacian method (Figure 4; Appendices D and E). For the BHT, the small-scale R ¼ 5 m had negligible improvement of the accuracy over the combined result of R ¼ 19 m and R ¼ 49 m (Appendix E), though it added some very narrow channels to the result (Appendix D). The performance achieved by a single-scale R ¼ 19 m approximates the merged Laplacian result of four scales (Appendix E). This agrees with the observation in Figure 4, that sometimes a single scale BHT may suffice with a carefully chosen SE that reflects the widest channel width. The supplement of R ¼ 49 m mainly includes the connectivity of a few channels wider than 40 m (Appendix D). Given the strong performance of the BHT method in all test sites above and the straightforward settings of parameters, we propose that BHT is more suitable for LRLG channel extractions than the Laplacian method.
Many studies of floodplain channel processes focused on the wide perennial channels visible in optical remote sensing of large floodplains (e.g., Trigg et al., 2012), but only a handful have focused on the small, intermittent channels, and channel processes as reported here (Czuba et al., 2019;Pinel et al., 2020). These studies found that incorporating the subtle channels into numerical simulations substantially improves the accuracy of model products. On the other hand, few have addressed inundation and sedimentation processes due to through-bank fluxes (e.g., Czuba et al., 2019;Dunne et al., 1998). We contend that part of the problem is that through-bank channel networks typically are not well inventoried or characterized in floodplain analyses. Of course, water conveyance through these channel systems becomes less relevant when the free surface elevation is higher than the local topography.
Accurate extraction of surface water conveyance features is important because floodplain inundation from the main river is not necessarily limited to times when the river stage is above bankfull (Park & Latrubesse, 2017;Rowland et al., 2009). Indeed, opportunities for inundation at high but below bankfull conditions occur with the aid of the channel networks presented here. Therefore, a quantitative and accurate approach is needed to highlight through-bank channel architecture to address questions related to, for example, ecosystem hydraulic connectivity (e.g., Miranda, 2005) or floodplain material exchange processes (Dunne et al., 1998).
At CNP the through-bank channels are distinct because they extend from the main river as well-defined conduits with up to 3 m of incision at the mouth, and often they form their own levees ( Figure 6). Through-bank inundation has been reported in field studies (Kaase & Kupfer, 2016) and shown in 2-D simulations (Kupfer et al., 2015;Meitzen, 2011), but detailed analyses of through-bank fluxes and frequency have not been reported. The Cely map and examination of the extraction results together helped identify 36 through-bank sites where channels cut through the bank and are of variable length but may extend across 10.1029/2020WR027603 Water Resources Research the floodplain (Figure 5d). Also, the through-bank channel mouths have a relatively uniform spacing of 1 km along the 35 km river reach and they typically occur at the outer bends of the main river.
Overall, the BHT approach shows that the total channel length is at least 298 km (match + omission + commission, Table 1) in the 69 km 2 high quality lidar area, giving an average drainage density (Dd) of 4.3 km/ km 2 . This value is far higher than Dd measurements published for other floodplains. For example, in the Middle Amazon River Dd is up to 2.0 km/km 2 based on 30 m Landsat imagery (Trigg et al., 2012). Likewise, the Middle Paraná River has a Dd of up to 3.2 km/km 2 (Drago, 2007), and for the Cosumnes River, CA, Dd was 2.6 km/km 2 (Florsheim & Mount, 2003). One important reason that we report such a high Dd is the use of the through-canopy, high-resolution lidar DEM whereas data sources for the other studies were much coarser or older. Given the increasing prevalence of high-quality lidar data on floodplains, we contend that most subtle channels can be identified by the BHT method.
A review of the literature reveals that there are several studies focused on extracting surface features from low-relief landscapes. However, upon closer inspection of the respective study areas it becomes clear that the term "low relief" means different things to different authors. In other words, "low relief" is implicitly understood, but not consistently defined. For example, Qin et al. (2012) described a 4 km 2 site with 60 m total relief and 2.4°gradient (4.2 × 10 −2 m/m) as "low relief". A common misconception with the use of "low relief" is that researchers tend to define the term on the local scale, and then even with different definitions of local scale (e.g., Fagan & Nanson, 2004;Passalacqua et al., 2012). In fact, if one assumes that "local relief" is defined as the total elevation range within a 5 m radius circle, it becomes apparent that the local relief is less than 3 m in CNP and in smooth areas with moderate slopes up to 16.7°. However, as the area included for relief analyses expands, the changing rates of relief will vary between the upland and the floodplain. This is partly due to the regional gradient, but largely an artifact of the scale of observation. The difference between high relief and low relief is ultimately revealed by analyses applied to larger observational scales, not the local scale. Therefore, descriptions of relief should be accompanied by quantitative assessments of relief at different scales. For our study site the landscape has a total relief of 10.5 m over a floodplain of 69 km 2 , and reduces to 1.9 m within the 1st to 99th percentile after detrending the DEM. In order to resolve uncertainty regarding use of the term "low-relief" we propose that researchers report the area over which relief was assessed, and the value of relief within that area, across the entire study site. For the Congaree River floodplain the landscape is low-relief based on the observation of less than 3 m of mean relief with 500 m radius circles that cover the entire 93 km 2 study site.

Conclusions
The extraction of terrain features from low-relief and low-gradient landscapes is challenging because the features of interest typically are very subtle, especially in coastal plain floodplains. The bottom hat transform (BHT) approach to channel extraction in such terrain consistently outperformed approaches based on Laplacian of the DEM and cumulative contributing area. Moreover, the BHT is advantageous in that DEM preprocessing is not required for channel extractions. The BHT extraction results both replicated and augmented the channel reference map used as ground truth. We determined that the area has 298 km of channels in 69 km 2 , giving a drainage density of 4.3 km/km 2 . This high drainage density can be expected to facilitate floodplain inundation during high, but below bankfull flows via through-bank channels. Moreover, a literature review of feature extraction from "low relief" landscapes reveals that the term "low relief" is not well defined. We propose that a low-relief landscape be defined as one having less than 3 m of relief within a 500 m radius.

Appendix A
Appendix B

Appendix C: A Brief Introduction to Morphological Image Processing With an Example
Originally, the morphological operations were defined on binary images. Dilation expands the boundary of an image feature and erosion shrinks the boundary of the image feature, both by the SE with respect to its origin. The operations can expand to grayscale images: Here (x,y) is a pixel of a grayscale image f and (s,t) is a pixel of a symmetrical SE with respect to its origin. The above equations also apply to binary images because the target features have the value of 1 and the background is 0. Therefore, dilation and erosion can be easily implemented by focal statistical tools in GIS software ( Figure C1a-C1b). The higher-level operator closing is defined as dilation followed by erosion. Closing can fill in gaps in binary images or fill dark areas (depressions) in grayscale images, and only fill areas narrower than the SE ( Figure C1c). In contrast, opening is defined as erosion followed by dilation. Opening removes image components in binary images or removes bright areas (peaks) in grayscale images, both narrower than the SE. BHT is defined as the difference between the closing of the DEM and the original DEM ( Figure C1d). See Chapter 9 of Gonzalez and Woods (2018) for more details of morphological image processing.
In the thresholding of BHT, the optimal constant c was determined by trial and error in the range of 0.0 to 0.5 m given the overall low relief. An increment of 0.05 m was chosen. Effort was made so that shallow channels are extracted but not necessarily connected to adjacent depressions ( Figure C1e). The optimal c should highlight most channels. Small holes and gaps less than 5 m were accepted because they can be filled through Figure A1. Map of oxbows by Shelley and Cohen (2007). Figure B1. Map of Congaree National Park created by John E. Cely.

10.1029/2020WR027603
Water Resources Research morphological filtering ( Figure C1f). For morphological filtering, we applied an opening followed by a closing both with R = 1 m, and then followed by a closing with R = 2 m. The last step was not applied to the BHT R = 5 m because it connected areas that are not necessarily channels. The binary raster was then converted to polygons. Holes in the polygons of up to 200 m 2 and with less than 50% ground return area were considered as channel (wet) areas, and then were filled.

Water Resources Research
For the SVM classification, 50 channel polygons and 50 nonchannel polygons were selected as training samples, for each scale. Similar samples were used among different scales when possible. The two geometric features used for the classification include polygon length to width ratio and the area ratio between the polygon and its minimal bounding circle. We improved the equation for calculating length of polygons by Cazorzi et al. (2013) from four scanning directions to eight, to determine the channel length and the polygon length to width ratio. Specifically, C-SVM (Chang & Lin, 2011) was used for the classification. The optimal parameter C in C-SVM (not to be confused with c in BHT thresholding) was determined with the fivefold cross validation of the training samples. The default value of C is 1 and was multiplied or divided by 10 to train the model, and to calculate the accuracy of cross validation. The C value with the highest accuracy was chosen as the optimal parameter. Then the classification model was created with C and all (100) training samples. Training accuracies of the 100 samples were reported in Table E1. The classification model was then applied to all polygons. The outcome of the classification is either 1 (channel) or 0 (nonchannel). The extracted lengths and accuracies were based on the198 km of reference channels (Table E1).

Data Availability Statements
The lidar data set is available at NOAA Digital Coast website (https://coast.noaa.gov/htdata/lidar1_z/ geoid12a/data/4815/). Cely Map is accessible in Appendix B and a hard copy with higher resolution can be purchased from Congaree National Park (https://www.nps.gov/cong/contacts.htm). The oxbow map is accessible in Appendix A and can be obtained from South Carolina Geological Survey (https://ngmdb. usgs.gov/Prodesc/proddesc_94282.htm). Four Pictometry images collected in 2009 were used to help digitize Tom's Creek. The scenes are SCRICH054070OrthoSectorTile3, SCRICH055069OrthoSectorTile4, SCRICH055070OrthoSectorTile1, and SCRICH055070OrthoSectorTile3. The images can be purchased from the vendor EagleView (https://www.eagleview.com/product/imagery/). They may also be freely obtained by institutions through a license agreement with vendor. Note. The threshold of 1 km channel length was not applied, different from Tables 1-4. Therefore, more commission channels are included in those results; see Appendix D for example.