How concave are river channels?

For over a century, geomorphologists have attempted to unravel information about landscape evolution, and processes that drive it, using river profiles. Many studies have combined new topographic datasets with theoretical models of channel incision to infer erosion rates, identify rock types with different resistance to erosion, and detect potential regions of tectonic activity. The most common metric used to analyse river profile geometry is channel steepness, or k_s. However, the calculation of channel steepness requires the normalisation of channel gradient by drainage area. This normalisation requires a power law exponent that is referred to as the channel concavity index. Despite the concavity index being crucial in determining channel steepness, it is challenging to constrain. In this contribution, we compare both slope–area methods for calculating the concavity index and methods based on integrating drainage area along the length of the channel, using so-called “chi” (χ) analysis. We present a new χ-based method which directly compares χ values of tributary nodes to those on the main stem; this method allows us to constrain the concavity index in transient landscapes without assuming a linear relationship between χ and elevation. Patterns of the concavity index have been linked to the ratio of the area and slope exponents of the stream power incision model (m∕n); we therefore construct simple numerical models obeying detachment-limited stream power and test the different methods against simulations with imposed m and n. We find that χ-based methods are better than slope–area methods at reproducing imposed m∕n ratios when our numerical landscapes are subject to either transient uplift or spatially varying uplift and fluvial erodibility. We also test our methods on several real landscapes, including sites with both lithological and structural heterogeneity, to provide examples of the methods' performance and limitations. These methods are made available in a new software package so that other workers can explore how the concavity index varies across diverse landscapes, with the aim to improve our understanding of the physics behind bedrock channel incision.

1 Introduction

Geomorphologists have been interested in understanding controls on the steepness of river channels for centuries. In his seminal “Report on the Henry Mountains”, Gilbert (1877) remarked that “We have already seen that erosion is favoured by declivity. Where the declivity is great the agents of erosion are powerful; where it is small they are weak; where there is no declivity they are powerless” (p. 114). Following Gilbert's pioneering observations of landscape form, many authors have attempted to quantify how topographic gradients (or declivities) relate to erosion rates. Landscape erosion rates are thought to respond to tectonic uplift (Hack, 1960). Therefore, extracting erosion rate proxies from topographic data provides novel opportunities for identifying regions of tectonic activity (Cyr et al., 2010; Lague and Davy, 2003; Seeber and Gornitz, 1983; Snyder et al., 2000; Wobus et al., 2006a), and may even be able to highlight potentially active faults (Kirby and Whipple, 2012). Analysing channel networks is particularly important for detecting the signature of external forcings from the shape of the topography, as fluvial networks set the boundary conditions for their adjacent hillslopes, therefore acting as the mechanism by which climatic and tectonic signals are transmitted to the rest of the landscape (Burbank et al., 1996; Hurst et al., 2013; Whipple, 2004; Whipple and Tucker, 1999).

Channels do not yield such information easily, however. Any observer of rivers or mountains will note that headwater channels tend to be steeper than channels downstream. Declining gradients along the length of the channel lead to river longitudinal profiles that tend to be concave up. Therefore, the gradient of a channel cannot be related to erosion rates in isolation; some normalising procedure must be performed. Over a century ago, Shaler (1899) postulated that as channels gain drainage area their slopes would decline, hindering their ability to erode. Authors such as Hack (1957), Morisawa (1962), and Flint (1974) expanded upon this idea in the early 20th century by quantifying the relationship between slope and drainage area, often used as a proxy for discharge. Flint (1974) found that channel gradient appeared to systematically decline downstream in a trend that could be described by a power law:

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

where θ is referred to as the concavity index since it describes how concave a profile is; the higher the value, the more rapidly a channel's gradient decreases downstream. Note that the term “concavity” is sometimes applied to river long profiles (Conway et al., 2015; Demoulin, 1998; Roy and Sinha, 2018), so most studies refer to θ, derived from slope–area data, as the “concavity index” (Kirby and Whipple, 2012). The term k_s is called the steepness index, as it sets the overall gradient of the channel. If we take the logarithm of both sides of Eq. (1), we find a linear relationship in log[S]–log[A] space with a slope of θ and an intercept (the value of log[S], where log[A] = 0) of log[k_s]:

$$\begin{array}{}\text{(2)}& \log \left[S\right]=-\mathit{\theta }\log \left[A\right]+\log \left[k_{\mathrm{s}}\right].\end{array}$$

A number of studies (DiBiase et al., 2010; Harel et al., 2016; Mandal et al., 2015; Ouimet et al., 2009; Scherler et al., 2014) have demonstrated that k_s is positively correlated with erosion rate, mirroring the predictions of Gilbert (1877) over a century earlier. Many authors have used channel steepness to examine fluvial response to climate, lithology, and tectonics (Cyr et al., 2010; DiBiase et al., 2010; Flint, 1974; Harkins et al., 2007; Kirby and Whipple, 2001, 2012; Kobor and Roering, 2004; Lague and Davy, 2003; Snyder et al., 2000; Tarboton et al., 1989; Vanacker et al., 2015; Wobus et al., 2006a).

The noise inherent in S–A analysis prompted Leigh Royden and colleagues to develop a method that compares the elevations of channel profiles, rather than slope (Royden et al., 2000). We can modify the Royden et al. (2000) approach to integrate Eq. (1), since S = dz∕dx where z is elevation and x is distance along the channel (Whipple et al., 2017), resulting in

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

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

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

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

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

Equation (5) shows that steepness of the channel (k_s) is related to the slope of the transformed channel in χ–elevation space. In both Eqs. (2) and (5), the numerical value of k_s depends on the value chosen for the concavity index, θ. In order to compare the steepness of channels in basins of different sizes, a reference concavity index is typically chosen (θ_ref), which is then used to extract a normalised channel steepness from the data (Wobus et al., 2006a):

$$\begin{array}{}\text{(6)}& k_{\mathrm{sn},i}=A_i^{{\mathit{\theta }}_{\mathrm{ref}}}{S}_i,\end{array}$$

where the subscript i is a data point. This data point often represents channel quantities that are smoothed; see discussion below.

1.1 Choosing a concavity index to extract channel steepness

