Abstract

High-resolution digital elevation models (DEMs) offer opportunities for channel network extraction due to its representation of realistic topography. Channels are generally surrounded by well-defined banks that have a distinct signature in the contour lines. Contour curvature is one of the important topographic attributes usually used for channel head identification; however, the curvature at channel heads may vary considerably between and even within watersheds. Therefore, uncertainty exists in the extracted channel heads due to the specified curvature threshold. In this paper, the locations of channel heads in 14 small mountainous watersheds are obtained using a nonparametric method based on the shape of contour lines generated from DEMs with a spatial resolution of 1 m, and the channel head curvature is computed from the extracted channel heads. The spatial distributions of the channel head curvature in these 14 watersheds have been analyzed, and another two watersheds with field-mapped channel heads are selected for validation. The results indicate that: (1) the channel head curvature is sensitive to the local terrain and varies within and between watersheds; (2) the Gamma distribution effectively fits the spatial distribution of the channel head curvature in all the selected watersheds; and (3) constant threshold-based methods for channel head identification gain significant location errors even within a single watershed.

INTRODUCTION

A channel network is an important geomorphologic attribute and hydrologic feature to control a watershed-scale hydrologic response such as rainfall-runoff, flow path, and flood hydrograph. Channel heads are the beginning of the channel network where water and sediment begin to flow and mark a boundary between hillslope processes and fluvial processes (Julian et al. 2012; Bierman & Montgomery 2014). The locations of channel heads control the total length of channels in stream networks, which determine how rapidly water and sediment move through a watershed. Their positions provide a direct insight into the control on source areas, drainage density, and water travel times for the watershed (Dietrich & Dunne 1993; Gomi et al. 2002). Thus, locating channel heads correctly is crucial for hydrological modeling, the transport of water, sediment, nutrients, and carbon, landscape evolution, as well as hydrological or geomorphological processes within individual watersheds.

