A segmentation approach for the reproducible extraction and quantification of knickpoints from river long profiles

Changes in the steepness of river profiles or abrupt vertical steps (i.e. waterfalls) are thought to be indicative of changes in erosion rates, lithology or other factors that affect landscape evolution. These changes are referred to as knickpoints or knickzones and are pervasive in bedrock river systems. Such features are thought to reveal information about landscape evolution and patterns of erosion, and therefore their locations are often reported in the geomorphic literature. It is imperative that studies reporting knickpoints and knickzones use a reproducible method of quantifying their locations, as their number and spatial distribution play an important role in interpreting tectonically active landscapes. In this contribution we introduce a reproducible knickpoint and knickzone extraction algorithm that uses river profiles transformed by integrating drainage area along channel length (the so-called integral or χ method). The profile is then statistically segmented and the differing slopes and step changes in the elevations of these segments are used to identify knickpoints, knickzones and their relative magnitudes. The output locations of identified knickpoints and knickzones compare favourably with human mapping: we test the method on Santa Cruz Island, CA, using previously reported knickzones and also test the method against a new dataset from the Quadrilátero Ferrífero in Brazil. The algorithm allows for the extraction of varying knickpoint morphologies, including stepped, positive slope-break (concave upward) and negative slope-break knickpoints. We identify parameters that most affect the resulting knickpoint and knickzone locations and provide guidance for both usage and outputs of the method to produce reproducible knickpoint datasets.

1 Introduction

Landscapes are shaped by competition between crustal processes such as tectonic plate motion or dynamic topography and deposition or erosion at the Earth's surface. This competition, if unperturbed, tends toward a topographic steady state at which vertical motions are counterbalanced by erosion (e.g. Hack, 1960; Willett and Brandon, 2002). In unglaciated landscapes, the main driver of erosion is the river system (Ahnert, 1970), which incises the landscape to remove and transport material from uplands to active basins. The analysis of river long profiles has been a key method to interpret landscape evolution (e.g. Wobus et al., 2006), from the early recognition of graded rivers (e.g. Gilbert, 1877) to the generalised recognition that river profiles reflect varying erosion processes (e.g. Mackin, 1948; Hack, 1960; Howard, 1965; Howard et al., 1994; Dietrich et al., 2003; Kirby and Whipple, 2012).

In a river system, topographic steady state requires spatially stable rock uplift and climatic conditions over a long period of time (Willett and Brandon, 2002). In most landscapes, however, these conditions are unlikely (Baldwin et al., 2003). Many processes have been suggested to result in both spatial and temporal variations in uplift rate, such as varying tectonic stress (e.g. Kirby and Whipple, 2012), complex mantle processes inducing vertical motions (e.g. Faccenna and Becker, 2010; Braun, 2010), uplift driven by differential rock density (Braun et al., 2014) and base-level variations linked to eustatic variations (e.g. Powell, 1875; Lambeck and Chappell, 2001; Schumann et al., 2016). River systems affected by these processes respond by transmitting signals upstream through the channel network (e.g. Whipple et al., 1999; Royden and Perron, 2013), eventually driving drainage network reorganisation and resulting in additional transient signals (e.g. Mather, 2000; Castelltort et al., 2012; Willett et al., 2014; Whipple et al., 2017b; Mudd, 2017). Moreover, river profiles are also affected by intrinsic landscape properties, such as fracture density (e.g Whipple, 2002) or differential lithology (e.g. Stock and Montgomery, 1999; Forte et al., 2016), which can also lead to morphological adjustment of the channel (e.g. Kirby and Whipple, 2012). The most direct and widely observed expression of river adjustment to transient or intrinsic perturbations is a discrete change in river gradient, commonly referred to as a “knickpoint”.

Changes in channel gradient linked to different lithologies have been recognised in geomorphological studies for centuries. Lapparent (1896) suggested that these changes may represent “successive reaches” with different base levels, hypothesising that these reaches somehow migrate upstream. Davis (1889) recognised the tectonic genesis of some of these signals, describing how landscapes experience erosion cycles with periods of “rejuvenation” followed by periods of gradual adjustment and thus transience. However, these early studies did not name such morphologies as distinct entities. The term knickpoint was first introduced into the geomorphological literature by Knopf (1924), borrowing the word from chemical sciences to “denote an abrupt change in direction from a gentle concave curve to a curve that is convex upward” (p. 636).

Based on earlier observations on the topography and geology of the Appalachians (e.g. Barrell, 1920; Bascom, 1921), Knopf (1924) described a knickpoint as a migrating steepened boundary between two river reaches. She went on to state that the downstream reach should flow with a gradient determined by the present-day balance between uplift and erosion, and the upstream reach should flow with a gradient representing an older such balance. The recognition of knickpoints and their significance in transient landscapes has driven much research into interpreting topography (e.g. Wobus et al., 2006; Crosby and Whipple, 2006; Abbühl et al., 2011; Kirby and Whipple, 2012), as well as using river profiles to extract past uplift histories (e.g. Pritchard et al., 2009).

The diverse nature of knickpoint formation means that these features have been used to investigate many geomorphological problems. For example, retreat rates have been used to link knickpoints with tectonic events and faulting (e.g. Attal et al., 2008, 2011; Whittaker and Boulton, 2012) or climatically triggered base-level fall (e.g. Crosby and Whipple, 2006; Baynes et al., 2015; Neely et al., 2017). Although migrating knickpoints are commonly associated with base-level variations, Haviv et al. (2010) highlighted the role of differential lithologies in the retreat rates of vertical knickpoints within tectonically and climatically stable landscapes. Furthermore, Scheingross and Lamb (2016) and Scheingross et al. (2017) noted the importance of sediment supply and hydraulic conditions in waterfall retreat, providing a quantitative interpretation of the early observations of Lapparent (1896) on waterfall migration. Cook et al. (2013) observed an important correlation between knickpoint retreat and bedload transport, further highlighting the importance of sediment transport. Bishop and Goldrick (2010) demonstrated that considering the role of resistant lithologies is crucial when studying landscape evolution, as they can considerably slow down landscape response time to transient signals. Other studies have linked knickpoints directly to landscape characteristics such as heterogeneous lithology (e.g. Tucker and Slingerland, 1996; Stock and Montgomery, 1999; Kirby et al., 2003; Duvall, 2004). Recent analogue experiments on knickpoint retreat (e.g. Baynes et al., 2018) have highlighted the interconnectivity of all these processes and the need to consider both internal and external landscape characteristics.

These examples demonstrate the importance but also the diversity of transient and lithologic signals in landscapes and highlight the fact that different processes can generate remarkably similar channel morphology. It is therefore crucial to define knickpoints morphologically before drawing interpretations about their significance in terms of processes or genesis. In this contribution, we aim to provide a method for reproducibly and systematically extracting knickpoints within real landscapes based on river profile morphology.

1.1 Knickpoint morphology and detection

1.1.1 Morphological description

Knickpoints can be defined as discrete changes in river gradient (Whipple et al., 1999). Haviv et al. (2010) proposed two endmember knickpoints: break-in-slope knickpoints expressed by an abrupt change in river gradient and break-in-elevation knickpoints characterised by a step in the elevation, such as a waterfall, with similar gradients on both sides of the knickpoint. These knickpoints are now commonly referred to as slope-break knickpoints and vertical-step knickpoints (e.g. Kirby and Whipple, 2012; Neely et al., 2017). Kirby and Whipple (2012) suggest that although vertical-step knickpoints tend to be linked to discrete heterogeneities along the river profile (e.g. caused by geological boundaries), both morphologies can either be fixed or mobile and each style of knickpoint may be generated by a range of processes.