The choice of the reference concavity index is important in determining the relative k_sn values amongst different sections in the channel network, which we illustrate in Fig. 1. This figure depicts hypothetical slope–area data, which appear to lie along a linear trend in slope–area space. Choosing a reference concavity index based on a regression through these data will result in the entire channel network having similar values of k_sn. Based on the data in Fig. 1, there is no evidence that the correct concavity index is anything other than the one represented by the linear fit through the data. However, these hypothetical data are in fact based on numerical simulations, presented in Sect. 3, in which we simulated a higher uplift rate in the core of the mountain range. The correct concavity index is therefore lower than that indicated by the log[S]–log[A] data, and instead the data show a strong spatial trend in channel steepness (interpretation 2 in Fig. 1). The simplest interpretation based on log[S]–log[A] data alone would have been entirely incorrect. This situation is analogous to the one described by Kirby and Whipple (2001), where downstream reductions in uplift rates in the Siwalik Hills of India and Nepal resulted in elevated apparent concavities. Conversely, if a single reference concavity index is chosen in an area with changing concavity indices, then spurious patterns in k_sn may arise (Gasparini and Whipple, 2014). These examples highlight that selecting the correct concavity index is crucial if we are to correctly interpret channel steepness data.

Figure 1 Sketch illustrating the effect of choosing different reference concavities. The data can be well fitted with a single regression, suggesting that all parts of the channel network have similar values of k_sn (interpretation 1). However, if a lower θ_ref is chosen, the k_sn values will be systematically higher for channels at lower drainage area (interpretation 2). This sketch is based on data from a numerical simulation where the latter situation has been imposed via higher uplift rates in the core of the mountain range, showing the potential for incorrect concavities to be extracted from slope–area data alone.

Extracting a reliable reference concavity indices from slope–area data on real landscapes is challenging: topographic data can be noisy, leading to a wide range of channel gradients for small changes in drainage area. The branching nature of river networks also results in large discontinuities in drainage areas where tributaries meet, resulting in significant data gaps in S–A space (Fig. 2). Wobus et al. (2006a) made recommendations for preprocessing of slope–area data that are still used in many studies. First, the digital elevation model (DEM) is smoothed (although with improved DEMs this is now rare), then topographic gradient is measured over either a fixed reach length or a fixed drop in elevation (Wobus et al., 2006a, recommends the latter), and then the data are averaged in logarithmically spaced bins. More recently, authors have proposed alternative channel smoothing strategies (Aiken and Brierley, 2013; Schwanghart and Scherler, 2017); all these methods use some form of smoothing and averaging.

Figure 2 A typical slope–area plot. This example is from a basin near Xi'an, China, with an outlet at approximately 34^∘26^′23.9${}^{\prime \prime }$ N, 109^∘23^′13.4${}^{\prime \prime }$ E. The data are taken from only the trunk channel. The slope–area data typically contain gaps due to tributary junctions, as well as wide ranges in slope for the reaches between junctions due to topographic noise inherent in deriving slope values. The result is a high degree of scatter in the data. These data are produced by averaging slope values over a fixed vertical interval of 20 m.

In order to circumvent these problems with S–A analysis, many authors have since used the integral approach (Eq. 5) to analyse channel concavity indices. This method also aims to normalise river profiles for their drainage area, but rather than comparing slope to area, it integrates area along channel length (Perron and Royden, 2013; Royden et al., 2000). Perron and Royden (2013) showed that the concavity index could be extracted from a channel by selecting the value of θ_ref that results in the most linear channel profile in χ–elevation space. However, if there are sections of the channel with different k_sn values, this will hinder our ability to extract the θ_ref value, as it is not appropriate to fit a single line throughout the entire profile (Mudd et al., 2014). Mudd et al. (2014) introduced a method to statistically determine the most likely concavity index by computing the best-fit series of linear segments using an algorithm that balanced fit of the data against overparameterisation using the Akaike information criterion (Akaike, 1974). This method, however, requires a number of input parameters and also performs computationally expensive segmentation on the χ–elevation data prior to calculating the concavity index (Mudd et al., 2014).

Perron and Royden (2013) suggested a second independent means to calculate the concavity index, which does not assume linearity of the profiles in χ–elevation space and may therefore be used in transient landscapes. This method is instead based on searching for collinearity of tributaries with the main stem channel (Mudd et al., 2014; Perron and Royden, 2013), and has since been used as a basis for other techniques that aim to minimise some quantitative description of scatter between tributaries and the trunk channel Goren et al. (2014); Hergarten et al. (2016).

Although the collinearity test does not assume any linearity of profiles in χ–elevation space, it does rely on the assumption that points in the channel network with the same value of χ will have the same elevation. Perron and Royden (2013) noted that this will hold true if transient erosion signals propagate vertically through the network at a constant rate, which has been predicted by some theoretical models of fluvial incision (Wobus et al., 2006b). Royden and Perron (2013) went on to demonstrate that changes in erosion rates in channel networks would lead to distinct segments that migrate upstream, which they termed “slope patches”, where the local k_s reflects local erosion rate. However, here, we wish to avoid basing our formulation on any theoretical models of fluvial incision, as this introduces assumptions regarding bedrock erosion processes. Y. Wang et al. (2017) highlighted that slope–area analysis could be used to complement integral analysis where concavity indices might vary spatially, since slope–area is agnostic with regard to these incision processes, but also found that integral-based analyses have lower uncertainties than slope–area analysis. Therefore, here, we use simple geometric relationships of knickpoint propagation to relate the collinearity test to channel concavity indices without relying on theoretical models of fluvial incision, as set out in Sect. 1.2.

1.2 Connecting the concavity index to collinearity

Playfair (1802) noted that tributary valleys tended to join the principal valley at a common elevation, suggesting that, at their outlets, the tributary streams must lower at the same rate as the principal streams into which they drain. Therefore, any change in incision rate on the main stem channel will be transmitted to the upstream tributaries. Using simple geometric relationships, Niemann et al. (2001) showed that a knickpoint should migrate upstream with a horizontal celerity (Ce_h, in length per time) of

$$\begin{array}{}\text{(7)}& {\mathrm{Ce}}_{\mathrm{h}}={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{2}}-S_{\mathrm{1}}}}\mathrm{\Delta }E,\end{array}$$

where S₁ is the channel slope prior to disturbance, S₂ is the channel slope after disturbance (e.g. due to a change in incision rate E), and ΔE is the difference between the incision rate before and after disturbance (which can be equated to uplift rates U₁ and U₂ in units of length per time, $\mathrm{\Delta }E=U_{\mathrm{2}}-U_{\mathrm{1}}$). Wobus et al. (2006b) simply inserted Eq. (1) into Eq. (7) so that the horizontal celerity is simply a function of drainage area, assuming that the concavity index is independent of rock uplift rate:

$$\begin{array}{}\text{(8)}& {\mathrm{Ce}}_{\mathrm{h}}={\displaystyle \frac{U_{\mathrm{2}}-U_{\mathrm{1}}}{k_{\mathrm{s}\mathrm{2}}-k_{\mathrm{s}\mathrm{1}}}}{A}^{\mathit{\theta }}.\end{array}$$

Noting that vertical celerity is simply the horizontal celerity multiplied by the local slope after disturbance S₂, Wobus et al. (2006b) showed that the vertical celerity (Ce_v) was not a function of drainage area:

$$\begin{array}{}\text{(9)}& {\mathrm{Ce}}_{\mathrm{v}}={\displaystyle \frac{U_{\mathrm{2}}-U_{\mathrm{1}}}{k_{\mathrm{s}\mathrm{2}}-k_{\mathrm{s}\mathrm{1}}}}{k}_{\mathrm{s}\mathrm{2}}.\end{array}$$

Thus, if we assume spatially homogeneous uplift and constant erodibility (i.e. channels with the same slope and drainage area erode at the same rate), then the vertical celerity propagating up the principal stream and all tributaries will be a constant. Equation (9) is derived from purely geometric relationships, suggesting that collinearity can be used to estimate the concavity index without assuming any theoretical models of fluvial incision under the assumption that the incision process will result in local slope–area relationships reflecting Eq. (1).

2 Calculation of concavity indices using collinearity

In this study, we revisit commonly used methods for estimating the concavity index using both slope–area analysis and collinearity methods based on integral analysis. Our objective is to determine the strengths and weaknesses of established methods alongside several new methods developed for this study, as well as to quantify the uncertainties in the concavity index. We present these methods in an open-source software package that can be used to constrain concavity indices across multiple landscapes. This information may give insight into the physical processes responsible for channel incision into bedrock, which are as yet poorly understood.

2.1 Slope–area analysis

For slope–area analysis, in this paper, we forgo initial smoothing of the DEM and use a fixed elevation drop along a D8 drainage pathway implemented using the network extraction algorithm of Braun and Willett (2013). We calculate the best-fit concavity indices using two different methods: (i) concavity indices extracted from all slope–area (S–A) data (i.e. no logarithmic bins, every tributary) and (ii) concavity indices of contiguous channel profile segments with consistent S–A scaling within the log-binned S–A data of the trunk stream, calculated using the statistical segmentation algorithm described in Mudd et al. (2014). We report the different extracted concavity indices and their uncertainties in the results below.

Figure 3 Sketch illustrating the methodology of the χ method using all profile data. In panel (a), the χ profiles of both the trunk channel and a tributary are shown. We take the χ coordinate of the nodes on the tributary channel and then project them onto a linear fit of the trunk channel to determine the residuals between tributary and trunk channel. We do this for all nodes and for all concavity indices. For each concavity index, the residuals are then used to calculate a maximum likelihood estimator (MLE), which varies as a function of the concavity index (panel b). The highest value of MLE is used to select the likely θ.

2.2 Methods for calculating collinearity using integral analysis

Here, we present two new methods of identifying collinear tributaries in χ–elevation space in order to constrain the best-fit concavity indices from fluvial profiles. Rather than fitting segments to the profiles, which is computationally expensive, we directly compare all the elevation data of the tributaries in each drainage basin to the main stem. This is not completely straightforward, however: because the χ coordinate integrates area and channel distance, it is very unlikely that a pixel on a tributary channel shares a χ coordinate with any pixel on the main stem. Instead, for every tributary pixel, we compare the tributary elevation with an elevation on the main stem at the same χ computed with a linear fit between the two pixels with the nearest χ coordinates (Fig. 3). We then calculate a maximum likelihood estimator (MLE) for each tributary. The MLE is calculated with

