Articles | Volume 7, issue 1
Research article
18 Feb 2019
Research article |  | 18 Feb 2019

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

Boris Gailleton, Simon M. Mudd, Fiona J. Clubb, Daniel Peifer, and Martin D. Hurst

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. Hack1960; Willett and Brandon2002). In unglaciated landscapes, the main driver of erosion is the river system (Ahnert1970), 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. Gilbert1877) to the generalised recognition that river profiles reflect varying erosion processes (e.g. Mackin1948; Hack1960; Howard1965; Howard et al.1994; Dietrich et al.2003; Kirby and Whipple2012).

In a river system, topographic steady state requires spatially stable rock uplift and climatic conditions over a long period of time (Willett and Brandon2002). 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 Whipple2012), complex mantle processes inducing vertical motions (e.g. Faccenna and Becker2010; Braun2010), uplift driven by differential rock density (Braun et al.2014) and base-level variations linked to eustatic variations (e.g. Powell1875; Lambeck and Chappell2001; 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 Perron2013), eventually driving drainage network reorganisation and resulting in additional transient signals (e.g. Mather2000; Castelltort et al.2012; Willett et al.2014; Whipple et al.2017b; Mudd2017). Moreover, river profiles are also affected by intrinsic landscape properties, such as fracture density (e.g Whipple2002) or differential lithology (e.g. Stock and Montgomery1999; Forte et al.2016), which can also lead to morphological adjustment of the channel (e.g. Kirby and Whipple2012). 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. Barrell1920; Bascom1921), 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 Whipple2006; Abbühl et al.2011; Kirby and Whipple2012), 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 Boulton2012) or climatically triggered base-level fall (e.g. Crosby and Whipple2006; 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 Slingerland1996; Stock and Montgomery1999; Kirby et al.2003; Duvall2004). 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 Whipple2012; 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:

(1) S = k s A - θ ,

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 ks is called the steepness index, as it sets the overall gradient of the channel, and a number of authors have noted that ks 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 ks (Fig. 1a).

Figure 1Different 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=dz/dx where z is elevation and x is distance along the channel (e.g. Whipple et al.2017a), resulting in

(2) z ( x ) = z ( x b ) + k s A 0 θ x b x A 0 A ( x ) θ d x ,

where A0 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):

(3) χ = x b x A 0 A ( x ) θ d x .

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