As discussed in Goldrick and Bishop (2007) and Kirby and Whipple (2012), both morphologies can be detected using a slope–area plot (Fig. 1) or a slope–distance plot. It has long been observed that channel gradients vary systematically as a function of drainage area. For example, Gilbert (1877) stated the following: “In general we may say that, ceteris paribus, declivity bears an inverse relation to quantity of water” (p. 114). How do we then find anomalous channel gradients? In the mid-twentieth century, authors such as Hack (1957) and Morisawa (1962) found systematic, quantitative relationships between channel gradient and drainage area that are often used as a proxy for discharge. Morisawa (1962) and later Flint (1974) recognised that channel gradients often declined systematically downstream in a trend that could be described by a power law:

$$\begin{array}{}\text{(1)}& S=k_{\mathrm{s}}{A}^{-\mathit{\theta }},\end{array}$$

where θ is referred to as the concavity index since it describes how concave a profile is: the higher the value, the more rapidly a channel's gradient decreases downstream. The term k_s is called the steepness index, as it sets the overall gradient of the channel, and a number of authors have noted that k_s frequently scales with erosion rate in lithologically homogeneous landscapes (e.g. Ouimet et al., 2009; DiBiase et al., 2010; Scherler et al., 2014; Mandal et al., 2015; Harel et al., 2016). A knickpoint might manifest itself as an abrupt change in slope–area scaling and lead to local variations in k_s (Fig. 1a).

Figure 1 Different methods to detect knickpoints. (a) Cartoon showing how vertical-step and slope-break knickpoints appear in slope–area plots; adapted from Kirby and Whipple (2012). (b) A slope–area plot derived from SRTM 30 m resolution data in Romania; the catchment's outlet coordinates are 45.252842, 26.375697 (WGS84). Different colours represent different tributaries, small “+” symbols are individual data points and circles are logarithmically binned data. A single slope-break knickpoint can be interpreted but minor knickpoints are more difficult to extract. (c) The same basin represented in a χ-elevation plot using θ=0.15.

However, using slope–area data derived from digital elevation models (DEMs) suffers from noise in channel slopes, leading to scattering of gradient data, as discussed in Perron and Royden (2013). Wobus et al. (2006) proposed methods to reduce the effect of noise and extract trends from slope–area plots. These recommendations include regular sampling of elevations to extrapolate artefact-free contour lines or logarithmic binning by drainage area. Smoothing induces inexorable data loss and may result in difficulties detecting subtle but important features such as knickpoints (Fig. 1b).

Alternatively, we can integrate Eq. (1), since $S=\mathrm{d}z/\mathrm{d}x$ where z is elevation and x is distance along the channel (e.g. Whipple et al., 2017a), resulting in

$$\begin{array}{}\text{(2)}& z\left(x\right)=z\left(x_b\right)+\left({\displaystyle \frac{k_{\mathrm{s}}}{{A_{\mathrm{0}}}^{\mathit{\theta }}}}\right)\underset{x_b}{\overset{x}{\int }}{\left({\displaystyle \frac{A_{\mathrm{0}}}{A\left(x\right)}}\right)}^{\mathit{\theta }}\mathrm{d}x,\end{array}$$

where A₀ is a reference drainage area introduced to non-dimensionalise the area term within the integral in Eq. (2). We can then define a longitudinal coordinate, χ (Royden et al., 2000):

$$\begin{array}{}\text{(3)}& \mathit{\chi }=\underset{x_b}{\overset{x}{\int }}{\left({\displaystyle \frac{A_{\mathrm{0}}}{A\left(x\right)}}\right)}^{\mathit{\theta }}\mathrm{d}x.\end{array}$$

χ has dimensions of length and is defined such that at any point in the channel,

$$\begin{array}{}\text{(4)}& z\left(x\right)=z\left(x_b\right)+\left({\displaystyle \frac{k_{\mathrm{s}}}{{A_{\mathrm{0}}}^{\mathit{\theta }}}}\right)\mathit{\chi }.\end{array}$$

The χ approach to represent normalised long profiles (Eqs. 4 and 3) can serve as an alternative method to explore the slope–area relationship within a drainage network. The χ coordinate integrates information about drainage area, while requiring less smoothing and lumping than log(S)–log(A) plots (Fig. 1c). This approach has been widely used in recent studies (e.g. Perron and Royden, 2013; Mudd et al., 2014; Willett et al., 2014; Mouchené et al., 2017; Whipple et al., 2017b; Neely et al., 2017; Moodie et al., 2017).

1.1.2 Existing algorithms

Traditional knickpoint identification from DEMs relied upon user-based selection along river long profiles (e.g. Hayakawa and Oguchi, 2006; Wobus et al., 2006). Several computational methods have been proposed for extracting knickpoints from DEM-derived datasets. The first (semi-)automated methods taking advantage of digital topographic data used long-profile geometry to isolate knickpoints or knickzones. Hayakawa and Oguchi (2006) proposed a semi-automated extraction method based on decreasing gradient with increasing length. This method involved the use of ArcGIS and spreadsheet software to process the outputs for each river. Recognising the need for automated regional knickpoint mapping methods in geomorphological studies, Gonga-Saholiariliva et al. (2011) proposed an automated algorithm to map abrupt changes in river gradient using slope, profile and plan-view curvature. Gallen et al. (2013) used systematic changes in profile convexity over given thresholds (>20 m in elevation drop coupled with a slope threshold ≥0.1) to isolate knickpoints in fluvially dominated channels with the aim of reconstructing rejuvenation events, both climatically and tectonically driven, in the southern Appalachians. A similar method has been implemented in ArcGIS by Queiroz et al. (2015). More recently, Zahra et al. (2017) published an ArcGIS toolset (called KET) that automates and optimises the Hayakawa and Oguchi (2006) method. These methods are based on the direct use of channel elevation, gradient and curvature, and they are therefore susceptible to the previously described limitations related to noise. Furthermore, the Hayakawa and Oguchi (2006) method does not incorporate drainage area information, which is an important parameter to consider when studying knickpoints over large spatial scales or when interpreting the retreat rates of these features.

Another set of methods exploits the use of k_s from Eq. (1) (or k_sn when calculated using a fixed value of θ) to extract knickpoints from slope–area plots, as reviewed by Neely et al. (2017). These methods suffer from limitations linked to slope–area scattering, noise sensitivity and difficulty in precisely locating knickpoints because of the stepped nature of drainage area (increasing instantaneously downstream when a new tributary reaches the river channel). To ameliorate problems with noise and data scattering, Bennett et al. (2016) devised a method that first calculates k_sn on channel profiles smoothed using the algorithm of Schwanghart and Scherler (2014). This derives k_sn either from the regression of slope–area plots or using the first-order derivative of χ plots. The method selects a knickpoint for which the ratio between downstream and upstream k_sn, averaged with two 2 km long serial windows, exceeds a factor of 2.

Neely et al. (2017) developed an algorithm focused on knickzone detection (KZ-Picker). Knickzones are selected from normalised profiles (using the approach of Perron and Royden, 2013) by comparison with a reference profile calculated for a defined concavity index (θ in Eq. 1). This reference profile is a line in χ-elevation space between the outlet and headwaters of the channel, and knickzones are then defined based on the deviation of the χ profile from the reference. After initial detection, knickzones are quantified by their relief (elevation drop) and adjusted using several filters or lumping-window parameters. This method is well adapted to detect knickzones that are composed of a base and a lip separating a steepened reach. An example of output produced by this algorithm and compared to ours is presented in Sect. 5.4.

Another method for extracting knickpoints has recently been implemented using TopoToolbox (Schwanghart and Scherler, 2014). Although unpublished, the code is available and also aims to reproducibly extract knickpoint locations from river profiles. It selects knickpoints by creating reference channel profiles that are concave up and then selecting knickpoints for which the actual channels are the most different from the reference channels. Although not based on the slope–area relationship, this method is perhaps the closest algorithmic attempt to match the knickpoint definition of early workers (e.g. Knopf, 1924). A sensitivity parameter defines the number of iterations and indirectly the number of knickpoints detected. After knickpoint extraction, a value is attributed to each identified knickpoint quantifying the divergence of the long profile from the reference profile. We discuss the similarities and differences of this method compared to our method in Sect. 6.