$$\begin{array}{}\text{(10)}& \mathrm{MLE}=\prod _{i=\mathrm{1}}^N\exp [-{\displaystyle \frac{r_i^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}],\end{array}$$

where N is the number of nodes in the tributary, r_i is the calculated residual between the elevation of tributary node i and the linear regression of elevation on the main stem, and σ is a scaling factor. If r_i is 0 for all nodes, then MLE = 1 (i.e. MLE varies between 0 and 1, with 1 being the maximum possible likelihood).

For a given drainage basin, we can multiply the MLE for each tributary to get the total MLE for the basin, and we can do this for a range of concavity indices to calculate the most likely value of θ. Because Eq. (10) is a product of negative exponentials, the value of the MLE will decrease as N increases, and in large datasets this results in MLE values below the smallest number that can be computed, meaning that in large datasets MLE values can often be reported as zero. To counter this effect, we increase σ until all tributaries have non-zero MLE values. As σ is simply a scaling factor, this does not affect which concavity index is calculated as the most likely value once all tributaries have non-zero MLEs (see the Supplement).

There are two disadvantages to using Eq. (10) on all points in the channel network. Firstly, because the MLE is calculated as a product of exponential functions, each data point will reduce the MLE, and so tributaries will influence MLE in proportion to their length. Secondly, because we use all data, we cannot estimate uncertainty when computing the most likely concavity index. Therefore, we apply a second method to the χ–elevation data that mitigates these two shortcomings by bootstrap sampling the data. This method evaluates a fixed number of discrete points on each tributary, but the points are selected randomly and this random selection is done iteratively, building up a population of MLE values for each concavity index.

For each iteration of the bootstrap method, we create a template of points in χ space, measured from the confluence of each tributary from the trunk channel (Fig. 4). We start by selecting a maximum value of χ upstream of the tributary junction, and then separate this space into N_BS nodes. We create evenly spaced bins between the maximum value of χ in the template, and then in each iteration randomly select one point in each bin. Using this template on each tributary, we calculate the residuals between the tributary and the trunk channel using Eq. (10). If, for a given tributary, a point in the template is located beyond the end of the tributary, then the point is excluded from the calculation of MLE. Figure 4 provides a schematic visualisation of this method.

Figure 4 Sketch showing how we compute residuals for our χ bootstrap method of determining the maximum likelihood estimator (MLE) of θ, and then use the uncertainty in MLE values to compute the uncertainty in θ.

We repeat these calculations over many iterations, and for each concavity index we compute the median MLE, the minimum and maximum MLEs, and the first and third quartile MLEs. We approximate the uncertainty range by first taking the most likely concavity index (having highest median MLE value amongst all concavity indices tested). We then find the span of concavity indices whose third quartile MLE values exceed the first quartile MLE value of the most likely concavity indices (Fig. 4).

One complication of using collinearity to calculate the most likely concavity index is that occasionally one may find a hanging tributary (Crosby et al., 2007; Wobus et al., 2006b), which could occur for a variety of reasons, such as the presence of geologic structures or lithologic variability. A hanging tributary can skew the overall MLE values in a basin, so in each basin we test the MLE and RMSE values in each tributary for outliers and iteratively remove these outlying tributaries, testing for the most likely concavity index on each iteration. However, we find that eliminating outlying tributaries has a minimal effect on the calculated concavity index. The other primary complication is that one must assume a concavity index prior to performing the χ transformation (Eq. 4), and thus slope–area analysis may be more suited to detecting changes in concavity within basins (Y. Wang et al., 2017). We suggest here an alternative approach of calculating concavity using χ methods in many small basins to look for any systematic changes. Before we can perform such analyses, however, we must constrain our confidence in estimates of the concavity index.

We also implement a disorder statistic (Goren et al., 2014; Hergarten et al., 2016; Shelef et al., 2018) that aims to quantify differences in the χ–elevation patterns between tributaries and the trunk channel. Here, we follow the method of Hergarten et al. (2016). The disorder statistic is calculated by first taking the χ–elevation pairs of every point in the channel network, ordered by increasing elevation. We calculate the sum:

$$\begin{array}{}\text{(11)}& S=\sum _{i=\mathrm{1}}^N|{\mathit{\chi }}_{\mathrm{s},i+\mathrm{1}}-{\mathit{\chi }}_{\mathrm{s},i}|,\end{array}$$

where the subscript s, i represents the ith χ coordinate that has been sorted by its elevation. The sum, S, is minimal if elevation and χ are related monotonically. However, it scales with the absolute values of χ, which are sensitive to the concavity index (see Eq. 4), so following Hergarten et al. (2016), we scale the disorder metric, D, by the maximum value of χ in the tributary network (χ_max):

$$\begin{array}{}\text{(12)}& D={\displaystyle \frac{\mathrm{1}}{{\mathit{\chi }}_{max}}}\left(\sum _{i=\mathrm{1}}^N\right|{\mathit{\chi }}_{\mathrm{s},i+\mathrm{1}}-{\mathit{\chi }}_{\mathrm{s},i}|-{\mathit{\chi }}_{max}).\end{array}$$

The disorder metric relies on the use of all the data in a tributary network, meaning that only one value of D can be calculated for each basin. Therefore, we cannot estimate the uncertainty in concavity index using this statistic alone. Furthermore, the random sampling approach we take with the previous χ methods is not appropriate, as skipping nodes in the χ–elevation sequence will lead to large values of S and substantially increase the disorder metric. We therefore employ a bootstrap approach based on the analysis of entire tributaries within each basin. First, we find every combination of three tributaries plus the trunk stream in the basin. For each combination, we then iterate through a range of concavity indices and calculate the disorder metric. This allows us to find the concavity index that minimises the disorder metric for each combination, resulting in a population of best-fit concavities, from which we calculate the median and interquartile range.

3 Testing on numerical landscapes

In real landscapes, we can only approximate the concavity index based on topography. However, we can create simulations where we fix a known concavity index and see if our methods reproduce this value. To do this, we rely on simple simulations driven by the general form of the stream power incision model, first proposed by Howard and Kerby (1983):

$$\begin{array}{}\text{(13)}& E=K{A}^m{S}^n,\end{array}$$

where E is the long-term fluvial incision rate, A is the upstream drainage area, S is the channel gradient, K is the erodibility coefficient, which is a measure of the efficiency of the incision process, and m and n are exponents. A number of variations of this equation are possible: some authors have proposed, for example, modifications that involve erosion thresholds (Tucker and Bras, 2000) or modulation by sediment fluxes (Sklar and Dietrich, 1998). However, Gasparini and Brandon (2011) showed that many of the modified versions of Eq. (13) could be captured simply by modifying the exponents m and n.

We have chosen this model because it can be related to the concavity index and therefore can be used to test the different methods under idealised conditions. We can relate the stream power incision model to Eq. (1) by rearranging Eq. (13) to solve for channel slope and relating it to local erosion rate, E:

$$\begin{array}{}\text{(14)}& S=({\displaystyle \frac{E}{K}}{)}^{\mathrm{1}/n}{A}^{-m/n}.\end{array}$$

Comparing Eqs. (1) and (14) reveals that the ratio between area and slope exponents in the stream power incision model, m∕n, is therefore equivalent to θ from Eq. (1). The channel steepness index, k_s, is related to erosion rate by

$$\begin{array}{}\text{(15)}& k_{\mathrm{s}}=({\displaystyle \frac{E}{K}}{)}^{\mathrm{1}/n}.\end{array}$$

The stream power incision model also makes predictions about how tectonic uplift can be translated into local erosion rates (Whipple and Tucker, 1999), and the predicted relationship between the channel steepness index and uplift has been exploited by a number of studies to identify areas of tectonic activity (Kirby and Whipple, 2012; Kirby et al., 2003; Wobus et al., 2006a). Furthermore, many workers have used the framework of the stream power incision model to extract uplift histories (Fox et al., 2014; Goren et al., 2014; Pritchard et al., 2009; Roberts and White, 2010). However, the ability of these studies to extract information from channel profiles is dependent on the both the m∕n ratio, equivalent to θ, and the slope exponent, n, which are key unknowns within these theoretical models of fluvial incision. The m∕n ratio is frequently assumed to be equal to 0.5, with n assumed to be unity, despite recent compilations of data from multiple landscapes showing that this may not be the case (Clubb et al., 2016; Harel et al., 2016; Lague, 2014), and numerical modelling studies showing that m∕n = 0.5 leads to unrealistic relief structures (Kwang and Parker, 2017).

To test the relative efficacy of our methods for calculating concavity indices, we first run each method on a series of numerically simulated landscapes in which the m∕n ratio is prescribed. We employ a simple numerical model, following Mudd (2016), where channel incision occurs based on Eq. (13). For computational efficiency, we do not include any other processes (e.g. hillslope diffusion) within our model. The elevation of the model surface therefore evolves over time according to

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial z}{\partial t}}=U-K{A}^m{S}^n,\end{array}$$