(4) z ( x ) = z ( x b ) + k s A 0 θ χ .

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 Royden2013; 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 Oguchi2006; 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 ks from Eq. (1) (or ksn 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 ksn on channel profiles smoothed using the algorithm of Schwanghart and Scherler (2014). This derives ksn 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 ksn, 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 Royden2013) 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 Scherler2014). 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. Knopf1924). 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 Oguchi2006) or MATLAB (e.g. Schwanghart and Scherler2014; 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 Scherler2014; 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 2Flowchart 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 Scherler2017).

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 (Wiener1949) 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 ksn 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 Whipple2012; Gasparini and Whipple2014; 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 Royden2013; 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 ksn, which is the gradient of the χ-elevation profile when we set A0=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) (Akaike1974; Hurvich and Tsai1989) 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 ks (or ksn 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 ks 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 nMC. 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 ntg (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 nsk such that the mean “skip” is nsk. We explore the sensitivity of the method to these parameters in the discussion.

The final ksn 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):

(5) z seg i = M χ i × χ i + b χ i ,

where i is the given node, zseg 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:

(6) M χ = E K × A 0 m 1 / n .

We note here that Mχ is the same as ksn if χ is calculated using A0=1 m2.

Figure 3Extraction of normalised channel steepness (ksn) 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 ksn 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 (Condat2013) is then applied on the signal to extract the main variations in ksn.


2.3 Knickpoint extraction from ksn 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 (ksn) 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 (Condat2013) to flatten the piecewise ksn patches and discretise all potential change points. This algorithm identifies where ksn and elevation are statistically varying the most within any transition zones. It also combines segments that have very small changes in ksn 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:

(7) minimise x N 1 2 k = 1 N y [ k ] - x [ k ] 2 + λ k = 1 N - 1 x [ k + 1 ] - x [ k ] ,

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 y1,y2,y3,yN, in this case ksn ordered by ascending χ within each river, x the denoised signal x1,x2,x3,xN, referred to as denoised ksn, and λ is a regularisation parameter (Condat2013). 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 λ+ 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 ksn data. These represent first-order knickpoints that we quantify by their change in denoised ksn, which we call Δksn. Δksn 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 ksn data can still contain closely clustered steps in ksn 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 ksn data are within the combining window and both have the same sign of Δksn, 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 rcomb 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 4Knickpoint extraction from the denoised ksn profiles. (a) The first step extracts all variations of ksn, quantifying each with Δksn, which we call the “raw” knickpoint dataset. Negative and positive changes represent decreases or increases in ksn, respectively. (b) After the detection of changes in ksn, knickpoints are combined. All knickpoints within a node window will be combined, summing their values (i.e. a sum of Δksn). 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 ksn 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 ksn values despite a jump in zseg. 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 zseg using Eq. (5) in order to isolate the main jumps in profile elevation. We differentiate this value along the river nodes (Δzseg) to detrend the elevation signal and focus on the stepped variations. For each node in the channel, the mean and standard deviation of Δzseg is calculated within a window of surrounding nodes; the window width in nodes is called rW. The nodes within the first and last half-windows are calculated using the first and last window, respectively. Δzseg 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 Δzseg 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 Δzseg as a quantitative measure of the size of each vertical-step knickpoint.

Figure 5Extraction of knickpoints from the segmented elevation (Eq. 5). (a) Expression of a vertical-step knickpoint in a χzseg 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 ksn or Δzseg, 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

(8) s = TP TP + FN ,

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

(9) r = TP TP + FP ,

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

(10) q = TP TP + FP + FN .

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 Δksn 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. Dorr1969; 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 (

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−1 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 ksn, 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. Dorr1969; Chemale et al.1994; Alkmim and Marshak1998). 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 Bezerra2018). 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 Chamberlin1915; James1933; Varajão1991; Salgado et al.2008; Peifer Bezerra2018). An alternative hypothesis is that the relief of the QF results from a complicated history of geographic cycles interrupted by epeirogenic uplift (e.g. King1956; Dorr1969; Barbosa1980).

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 Bezerra2018). 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 ksn|>0.8 and Δzseg>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 1Parameter 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.

Download Print Version | Download XLSX

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 2Accuracy metrics for calibration site I (Smugglers catchment, California, USA).

Download Print Version | Download XLSX

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 ksn|>0.8 and Δzseg>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 3Accuracy metrics for calibration site II (Ribeirão Caraça basin, Caraça Range, QF, Brazil).

Download Print Version | Download XLSX

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 ksn because it determines the values of χ (Eq. 6), and a higher concavity index will produce higher ksn 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 ksn 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 zscore (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 ksn, 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 ksn 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. Δksn 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 ntg) and the average number of nodes skipped (noted nsk). We therefore ran sensitivity analyses on these parameters testing every combination for the following ranges of values: from 5 to 120 ntg and values of 1 to 4 for the nsk parameter. Our results show that both of these parameters affect the segment lengths. Increasing either the number of ntg or the nsk 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 (nMC) processed for each segment from 5 to 500 and find that the results become insensitive to nMC when nMC>50.

The results of the vertical-step knickpoint detection can change with the size of the moving window that detects sudden changes in zseg compared to neighbouring nodes (Sect. 2). We tested the following combination of parameters for vertical-step knickpoint detection: rW 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 rW 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 6Tσ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 ksn. 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 ntg and nsk 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 ksn 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 Δksn knickpoints. If high values for the ntg and nsk parameters are used, the resulting knickpoint dataset will be sparser but will not necessarily detect local changes in ksn 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 ntg=80 and nsk=1.

Once segmentation is performed, we use the TVD routines to isolate changes in ksn, which require an additional parameter (λ) to control the degree of denoising (Eq. 7). As the relative magnitude of ksn is controlled by the θ value, we also determine the λ value for each value of θ that best isolates changes in ksn 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 ksn and denoised ksn 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 Scherler2017).

Figure 6The test location on Santa Cruz Island, CA, USA. (a) Map of channel network extracted with the Pelletier method (Pelletier2013) and coloured by ksn 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 rcomb 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 (rW) 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 nMC (>50) the results from one run to the next are nearly identical.