1.1.3 Motivation for a new method

Despite the large number of past approaches to selecting knickpoints, we have developed a new method because (i) many authors still select knickpoints based on qualitative interpretation of channel long profiles or slope–area data and we desired an open-source, reproducible method that has no reliance on proprietary software such as ArcGIS (e.g. Hayakawa and Oguchi, 2006) or MATLAB (e.g. Schwanghart and Scherler, 2014; Neely et al., 2017). (ii) Channel erosion is modelled to scale with discharge, and therefore we wished to use a method that includes discharge (or its proxy drainage area). (iii) Existing slope–area approaches make it difficult to pinpoint knickpoint locations (Fig. 1), and therefore we choose to use a χ-based approach. (iv) We wished to develop a method that not only selected knickpoint locations but also included metrics of changes in normalised channel steepness, as that metric is frequently used in tectonic geomorphology, and (v) we aimed to create a method allowing for differentiation between different knickpoint morphologies (e.g. slope break vs. vertical step).

Although the newest methods (Schwanghart and Scherler, 2014; Neely et al., 2017) meet a subset of these criteria, they both only describe a specific morphology of a knickpoint and/or knickzone and use indirect methods to quantify their magnitude (e.g. derived from comparison with a reference profile). Our aim here is to provide a method that selects locations, styles (e.g. vertical step, slope break) and magnitudes (e.g. main features or secondary ones) of knickpoints and knickzones that is free of manual selection in order to complement these existing methods that are more focused on identifying locations of a particular style of knickpoint and knickzone (e.g. waterfall).

We provide comparisons with two existing methods in Sect. 5.4. These have been chosen for the following reasons: (i) the knickpoint-extracting algorithms are open source (with the limitation of MATLAB licences), (ii) the methods are objective, reproducible and provide a quantification of knickpoint magnitude in order to compare it with ours, and (iii) Schwanghart and Scherler (2014) is purely based on channel morphology, while Neely et al. (2017) use the slope–area relationship and χ, thus providing a reasonable comparison of our algorithm with the range of existing methods.

2 Methods

An overview of our knickpoint identification method can be found in Fig. 2.

Figure 2 Flowchart of the knickpoint detection algorithm.

2.1 DEM preprocessing and river network extraction

Firstly, we fill the DEM using the filling algorithm of Wang and Liu (2006) to make sure that each cell has a flow direction and to avoid internal basins generated by DEM noise (e.g. Barnes et al., 2014). This approach is suitable for cases in which no feature is spuriously damming the DEM. Spurious damming can occur when vegetation, bridges or other features lead to high elevations over the channel when in fact the channel sits at a lower elevation. The filling process will create flat surfaces behind such spurious dams and will therefore hinder channel profile analysis.

If features that lead to spurious damming are present, we give users the option to use a breaching or carving algorithm. This excavates through spurious dams to avoid overfilling. The depression-breaching algorithm in our code is that created by Lindsay (2016) and adapted from Barnes (2016) within our method. It is also possible to supply the algorithm with preprocessed DEMs (e.g. Schwanghart and Scherler, 2017).

From the preprocessed, carved or filled DEM, we provide several methods of extracting the river network, including the DrEICH method (Clubb et al., 2014), a curvature method proposed by Pelletier (2013) and a method that uses a Wiener filter (Wiener, 1949) that combines elements of the methods of Pelletier (2013) and Passalacqua et al. (2010) first implemented by Grieve et al. (2016) and Clubb et al. (2017). Grieve et al. (2016) found this latter method least sensitive to DEM resolution. Finally, we include extraction based on a drainage area threshold more suitable for low-resolution DEMs (e.g. SRTM, ASTER) or large-scale studies in which the location of channel heads is less important. We also ensure during the preprocessing that no catchments are beheaded by the edge of the DEM, as the χ coordinate is a function of drainage area and therefore incomplete basins will have incorrect χ values.

2.2 k_sn extraction

Following channel extraction, we then calculate the χ coordinate for the resulting network. A key parameter that must be constrained prior to calculation of χ is the concavity index (θ). Changing the concavity index significantly affects values of the χ coordinate (e.g. Kirby and Whipple, 2012; Gasparini and Whipple, 2014; Mudd et al., 2018a) and therefore subsequent knickpoint extraction. We select the concavity index using a method developed by Mudd et al. (2018a). This method calculates the χ coordinates for a range of concavities within each watershed and determines the most likely concavity index by directly comparing the collinearity of points on each tributary with the trunk channel (Perron and Royden, 2013; Mudd et al., 2018a). This approach does not assume linearity in χ-elevation space and is therefore applicable in transient landscapes (Mudd et al., 2018a).

Once we determine θ values for each basin, we calculate χ and then use χ-elevation profiles to determine changes in k_sn, which is the gradient of the χ-elevation profile when we set A₀=1 (see Eq. 2). Theoretical work by Royden and Perron (2013) suggested that in eroding landscapes changes in erosion rates would be represented by changes in χ-elevation gradient between segments of channels that would be linear in χ-elevation space, which Royden and Perron (2013) called slope patches. Mudd et al. (2014) devised a statistical method that identified the most likely linear segments in χ-elevation space. This technique searched all possible combinations of channel pixels and used the corrected Akaike information criterion (AIC) (Akaike, 1974; Hurvich and Tsai, 1989) to balance the goodness of fit of linear segments against overfitting the data. Here we use this same algorithm to search for breaks in slope within the profile corresponding to knickpoint locations.

Knickpoints will manifest themselves as changes in the slope of these patches equivalent to the slope-break knickpoints of Kirby and Whipple (2012), whereas knickzones will be represented by patches with locally high gradients. That is, knickpoints and knickzones result in either changes in or locally high values of k_s (or k_sn if calculated with a fixed concavity index). The segmentation algorithm casts the profile as a series of linear segments, and each segment has a gradient and an intercept. The gradient reflects k_s of the segment and the intercept can be used to detect vertical-step knickpoints, as it detects elevation jumps between adjacent slope patches.

The method developed by Mudd et al. (2014) subsamples underlying topographic data iteratively: on each iteration nodes from the channel network are chosen randomly and segmentation is applied to this subset of nodes. The number of iterations is called n_MC. This iterative approach was taken because it significantly reduces the sensitivity of the results to user parameters (Mudd et al., 2014). The computational expense of the segmentation scales highly non-linearly with the number of nodes, so channel profiles are broken into subsections of length n_tg (called the “target nodes” in Mudd et al., 2014). The sampling of the underlying data on each iteration is random: after sample nodes are “skipped” randomly, the number of nodes skipped varies with a uniform distribution from zero to twice a parameter n_sk such that the mean “skip” is n_sk. We explore the sensitivity of the method to these parameters in the discussion.

The final k_sn values are an average of many iterations using different channel profiles subsampled from the raw data, as are intercepts of local segments. These averaged values are used to build segmented elevation. Each node then represents an average of the best-fit segments for every iteration of the segmentation routine (Fig. 3a):

$$\begin{array}{}\text{(5)}& z_{\mathrm{seg}i}=M_{\mathit{\chi }i}\times {\mathit{\chi }}_i+b_{\mathit{\chi }i},\end{array}$$

where i is the given node, z_seg its elevation on the segment, M_χ the average gradient of the segments and b_χ the averaged intercept of the segments. M_χ can be expressed with the following equation:

$$\begin{array}{}\text{(6)}& M_{\mathit{\chi }}={\left({\displaystyle \frac{E}{K\times A_{\mathrm{0}}^m}}\right)}^{\mathrm{1}/n}.\end{array}$$