where U is the uplift rate. Fluvial incision is solved using the algorithm of Braun and Willett (2013), where the drainage area is computed using the D8 flow direction algorithm to improve speed of computation and the topographic gradient is calculated in the direction of steepest descent. In our model, we perform a direct numerical solution of Eq. (16) where n=1 and use Newton–Raphson iteration where n≠1. These simulations are performed using the MuddPILE numerical model (Mudd et al., 2017), first used by Mudd (2016). We set the north and south boundaries of the model domain to fixed elevations, whereas the east and west boundaries are periodic. Our model domain is 30 km in the X direction and 15 km in the Y direction, with a grid resolution of 30 m. This allows us to test the methods of estimating m∕n on several drainage basins in each model domain, and at a resolution comparable to that of globally available DEMs.

3.1 Transient landscapes

In order to test the methods' ability to identify the correct m∕n value, we ran a series of numerical experiments with varying m∕n ratios: $m/n=\mathrm{0.5}$, $m/n=\mathrm{0.35}$, and $m/n=\mathrm{0.65}$. For each ratio, we also performed simulations with varying values of n, as the n exponent has been shown to impact the celerity of transient knickpoint propagation through the channel network (Royden and Perron, 2013). Crucially, Royden and Perron (2013) showed that when n is not unity, upstream propagating knickpoints will erase information about past base level changes encoded in the channel profiles. This may cloud selection of the correct m∕n ratio, but Lague (2014) and Harel et al. (2016) have suggested many, if not most, natural landscapes have evidence for an n exponent that is not unity. Therefore, we ran simulations with n=1, n=2, n=1.5, and n=0.66 for each m∕n ratio, varying m accordingly (see the Supplement for details of each model run).

Figure 5 Shaded relief plots of the model runs with temporally varying uplift, with drainage basins plotted by the best-fit θ predicted from the χ bootstrap analysis (column a), and slope–area analysis (column b). Each row represents a model run with a different m∕n ratio. The basins are coloured by the predicted θ, where darker colours indicate a higher concavity index. The extracted channel network for each basin is shown in blue.

We initialised the model runs using a low relief surface that is created using the diamond–square algorithm (Fournier et al., 1982). We found this approach resulted in drainage networks that contained more topological complexity than those initiated from simple sloping or parabolic surfaces. Our aim was to test the ability of each method to extract the correct m∕n ratio without assuming that the landscapes were in steady state; therefore, each simulation was forced with varying uplift through time to ensure that the channel networks were transient.

Each model was run with a baseline uplift rate of 0.5 mm yr⁻¹, which was increased by a factor of 4 for a period of 15 000 years, then decreased back to the baseline for another 15 000 years. For the runs with n=2, the cycles were set to 10 000 years, which was necessary to preserve evidence of transience, as knickpoints propagate more rapidly through the channel network as n increases. Relief is very sensitive to model parameters, and we found in numerical experiments that basin geometry was sensitive to relief, mirroring the results of Perron et al. (2008). We wanted modelled landscapes to have comparable relief and similar basin geometry across our simulations to ensure similar landscape configurations for different values of m, n, and m∕n. We therefore calculated the χ coordinate and solved Eq. (5) to find the K value for each modelled landscape that produced a relief of 200 m at the location with the greatest χ value given an uplift rate of 0.5 mm yr⁻¹.

We analysed these model runs using each of the methods of estimating the best-fit m∕n outlined in Sect. 2.2. We extracted a channel network from each model domain using a contributing area threshold of 9×10⁵ m². We performed a sensitivity analysis of the methods to this contributing area threshold (see the Supplement) and found that the estimated best-fit m∕n ratios were insensitive to the value of the threshold.

Drainage basins were selected by setting a minimum and maximum basin area, 9×10⁶ and 4.5×10⁷ m², respectively; these values were chosen so extracted basins represented a good balance between the number of extracted basins and the number of tributaries in each basin. Nested basins were removed, as were basins that bordered the edge of the model domain. We exclude basins on the domain boundaries as the calculation of the χ coordinate for the integral profile analysis is dependent on drainage area, which may not be realistic at the edge of the domain. Elimination of basins on the edge of the DEM is essential for real landscapes, as a basin beheaded by raster clipping will have incorrect χ values, and we wanted to ensure both simulations and analyses on real basins used the same extraction algorithms. For each basin, we identified the best-fit concavity index predicted in five ways (as described in the methods section): (i) regression of all χ–elevation data; (ii) use of χ–elevation data processed by our method of sampling points with the bootstrap method; (iii) minimisation of the disorder metric from χ–elevation data using a similar technique to Hergarten et al. (2016); (iv) regression of all slope–area data; and (v) regressions through slope–area data for individual segments of the main stem. For all but the final method, the analyses use all tributaries in the basins.

Figure 6 Plots showing the predicted best-fit θ for each basin and each method for $m/n=\mathrm{0.5}$, where n = 1, n = 2, n = 1.5, and n = 0.66. The χ methods are shown in shades of red and the slope–area methods are shown in shades of blue.

Figure 5 shows the spatial distribution of the predicted m∕n ratio for a series of basins from these cyclic model runs, where m∕n = 0.35, 0.5, and 0.65, and n=1. We also plot the m∕n ratio predicted for each basin from all methods with varying values of n, an example of which is shown in Fig. 6. Our modelling results show that for each value of m∕n ratio tested, the method using all χ data identifies the correct ratio for every basin in the model domain. The bootstrap approach provides an estimate of the error on the best-fit m∕n ratio for each basin: Fig. 6 shows that there is no error on the predicted m∕n ratio, meaning that an identical m∕n ratio is predicted with each iteration of the bootstrap approach. The slope–area methods, in contrast, show more variation in the predicted m∕n ratio for each value of m∕n and n tested (Figs. 5 and 6). Furthermore, the segmented slope–area data show a higher uncertainty in the predicted m∕n ratio compared to the other methods. The results of the model runs for all values of m∕n and n are presented in the Supplement.

3.2 Spatially heterogeneous landscapes

Alongside these temporally transient scenarios, we also wished to test the ability of each method to identify the correct m∕n ratio in spatially heterogeneous landscapes, simulating the majority of real sites where lithology, climate, or uplift are generally non-uniform. Therefore, we performed additional runs where m∕n = 0.5, n=1, but U and K varied in space. We generated the model domains using the same diamond–square initial condition as the spatially homogeneous runs. For the run with spatially varying K, we calculate the steady-state value of K required to produce a surface with a relief of 400 m and an uplift rate of 1 mm yr⁻¹ using the same method as for the previous runs. From this baseline value of K, we calculated a maximum K value which is 5 times that of the baseline. We then created 10 “patches” within the initial model domain where K was assigned randomly between the baseline and the maximum. These are rectangular in shape with K values that taper to the baseline K over 10 pixels. We acknowledge this pattern is not very realistic but emphasise that the aim is not to recreate real landscapes but rather to confuse the algorithms for quantifying the concavity index. This allows us to test if they can still detect modelled concavity indices even if we violate some of the assumptions implicit in our theoretical framework.