Figure 7Knickpoint 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 8Knickpoint 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 ksn. 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 9Sensitivity of the knickpoint extraction to the concavity index (θ). As different values of θ result in different values of ksn, we use a normalised zscore (i.e. the difference to the mean normalised by the standard deviation) to compare the overall spread of Δksn. The plot shows probability distributions of the zscore of Δksn 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 Whipple2006). 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 ksn 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, A0 from Eq. (3) should be set such that the value of the χ coordinate is the same order of magnitude as the elevation. However, if A0≠1, then the gradient of the segment corresponds to Mχ in Eq. (6) rather than to ksn. 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 A0 varying from 1 to 500 (Fig. 10). This sensitivity analysis shows that, as A0 is increased, the extreme values of Δksn 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 Δksn and Δzseg) cannot be compared if calculated with different A0 from Eq. (3). However, the location of the isolated main knickpoints can still be compared.

Figure 10The effect of varying A0 on knickpoint extraction (Eq. 3). The reference area (A0) will affect knickpoint magnitude and can be increased to reduce exaggerations in χ-elevation gradients. Changing A0 does not affect the relative order of knickpoints: the largest knickpoints remain the largest for all values of A0. Increasing A0, however, reduces the spread in the zscore of the changes in channel steepness. This value has to be set only if necessary (e.g. if the high-gradient effect is important): A0≠1 implies that the magnitude is not Δksn but ΔMχ from Eq. (6). Moreover, overestimating A0 can mask knickpoints that would be detected with A0=1 m2. 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 A0 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 A0 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 ksn increases downstream (positive slope-break knickpoints) or (ii) identifying locations where ksn 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 Whipple2012). 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. Knopf1924) can be isolated by only displaying the knickpoints with a drop in Δksn (Fig. 11b). In order to examine steepened reaches or knickzones, we can also isolate locations where Δksn increases. Finally, waterfall detection can be achieved, if the resolution of the DEM allows it, by focusing on locations with a jump in zseg. 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 11Comparison 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 ksn 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 Whipple2006). 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 Scherler2014) 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 Scherler2014), 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 Δksn 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 (Condat2013). 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 ksn 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 (Condat2013). 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.

Code and data availability

Code used for analysis is located in the LSDTopoTools GitHub repository at (Mudd et al.2018b), and scripts for visualising the results can be found at (Mudd et al.2019a). We have also provided documentation detailing how to install and run the software, which can be found at (Mudd et al.2019b). As part of the Supplement we have also provided example parameter files which can be used to reproduce the results of all analyses performed in this study.


The supplement related to this article is available online at:

Author contributions

BG designed the study with contributions from SMM, FJC, DP and MDH . BG designed the algorithms and wrote the code with contributions from SMM and FJC. BG and DP ran the analysis on test sites. BG wrote the paper with contributions from SMM, FJC, DP and MDH.

Competing interests

The authors declare that they have no conflict of interest.


We thank Emma Graf for testing the software. Boris Gailleton was funded by European Union initial training grant 674899 – SUBITOP. Simon M. Mudd was supported by NERC grant NE/J009970/1, Fiona J. Clubb was supported by a Geo.X fellowship and Daniel Peifer had support from the Coordination for the Improvement of Higher Education Personnel (CAPES) under a Science without Borders fellowship BEX 12000/13-2. We thank the German Aerospace Center (DLR) for granting access to TanDEM-X data as part of the project DEM_GEOL1345. We would like to thank Giulia Sofia, Wolgang Schwanghart, Stefan Hergarten and an anonymous reviewer for their helpful comments and suggestions.

Edited by: Giulia Sofia
Reviewed by: Wolfgang Schwanghart, Stefan Hergarten, and one anonymous referee


Abbühl, L. M., Norton, K. P., Jansen, J. D., Schlunegger, F., Aldahan, A., and Possnert, G.: Erosion rates and mechanisms of knickzone retreat inferred from 10Be measured across strong climate gradients on the northern and central Andes Western Escarpment, Earth Surf. Proc. Land., 36, 1464–1473,, 2011. a

Ahnert, F.: Functional relationships between denudation, relief, and uplift in large, mid-latitude drainage basins, Am. J. Sci., 268, 243–263,, 1970. a

Akaike, H.: A New Look at the Statistical Model Identification, IEEE T. Automat. Contr., 19, 716–723,, 1974. a

Alkmim, F. F. and Marshak, S.: Transamazonian Orogeny in the Southern São Francisco Craton Region, Minas Gerais, Brazil: evidence for Paleoproterozoic collision and collapse in the Quadrilátero Ferrifero, Precambrian Res., 90, 29–58,, 1998. a