We note here that M_χ is the same as k_sn if χ is calculated using A₀=1 m².

Figure 3 Extraction of normalised channel steepness (k_sn) from a river profile. (a) Example of best-fit segmentation (Mudd et al., 2014) where + symbols are individual data points and the coloured lines are the segments. (b) The associated plot of k_sn plotted as a function of χ coordinate. The segmentation output results in some noise due to iterative sampling of the channel network (+ symbols). A total variation denoising filter (Condat, 2013) is then applied on the signal to extract the main variations in k_sn.

2.3 Knickpoint extraction from k_sn data

2.3.1 Change point detection

Change point detection is a common technique used within many fields (e.g. time series analysis) and a number of statistical tools have been developed to identify change points, reviewed and described by Truong et al. (2018). In our case, the signal (k_sn) is by definition piecewise stationary, and abrupt changes occur between each segment (i.e. knickpoints). Change point detection algorithms aim to estimate and isolate the exact location of these boundaries between stationary patches. Method choice depends on the nature of the original dataset (e.g. noise intensity) and the number of changes we aim to extract (e.g. predetermined or unknown). In our case, although the segmentation algorithm of Mudd et al. (2014) can result in very sharp segment boundaries, in many cases the transitions between segments is fuzzy. We therefore have an unknown number of change points to detect from a variably noisy signal. We therefore choose to use a signal processing filter (Condat, 2013) to flatten the piecewise k_sn patches and discretise all potential change points. This algorithm identifies where k_sn and elevation are statistically varying the most within any transition zones. It also combines segments that have very small changes in k_sn relative to the noise in the data (Fig. 3b).

We denoise the data using a one-dimensional total variation denoising (TVD) filter, a signal processing filter adapted from an optimised algorithm by Condat (2013) solving the following equation:

$$\begin{array}{}\text{(7)}& \underset{x\in {\mathrm{\Re }}^N}{\mathrm{minimise}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{k=\mathrm{1}}^N{\left|y\left[k\right]-x\left[k\right]\right|}^{\mathrm{2}}+\mathit{\lambda }\sum _{k=\mathrm{1}}^{N-\mathrm{1}}\left|x[k+\mathrm{1}]-x\left[k\right]\right|,\end{array}$$

where N represents the number of samples (nodes) per population (in this case a river channel from source to the next higher-order stream or the outlet), y represents the raw signal $y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}},\mathrm{\dots }{y}_N$, in this case k_sn ordered by ascending χ within each river, x the denoised signal $x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},\mathrm{\dots }{x}_N$, referred to as denoised k_sn, and λ is a regularisation parameter (Condat, 2013). This method minimises variations, whereby the parameter λ must be real and greater than zero. Greater λ values result in less variation in the processed signal, and $\mathit{\lambda }\to +\mathrm{\infty }$ results in no variation in the processed signal whatsoever. The selection and sensitivity of this parameter are discussed in Sect. 5.1.

After denoising the data, our method then iterates through all nodes in each channel and identifies change points as any variation in the denoised k_sn data. These represent first-order knickpoints that we quantify by their change in denoised k_sn, which we call Δk_sn. Δk_sn is a quantitative measure of the magnitude of the slope-break component of the knickpoint (Fig. 4a). We refer to change points as knickpoints in the rest of the paper.

2.3.2 Combining knickpoints

Denoised k_sn data can still contain closely clustered steps in k_sn values, which may in fact represent a single knickpoint. We therefore use an algorithm to determine which of these clusters can be combined. Iterating through each river, the algorithm tests the neighbouring nodes of each raw knickpoint in a window that we call the “combining window”. If two knickpoints in the denoised k_sn data are within the combining window and both have the same sign of Δk_sn, the two knickpoints are merged and their magnitude summed. This process is repeated using newly merged knickpoints until no nodes are within the combining window or until a change in knickpoint sign (Fig. 4b). The combined knickpoint is then centred between the combined nodes. The width of the combining window (which we denote r_comb and is defined by a number of nodes rather than a flow distance) is a user-defined parameter, the selection of which we address in Sect. 5.1.

Figure 4 Knickpoint extraction from the denoised k_sn profiles. (a) The first step extracts all variations of k_sn, quantifying each with Δk_sn, which we call the “raw” knickpoint dataset. Negative and positive changes represent decreases or increases in k_sn, respectively. (b) After the detection of changes in k_sn, knickpoints are combined. All knickpoints within a node window will be combined, summing their values (i.e. a sum of Δk_sn). This process is repeated as long as the subsequent raw knickpoint is within a node window and as long as the polarity (i.e. if it is negative or positive) does not change.

2.3.3 Vertical-step knickpoint detection

Small variations between segments with similar k_sn values may be ignored by denoising, which may seem trivial if the aim is to isolate the main variations in channel steepness. However, this may lead to vertical-step knickpoints being missed if channel segments above and below the vertical-step knickpoint have similar k_sn values despite a jump in z_seg. We therefore use a second approach to extract knickpoints, allowing us to identify both slope-break and vertical-step knickpoints.

The algorithm calculates changes in z_seg using Eq. (5) in order to isolate the main jumps in profile elevation. We differentiate this value along the river nodes (Δz_seg) to detrend the elevation signal and focus on the stepped variations. For each node in the channel, the mean and standard deviation of Δz_seg is calculated within a window of surrounding nodes; the window width in nodes is called r_W. The nodes within the first and last half-windows are calculated using the first and last window, respectively. Δz_seg is then compared to the standard deviation of the nodes within the corresponding window multiplied by a coefficient (which we call T_σ), and the node is selected as a vertical-step knickpoint if Δz_seg is greater (Fig. 5b). This approach ensures that the selected vertical-step knickpoints show an anomalous increase in elevation. The selection of the window width and the coefficient is discussed in Sect. 5.1. We can then use Δz_seg as a quantitative measure of the size of each vertical-step knickpoint.

Figure 5 Extraction of knickpoints from the segmented elevation (Eq. 5). (a) Expression of a vertical-step knickpoint in a χ−z_seg profile compared to a slope-break knickpoint. (b) Representation of the identification window and the corresponding standard deviation around the reference node (in red). μ is the mean and T_σ the coefficient applied to the standard deviation. This process is repeated for each node. Reference nodes outside their own window are considered to be outliers.

2.4 Accuracy metrics

The accuracy of the method is assessed using a true positive (TP), false positive (FP) and false negative (FN) approach. This comparison method is often use to test algorithm performances on point data, such as channel heads (e.g. Orlandini et al., 2011; Clubb et al., 2014) or knickzone locations (e.g. Neely et al., 2017). We test the algorithm with these accuracy metrics using two sites where locations of hand-picked knickpoints based on field observations and river profiles are available. Knickpoints were identified at Santa Cruz Island (California, USA) by Neely et al. (2017), and we introduce a new dataset in the Quadrilátero Ferrífero, Minas Gerais, Brazil.

We define as TP a reference knickpoint detected by the algorithm, as FP a knickpoint detected by the algorithm that is not a reference knickpoint, and as FN reference knickpoints not detected by the algorithm. Neely et al. (2017) propose a fourth kind of prediction called “mixed” to assess the knickzone base and lip detection, whereby only one of the two knickzone boundaries is detected. We chose not to use this approach as we define a knickpoint as a point location showing an increase or decrease in k_sn or Δz_seg, which is more applicable to varying knickpoint morphologies. The definition of the different knickpoint predictions allows for the calculation of sensitivity, s, reliability, r, and metrics. We also add an overall quality metric, q, described in Heipke et al. (1997). The sensitivity can be expressed as

$$\begin{array}{}\text{(8)}& s={\displaystyle \frac{\sum \mathrm{TP}}{\sum \mathrm{TP}+\sum \mathrm{FN}}},\end{array}$$