Figure 7 Results of the model runs with spatially varying erodibility (K, a, c, e) and uplift (U, b, d, f). The rectangular patches of low relief are areas of high erodibility in (a), (c), and (e). Panels (a)–(d) show the spatial pattern of predicted θ from the χ bootstrap analysis and the slope–area analysis, where the basins are coloured by θ (darker colours indicate higher concavity index). Panels (e) and (f) show density plots of the distribution of θ for each method, where the dashed line marks the correct θ = 0.5.

For the spatially varying uplift run, we varied uplift in the N–S direction by modelling it as a half sine wave:

$$\begin{array}{}\text{(17)}& U=U_{\mathrm{A}}\sin \left(\right(\mathit{\pi }y)/L)+U_{min},\end{array}$$

where y is the northing coordinate and L is the total length of the model domain in the y direction, U_A is an uplift amplitude, set to 0.2 mm yr⁻¹, and U_min is a minimum uplift, expressed at the north and south boundaries, of 0.2 mm yr⁻¹. Both scenarios, with spatially varying erodibility and uplift, were run to approximately steady state; the maximum elevation change between 15 000-year printing intervals was less than a millimetre.

Inherent in each collinearity-based method of quantifying the most likely m∕n ratio is the assumption that U and K do not vary in space (Perron and Royden, 2013); our spatially heterogeneous experiments therefore violate basic assumptions of the integral method. These conditions, however, are likely true in virtually all natural landscapes. Therefore, our aim here was to test if we could recover m∕n ratios from numerical landscapes that are more similar to real landscapes than those with spatially homogeneous U and K.

Figure 7 shows the distribution of predicted m∕n ratios for the runs with spatially varying K and U from both the integral bootstrap approach and the slope–area method. In comparison to our model runs where K and U were uniform, each method performs worse at identifying the correct m∕n ratio of 0.5. However, in both model runs, the integral methods identified the correct ratio in a higher proportion of the drainage basins than the slope–area methods. Furthermore, the distribution of m∕n predicted by the integral methods reaches a peak at the correct m∕n ratios of 0.5, suggesting that even in spatially heterogeneous landscapes the methods can still be applied. Our run with the random distribution of erodibility patches shows that the correct calculation of the m∕n ratio is highly dependent on the spatial continuity of K; in basins contained within a single patch (e.g. basins 4, 5, and 6), the integral profile method correctly identified the m∕n ratios. Figure 8 shows example χ–elevation plots at varying m∕n ratios for basin 2, which encompasses several patches with varying K values. Within this basin, tributaries that drain a patch with the same K value are still collinear in χ–elevation space. Based on these results, we suggest that, in real landscapes, monolithologic catchments should be analysed wherever possible in order to select an appropriate concavity index.

Figure 8 Example χ–elevation plots for the model run with spatially varying erodibility, where points are coloured by K. The m∕n increases in each plot from 0.2 to 0.9. Tributaries with the same K value are collinear in χ–elevation space.

4 Constraining concavity indices in real landscapes

Our numerical modelling results suggest that the integral profile analysis is most successful in identifying the correct concavity index out of the entire range of m∕n and n values tested. However, these modelling scenarios cannot capture the range of complex tectonic, lithologic, and climatic influences present in nature. Therefore, we repeat our analyses on a range of different landscapes with varying climates, relief structures, and lithologies, to provide some examples of the variation in the concavity index predicted using each method. For each field site, topographic data were obtained from OpenTopography, using the seamless DEM generated from NASA's Shuttle Radar Topography Mission (SRTM) at a grid resolution of 30 m. The Supplement contains metadata for each site so readers can extract the same topographic data used here.

Figure 9 Exploration of the most likely concavity indices in the Loess Plateau, China, Universal Transverse Mercator (UTM) zone 49 N. Basins with the most likely concavity index determined by the disorder method are displayed in panel (a); the basin number is followed by the most likely concavity index in the basin labels. The probability density of best-fit concavity index for all the basins (i.e. not the uncertainty within individual basins but rather the probability distribution of the best-fit concavity indices of all the basins) is displayed in panel (b). In basin 1, the most likely concavity index determined by the bootstrap and disorder methods is 0.45 and the χ–elevation plot for this concavity index is shown in panel (c).

4.1 An example of a relatively uniform landscape: Loess Plateau, China

In order to demonstrate the ability of the methods to constrain concavity indices in a relatively homogeneous landscape, we first analyse the Loess Plateau in northern China. The channels of the Loess Plateau are incising into wind-blown sediments that drape an extensive area of over 400 000 km² (Zhang, 1980) and can exceed 300 m thickness (Fu et al., 2017). The plateau is underlain by the Ordos Block, a succession of non-marine Mesozoic sediments which has undergone stable uplift since the Miocene (B. Wang et al., 2017; Yueqiao et al., 2003). Although there have been both recent (Wang et al., 2016) and historic (Wang et al., 2006) changes in sediment discharge from the plateau, the friable substrate means that channel networks and channel profiles might be expected to adjust quickly to perturbations in erosion rate. Indeed, Willett et al. (2014) suggested, based on differences in the χ coordinate across drainage divides, that the channel networks in large portions of the plateau are geomorphically stable. The stable tectonic setting and homogeneous, weak substrate of the Loess Plateau make an ideal natural laboratory for testing our methods on relatively homogeneous channel profiles.

We ran each of the methods on an area of the Loess Plateau approximately 11 000 km² in size near Yan'an, in the Chinese Shaanxi province (Fig. 9a). We find relatively good agreement between both the χ and slope–area methods of estimating the most likely concavity index. Figure 9b shows the probability distribution of concavities determined from the population of the most likely concavities from each basin (i.e. it does not include underlying uncertainty in each basin), but the peaks of these curves lie at a θ≈0.45 using both the disorder method and the bootstrap method, 0.4 using all slope–area data, and at approximately 0.5 using the all χ data method. This level of agreement gives the worker some confidence that channel steepness analyses in this area of the Loess Plateau using reference concavities between 0.4 and 0.5 should give an accurate representation of the relative steepness of the channels.