Attal, M., Tucker, G. E., Whittaker, A. C., Cowie, P. A., and Roberts, G. P.: Modelling fluvial incision and transient landscape evolution: Influence of dynamic Channel adjustment, J. Geophys. Res.-Earth, 113, F03013,, 2008. a

Attal, M., Cowie, P. A., Whittaker, A. C., Hobley, D., Tucker, G. E., and Roberts, G. P.: Testing fluvial erosion models using the transient response of bedrock rivers to tectonic forcing in the Apennines, Italy, J. Geophys. Res.-Earth, 116, F02005,, 2011. a

Baldwin, J. A., Whipple, K. X., and Tucker, G. E.: Implications of the shear stress river incision model for the timescale of postorogenic decay of topography, J. Geophys. Res.-Sol. Ea., 108, B3,, 2003. a

Barbosa, G.: Superfícies de erosão no Quadrilátero Ferrífero, Minas Gerais, Braz. J. Geol., 10, 89–101, 1980. a

Barnes, R.: RichDEM: Terrain Analysis Software, available at: (last access: 8 August 2018), 2016. a

Barnes, R., Lehman, C., and Mulla, D.: Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models, Comput. Geosci., 62, 117–127,, 2014. a

Barrell, J.: The Piedmont terraces of the Northern Appalachians, Am. J. Sci., 49, 327–362,, 1920. a

Bascom, F.: Cycles of Erosion in the Piedmont Province of Pennsylvania, J. Geol. 29, 540–559,, 1921. a

Baynes, E. R. C., Attal, M., Niedermann, S., Kirstein, L. A., Dugmore, A. J., and Naylor, M.: Erosion during extreme flood events dominates Holocene canyon evolution in northeast Iceland, P. Natl. Acad. Sci. USA, 112, 2355–2360,, 2015. a

Baynes, E. R. C., Lague, D., Attal, M., Gangloff, A., Kirstein, L. A., and Dugmore, A. J.: River self-organisation inhibits discharge control on waterfall migration, Sci. Rep., 8, 2444,, 2018. a

Bennett, G. L., Miller, S. R., Roering, J. J., and Schmidt, D. A.: Landslides, threshold slopes, and the survival of relict terrain in the wake of the Mendocino Triple Junction, Geology, 44, 363–366,, 2016. a

Bishop, P. and Goldrick, G.: Lithology and the evolution of bedrock rivers in post-orogenic settings: constraints from the high-elevation passive continental margin of SE Australia, Geol. Soc. Spec. Publ., 346, 267–287,, 2010. a

Braun, J.: The many surface expressions of mantle dynamics, Nat. Geosci., 3, 825–833,, 2010. a

Braun, J., Simon-Labric, T., Murray, K. E., and Reiners, P. W.: Topographic relief driven by variations in surface rock density, Nat. Geosci., 7, 534–540,, 2014. a

Castelltort, S., Goren, L., Willett, S. D., Champagnac, J.-D., Herman, F., and Braun, J.: River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain, Nat. Geosci., 5, 744–748,, 2012. a

Chemale, F., Rosière, C. A., and Endo, I.: The tectonic evolution of the Quadrilátero Ferrífero, Minas Gerais, Brazil, Precambrian Res., 65, 25–54,, 1994. a

Clubb, F. J., Mudd, S. M., Milodowski, D. T., Hurst, M. D., and Slater, L. J.: Objective extraction of channel heads from high-resolution topographic data, Water Resour. Res., 50, 4283–4304,, 2014. a, b, c

Clubb, F. J., Mudd, S. M., Milodowski, D. T., Valters, D. A., Slater, L. J., Hurst, M. D., and Limaye, A. B.: Geomorphometric delineation of floodplains and terraces from objectively defined topographic thresholds, Earth Surf. Dynam., 5, 369–385,, 2017. a

Condat, L.: A Direct Algorithm for 1D Total Variation Denoising, IEEE Signal Proc. Let., 20, 1054–1057,, 2013. a, b, c, d, e, f

Cook, K. L., Turowski, J. M., and Hovius, N.: A demonstration of the importance of bedload transport for fluvial bedrock erosion and knickpoint propagation, Earth Surf. Proc. Land., 38, 683–695,, 2013. a