where ∑TP and ∑FN are the sum of TP and FN. This metric measures the method's ability to detect a knickpoint that a user would have manually picked. s=1 implies the detection of all the locations of reference knickpoints. The reliability can be expressed as

$$\begin{array}{}\text{(9)}& r={\displaystyle \frac{\sum \mathrm{TP}}{\sum \mathrm{TP}+\sum \mathrm{FP}}},\end{array}$$

where ∑TP and ∑FP are the sum of TP and FP. This metric measures the occurrences of the method identifying knickpoints that a user would not have picked. The overall quality metric can be expressed as

$$\begin{array}{}\text{(10)}& q={\displaystyle \frac{\sum \mathrm{TP}}{\sum \mathrm{TP}+\sum \mathrm{FP}+\sum \mathrm{FN}}}.\end{array}$$

A q value of unity implies perfect agreement between algorithmically and hand-picked knickpoints. We focus on these metrics instead of the knickpoint magnitude, as it is more difficult to predict and is dependent on many parameters within the extraction of the Δk_sn values.

3 Test locations

In order to test the performance of our method, we extract knickpoints from two field sites with independently mapped knickpoint and knickzone locations. The first of these sites is Smugglers Basin on Santa Cruz Island (California, US), where knickpoints and knickzones were mapped by Neely et al. (2017) using a combination of fieldwork and supervised selection from river long profiles. Smugglers Basin is undergoing transient adjustment to climatic and tectonic signals (Neely et al., 2017). The second field site is located in the Quadrilátero Ferrífero (Minas Gerais, Brazil), where we present a new dataset of extracted knickpoint and knickzone locations from field observations and river profiles. Quadrilátero Ferrífero represents a more stable site in terms of climate and tectonics (e.g. Dorr, 1969; Salgado et al., 2008), and therefore knickpoints in this landscape have been linked instead to changes in lithology.

3.1 Santa Cruz Island, USA

The first calibration test site is the headwaters of the Smugglers Cove catchment, located in the SE of Santa Cruz Island, the largest of the California Channel Islands (California, USA; Fig. 6a). Lidar data at 1 m resolution are available in the basin via the 2010 US Geological Survey Channel Islands Lidar Collection, available from OpenTopography (opentopography.org).

The basin has a total relief of approximately 550 m and drains to the Pacific Ocean. Previous work has estimated uplift rates of ≈1 mm yr⁻¹ using dated terraces and fault activity (e.g. Pinter et al., 1998; Muhs et al., 2014), and the site has experienced regional sea-level variations (e.g. Schumann et al., 2016; Pinter et al., 2018). This, along with bedrock heterogeneity, has led to numerous knickzones in the catchment which have been mapped and tested against a previous knickzone extraction algorithm by Neely et al. (2017). A total of 18 knickzone bases and lips have been reported based on topographic expression and field observations across the whole catchment. As the Neely et al. (2017) algorithm is targeted specifically at knickzones, we compare the mapped knickzone bases and lips with those picked by our algorithm. Knickzone bases and lips are the equivalent of increases and decreases in k_sn, respectively.

We extracted channel heads using a curvature-based method of channel extraction, following Pelletier (2013) and Grieve et al. (2016). This method has an estimated accuracy of ≈10 m horizontally along drainage paths (Clubb et al., 2014). Before extracting channel steepness, we calculated the best-fit concavity index for the basin by maximising collinearity between the main-stem channel and the tributaries in χ-elevation space using the bootstrapping method of Mudd et al. (2018a): the best-fit θ at the site is 0.25.

3.2 Quadrilátero Ferrífero, Minas Gerais, Brazil

The second calibration test site is located in the eastern part of the Quadrilátero Ferrífero (QF, Brazil), in a basin draining the Caraça Range (Fig. 8). The QF is an area of relatively high elevation in southeastern Brazil, and the Caraça Range is its most pronounced topographic feature with a maximum elevation of ≈2100 m and maximum relief of ≈1500 m. Tectonic activity is thought to have ceased by ≈500 Ma (e.g. Dorr, 1969; Chemale et al., 1994; Alkmim and Marshak, 1998). Upstream areas are primarily underlain by resistant rocks (e.g. quartzites and banded iron formations), whereas less resistant rocks often underlie downstream areas (e.g. schists and phyllites). The association of mountainous topography and long-term tectonic stability have led to controversy in the post-orogenic evolution of the QF (Peifer Bezerra, 2018). The most accepted hypothesis is that differential denudation of lithologies with different resistance to denudation has led to a geomorphic differentiation whereby the uplands, underlain by strong rocks, are high because they have been denuded less and more slowly than their surroundings (e.g. Harder and Chamberlin, 1915; James, 1933; Varajão, 1991; Salgado et al., 2008; Peifer Bezerra, 2018). An alternative hypothesis is that the relief of the QF results from a complicated history of geographic cycles interrupted by epeirogenic uplift (e.g. King, 1956; Dorr, 1969; Barbosa, 1980).

Knickpoints are common features in the rivers flowing away from the Caraça Range (Fig. 8). These rivers have headwaters at high elevations (≈2000 m), and their long profiles display many convexities associated with substantial elevation drops (up to 1.4 km of descent over ≈15 km of downstream distance) and steep channel and hillslope gradients. These rivers flow over quartzite terrains, transitioning in their distal part to schists (see Supplement Sect. S5.2). The origin of these knickpoints is unresolved, being possibly the result of spatial variations in rock resistance or alternatively resulting from transient uplift signals that have failed to progress beyond quartzite units (Peifer Bezerra, 2018). We used a TanDEM-X DEM with 12 m resolution to extract knickpoints from the QF. Before extracting channel steepness, we estimated the best-fit concavity index as 0.15 using the methods presented in Mudd et al. (2018a).

4 Results

4.1 Performance at Santa Cruz Island

We carried out knickpoint extraction on Santa Cruz Island (Fig. 6b) initially with parameters detailed in Table 1; the full parameter file is available in the Supplement. As explained in Sect. 2, extraction prior to post-processing thinning generates a dense dataset of knickpoints both within and outside knickzones identified by the calibration dataset (see Supplement Sect. S5.1). Therefore, we apply a threshold approach to thin the dataset by removing small knickpoints. We set cut-off values of |Δk_sn|>0.8 and Δz_seg>2.1, whereby knickpoints smaller than these thresholds are ignored. These values are set for this case study with the specific aim of isolating the main knickpoints while matching with the calibration dataset. This approach is fully reproducible and does not involve manual picking of knickpoints.

Table 1 Parameter values used for the two field sites. Differences in parameter values between the two sites are due to differing DEM resolution (1 m for Santa Cruz Island and 12 m for the Ribeirão Caraça). Sensitivity to these parameters is described in Sect. 4.3. Note that although the parameter values have been carefully optimised for knickpoint analysis, we suggest the values below as defaults for each of these two data resolutions in order to allow for a rapid initial knickpoint extraction for other landscapes.

Our thinning procedure reduced the number of slope-break knickpoints from 398 to 160 and the number of vertical-step knickpoints from 40 to 17. This is a relatively high number of knickpoints compared to the calibration bases and lips (18 pairs). However, this disparity can partly be explained by the differences in methods: our algorithm details discrete changes in channel morphology, whereas the calibration knickzones are identified over longer channel reaches. Therefore, one mapped knickzone may contain several algorithmically identified knickpoints.

Neely et al. (2017) propose an error radius of 50 m around each base and lip in order to test the performance of their algorithm: we used the same approach when comparing our extracted knickpoints to the calibration data (Fig. 6b). A TP is determined as any knickpoint within the calibration knickzone or the corresponding 50 m radius. An FP is determined as any knickpoint which does not lie within this radius, and an FN is determined as a base or a lip which is not identified by our algorithm. The reliability, sensitivity and overall quality metrics are presented in Table 2. High sensitivity (s=0.93) but lower reliability (r=0.53) and overall quality (q=0.51) suggest that the algorithm detects the bulk of human-selected knickpoints, but also a significant amount of other knickpoint features. The implications of these results are discussed below.