As well as determining the best-fit concavity index for the landscape as a whole, we can also examine the channel networks in individual basins: Fig. 9c shows the χ–elevation profiles for an example basin. In this basin, the tributaries are well aligned with the trunk channel at the most likely θ of 0.45, both using all the χ data and with the bootstrap approach. In our explorations of different landscapes, the Loess Plateau is the landscape that most resembles the idealised landscapes that we find in our model simulations. The Loess Plateau is notable for the homogeneity of its substrate over a large area; most locations on Earth are not as homogeneous.

4.2 An example of lithologic variability: Waldport, Oregon, USA

Many studies analysing the steepness of channel profiles are focused in areas where external factors, such as lithology or tectonics, are not uniform. Here, we select an example of a landscape with two dominant lithologic types in a location along the Oregon coast near the town of Waldport, Oregon (Fig. 10). The Oregon Coast Range is dominated by the Tyee Formation, made up primarily of turbidites deposited during the Eocene (Heller et al., 1987). In addition to these sedimentary units, our selected landscape also contains the Yachats Basalt, which erupted mostly as sub-areal flows between 3 and 9 m in thickness during the late Eocene (Davis et al., 1995). Erosion rates inferred from ¹⁰Be concentrations in stream sediments are between 0.11 to 0.14 mm yr⁻¹ (Bierman et al., 2001; Heimsath et al., 2001), similar to rock uplift rates of 0.05–0.35 mm yr⁻¹ inferred from marine terraces (Kelsey et al., 1994). Short-term erosion rates derived from stream sediments fall into the range of 0.07 to 0.18 mm yr⁻¹ (Wheatcroft and Sommerfield, 2005), leading a number of authors to suggest that the Coast Range is in topographic steady state, where uplift is balanced by erosion (Reneau and Dietrich, 1991). Thus, our site contains a clear lithologic contrast but has been selected to minimise spatial variations in uplift or erosion rates.

Figure 10 Exploration of the most likely concavity index near Waldport, Oregon, UTM zone 10 N. Basin numbers and the underlying lithology is displayed in panel (a). The most likely concavity index determined by the bootstrap method as a function of the percent of each basin in the different lithologies is shown in panel (b). Panel (c) shows the χ–elevation plot for a basin that has two bedrock types; the channel pixels are coloured by lithology. The plot uses the typical concavity index for basalt (0.7).

We find that whereas basins developed on basalt have a relatively uniform concavity index of approximately 0.7, the most likely concavity indices in the sandstone show considerably more scatter (Fig. 10b), with a lower average θ. We present these data as an example of spatially varying concavity indices as a function of lithology. This is consistent with results of VanLaningham et al. (2006), who found high concavities in volcanic rocks around Waldport but lower elsewhere, and found high values of the concavity index in sedimentary rock but with a higher degree of scatter along the Oregon Coast Range.

Figure 11 Basins analysed near the Gulf of Evia, Greece, UTM zone 34 N, that interact with active normal faults previously studied by Whittaker and Walker (2015).

Whipple and Tucker (1999) suggested that the concavity index is controlled primarily by discharge–drainage area and channel width–drainage area relationships, which may be influenced by lithology, but other authors have found systematic variations in concavity indices with lithology (Duvall et al., 2004; Lima and Flores, 2017; VanLaningham et al., 2006). Lima and Flores (2017) suggested that the thickness of basalt flows could influence concavity indices, with different knickpoint propagation mechanisms in massive versus thinly bedded flows. Duvall et al. (2004) suggested that having hard rocks in headwaters and weak below might influence concavity indices, which may be tested by comparing concavity indices in both monolithologic basins and basins with mixed lithology. The χ profiles in basin 17 (Fig. 10c) are notable because this basin features two bedrock types: basalt in the lower reaches and sandstone in the headwaters. If the selected concavity index is too high, tributaries will fall below the trunk channel in χ–elevation space. In Fig. 10c, θ is chosen to reflect the typical value of the basalt basins, and tributary channels in the sandstone fall below the trunk channel, meaning that changes in the concavity index can be seen within basins. We therefore suggest that workers must be cautious when using a reference concavity index in determining channel steepness indices in basins with heterogeneous lithology.

4.3 An example of a tectonically active site: Gulf of Evia, Greece

The steepness of channel profiles and presence of steepened reaches (knickpoints) in tectonically active areas can reveal spatial patterns in the distribution of erosion and/or uplift (Densmore et al., 2007; DiBiase et al., 2010; Vanacker et al., 2015) and have the potential to allow identification of active faults (Kirby and Whipple, 2012). However, these systematic spatial patterns in channel steepness may challenge our ability to constrain the concavity index. Our third example is in a tectonically active landscape where we have found spatial variations in the most likely concavity index between catchments proximal to active normal faults. We explore a series of basins draining across faults in the Sperchios Basin, Gulf of Evia, Greece (Fig. 11), predominantly cut into clastic sediments (Eliet and Gawthorpe, 1995). Previous work (Whittaker and Walker, 2015) has demonstrated that catchment morphology reflects interaction with these faults. The rivers are typically characterised by convex longitudinal profiles that commonly have two knickpoints. The upper set of knickpoints is attributed to the initiation of faulting and the resulting growth of topography. The lower set of knickpoints is interpreted as the result of subsequent increase (3–5 ×) in throw rate due to fault linkage (Whittaker and Walker, 2015). The elevations of each group of knickpoints both scale with footwall relief, suggesting that fault throw rates scale with fault segment length. The Gulf of Evia therefore represents a natural experiment where uplift and erosion rates are expected to vary both temporally and spatially.

Steep, smaller catchments tend to drain across the footwalls of these faults, whilst larger catchments drain the landscape behind the faults, through the relay zones between fault segments. We derived the best-fit concavity index for each catchment following each of the five methods (Fig. 12). Given the presence of knickpoints along the river profiles, it is not appropriate to derive the concavity index by linear regression of all log[S]–log[A] data. We find that the concavity index estimated from segmented slope–area analysis is highly variable between catchments (Fig. 12, inset), with a tendency toward abnormally large values, exceeding the upper range of values typically predicted by incision models (Whipple and Tucker, 1999). Values of θ derived using the χ methods are predicted to be relatively low, typically 0.1–0.6 (Fig. 12), and whilst the χ methods do not agree perfectly, they do covary, and are for the most part within uncertainty of each other (with the exception of basins 1 and 20).

