Spatial variation of channel head curvature in small mountainous watersheds

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. The Study area and data source section introduces the study sites and the topographic data used for extracting channel heads. The Methodology section 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, Discussion and Conclusions are presented in the following sections.

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 km². 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.

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.

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.

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.

In this part, the Extracted channel heads section shows the extracted channel heads and their curvature values based on two typical watersheds, and the Statistics of the channel head curvature, Distribution evaluation and analysis, and Fitted Gamma distributions sections focus on the statistical analysis and distribution of extracted channel head curvatures within all the 14 watersheds.

Channel networks are identified based on the nonparametric method developed by Hooshyar et al. (2016), as briefly described in the Channel head identification section. 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.

As shown in Figure 4, the channel head curvature varies from 0.023 to 0.115 m⁻¹ in Figure 4(a) and from 0.026 to 0.258 m⁻¹ in Figure 4(b). The curvature calculation is essential for channel network extraction, and Figure 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 Figure 5(b) and 5(d), respectively.

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/km² with a mean value of 7.48 km/km², 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⁻¹. 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.

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.

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.

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.

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.

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.

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.

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⁻¹ to locate channel heads. Lashermes et al. (2007) developed a method to determine the curvature threshold, denoted by approximately 0.025 m²/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⁻¹, 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.

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⁻¹ (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.

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.

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.

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