Table 2 Accuracy metrics for calibration site I (Smugglers catchment, California, USA).

4.2 Performance at Quadrilátero Ferrífero, Minas Gerais, Brazil

The application of our method in the Ribeirão Caraça basin resulted in a dense dataset of knickpoints (n=252); see Table 1 for parameter values and the Supplement for the full parameter file. To thin this dataset, we removed knickpoints with attributes lower than the cut-off values of |Δk_sn|>0.8 and Δz_seg>2.1 for the slope-break and vertical-step knickpoints, respectively. This filtering procedure decreased the number of slope-break knickpoints from 252 to 108, whereas the number of vertical-step knickpoints diminished from 44 to 23. We tested the performance of our method compared to human-selected knickpoints for the Ribeirão Caraça basin using the metrics TP, FP and FN (Table 3). We used the same error radius as was used on Santa Cruz Island for consistency. These metrics (see Sect. 2.4) indicate that the sensitivity of our method is high for the Ribeirão Caraça basin (s=0.89), and thus the bulk of human-selected knickpoints are captured by our algorithm. On the other hand, the reliability (r=0.60) and the overall quality (q=0.56) are lower because the number of false positives is high, indicating that our algorithm determines a relatively high number of knickpoints compared to human selection. In summary, our algorithm captures knickpoints that are visually selected for the Ribeirão Caraça basin and many knickpoints that are not recognised by traditional field mapping of knickpoints, but are morphologically similar as defined by our algorithm.

Table 3 Accuracy metrics for calibration site II (Ribeirão Caraça basin, Caraça Range, QF, Brazil).

4.3 Sensitivity to algorithm parameters

One important parameter in our method of knickpoint detection is the concavity index (θ). The concavity index controls the magnitude of k_sn because it determines the values of χ (Eq. 6), and a higher concavity index will produce higher k_sn values for the same channel. We ran the algorithm on Santa Cruz Island for θ values ranging from 0.05 to 0.95 in steps of 0.05.

Because the value of θ affects the k_sn order of magnitude, λ must be adapted to keep denoising the signal. We therefore tested a wide range of λ values for each θ value. From these tests (see Supplement Sect. S4.1) we determined default λ values appropriate for a range of θ values. These default values are implemented internally in the code, but can be modified if needed. The sensitivities of knickpoint locations to θ using default λ values are presented in Fig. 9. This analysis shows that the general spread of the data, represented by its z_score (difference between the data point and the mean normalised by the standard deviation), is not significantly impacted by different θ values. However, the relative magnitude of each knickpoint, measured by changes in k_sn, depends on the chosen value of θ. Therefore, if the intention of the user is to find the spatial distribution of the largest knickpoints then it is essential that θ is picked with care (see Supplement Sect. S4.2 for more illustrations).

Because k_sn values are sensitive to the value of the concavity index, θ, it is important to note that basins with different θ values should be analysed separately to isolate knickpoint locations. Δk_sn values are therefore also dependent on the value of θ, so the relative magnitudes of knickpoints and knickzones should only be compared amongst basins with the same θ value. On the other hand, the locations of knickpoints and knickzones are relatively insensitive to θ so the method can be used to determine the spatial distribution of knickpoints across large areas even in the event that the concavity index may vary spatially.

The extraction of channel steepness will also be influenced by parameters in the segment fitting algorithm (Mudd et al., 2014): the number of target nodes (noted n_tg) and the average number of nodes skipped (noted n_sk). We therefore ran sensitivity analyses on these parameters testing every combination for the following ranges of values: from 5 to 120 n_tg and values of 1 to 4 for the n_sk parameter. Our results show that both of these parameters affect the segment lengths. Increasing either the number of n_tg or the n_sk parameter leads to longer segments (see Supplement Sect. S4.3 for more details). This affects the number of knickpoints detected. We also tested the number of Monte Carlo iterations (n_MC) processed for each segment from 5 to 500 and find that the results become insensitive to n_MC when n_MC>50.

The results of the vertical-step knickpoint detection can change with the size of the moving window that detects sudden changes in z_seg compared to neighbouring nodes (Sect. 2). We tested the following combination of parameters for vertical-step knickpoint detection: r_W from 10 to 200 nodes over intervals of 10 nodes and T_σ from 5 to 10 over intervals of 0.5. Our results show that the extraction is insensitive to r_W above a threshold minimum value of around 80 in our case. Below this value, the algorithm begins to identify steep channels as a succession of steps and will detect each node in the steep section as a knickpoint. We find that the number of extracted knickpoints becomes much higher if T_σ<6, whereas T_σ>8 results in very few knickpoints being detected. We therefore suggest selecting a value of $\mathrm{6}\le T_{\mathit{\sigma }}\le \mathrm{8}$.

The resolution of the DEM may also affect the location of extracted knickpoints and knickzones. We conducted a sensitivity analysis on raster resolution by resampling the original 1 m lidar-derived DEM into coarser grids to represent commonly available resolutions of 5 m (e.g. NED or NetMap), 10 m (e.g. NED or TanDEMX) and 30 m (ASTER or SRTM). Our results (see Supplement Sect. S4.7) show a decreasing number of detected knickpoints at coarser grid resolutions. This is directly linked to the amount of nodes in each river profile: as the resolution decreases, the number of nodes per river also decreases, meaning that fewer segments are used to extract k_sn. Therefore, fewer knickpoints are detected as knickpoints tend to be located near the segment boundaries. Furthermore, with lower-resolution grids the knickpoints that are detected tend to represent larger-scale variations in the channel profile. Vertical-step knickpoints also tend to be identified as steepened reaches rather than purely vertical regions of the channel profile, as the grid resolution prohibits the identification of small waterfalls. In order to show an overview of the algorithm performance in different field sites and DEM datasets, we extracted knickpoints from an additional test site using a 30 m DEM derived from SRTM (Supplement, Fig. S21).

5 Discussion

5.1 Selecting parameter values

Ideally our method for knickpoint detection could proceed without any human supervision. Due the method's sensitivity to grid resolution, roughness and the intrinsically heterogeneous nature of landscapes, the method does, however, retain some user-defined parameters. The sensitivity analysis performed on the Santa Cruz Island data (Sect. 4.3) indicates which of these must be selected with care.

We found that changing the concavity index does not change the location of the knickpoints substantially, but it does control their relative magnitude (Sect. 4.3), and therefore if the user is interested in knickpoint magnitude then θ should be selected carefully (e.g. Mudd et al., 2018a). The parameters linked to segmenting the χ-elevation profiles (Mudd et al., 2014) that affect results are the n_tg and n_sk parameters (Sect. 4.3). Increasing both of these increases the length of the segments, whereby setting these parameters to smaller values result in a large number of detected changes in k_sn which must thereafter be thinned. The one potential advantage of smaller segments is that more vertical-step knickpoints can be detected (i.e. waterfalls). Smaller segments also affect the relative values of knickpoint magnitude because short, steep reaches can be extracted and will generate high-magnitude Δk_sn knickpoints. If high values for the n_tg and n_sk parameters are used, the resulting knickpoint dataset will be sparser but will not necessarily detect local changes in k_sn due to local layers of hard rock or a change in erosion processes, for example. Larger segments are also less sensitive to topographic noise. After running sensitivity analyses, we recommend default parameters of n_tg=80 and n_sk=1.