Crosby, B. T. and Whipple, K. X.: Knickpoint initiation and distribution within fluvial networks: 236 waterfalls in the Waipaoa River, North Island, New Zealand, Geomorphology, 82, 16–38,, 2006. a, b, c, d

Davis, W. M.: The Geographical Cycle, Geogr. J., 14, 481–504,, 1889. a

DiBiase, R. A., Whipple, K. X., Heimsath, A. M., and Ouimet, W. B.: Landscape form and millennial erosion rates in the San Gabriel Mountains, CA, Earth Planet. Sci. Lett., 289, 134–144,, 2010. a

Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Stock, J. D., Heimsath, A. M., and Roering, J. J.: Geomorphic transport laws for predicting landscape form and dynamics, Geophysical Monograph Series, American Geohpysical Union, 135, 103–132,, 2003. a

Dorr, J. V. N.: Physiographic, stratigraphic, and structural development of the Quadrilatero Ferrifero, Minas Gerais, Brazil, Geological Survey Professional Paper, p. 110, U.S. Government Printing Office, Washington, D.C., Minas Gerais, Brazil, 1969. a, b, c

Duvall, A.: Tectonic and lithologic controls on bedrock channel profiles and processes in coastal California, J. Geophys. Res., 109, F03002,, 2004. a

Faccenna, C. and Becker, T. W.: Shaping mobile belts by small-scale convection, Nature, 465, 602–605,, 2010. a

Flint, J. J.: Stream gradient as a function of order, magnitude, and discharge, Water Resour. Res., 10, 969–973,, 1974. a

Forte, A. M., Yanites, B. J., and Whipple, K. X.: Complexities of landscape evolution during incision through layered stratigraphy with contrasts in rock strength, Earth Surf. Proc. Land., 41, 1736–1757,, 2016. a

Gallen, S. F., Wegmann, K. W., and Bohnenstiehl, D. W. R.: Miocene rejuvenation of topographic relief in the southern Appalachians, GSA Today, 23, 4–10,, 2013. a

Gasparini, N. M. and Whipple, K. X.: Diagnosing climatic and tectonic controls on topography: Eastern flank of the northern Bolivian Andes, Lithosphere, 6, 230–250,, 2014. a

Gilbert, G.: Geology of the Henry Mountains, USGS Unnumbered Series, Government Printing Office, Washington, D.C., USA, 1877. a, b

Goldrick, G. and Bishop, P.: Regional analysis of bedrock stream long profiles: Evaluation of Hack's SL form, and formulation and assessment of an alternative (the DS form), Earth Surf. Proc. Land., 32, 649–671,, 2007. a

Gonga-Saholiariliva, N., Gunnell, Y., Harbor, D., and Mering, C.: An automated method for producing synoptic regional maps of river gradient variation: Procedure, accuracy tests, and comparison with other knickpoint mapping methods, Geomorphology, 134, 394–407,, 2011. a

Grieve, S. W. D., Mudd, S. M., Milodowski, D. T., Clubb, F. J., and Furbish, D. J.: How does grid-resolution modulate the topographic expression of geomorphic processes?, Earth Surf. Dynam., 4, 627–653,, 2016. a, b, c

Hack, J.: Studies of longitudinal profiles in Virginia and Maryland, U.S. Geological Survey Professional Paper 294-B, United States Government Printing Office, Washington, D.C., USA, 1957. a

Hack, J. T.: Interpretation of erosional topography in humid temperate regions, Am. J. Sci., 258, 80–97, 1960. a, b

Harder, E. C. and Chamberlin, R. T.: The Geology of Central Minas Gerais, Brazil: Part II, J. Geol., 23, 385–424,, 1915. a

Harel, M. A., Mudd, S. M., and Attal, M.: Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184–196,, 2016. a

Haviv, I., Enzel, Y., Whipple, K. X., Zilberman, E., Matmon, A., Stone, J., and Fifield, K. L.: Evolution of vertical knickpoints (waterfalls) with resistant caprock: Insights from numerical modeling, J. Geophys. Res.-Earth, 115, F03028,, 2010. a, b, c

Hayakawa, Y. S. and Oguchi, T.: DEM-based identification of fluvial knickzones and its application to Japanese mountain rivers, Geomorphology, 78, 90–106,, 2006. a, b, c, d, e

Heipke, C., Mayer, H., Wiedemann, C., and Jamet, O.: Evaluation of Automatic Road Extraction, Int. Arch. Photogramm., 32, 151–160, 1997. a