The extraction of channel network has been of great interest to geomorphologists and hydrologists in the past decades. Several traditional methods for automatic channel network extraction have been developed in the previous studies. These methods usually need to specify a constant threshold for channel head identification in a watershed, and the threshold is usually provided for merely contributing area A (O'Callaghan & Mark 1984; Band 1986; Tarboton et al. 1991), slope–area relationship (Montgomery & Dietrich 1988; Willgoose et al. 1991; Dietrich et al. 1992; Dietrich et al. 1993; Ijjasz-Vasquez & Bras 1995), and Strahler order (Horton 1945; Strahler 1952; Strahler 1957) of surface flow paths extracted from the digital elevation model (DEM) (Peckham 1995). Montgomery & Dietrich (1988) suggested that the source area and slope downstream of the channel head indicated a strong inverse correlation and that drier regions tended to have larger drainage areas for the same gradient. The distance errors between extracted channel heads using threshold-based methods and the mapped channel heads through field investigations are significant in many studied watersheds (Tarolli & Dalla Fontana 2009; Orlandini et al. 2011; Clubb et al. 2014).

The rapid development of the airborne light detection and ranging (LiDAR) technique with its great advantage of obtaining high-resolution DEMs has brought new opportunities for using direct topographic attributes such as curvature (Lashermes et al. 2007; Sofia et al. 2011; Pelletier 2013), openness (Sofia et al. 2011), slope direction (Lashermes et al. 2007), and wet channel segments (Hooshyar et al. 2015) to extract the channel network. These methods are more complex and computationally expensive compared with traditional methods using a drainage area or slope–area threshold (Clubb et al. 2014; Hooshyar et al. 2016). However, topographic attributes obtained from high-resolution DEMs represent a detailed description of surface geometric properties and terrain convergences, which largely improve the accuracy or precision of channel network extraction (Tarolli & Dalla Fontana 2009; Sofia et al. 2011). Particularly, the curvature is often used to describe the physical characteristics of a drainage basin, which assists in understanding erosion and runoff processes and has been widely used for identifying the locations of channel heads (Passalacqua et al. 2010; Pelletier 2013). The curvature-based method for channel head identification originates from the traditional contour crenulation method (Strahler 1953; Marie 1957), in which the channel begins when the landscape is convergent and the convergence is measured by the curvature. Similar to the methods based on the upslope area threshold and the slope–area threshold, the channel head curvature threshold needs to be specified for the curvature-based method.

The channel head curvature may vary spatially from watershed to watershed. Moreover, the magnitude of the curvature is dependent on the filtering method used for smoothing DEM and the resolution of the DEM (Tarolli & Dalla Fontana 2009). Passalacqua et al. (2010) developed a geometric framework for channel network extraction using nonlinear geometric filtering (scale-adaptive filtering) and showed the impact of filtering methods on the calculation of curvature from DEMs through analysis and comparison between the Gaussian filter and the Perona–Malik filter. Grieve et al. (2016) explored the relationship between the curvature and the grid resolution as well as the estimation of the hillslope sediment transport coefficient and concluded that curvature values are generally insensitive to a grid resolution, particularly in landscapes with broad hilltops and valleys; however, curvature distributions changed and became obviously condensed around the mean. Therefore, the curvature threshold may vary with watersheds, the spatial resolution of DEM, and the filtering methods used for smoothing DEM. For that purpose, the channel head curvature threshold is identified from the curvature statistics of individual watershed (Lashermes et al. 2007; Tarolli & Dalla Fontana 2009). Although the thresholds identified from curvature statistics vary among watersheds, they are usually regarded as constant when locating the channel heads within a single watershed. However, the channel head curvature in a watershed may vary from tributary to tributary (Clubb et al. 2014), due to the significant difference in local topography, soil property, and land cover.

Due to the potential limitation of the constant channel head curvature threshold, the following questions emerge: (1) How much spatial variability does channel head curvature have within a watershed? (2) What is the appropriate distribution function for describing the spatial distribution of the channel head curvature within a watershed? (3) How does the spatial distribution vary among watersheds? To address these questions, the curvature of channel heads needs to be quantified. At present day, mapped channel heads can be obtained from the field investigation. For example, Clubb et al. (2014) compiled a series of channel head data in the field, both from their own field mapping and from the work of Julian et al. (2012), which provided quite comprehensive data for channel head identification. However, most channel heads, in general, are still unmapped and some of channel heads cannot be accessed through investigation because of restriction by the surroundings. Therefore, the observations are relatively rare for distribution analysis within and between watersheds due to the limited datasets of field-mapped channel heads. In this study, the channel heads are identified from high-resolution topographic data using a nonparametric method (i.e., a curvature threshold for channel head is not required), which was developed by Hooshyar et al. (2016) and demonstrated to be well performed in terms of identifying small channels; moreover, the average and standard deviation of distance errors between mapped and predicted channel heads were better or more comparable to existing methods, such as DrEICH proposed by Clubb et al. (2014) and Pelletier (2013). In this method, the transitions from unchannelized to channelized sections are defined as channel heads, which are extracted using the shape of contours which reflect the entire cross-sectional profile including possible banks. Therefore, this method is conducted based on the fact that channels and valleys have a fundamental geomorphologic difference, and channel head identification is performed independently for each tributary. Due to its physical basis, this method is robust and ensures the accuracy of channel head extraction.

In this paper, the channel heads in 14 small mountainous watersheds are extracted using the nonparametric method developed by Hooshyar et al. (2016), the curvature for each channel head is computed, and the distribution of the channel head curvature within each watershed is analyzed based on L-moment diagram. The rest of this paper is organized as follows: Section ‘Study Area and Data Source’ introduces the study sites and the topographic data used for extracting channel heads. Section ‘Methodology’ briefly describes the method for calculating curvature, extracting channel heads, and statistical analysis including the L-moment diagram and the probability plot correlation coefficient. Results, Discussions and Conclusions are presented in the following sections.

STUDY AREA AND DATA SOURCE

All the study watersheds are in the mountainous area in the United States as shown in Figure 1. The base-map, with a limited resolution, is added from the online database of ARCGIS. These 14 watersheds are selected based on the availability of LiDAR data and the fact that human activities such as urbanization are minimal in these watersheds. The name, location, drainage area, and climate aridity index (AI) for each watershed are shown in Table 1. The drainage area of these watersheds varies from 5.8 to 21.96 km2. The climate AI, defined as the ratio of mean annual precipitation to mean annual potential evapotranspiration, is obtained from the CGIAR-CSI website (http://srtm.csi.cgiar.org/) and varies from 0.52 to 2.25. The mean annual precipitation and temperature data from 1981 to 2010 are downloaded from the PRISM website (http://prism.oregonstate.edu/normals/), and the mean annual precipitation and temperature vary in the range of 724.23–2,436.84 mm and 2.72–24.33 °C, respectively. The patterns of parent material have a great effect on channel development; in these studied areas, there are mainly three types of parent materials, namely colluvium, residuum, and alluvium. In detail, the parent material of most watersheds is residuum or colluvium. Colluvium mainly derived from sandstone, siltstone, or granite exists widely in Lake Creek and Mica Creek; and residuum weathered from igneous, sandstone, shale, granite, or granodiorite exists widely in Big Creek, Mud Creek, Wheeler Ridge, Cottonwood Creek, and Little San Gorgonio Peak. The alluvium only exists in Waniha River in combination with colluvium. Hackmans Falls is glaciofluvial deposits and the remaining watersheds are residuum and colluvium.

Table 1

Location, drainage area, AI, precipitation, and temperature for study watersheds

Watershed index Watershed name County and state Drainage area (km2Aridity index Annual precipitation (mm) Annual temperature (°C) 
Lake Creek Lane County, OR 20.61 1.67 1,481.75 11.29 
Big Creek Monterey County, CA 14.40 0.52 724.23 14.90 
Elder Creek Mendocino County, CA 13.48 1.33 1,991.09 12.91 
Kenny Creek Mendocino County, CA 9.12 1.44 2,436.84 12.87 
Mud Creek Mendocino County, CA 13.27 1.38 1,652.21 13.62 
Waniha River Kauai County, HI 21.96 2.25 941.07 24.33 
Moscow Mountain Latah County, ID 16.81 0.72 1,080.72 7.44 
Hackmans Falls Sierra County, CA 8.28 1.31 1,880.05 7.92 
Lake Bonneville Tooele County, UT 7.65 0.61 730.24 6.96 
10 North Beaver Creek Denver County, CO 7.99 0.84 763.30 2.72 
11 Mica Creek Shoshone County, ID 11.73 1.16 1,495.29 6.02 
12 Wheeler Ridge Kern County, CA 5.80 2.10 1,765.29 8.70 
13 Cottonwood Creek Tuolume County, CA 6.75 1.03 1,109.57 11.78 
14 Little San Gorgonio Peak Riverside County, CA 6.45 0.74 799.67 9.21 
Watershed index Watershed name County and state Drainage area (km2Aridity index Annual precipitation (mm) Annual temperature (°C) 
Lake Creek Lane County, OR 20.61 1.67 1,481.75 11.29 
Big Creek Monterey County, CA 14.40 0.52 724.23 14.90 
Elder Creek Mendocino County, CA 13.48 1.33 1,991.09 12.91 
Kenny Creek Mendocino County, CA 9.12 1.44 2,436.84 12.87 
Mud Creek Mendocino County, CA 13.27 1.38 1,652.21 13.62 
Waniha River Kauai County, HI 21.96 2.25 941.07 24.33 
Moscow Mountain Latah County, ID 16.81 0.72 1,080.72 7.44 
Hackmans Falls Sierra County, CA 8.28 1.31 1,880.05 7.92 
Lake Bonneville Tooele County, UT 7.65 0.61 730.24 6.96 
10 North Beaver Creek Denver County, CO 7.99 0.84 763.30 2.72 
11 Mica Creek Shoshone County, ID 11.73 1.16 1,495.29 6.02 
12 Wheeler Ridge Kern County, CA 5.80 2.10 1,765.29 8.70 
13 Cottonwood Creek Tuolume County, CA 6.75 1.03 1,109.57 11.78 
14 Little San Gorgonio Peak Riverside County, CA 6.45 0.74 799.67 9.21 
Figure 1

Map for the location of study watersheds.

Figure 1

Map for the location of study watersheds.

Airborne LiDAR has become an important technique to obtain topographic data with a submeter and a finer spatial resolution, and plays an important role in the topographic and hydrologic analysis. In many studies, LiDAR data have been used to extract channel networks (Lashermes et al. 2007; Orlandini & Moretti 2009a, 2019b; Passalacqua et al. 2010; Orlandini et al. 2011; Sofia et al. 2011; Pelletier 2013; Clubb et al. 2014), the junction angle of channels (Hooshyar et al. 2017), and topographic depressions (Le & Kumar 2014). LiDAR data used in this study are obtained from the website of the Open Topography data (http://www.opentopography.org/). The ground returns of LiDAR, in which vegetation and buildings are filtered out, are used in this study. The DEMs with a resolution of 1 m and the projection coordination system (NAD_1983_UTM) are generated from the ground returns using the interpolation (spline) function of ArcGIS.

METHODOLOGY

Curvature calculation

The curvature is one of the important topographic attributes to display the structure of terrain and surface (Thompson et al. 2001; Schmidt et al. 2003) and is a critical variable in the erosional and depositional processes that drain water and the sediment. Terrain surface curvature includes a profile curvature in the direction of steepest slope (also called the gradient direction) and contour curvature perpendicular to the gradient direction. Profile curvature represents the change of slope in the gradient direction, whereas contour curvature indicates the change of aspect angle and reflects the convergence/divergence of the surface. Positive contour curvature implies flow convergence and indicates potential valley. Therefore, contour curvature is focused on this paper and is computed as follows (Mitášová & Hofierka 1993): 
formula
(1)
where z is elevation; x and y are the coordinate values in DEMs; zx (zy) represents the first derivative of elevation (z) with respect to x (y); zxx and zyy are the second derivatives; zxy is the first derivative of zx with respect to y. Hooshyar et al. (2016) adopted a fourth-order polynomial function shown as Equation (2) to fit each 3 × 3 moving window of DEM in order to robustly calculate contour curvature for each pixel (Zevenbergen & Thorne 1987): 
formula
(2)
Based on the coefficients of the polynomial function in Equation (2), Equation (1) can be transformed into the following equation for computing contour curvature: 
formula
(3)

Channel head identification

As previously mentioned, the channel head is identified based on the shape of contours (Hooshyar et al. 2016), and a spatially constant curvature threshold is not required in this method. Figure 2 illustrates the flowchart of channel head identification, mainly including six steps: (1) filtering DEM using the Perona–Malik Filter (Passalacqua et al. 2010); (2) calculating the curvature based on Equation (3); (3) removing sinks and calculating the curvature-adjusted flow direction given the flow direction from the traditional D (Tarboton 1997); (4) extracting valley and ridge skeletons; (5) extracting valley network by thinning and connecting the skeleton; and (6) identifying channel heads using k-means clustering and extracting the channel network. k-means clustering (Mac Queen 1966) is a widely used unsupervised classification algorithm to partition samples into k-clusters by minimizing the overall distance between samples and the centroid of their relevant clusters: 
formula
(4)
where k is the total number of clusters; Ni is the number of samples in cluster i; and pji is the jth sample of cluster i and is the corresponding centroid. In this paper, the contour lines within each tributary are divided into two clusters (channelized and unchannelized) based on their shapes. Each contour is processed as a 2 × M matrix consisting of M equally distanced points. The dissimilarity of two contour lines (e.g. L1 and L2) is defined as the minimum root mean square difference of one's best orientation. The best orientation of L1 to L2, denoted as ‘L1’, is obtained by translating, rotating, and scaling L1. Therefore, the dissimilarity between L1 and L2 is qualified by minimizing and then k-means clustering algorithm minimizes the dissimilarity index of contour lines in each cluster. Finally, the clustering result is a transition contour between channelized and unchannelized contour lines.
Figure 2

Flowchart of the procedure of channel head identification.

Figure 2

Flowchart of the procedure of channel head identification.

Figure 3 shows how channel heads are located. Figure 3(b) shows a zoomed-in area within a single valley tributary where the contour lines are classified into unchannelized and channelized sections by the means of k-means clustering based on the similarity of their shapes, and channel head is determined to be a member of the channelized cluster with the maximum elevation. The details on the procedure of valley and channel extraction are referred to Hooshyar et al. (2016). The curvature value for each pixel within the extracted channel network including channel head is obtained from this nonparametric method.

Figure 3

Map of the k-means clustering method showing how channel heads are located: (a) extracted valley and channel network of a representative section in Big Creek and (b) a zoomed-in area within a single valley tributary.

Figure 3

Map of the k-means clustering method showing how channel heads are located: (a) extracted valley and channel network of a representative section in Big Creek and (b) a zoomed-in area within a single valley tributary.

Distribution evaluation

The appropriate distribution function for describing the spatial variation of the channel head curvature is identified for the study watersheds. Following the concept of probability weight moment (Greenwood et al. 1979), Hosking (1990) proposed a regional frequency analysis method called L-moments, which is defined by a certain linear combination of sequence values in ascending order. Compared with the conventional moment method, the L-moment method is less affected by errors of individual points in samples (Vogel & Fennessey 1993; Vogel et al. 1994; Sankarasubramanian & Srinivasan 1999). The L-moment method has been widely accepted for evaluating the applicability of alternative distributions to data due to its reliability and accuracy. The details of L-moments are referred to Hosking & Wallis (1997).

L-moment ratio diagrams provide convenient and direct views on the comparison between the extracted channel head curvature within watersheds and theoretical statistical distributions. L-Cv vs. L-Skew is a commonly used type in L-moment diagrams, which enable us to compare the goodness of fit of a range of two parameters (one-parameter distribution is a special case of two-parameter distribution).

Among various theoretical distribution functions, the Gamma (G2) distribution is the most commonly advocated distribution for hydrological frequency analysis (Thom 1951; Buishand 1978; Geng et al. 1986; Duan et al. 1995; Lei et al. 2018). Gamma is bounded on the left at zero, is positively skewed and, most importantly, offers the possibility of describing a wide range of distribution shapes (Husak et al. 2010). For example, Thom (1951) suggested the Gamma (G2) distribution function for wet-day rainfall series; Geng et al. (1986) provided a review of the literature supporting the use of the G2 distribution for modeling wet-day rainfall; Karl & Knight (1998) utilized the G2 distribution to interpolate missing rainfall observations. In addition to G2 distribution, lognormal (LN2), Weibull (W2), and generalized Pareto (GP2) are also positively skewed and widely applied for statistical analysis (Saf 2009; Burgueño et al. 2010; Mosaffaie 2015). Therefore, in this paper, G2, LN2, W2, and GP2 distribution functions are selected as candidates for the spatial distribution of channel head curvature.

The probability plot correlation coefficient (PPCC) test is a simple but powerful goodness-of-fit test developed by Filliben (1975). Let denote the observed values and x(i) represents the ith largest value in a sample. The test uses the correlation r between the ordered observations xi and the corresponding fitted quantiles determined by plotting positions qi for each xi in which G(x) is the distribution of x. The PPCC statistic has a maximum value of 1 and values for r near 1.0 suggest that the observations could have been drawn from the fitted distribution. Essentially, r measures the linearity of the probability plot, providing a quantitative assessment of fit. If denotes the average value of the observations and denotes the average value of the fitted quantiles, then 
formula
(5)

RESULTS

In this part, Section ‘Extracted Channel Heads’ shows the extracted channel heads and their curvature values based on two typical watersheds, and Sections ‘Statistics of Channel Head Curvature, Distribution Evaluation and Analysis, and Fitted Gamma Distributions’ focus on the statistical analysis and distribution of extracted channel head curvatures within all the 14 watersheds.

Extracted channel heads

Channel networks are identified based on the nonparametric method developed by Hooshyar et al. (2016), as briefly described in the section ‘Channel Head Identification’. Two typical watersheds, Big Creek and Waniha River with the lowest and highest climate AI, respectively, are selected to show the extracted channel heads and relevant curvature values, and part of Big Creek and Waniha River watersheds are shown clearly in Figure 4.

Figure 4

Channel head curvature (m−1) within the part of watersheds: (a) Big Creek and (b) Waniha River.

Figure 4

Channel head curvature (m−1) within the part of watersheds: (a) Big Creek and (b) Waniha River.

As shown in Figure 4, the channel head curvature varies from 0.023 to 0.115 m−1 in Figure 4(a) and from 0.026 to 0.258 m−1 in Figure 4(b). The curvature calculation is essential for channel network extraction, and Figures 5(a) and 5(c) show the calculated curvature in Big Creek and Waniha River. It is obvious that those curvature pixels with spatially continuous high values are likely to be channelized. The extracted channel networks for Big Creek and Waniha River are also shown in Figures 5(b) and 5(d), respectively.

Figure 5

The calculated curvature and extracted channel network within Big Creek and Waniha River: (a–b) Big Creek and (c–d) Waniha River.

Figure 5

The calculated curvature and extracted channel network within Big Creek and Waniha River: (a–b) Big Creek and (c–d) Waniha River.

Statistics of the channel head curvature

Table 2 shows the drainage density, the total number of extracted channel heads, and the mean, the coefficient of variance (CV), and the skewness of channel head curvature values in each watershed. The total number of channel heads extracted within these watersheds varies from a minimum of 89 to a maximum of 1,722 with a mean value close to 710. Moreover, the drainage density defined as the ratio between the total length of extracted channels and the basin area varies in the range of 2.00–14.50 km/km2 with a mean value of 7.48 km/km2, which indicates a distinct difference in the channel network development among these selected small mountainous watersheds. The mean value of the channel head curvature varies across the study watersheds, ranging from 0.031 to 0.089 m−1. The CVs for the watersheds are all over 0.55, and the largest value is 0.76, which indicates the spatial variation of the channel head curvature within each watershed. For each watershed, the channel head curvature shows strong skewness, which is positively biased ranging from 0.83 to 1.80. Based on the results and analysis, it is evident that the channel head curvature varies considerably within and between these study watersheds. Therefore, it is uncertain how to define a spatially constant curvature threshold for channel head identification.

Table 2

Drainage density, total number of extracted channel heads, and statistical summary of the channel head curvature

Watershed index Drainage density (km/km2Number of channel heads Mean (m−1CV Skewness 
6.66 1,050 0.049 0.64 1.36 
8.48 1,159 0.069 0.59 1.16 
6.72 776 0.048 0.65 1.29 
9.87 685 0.062 0.61 0.95 
11.59 1,499 0.053 0.66 1.28 
7.44 1,722 0.086 0.59 1.10 
2.00 142 0.034 0.60 0.83 
5.67 410 0.037 0.75 1.73 
5.22 346 0.040 0.68 1.72 
10 5.27 300 0.031 0.61 1.23 
11 2.53 89 0.035 0.58 1.60 
12 9.89 489 0.037 0.72 1.46 
13 8.94 379 0.033 0.76 1.80 
14 14.50 867 0.089 0.73 1.59 
Watershed index Drainage density (km/km2Number of channel heads Mean (m−1CV Skewness 
6.66 1,050 0.049 0.64 1.36 
8.48 1,159 0.069 0.59 1.16 
6.72 776 0.048 0.65 1.29 
9.87 685 0.062 0.61 0.95 
11.59 1,499 0.053 0.66 1.28 
7.44 1,722 0.086 0.59 1.10 
2.00 142 0.034 0.60 0.83 
5.67 410 0.037 0.75 1.73 
5.22 346 0.040 0.68 1.72 
10 5.27 300 0.031 0.61 1.23 
11 2.53 89 0.035 0.58 1.60 
12 9.89 489 0.037 0.72 1.46 
13 8.94 379 0.033 0.76 1.80 
14 14.50 867 0.089 0.73 1.59 

Distribution evaluation and analysis

Figure 6 displays empirical and theoretical relationships between L-Cv and L-Skew for the channel head curvature distribution in individual watersheds. The curves in Figure 6 represent the theoretical distributions indicated, and each yellow point represents the empirical relationship between L-Cv and L-Skew in a single watershed. By comparing the data points with the theoretical distribution curves, we can see the degree to which the statistical characteristics of the extracted channel head curvature match those of the theoretical distributions. Figure 6 shows that these data points mainly locate between W2 and LN2 curves and fall primarily on the G2 curve, which suggests that G2 is obviously better than W2 and LN2 since more points are on or near the G2 curve.

Figure 6

L-Cv vs. L-Skew L-moment ratio diagram. L-moments of each watershed with dotted points.

Figure 6

L-Cv vs. L-Skew L-moment ratio diagram. L-moments of each watershed with dotted points.

G2 is identified as the potential distribution for the channel head curvature using the L-moment diagrams. To evaluate and analyze the distribution of the extracted channel head curvature, PPCC is further selected as a quantitative method for evaluating the goodness of fit of theoretical distributions. PPCC has been shown to be a powerful statistical method for evaluating alternative distributional hypotheses (Stedinger 1993). In this study, a relatively higher PPCC value (approximately 1.0) indicates better performance of G2 distribution for capturing the spatial variation of the channel head curvature in each watershed. It should be emphasized that PPCC values for the G2 distribution of the channel head curvature in the 14 study sites are all above 0.990 with a median value of 0.997. It can be concluded that G2 performs quite well for channel head curvature distribution within these study sites and, therefore, G2 is selected for fitting the channel head curvature.

Fitted Gamma distributions

The probability density function of Gamma distribution is given by: 
formula
(6)
where θ is the scale parameter, k is the shape parameter, and is the Gamma function evaluated at k. Figure 7 shows the histograms and fitted curves of Gamma distribution of the channel head curvature within each selected watershed. The histograms display positive (i.e., left) skewness for all the watersheds and the majority of channel head curvature values are less than 0.3 m−1 in most watersheds, and the peak value occurs in the range from 0.01 to 0.1 m−1. Obviously, the Gamma distribution fits well with the extracted channel head curvature within each watershed. The estimated parameters θ and k vary in the range of 0.011–0.067 and 1.956–3.056, respectively, which indicates the spatial variation of the channel head curvature across selected study watersheds. Generally, the channel head curvature values have better agreement with the fitted Gamma curve in the climbing stage, while the values during the falling stage are overestimated in almost all the selected watersheds. For those catchments with a larger drainage area, such as Lake Creek, Big Creek, Elder Creek, Waniha River, and Moscow Mountain, the channel head curvature values are significantly closer to the fitted curves probably due to the fact that large watersheds tend to have more channel heads. It should be noted that the fitted curves in Elder Creek, Kenny Creek, Mud Creek, Moscow Mountain, and Mica Creek are sharper and thinner compared to the rest and the mean annual precipitation in most of these study sites are above 1,500 as shown in Table 1.
Figure 7

Histograms and fitted curves of Gamma distribution of the channel head curvature within each selected watershed using a 0.01 m−1 bin size.

Figure 7

Histograms and fitted curves of Gamma distribution of the channel head curvature within each selected watershed using a 0.01 m−1 bin size.

DISCUSSIONS

Comparison with field-mapped channel heads

Since channel head identification is sensitive to the method used for channel extraction and may bring uncertainties in analyzing the distribution of the channel head curvature, here we evaluate the distribution of field-mapped channel head curvature. Clubb et al. (2014) have compared three state-of-the-art methods of channel head prediction with field-mapped channel heads in three catchments: Feather River in California, Mid Bailey Run, and Indian Creek in Ohio, and the location of mapped channel heads are freely open and available online (https://datashare.is.ed.ac.uk/handle/10283/524). For mapped channel head curvature extraction, bare earth DEMs with a resolution of 2.5 feet for Indian Creek and Mid Bailey Run were obtained from the website (http://ogrip.oit.ohio.gov/ProjectsInitiatives/OSIPDataDownloads.aspx) and were resampled to the UTM Zone 17N with a resolution of 1 m by the means of cubic interpolation. Since relatively few mapped channel heads (approximately 15) were collected in Feather River, we focused on the other two catchments, namely Mid Bailey Run and Indian Creek. Figure 8 illustrates the locations of the extracted and mapped channel heads in Indian Creek and Mid Bailey Run. It should be noted that not all the channel heads within the watersheds are located through field investigation; therefore, comparison here focuses on those tributaries where mapped channel heads are available. It can be easily seen that there are only slight distance errors between the extracted and mapped channel heads, which clearly indicate the reliability of the method used for channel head identification in this study.

Figure 8

The extracted and mapped channel heads: (a) Indian Creek and (b) Mid Bailey Run.

Figure 8

The extracted and mapped channel heads: (a) Indian Creek and (b) Mid Bailey Run.

The curvature values of the mapped channel heads are obtained from the corresponding pixels according to mapped channel heads' locations. Table 3 shows the mean, CV, and skewness of the mapped channel head curvature within these two watersheds. The mapped channel head curvature also varies in Indian Creek and Mid Bailey Run with the CV of 0.567 and 0.626, respectively. The skewness is positively biased with 1.061 and 0.488 in Indian Creek and Mid Bailey Run, respectively. It is consistent with the previous 14 watersheds that the channel head curvature indicates a significant positive skewness distribution.

Table 3

Statistical summary of the mapped channel head curvature in Indian Creek and Mid Bailey Run

Watershed Mapped channel heads Maximum (m−1Minimum (m−1Mean (m−1CV Skewness 
Indian Creek 36 0.240 0.013 0.082 0.567 1.061 
Mid Bailey Run 53 0.317 0.005 0.131 0.626 0.488 
Watershed Mapped channel heads Maximum (m−1Minimum (m−1Mean (m−1CV Skewness 
Indian Creek 36 0.240 0.013 0.082 0.567 1.061 
Mid Bailey Run 53 0.317 0.005 0.131 0.626 0.488 

We also examine the Gamma distribution of the mapped channel head curvature using one sample Kolmogorov–Smirnov test (K–S test) in MATLAB based on the null hypothesis that the Gamma distribution is fitted, and the results indicate that the K–S test accepts the null hypothesis at the 5% significance level with the p-value of 0.82 and 0.58 in Indian Creek and Mid Bailey Run, respectively.

Sensitivity analysis of the distribution

Hooshyar et al. (2016) utilized the nonparametric method in three study sites and concluded that the proposed method had equal or better performance compared with existing methods in some indices, such as distance error, reliability, and sensitivity. However, as mentioned before, methods proposed for channel head extraction may move the heads upstream or downstream, and thus distance error between extracted and mapped or realistic channel heads exists. Therefore, in this case, we also discuss this source of error by moving the extracted channel heads downstream by 5 and 10 m and then recalculate all the distributions. Figure 9 shows the L-moment ratio and the empirical distribution that the extracted curvature values approximately fitted. It can be easily seen that most points are on or near the G2 curve, while the rest (two or three out of 14) are a little far away, which indicates the applicability and robustness of the G2 distribution in capturing the spatial variation of the channel head curvature.

Figure 9

L-Cv vs. L-Skew L-moment ratio diagram. L-moments of each watershed with dotted points: (a) downstream 5 m and (b) downstream 10 m.

Figure 9

L-Cv vs. L-Skew L-moment ratio diagram. L-moments of each watershed with dotted points: (a) downstream 5 m and (b) downstream 10 m.

Constant curvature threshold analysis

In general, varying curvature thresholds may provide more reliable predictions of channel head across catchments having different topography; however, it is difficult to determine the varied thresholds for a watershed. And, the constant curvature threshold-based methods have been often used for identifying the location of channel heads because of its simplicity and efficiency in practice. For example, Pelletier (2013) adopted a curvature threshold of 0.1 m−1 to locate channel heads. Lashermes et al. (2007) developed a method to determine the curvature threshold, denoted by approximately 0.025 m2/m, for valley identification within individual watersheds by comparing the distribution of curvature of all grids with a Gaussian distribution through the quantile–quantile plot. We computed the curvature thresholds for each watershed using the method developed by Lashermes et al. (2007). The identified curvature threshold within these 14 study watersheds is shown in Table 4. Meanwhile, the mode of the extracted channel head curvature for each watershed derived from the method proposed by Hooshyar et al. (2016) is also listed in Table 4 for comparison. The mode of curvature means the value that appears most frequently in a series of channel head curvatures within a catchment. In this study, those channel head curvature values located within the highest bin of the histogram (Figure 7) were averaged as threshold. It is evident that the above-identified curvature threshold is relatively smaller or approximate to the mode of the extracted channel head curvature. Moreover, the identified curvature thresholds vary in the range of 0.008–0.032 m−1, which indicate again the spatial variation of the channel head curvature within different watersheds. Among the watersheds selected in this study, Big Creek and Waniha River are characterized by the lowest and highest AI, respectively, representing the typical arid and humid watersheds. The mode of the extracted channel head curvature in Waniha River is significantly larger than that in Big Creek, which is consistent with the calculated mean value as shown in Table 2.

Table 4

Mode of the extracted channel head curvature and curvature thresholds identified by the method proposed by Lashermes et al. (2007) within 14 study watersheds

Watershed index Watershed name Mode of curvature (m−1Identified curvature threshold (m−1
Lake Creek 0.025 0.014 
Big Creek 0.035 0.030 
Elder Creek 0.035 0.015 
Kenny Creek 0.054 0.014 
Mud Creek 0.045 0.015 
Waniha River 0.065 0.032 
Moscow Mountain 0.025 0.006 
Hackmans Falls 0.025 0.012 
Lake Bonneville 0.025 0.013 
10 North Beaver Creek 0.025 0.012 
11 Mica Creek 0.026 0.005 
12 Wheeler Ridge 0.035 0.008 
13 Cottonwood Creek 0.015 0.010 
14 Little San Gorgonio Peak 0.036 0.025 
Watershed index Watershed name Mode of curvature (m−1Identified curvature threshold (m−1
Lake Creek 0.025 0.014 
Big Creek 0.035 0.030 
Elder Creek 0.035 0.015 
Kenny Creek 0.054 0.014 
Mud Creek 0.045 0.015 
Waniha River 0.065 0.032 
Moscow Mountain 0.025 0.006 
Hackmans Falls 0.025 0.012 
Lake Bonneville 0.025 0.013 
10 North Beaver Creek 0.025 0.012 
11 Mica Creek 0.026 0.005 
12 Wheeler Ridge 0.035 0.008 
13 Cottonwood Creek 0.015 0.010 
14 Little San Gorgonio Peak 0.036 0.025 

We further discuss the performance of using the mode of the extracted channel head curvature as constant threshold to locate channel heads within these 14 watersheds. Table 5 shows the number of channel heads with the curvature value in the range of mode ± 0.01 m−1 (a bin size), in which the channel heads can be regarded as correctly identified by a constant threshold. The related percentage over the total number of channel heads is also computed and shown in Table 5, which shows that the percentage varies in the range of 17.6–50.0% and indicates that the location errors significantly exist when using the threshold-based method for channel head identification.

Table 5

Total number and responding percentage of channel heads with the curvature value in the range of mode ±0.01 m−1

Watershed index Number of channel heads with curvature within mode ±0.01 Percentage 
321 30.6 
244 21.1 
249 32.1 
167 24.4 
419 28.0 
343 19.9 
58 40.8 
161 39.3 
150 43.4 
10 150 50.0 
11 39 43.8 
12 162 33.1 
13 171 45.1 
14 153 17.6 
Watershed index Number of channel heads with curvature within mode ±0.01 Percentage 
321 30.6 
244 21.1 
249 32.1 
167 24.4 
419 28.0 
343 19.9 
58 40.8 
161 39.3 
150 43.4 
10 150 50.0 
11 39 43.8 
12 162 33.1 
13 171 45.1 
14 153 17.6 

Further analysis on the distribution of channel head curvature

The nonparametric method proposed by Hooshyar et al. (2016) has its own definite advantage in the way that the channel head prediction is performed separately within the individual valley tributary with no spatially constant curvature threshold, which makes it able to capture the spatial variation of the channel head curvature threshold through distribution analysis. Based on previous studies, we utilized this nonparametric method to extract the channel head curvature within 14 watersheds and analyzed a Gamma distribution of the channel head curvature, which can offer an overall view of the channel head curvature within and between watersheds. Indeed, there are also some limitations in this study. The curvature is the most commonly used topographic feature for channel head identification. The paper focuses on the distribution of the channel head curvature within and between watersheds, which occupies and underpins the whole study. However, other topographic features such as area, slope, and width have also been used to locate channel heads, and they were spatially constant thresholds and examined spatially varied within each watershed. Therefore, more effort should be exerted on these topographic features from the view of comprehensive analysis, not just be limited to a single attribute, in a future study. In addition, the distribution analysis and comparison of channel head curvatures in this paper are based on a specified grid resolution; however, the grid resolution has a great influence on the calculation of curvature from DEMs, and further research on these topics could improve understanding of spatial variation of the channel head curvature across a wider range of grids.

CONCLUSIONS

Due to the complexity associated with channel head location, the purpose of capturing the spatial variation of the channel head curvature in this study is very practical. Based on a nonparametric method, the channel heads are extracted from high-resolution digital elevation data at submeter obtained from LiDAR in 14 mountainous watersheds in the United States. Previous research on channel head identification is mainly concentrated on the spatially constant threshold, and this study built on previous research by analyzing the spatial variation of curvature within individual watersheds.

For these 14 watersheds selected in this study, the extracted channel head curvature has been shown to vary within and between watersheds. Through L-moment analysis, we concluded that the channel head curvature can be well represented by the Gamma distribution. Meanwhile, we also tested the curvature extraction and distribution analysis in Indian Creek and Mid Bailey Run where field-mapped channel heads are available; the comparisons between the extracted and mapped channel heads clearly indicate that our statistical analysis based on the extracted channel head curvature is reliable.

We discussed the performance of using curvature threshold to locate channel head within and between watersheds, and it can be concluded that different curvature thresholds are needed for different watersheds, particularly, a constant curvature threshold is not suitable for channel head identification ignoring the spatial variation of the channel head curvature between watersheds. Constant threshold-based methods for channel head identification undoubtedly gain significant distance errors even within a single watershed.

As with any new attempt, definitive conclusions regarding channel head curvature distribution and its spatial variation analyzed in this study must be based on the substantive test using large samples of watersheds representing a wide range of geographical and hydrogeological conditions. However, the paper is just an attempt for the channel head curvature distribution and constrained by the availability of local datasets. Work covering more study sites and future research aiming at relating the channel head curvature distribution to hydrological and geomorphic information, such as rainfall, evaporation, and land cover, can contribute to improving the reliability and accuracy of channel head identification.

ACKNOWLEDGEMENTS

This work is partially supported by the National Natural Science Foundation of China (No. 91647201, 91547116, and 51709033).

REFERENCES

REFERENCES
Bierman
P. R.
&
Montgomery
D. R.
2014
Key Concepts in Geomorphology
.
W.H. Freeman
,
New York
.
Buishand
T. A.
1978
Some remarks on the use of daily rainfall models
.
Journal of Hydrology
36
(
3
),
295
308
.
Burgueño
A.
,
Martínez
M. D.
,
Lana
X.
&
Serra
C.
2010
Statistical distributions of the daily rainfall regime in Catalonia (Northeastern Spain) for the years 1950–2000
.
International Journal of Climatology
25
(
10
),
1381
1403
.
Clubb
F. J.
,
Mudd
S. M.
,
Milodowski
D. T.
,
Hurst
M. D.
&
Slater
L. J.
2014
Objective extraction of channel heads from high-resolution topographic data
.
Water Resources Research
50
(
5
),
4283
4304
.
Dietrich
W. E.
,
Dunne
T.
1993
The channel head
. In:
Channel Network Hydrology
(
Beven
K.
&
Kirkby
M. J.
, eds).
John Wiley and Sons
,
New York
, pp.
175
219
.
Dietrich
W. E.
,
Wilson
C. J.
,
Montgomery
D. R.
,
McKean
J.
&
Bauer
R.
1992
Erosion thresholds and land surface morphology
.
Geology
20
(
8
),
675
679
.
Dietrich
W. E.
,
Wilson
C. J.
,
Montgomery
D. R.
&
McKean
J.
1993
Analysis of erosion thresholds, channel networks, and landscape morphology using a digital terrain model
.
Journal of Geology
101
(
2
),
259
278
.
Duan
J.
,
Sikka
A. K.
&
Grant
G. E.
1995
A comparison of stochastic models for generating daily precipitation at the H J Andrews experimental forest
.
Northwest Science
69
(
4
),
318
329
.
Geng
S.
,
Vries
F. W. T. P. D.
&
Supit
I.
1986
A simple method for generating daily rainfall data
.
Agricultural & Forest Meteorology
36
(
4
),
363
376
.
Greenwood
J. A.
,
Landwehr
J. M.
,
Matalas
N. C.
&
Wallis
J. R.
1979
Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form
.
Water Resources Research
15
(
5
),
1049
1054
.
Grieve
S. W. D.
,
Mudd
S. M.
,
Milodowski
D. T.
,
Clubb
F. J.
&
Furbish
D. J.
2016
How does grid-resolution modulate the topographic expression of geomorphic processes?
Earth Surface Dynamics
4
(
3
),
1
48
.
Hooshyar
M.
,
Kim
S.
,
Wang
D.
&
Medeiros
S. C.
2015
Wet channel network extraction by integrating LiDAR intensity and elevation data
.
Water Resources Research
51
(
12
),
10029
10046
.
Hooshyar
M.
,
Wang
D.
,
Kim
S.
,
Medeiros
S. C.
&
Hagen
S. C.
2016
Valley and channel networks extraction based on local topographic curvature and K-means clustering of contours
.
Water Resources Research
52
(
10
),
8081
8102
.
Hooshyar
M.
,
Singh
A.
&
Wang
D.
2017
Hydrologic controls on junction angle of river networks
.
Water Resources Research
53
(
5
),
4073
4083
.
Hosking
J. R. M.
1990
L-moments: analysis and estimation of distributions using linear combinations of order statistics
.
Journal of the Royal Statistical Society Series B (Methodological)
52
(
1
),
105
124
.
Hosking
J. R. M.
&
Wallis
J. R.
1997
Regional Frequency Analysis
.
Cambridge University Press
,
Cambridge
.
Husak
G. J.
,
Michaelsen
J.
&
Funk
C.
2010
Use of the gamma distribution to represent monthly rainfall in Africa for drought monitoring applications
.
International Journal of Climatology
27
(
7
),
935
944
.
Karl
T. R.
&
Knight
R. W.
1998
Secular trends of precipitation amount, frequency, and intensity in the United States
.
Bulletin of the American Meteorological Society
79
(
2
),
231
241
.
Lashermes
B.
,
Foufoula-Georgiou
E.
&
Dietrich
W. E.
2007
Channel network extraction from high resolution topography using wavelets
.
Geophysical Research Letters
34
(
23
),
102
120
.
Le
P. V. V.
&
Kumar
P.
2014
Power law scaling of topographic depressions and their hydrologic connectivity
.
Geophysical Research Letters
41
(
5
),
1553
1559
.
Lei
Y.
,
Hanson
L. S.
,
Ding
P.
,
Wang
D.
&
Vogel
R. M.
2018
The probability distribution of daily precipitation at the point and catchment scales in the United States
.
Hydrology & Earth System Sciences
22
(
2
),
6519
6531
.
MacQueen
J.
1966
Some methods for the classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume I. University of California Press, Berkeley, CA, pp. 281–291
.
Marie
M.
1957
Accuracy of determination of stream lengths from topographic maps
.
Eos Transactions American Geophysical Union
38
(
1
),
86
.
Montgomery
D. R.
&
Dietrich
W. E.
1988
Where do channels begin?
Nature
336
,
232
.
O'Callaghan
J. F.
&
Mark
D. M.
1984
The extraction of drainage networks from digital elevation data
.
Computer Vision, Graphics, and Image Processing
28
(
3
),
323
344
.
Orlandini
S.
&
Moretti
G.
2009a
Comment on “Global search algorithm for nondispersive flow path extraction” by Kyungrock Paik
.
Journal of Geophysical Research Earth Surface
114
(
F4
),
F04004
.
Orlandini
S.
&
Moretti
G.
2009b
Determination of surface flow paths from gridded elevation data
.
Water Resources Research
45
(
3
),
150
164
.
Orlandini
S.
,
Tarolli
P.
,
Moretti
G.
&
Fontana
G. D.
2011
On the prediction of channel heads in a complex alpine terrain using gridded elevation data
.
Water Resources Research
47
(
2
),
247
255
.
Passalacqua
P.
,
Trung
T. D.
,
Foufoula-Georgiou
E.
,
Sapiro
G.
&
Dietrich
W. E.
2010
A geometric framework for channel network extraction from LiDAR: nonlinear diffusion and geodesic paths
.
Journal of Geophysical Research Earth Surface
115
(
F1
),
F01002
.
Peckham
S. D.
1995
Self-similarity in the three-dimensional geometry and dynamics of large river basins
.
University of Colorado
.
Sankarasubramanian
A.
&
Srinivasan
K.
1999
Investigation and comparison of sampling properties of L-moments and conventional moments
.
Journal of Hydrology
218
(
1
),
13
34
.
Schmidt
J.
,
Evans
I. S.
&
Brinkmann
J.
2003
Comparison of polynomial models for land surface curvature calculation
.
International Journal of Geographical Information Science
17
(
8
),
797
814
.
Stedinger
J. R.
1993
Frequency analysis of extreme events
.
Handbook of Hydrology
.
McGraw-Hill
,
New York
, p.
18
.
Strahler
A. N.
1952
Hypsometric (area-altitude) analysis of erosional topography
.
Geological Society of America Bulletin
23
(
6
),
1117
1142
.
Strahler
A. N.
1953
Revision of Horons’ quantitative factors in erosional terrain
.
Eos Transactions American Geophysical Union
34
,
345
365
.
Strahler
A. N.
1957
Quantitative analysis of watershed geomorphology
.
Eos Transactions American Geophysical Union
38
(
6
),
913
920
.
Tarboton
D. G.
,
Bras
R. L.
&
Rodriguez-Iturbe
I.
1991
On the extraction of channel networks from digital elevation data
.
Hydrological Processes
5
(
1
),
81
100
.
Thom
H. C. S.
1951
A frequency distribution for precipitation
.
Bulletin of the American Meteorological Society
32
(
10
),
387
.
Vogel
R. M.
&
Fennessey
N. M.
1993
L moment diagrams should replace product moment diagrams
.
Water Resources Research
29
(
6
),
1745
1752
.
Vogel
R. M.
,
Member
J.
&
Fennessey
N. M.
1994
Flow-duration curves. I: new interpretation and confidence intervals
.
Journal of Water Resources Planning and Management
120
(
4
),
485
504
.
Willgoose
G.
,
Bras
R. L.
&
Rodriguez-Iturbe
I.
1991
A coupled channel network growth and hillslope evolution model: 1. Theory
.
Water Resources Research
27
(
7
),
1671
1684
.
Zevenbergen
L. W.
&
Thorne
C.
1987
Quantitative analysis of land surface topography
.
Earth Surface Processes & Landforms
12
(
1
),
47
56
.