Once segmentation is performed, we use the TVD routines to isolate changes in k_sn, which require an additional parameter (λ) to control the degree of denoising (Eq. 7). As the relative magnitude of k_sn is controlled by the θ value, we also determine the λ value for each value of θ that best isolates changes in k_sn based on our sensitivity analysis (Sect. 4.3). However, some landscapes that are either very gentle or steep may require changes to the λ value: low-relief landscapes may require a smaller λ value, whereas the opposite is true for steep landscapes. The user can check the efficacy of the selected λ value by plotting k_sn and denoised k_sn against χ or the flow distance. Guidance on the selection of λ is described in greater detail in the Supplement Sect. S4.1.

We also explored the possibility of using the TVD routine to denoise the river profile before extracting knickpoints in order to avoid dependency on the θ parameter. We applied the denoising routine on Δelevation in order to reduce the amount of variation. The intensity λ of denoising has to be manually selected and controls the amount of change from original data. Results from these tests are available in the Supplement (Figs. S18–S20). We found that additional denoising is still required during the Monte Carlo segment determination of Mudd et al. (2014). We suggest that prior smoothing of river profiles needs to be carefully considered, as it unavoidably leads to some modification of the existing profile. Users of our software may, if they wish, apply a technique for denoising river profiles prior to applying our method (e.g. Schwanghart and Scherler, 2017).

Figure 6 The test location on Santa Cruz Island, CA, USA. (a) Map of channel network extracted with the Pelletier method (Pelletier, 2013) and coloured by k_sn value calculated with Mudd et al. (2014). (b) Extracted knickpoints plotted after thinning the dataset as described in Sect. 4.1. The purple and green circles respectively represent the calibration knickzone bases and lips with the 50 m radius used for assessing algorithm performances. Stars and associated numbers are source numbers, which can be compared to Fig. 7. Topographic data are 1 m precision lidar DEM (see Supplement Sect. S1 for metadata) reprojected in WGS84 UTM zone 11N.

The width of the combining window can also be an important factor. As explained in Sect. 2, segment boundaries can still be fuzzy after the denoising process, generating successions of low-magnitude slope-break knickpoints. The combining window solves this issue by merging adjacent knickpoints within a certain radius. However, underestimating r_comb could result in retaining some of these low-magnitude knickpoints. Overestimating its size would possibly result in shifted knickpoint locations and misrepresentation of their magnitude if unrelated knickpoints are merged. In the case in which the DEM resolution is high enough to represent a close succession of knickpoints, we recommend carefully choosing a combining window smaller than the spacing between these features in order to avoid merging them.

Vertical-step knickpoint detection is controlled by two parameters: the window radius (r_W) and the standard deviation threshold for detecting anomalies (T_σ). Section 4.3 details the combined sensitivity analysis on these parameters and allows us to determine a set of values suitable for this analysis. However, if the user's specific aim is to detect vertical-step knickpoints (assuming that the DEM precision allows it), we recommend that users precisely constrain the standard deviation coefficient, the window size and the segment size in order to make sure that vertical-step knickpoints are extracted rather than slope break.

Although parameters in the method may be tuned and therefore the method can be supervised, it is reproducible. Workers using the method can report on the parameter values used and others can use these to reproduce the original results. One advantage of these adjustable parameters is that users can visually inspect outputs and change parameters such that the algorithm selects “obvious” knickpoints. However, we emphasise that this is not hand-picking of knickpoints: the algorithm output is a dense dataset of knickpoints. While sorting the dataset, once a threshold or statistical criterion is selected, all knickpoints and knickzones matching the selection are chosen. This means that one cannot eliminate knickpoints that qualitatively appear to be in the “wrong” place. As highlighted in Fig. 7, human-selected knickpoints and knickzones frequently produce biased knickpoint datasets that both include and exclude knickpoints and knickzones that have the same magnitude. We note that because the segmentation algorithm uses a Monte Carlo sampling routine (Mudd et al., 2014) there may be minor differences in results between two users, but by using a reasonable n_MC (>50) the results from one run to the next are nearly identical.

Figure 7 Knickpoint extraction for Santa Cruz Island, CA, USA, shown for the channel long profiles. These are the same knickpoints depicted in Fig. 6b. The stars and associated numbers correspond to the source numbers, and green and mauve circles correspond to the lips and bases of mapped knickpoints from Neely et al. (2017).

Figure 8 Knickpoint extraction on the Ribeirão Caraça basin (Caraça Range, QF, Brazil). (a) Map of knickpoints extracted with the algorithm after thinning the dataset as described in Sect. 4.2. Most of the calibration knickpoints are expressed by a succession of knickpoints detailing along-channel increases and/or decreases in k_sn. Streams depicted in panel (b) are shown as thick blue lines. (b) Longitudinal profile of the trunk stream (the Ribeirão Caraça river) highlighting the performance of the algorithm in picking along-channel breaks in steepness. (c) Example of a known waterfall (i.e. waterfall with a name) in the field: the Cascatinha waterfall. This waterfall features an elevation break of 40 m. Other known waterfalls include the Cascatona, Bocaina, Brumadinho and Quebra-ossos waterfalls.

Figure 9 Sensitivity of the knickpoint extraction to the concavity index (θ). As different values of θ result in different values of k_sn, we use a normalised z_score (i.e. the difference to the mean normalised by the standard deviation) to compare the overall spread of Δk_sn. The plot shows probability distributions of the z_score of Δk_sn represented by violin plots calculated with a kernel density estimation (bandwidth 0.20). The outliers and their relative magnitudes are affected by this parameter, whereas the general data distribution remains similar. The “min” and “max” stated above and below the violin plots respectively represents the minimum and maximum ΔM_χ for each run.

5.2 Quantification and selection of knickpoints

The aim of extracting knickpoints is mainly to link knickpoint location and magnitude to a specific event resulting in landscape transience (e.g. Crosby and Whipple, 2006). Therefore, an important step is to isolate the most significant knickpoint features from the dense raw dataset in order to interpret landscape evolution, which can be done using knickpoint magnitude. Knickpoint magnitude may be affected by the calculation of k_sn using the gradient of segments in χ-elevation space. Depending on the relief, particularly with a high value of θ, the absolute values of χ coordinates and associated elevation can differ by an order of magnitude. If the values of χ are low compared to the values for elevation, any changes in elevation at a knickpoint will result in a much higher segment gradient than if the χ values are of a similar magnitude as the elevation. This can result in the exaggeration of knickpoint magnitude in high-relief landscapes, for example, where it is more likely that χ values will be lower than the elevation values, eventually resulting in a bias during the sorting. We therefore suggest that, in such cases, A₀ from Eq. (3) should be set such that the value of the χ coordinate is the same order of magnitude as the elevation. However, if A₀≠1, then the gradient of the segment corresponds to M_χ in Eq. (6) rather than to k_sn. We wish to emphasise that this does not change the relative ordering between knickpoints. We illustrate this relationship by running a simple sensitivity analysis on the Santa Cruz Island dataset, with a range of A₀ varying from 1 to 500 (Fig. 10). This sensitivity analysis shows that, as A₀ is increased, the extreme values of Δk_sn within the dataset are reduced so that the effect of low absolute χ values on the gradient calculation is diminished. As for θ (see Sect. 4.3), knickpoint absolute magnitude (i.e. the direct value of Δk_sn and Δz_seg) cannot be compared if calculated with different A₀ from Eq. (3). However, the location of the isolated main knickpoints can still be compared.

Figure 10 The effect of varying A₀ on knickpoint extraction (Eq. 3). The reference area (A₀) will affect knickpoint magnitude and can be increased to reduce exaggerations in χ-elevation gradients. Changing A₀ does not affect the relative order of knickpoints: the largest knickpoints remain the largest for all values of A₀. Increasing A₀, however, reduces the spread in the z_score of the changes in channel steepness. This value has to be set only if necessary (e.g. if the high-gradient effect is important): A₀≠1 implies that the magnitude is not Δk_sn but ΔM_χ from Eq. (6). Moreover, overestimating A₀ can mask knickpoints that would be detected with A₀=1 m². The “min” and “max” stated above and below the violin plots represent the minimum and maximum ΔM_χ for each run.