Howard, A. D.: Geomorphological systems; equilibrium and dynamics, Am. J. Sci., 263, 302–312,, 1965. a

Howard, A. D., Dietrich, W. E., and Seidl, M. A.: Modeling fluvial erosion on regional to continental scales, J. Geophys. Res.-Sol. Ea., 99, 13971–13986,, 1994. a

Hurvich, C. M. and Tsai, C. L.: Regression and time series model selection in small samples, Biometrika, 76, 297–307,, 1989. a

James, P. E.: The surface configuration of southeastern Brazil, Ann. Assoc. Am. Geogr., 23, 165–193, 1933. a

King, L. C.: A geomorfologia do Brasil oriental, Rev. Bras. Geogr., 2, 147–265, 1956. a

Kirby, E. and Whipple, K. X.: Expression of active tectonics in erosional landscapes, J. Struct. Geol., 44, 54–75,, 2012. a, b, c, d, e, f, g, h, i, j, k

Kirby, E., Whipple, K. X., Tang, W., and Chen, Z.: Distribution of active rock uplift along the eastern margin of the Tibetan Plateau: Inferences from bedrock channel longitudinal profiles, J. Geophys. Res.-Sol. Ea., 108, B4,, 2003. a

Knopf, E. B.: Correlation of residual erosion surfaces in the eastern appalachian highlands, B. Geol. Soc. Am., 35, 633–668,, 1924. a, b, c, d

Lambeck, K. and Chappell, J.: Sea level change through the last glacial cycle., Science, 292, 679–86,, 2001. a

Lapparent, A. D.: Leçons de géographie physique, Masson et cie, éditeurs, Paris, France, 1896. a, b

Lindsay, J. B.: Efficient hybrid breaching-filling sink removal methods for flow path enforcement in digital elevation models, Hydrol. Proc., 30, 846–857,,, 2016. a

Mackin, J. H.: Concept of the graded river, GSA Bulletin, 59, 463–512,[463:cotgr];2, 1948. a

Mandal, S. K., Lupker, M., Burg, J.-P., Valla, P. G., Haghipour, N., and Christl, M.: Spatial variability of 10Be-derived erosion rates across the southern Peninsular Indian escarpment: A key to landscape evolution across passive margins, Earth Planet. Sci. Lett., 425, 154–167,, 2015. a

Mather, A. E.: Adjustment of a drainage network to capture induced base-level change: an example from the Sorbas Basin, SE Spain, Geomorphology, 34, 271–289,, 2000. a

Moodie, A. J., Pazzaglia, F. J., and Berti, C.: Exogenic forcing and autogenic processes on continental divide location and mobility, Basin Research, 30, 344–369,, 2017. a

Morisawa, M. E.: Quantitative Geomorphology of Some Watersheds in the Appalachian Plateau, GSA Bulletin, 73, 1025–1046,[1025:QGOSWI]2.0.CO;2, 1962. a, b

Mouchené, M., van der Beek, P., Mouthereau, F., and Carcaillet, J.: Controls on Quaternary incision of the Northern Pyrenean foreland: Chronological and geomorphological constraints from the Lannemezan megafan, SW France, Geomorphology, 281, 78–93,, 2017. a

Mudd, S. M.: Detection of transience in eroding landscapes, Earth Surf. Proc. Land., 42, 24–41,, 2017. a

Mudd, S. M., Attal, M., Milodowski, D. T., Grieve, S. W., and Valters, D. A.: A statistical framework to quantify spatial variation in channel gradients using the integral method of channel profile analysis, J. Geophys. Res.-Earth, 119, 138–152,, 2014. a, b, c, d, e, f, g, h, i, j, k, l, m

Mudd, S. M., Clubb, F. J., Gailleton, B., and Hurst, M. D.: How concave are river channels?, Earth Surf. Dynam., 6, 505–523,, 2018a. a, b, c, d, e, f, g

Mudd, S. M., Clubb, F. J., Gailleton, B., Hurst, M. D., Milodowski, D. T., and Valters, D. A.: The LSDTopoTools Chi Mapping Package (Version 1.11), Zenodo,, 2018b. a

Mudd, S. M., Clubb, F. J., Gailleton, B., Grieve, S. W. D., Valters, D. A., and Hurst, M. D.: LSDTopoTools Documentation (Version v2.0), Zenodo,, 2019a. a