Figure 12 The predicted best-fit θ determined using the χ methods (red points) and slope–area methods (blue points shown in inset). Basin numbers correspond to those plotted in Fig. 11.

Figure 13 Profile χ–elevation plots associated with best-fit θ for basin 7, a large catchment with many tributaries draining across a relay zone between normal fault segments (column a), and basin 10, a small, steep catchment draining directly across the footwall segment of a normal fault with few tributaries (column b). Note that MLE values depend on both the number of nodes and the vertical offset from the trunk channel.

A number of authors have suggested that in both highly transient and rapidly eroding landscapes processes other than fluvial plucking or abrasion become important in shaping the channel profile, such as debris flows and plunge pool erosion (DiBiase et al., 2015; Haviv et al., 2010; Scheingross and Lamb, 2017; Stock and Dietrich, 2003). Recent work has suggested that retreat of vertical waterfalls may result in similar concavities to fluvial incision processes operating in lower gradient settings (Shelef et al., 2018), whereas debris flows have been shown to lead to channels with low concavity indices (Stock and Dietrich, 2003). The lowest values of θ = 0.1 at the Evia site typically occur for the small, steep catchments draining across the footwalls of the fault segments (e.g. basin 10; Fig. 13), with higher θ values typical for catchments that do not cross faults or those that cross relay zones (e.g. basin 7; Fig. 13). Plots of χ–elevation such as in Fig. 13 demonstrate that there can be considerable variability in the morphology of tributaries as they respond to adjustment in the trunk channel.

Our aim here is not to provide a comprehensive examination of the topography and tectonic evolution of the Sperchios Basin (Whittaker and Walker, 2015) but to demonstrate the impact of tectonic transience on our ability to quantify the concavity index. Low concavity indices in steep, small catchments draining across the faults may reflect the contribution of debris flow processes to valley erosion at smaller drainage areas, which tends to lead to lower apparent concavity indices in the topography (Stock and Dietrich, 2003). Additionally, these catchments may in effect behave as fluvial hanging valleys (Wobus et al., 2006b). Concavity indices derived using the bootstrap points method are in all cases equal to or lower than values derived using all χ data. This is noteworthy because of the difference in how tributaries are weighted between the two techniques. Using all χ data, longer tributaries have more influence on the calculation of the most likely concavity index, whereas the bootstrap points method weights each tributary equally (since the same number of points are sampled on each tributary). Thus, if the steepness of the channels at low drainage area is influenced by debris flow processes (Stock and Dietrich, 2003), we would expect this to be more influential on the derived concavity index when using the bootstrap points method, resulting in lower θ values.

Finally, it is recognised that transient landscapes are likely settings for drainage network reorganisation (Willett et al., 2014). In the absence of lithologic variability, climate gradients and tectonic transience, gradients in χ in the channel network between adjacent drainage basins are predicted to indicate locations where drainage divides are migrating (toward the catchment with higher χ) and drainage network reorganisation is ongoing (Willett et al., 2014). On the other hand, numerical simulations suggest that spatial variability in uplift is more important than temporal gradients in uplift rates (Whipple et al., 2016). Rivers draining across normal fault systems are often routed through the relay zones between fault tips, where uplift rates are lowest, capturing and rerouting much of the drainage area above the footwall (Paton, 1992). In the Sperchios Basin, this has resulted in strong gradients in χ across topographic divides (Fig. 14), particularly between the large catchments draining the landscape behind the footwall (which have likely been gaining drainage area), and the short, steep catchments draining across the footwall (which have likely been truncated).

Figure 14 Spatial distribution of the χ coordinate in the channel network calculated using A₀ = 1 m²; θ = 0.45. Gradients in χ across topographic divides (black) can indicate planform disequilibrium such that the drainage network may be reorganising. Divides will tend to migrate from low values of χ towards high values in the channel network.

Our analysis of the topography in the Sperchios Basin, whilst not exhaustive, highlights that river profiles and the resulting concavities (and/or k_sn) derived from topography are not alone sufficient to interpret the history of landscape evolution but must be considered alongside other observational data and in the context of a process-based understanding of landscape evolution and tectonics.

5 Conclusions

For over a century, geomorphologists have sought to link the steepness of bedrock channels to erosion rates, but any attempt to do so requires some form of normalisation. This normalisation is required because in addition to topographic gradient, the relative efficacy of incision processes is thought to correlate with other landscape properties that are a function of drainage area, such as discharge or sediment flux. Theory developed over the last four decades suggest that the concavity index may be used to normalise channel gradient, and over the last two decades many authors have compared the steepness of channels normalised to a reference concavity index derived from slope–area data (Kirby and Whipple, 2001; Snyder et al., 2000). In recent years, an integral method of channel analysis has also been developed (Perron and Royden, 2013) that can complement slope–area analysis and via alignment of tributaries provide an independent test of the concavity index.

In this contribution, we have developed a suite of methods to quantify the most likely concavity index using both slope–area analysis and the integral method. In addition to traditional slope–area methods, we also present methods of analysing χ-transformed channel networks that do not require the profiles to be linear from source to outlet but constrain concavity-based collinearity of each tributary and the trunk channel. In a second method, we quantify uncertainty on the predicted value of θ using a subset of points on the tributary network that are randomly assigned within a bootstrap sampling framework. We also test a similar disorder metric that is a minimum when tributaries and trunk channel are most collinear. We test these methods against idealised, modelled landscapes that obey the stream power incision model but have been subject to transient uplift, as well as spatially varying uplift and erodibility, where the concavity index is imposed through the ratio of the exponents m and n.

We find that χ-based methods are best able to reproduce the concavity indices imposed on the model runs. We recommend users calculate the most likely concavity indices using the bootstrap and disorder methods as these provide estimates of uncertainty, although the disorder method is the most tightly constrained of the χ-based methods. The most likely concavities determined from χ-based methods on transient landscapes have low uncertainty because the transient models do not violate any assumptions underlying χ-based methods. The spatially variable model runs, where assumptions of the χ method are violated, still perform better than slope–area analysis in extracting the correct concavity index. This gives us some confidence that in real landscapes, where non-uniform uplift and spatially varying erodibility are likely pervasive, calculated concavities may still reveal useful information about the incision processes. Thus, we hope future workers can calculate reliable, reproducible concavity indices for many small basins in regions with spatially varying uplift, climate, or lithology to test if patterns in the concavity index can be linked to variations in these landscape properties.