Our sensitivity analyses suggest that two different approaches may be used to select knickpoints. The first of these is that a single θ and A₀ can be fixed for an entire landscape: the knickpoint magnitudes can directly be used to isolate the main knickpoint locations and relative importance. However, this approach may lead to some errors due to inevitable landscape heterogeneity over larger scales. The second approach is to calculate θ and A₀ values separately for individual basins, which allows knickpoints to be extracted with greater precision than if a single value is set for the entire landscape. However, this approach means that the knickpoint extraction has to be processed independently for each catchment, and only the location (e.g. latitude, longitude, elevation) is comparable between different catchments. Which approach is taken is dependent on the aims of each particular study and should be carefully considered on a case-by-case basis.

5.3 Knickpoint and knickzone morphology

Along with the calculation of knickpoint magnitude, our algorithm allows for the characterisation of knickpoint morphology. We can identify different knickpoint or knickzone types by (i) identifying locations where k_sn increases downstream (positive slope-break knickpoints) or (ii) identifying locations where k_sn decreases (negative slope-break knickpoints) and (iii) identifying locations where a sudden change in elevation occurs (vertical step knickpoints). This approach is suitable to identify the most common morphologies described in the literature (e.g. Haviv et al., 2010; Kirby and Whipple, 2012). However, we wish to emphasise that this algorithm can also be used to focus on one particular knickpoint morphology. For example, the classical convex-upwards knickpoint expression (e.g. Knopf, 1924) can be isolated by only displaying the knickpoints with a drop in Δk_sn (Fig. 11b). In order to examine steepened reaches or knickzones, we can also isolate locations where Δk_sn increases. Finally, waterfall detection can be achieved, if the resolution of the DEM allows it, by focusing on locations with a jump in z_seg. We provide all these different knickpoint types for the Smugglers catchment in the Supplement (Fig. S12).

5.4 Comparison with other knickpoint extraction techniques

For each of our two study sites, we have presented performance metrics of our method compared to knickpoints selected by humans. We find that our method has a high sensitivity, meaning that nearly all human-identified knickpoints were selected by the algorithm, but a lower reliability. This suggests that our algorithm also identifies many changes in channel steepness which are not selected as knickpoints through field mapping techniques. This raises the question of whether the algorithmic selection of knickpoints is more or less trustworthy than those selected by humans.

Figure 11 Comparison of results from the Smugglers catchment for our algorithm and the most recent similar ones. (a) Results for a single source from KZ-Picker (Neely et al., 2017) and our results. The results from Neely et al. (2017) are directly taken from their study to ensure objectivity. Only the slope-break knickpoints are displayed to make the comparison valid. (b) Basin-wide comparison between our algorithm outputs and the one recently implemented in Schwanghart and Scherler (2014) using a tolerance of 5. We only display the knickpoints showing a decrease in k_sn in order to provide a relevant comparison with the knickpoint morphology detected by Schwanghart and Scherler (2014). Differences in channel length are due to different methods for extracting channel heads between the two techniques.

Knickpoints identified for geomorphic studies should be reproducible, in that two workers should be able to select the same locations and magnitudes from the same river profile. This is challenging when mapping features in the field, as different workers may have different criteria for what constitutes a knickpoint. Furthermore, knickpoint selection should be objective: the same morphological criteria should be used to identify all features in the dataset. A common problem with field mapping by humans is that some specific features are picked in order to interpret a signal, whereas others with a similar morphology may be omitted. Our approach allows for the production of an objective dataset of knickpoint locations and magnitudes that can later be correlated by the user with process-based interpretations. Algorithmic extraction also allows for coverage of much larger areas compared to field mapping that can later be calibrated with additional data (e.g. Crosby and Whipple, 2006). As illustrated by our accuracy metrics, our algorithm produces datasets significantly denser than hand-picked knickpoints. However, it is possible to thin the number of knickpoints by applying threshold metric values selected based on statistical criteria and making the number of identified features similar to human-picked datasets. Such a process is objective in the sense that no hand selection is involved; only the morphology drives the thinning.

To provide a full assessment of our methods, we compare the output to that generated by two other algorithms as explained in Sect. 1.1.3: TopoToolbox (Schwanghart and Scherler, 2014) and KZ-Picker (Neely et al., 2017). Figure 11a expresses the differences between KZ-Picker and our algorithm for a single channel, whereby KZ-Picker identifies the main knickzone (in red) and quantifies its magnitude by the difference in elevation between the toe and lip of the knickzone. The purpose of the KZ-Picker is to find broad zones of steepened channels and is less granular than our method (e.g. Sect. 4.1). It is also not constructed to identify discrete vertical-step knickpoints. Because the raw output from our algorithm is denser than the KZ-Picker, the main knickpoints from our algorithm require more sorting based on their magnitudes, which results in extra steps to explore the data.

Figure 11b provides a basin-wide comparison of our outputs with those from TopoToolbox (Schwanghart and Scherler, 2014), with the tolerance parameter of the TopoToolbox method fixed to 5. In order to ensure that the comparison is valid we only compare it to our negative Δk_sn knickpoints, which should quantify similar features. The TopoToolbox method effectively identifies the main knickpoints expressed by the difference to an idealised profile that is concave up. However, reducing the tolerance parameter increases the number of knickpoints detected (e.g. 10: 12 knickpoints, 5: 44 knickpoints, 1: 343 and 0.1: 2234), meaning that the TopoToolbox method can result in a network of knickpoints that has a similar density to our method. However, the TopoToolbox method relies on profiles in elevation plotted against flow distance, so further processing is required to analyse changes in channel steepness using this method. Because the selection of knickpoints in this method is not normalised for drainage area, the largest knickpoints selected may not correspond to the largest changes in channel steepness. However, it has fewer parameters and is more computationally efficient than our method.

While the KZ-Picker and the TopoToolbox methods are well adapted for identifying specific types of knickpoint, neither allows for the separate identification and quantification of positive slope-break, negative slope-break and vertical-step knickpoints. Each method produces slightly different data products that can be used to interpret different components of the channel network, making these methods complementary.

Finally, we chose to build our change point detection method using the TVD routine (Condat, 2013). However, as explained in Sect. 2, alternative methods could be used. The algorithm therefore provides the raw data before the TVD routine, meaning that these data can be ingested by other change point detection techniques, e.g. the methods reviewed in Truong et al. (2018) and the associated open-source code.

6 Conclusions

We have developed a new method for extracting knickpoints and knickzones from topographic data. Our method extracts slope-break knickpoint locations using changes in channel steepness k_sn calculated by combining a statistical method for segmenting channels into reaches of different channel steepness (Mudd et al., 2014) and a recently introduced denoising technique (Condat, 2013). The method also identifies vertical-step knickpoints by searching for breaks in elevation between channel segments of similar channel steepness. Our algorithms provide a dense dataset of objectively extracted knickpoint locations, along with the relative magnitude of each knickpoint defined by either the change in channel steepness (for slope-break knickpoints) or the jump in elevation (for vertical-step knickpoints) to quantify knickpoint morphologies.

We tested our algorithm on two datasets for which knickpoints were independently field mapped and found that our method successfully extracted the human-identified knickpoints in the vast majority of cases. In general the method identifies more knickpoints compared to field mapping, as illustrated by our accuracy metrics, especially in the case of knickzones in which one broad steepened reach may result in multiple discrete segments in χ-elevation space. We provide tools for sorting and thinning the dense dataset in order to isolate the most significant breaks in the channel profile without involving any human-based selection. Resulting knickpoints can be compared with lithological, climatic or tectonic datasets. Our method therefore provides an objective, systematic and reproducible technique for quantifying knickpoints and knickzones, which can then be used to inform process-based interpretations of landscape evolution.