Mudd, S. M., Clubb, F. J., Gailleton, B., Valters, D. A., Hurst, M. D., and Grieve, S. W. D.: LSDMappingTools (Version v0.1), Zenodo,, 2019b. a

Muhs, D. R., Simmons, K. R., Schumann, R. R., Groves, L. T., DeVogel, S. B., Minor, S. A., and Laurel, D. A.: Coastal tectonics on the eastern margin of the Pacific Rim: Late Quaternary sea-level history and uplift rates, Channel Islands National Park, California, USA, Quaternary Sci. Rev., 105, 209–238,, 2014. a

Neely, A. B., Bookhagen, B., and Burbank, D. W.: An automated knickzone selection algorithm (KZ-Picker) to analyze transient landscapes: Calibration and validation, American Geophysical Union,, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t

Orlandini, S., Tarolli, P., Moretti, G., and Dalla Fontana, G.: On the prediction of channel heads in a complex alpine terrain using gridded elevation data, Water Resour. Res., 47, 2,, 2011. a

Ouimet, W. B., Whipple, K. X., and Granger, D. E.: Beyond threshold hillslopes: Channel adjustment to base-level fall in tectonically active mountain ranges, Geology, 37, 579–582,, 2009. a

Passalacqua, P., Do Trung, T., Foufoula-Georgiou, E., Sapiro, G., and Dietrich, W. E.: A geometric framework for channel network extraction from lidar: Nonlinear diffusion and geodesic paths, J. Geophys. Res.-Earth, 115, F01002,, 2010. a

Peifer Bezerra, D.: The pattern and style of landscape evolution in post-orogenic settings, available at: (last access: 15 December 2018), 2018. a, b, c

Pelletier, J. D.: A robust, two-parameter method for the extraction of drainage networks from high-resolution digital elevation models (DEMs): Evaluation using synthetic and real-world DEMs, Water Resour. Res., 49, 75–89,, 2013. a, b, c, d

Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Proc. Land., 38, 570–576,, 2013. a, b, c, d

Pinter, N., Lueddecke, S. B., Keller, E. A., and Simmons, K. R.: Late Quaternary slip on the Santa Cruz Island fault, California, Bull. Geol. Soc. Am., 110, 711–722,<0711:LQSOTS>2.3.CO;2, 1998. a

Pinter, N., Hardiman, M., Scott, A. C., and Anderson, R. S.: Discussion of “Fluvial system response to late Pleistocene-Holocene sea-level change on Santa Rosa Island, Channel Islands National Park, California” (Schumann et al., 2016. Geomorphology, 268: 322–340), Geomorphology, 301, 141–143,, 2018. a

Powell, J. W.: Exploration of the Colorado River of the West and its tributaries: Explored in 1869, 1870, 1871, and 1872, under the direction of the Secretary of the Smithsonian Institution, USGS Unnumbered Series, Government Printing Office, Washington, D.C., USA, 1875. a

Pritchard, D., Roberts, G. G., White, N. J., and Richardson, C. N.: Uplift histories from river profiles, Geophys. Res. Lett., 36, L24301,, 2009. a

Queiroz, G. L., Salamuni, E., and Nascimento, E. R.: Knickpoint finder: A software tool that improves neotectonic analysis, Comput. Geosci., 76, 80–87,, 2015. a

Royden, L. and Perron, J.: Solutions of the stream power equation and application to the evolution of river longitudinal profiles, J. Geophys. Res.-Earth, 118, 497–518,, 2013. a, b, c

Royden, L., Clark, M., and Whipple, K.: Evolution of river elevation profiles by bedrock incision: analytical solutions for transient river profiles related to changing uplift and precipitation rates, in: Eos, Transactions American Geophysical Union 81 (Fall Meeting Supplement Abstract T62F–09), vol. 81, 1–2, Fall Meeting Supplement, 15–19 December 2000, San Fransisco, California, 2000. a

Salgado, R., Augusto, A., Braucher, R., Chicarino, A., Coun, F., Fortes, A., and Chicarino, D.: Relief evolution of the Quadrilátero Ferrífero (Minas Gerais, Brazil) by means of 10Be cosmogenic nuclei, Z. Geomorphol., 52, 317–323, 2008. a, b

Scheingross, J. S. and Lamb, M. P.: Sediment transport through self-adjusting, bedrock-walled waterfall plunge pools, J. Geophys. Res.-Earth, 121, 939–963,, 2016. a

Scheingross, J. S., Lo, D. Y., and Lamb, M. P.: Self-formed waterfall plunge pools in homogeneous rock, Geophys. Res. Lett., 44, 200–208,, 2017. a

Scherler, D., Bookhagen, B., and Strecker, M. R.: Tectonic control on 10Be-derived erosion rates in the Garhwal Himalaya, India, J. Geophys. Res.-Earth, 119, 83–105,, 2014. a

Schumann, R. R., Pigati, J. S., and McGeehin, J. P.: Fluvial system response to late Pleistocene-Holocene sea-level change on Santa Rosa Island, Channel Islands National Park, California, Geomorphology, 268, 322–340,, 2016. a, b

Schwanghart, W. and Scherler, D.: Short Communication: TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences, Earth Surf. Dynam., 2, 1–7,, 2014. a, b, c, d, e, f, g, h, i

Schwanghart, W. and Scherler, D.: Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques, Earth Surf. Dynam., 5, 821–839,, 2017. a, b

Stock, J. D. and Montgomery, D. R.: Geologic constraints on bedrock river incision using the stream power law, J. Geophys. Res.-Sol. Ea., 104, 4983–4993,, 1999. a, b

Truong, C., Oudre, L., and Vayatis, N.: A review of change point detection methods, available at: (last access: 15 November 2018), 2018. a, b

Tucker, G. E. and Slingerland, R.: Predicting sediment flux from fold and thrust belts, Basin Res., 8, 329–349,, 1996. a

Varajão, G. V.: a Questão Da Correlação Das Superfícies De Erosão Do Quadrilátero Ferrífero , Minas Gerais, Rev. Bras. Geoci., 21, 138–145, 1991. a

Wang, L. and Liu, H.: An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling, Int. J. Geogr. Inf. Sci., 20, 193–213,, 2006.  a

Whipple, K. X.: Implications of sediment-flux-dependent river incision models for landscape evolution, J. Geophys. Res., 107, 2039,, 2002. a

Whipple, K. X., Kirby, E., and Brocklehurst, S. H.: Geomorphic limits to climate-induced increases in topographic relief, Nature, 401, 39–43,, 1999. a, b

Whipple, K. X., DiBiase, R. A., Ouimet, W. B., and Forte, A. M.: Preservation or piracy: Diagnosing low-relief, high-elevation surface formation mechanisms, Geology, 45, 91–94,, 2017a. a

Whipple, K. X., Forte, A. M., DiBiase, R. A., Gasparini, N. M., and Ouimet, W. B.: Timescales of landscape response to divide migration and drainage capture: Implications for the role of divide mobility in landscape evolution, J. Geophys. Res.-Earth, 122, 248–273,, 2017b. a, b

Whittaker, A. C. and Boulton, S. J.: Tectonic and climatic controls on knickpoint retreat rates and landscape response times, J. Geophys. Res.-Earth, 117, F2,, 2012. a

Wiener, N.: Extrapolation, interpolation, and smoothing of stationary time series: with engineering applications, Technology Press of the Massachusetts Institute of Technology, MIT, Massuchussetts, USA, 1949. a

Willett, S. D. and Brandon, M. T.: On steady state in mountain belts, Geology, 30, 175–178,<0175:OSSIMB>2.0.CO;2, 2002. a, b

Willett, S. D., McCoy, S. W., Taylor Perron, J., Goren, L., and Chen, C. Y.: Dynamic reorganization of River Basins, Science, 343, 1248765,, 2014. a, b

Wobus, C. W., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: procedurses, promise, and pitfalls, Geol. Soc. Am. Spec. P., 398, 55–74,, 2006. a, b, c, d

Zahra, T., Paudel, U., Hayakawa, Y. S., and Oguchi, T.: Knickzone Extraction Tool (KET) – A new ArcGIS toolset for automatic extraction of knickzones from a DEM based on multi-scale stream gradients, Open Geosci., 9, 73–88,, 2017. a

Short summary
The shape of landscapes is influenced by climate changes, faulting or the nature of the rocks under the surface. One of the most sensitive parts of the landscape to these changes is the river system that eventually adapts to such changes by adapting its slope, the most extreme example being a waterfall. We here present an algorithm that extracts changes in river slope over large areas from satellite data with the aim of investigating climatic, tectonic or geologic changes in the landscape.