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 3Sketch 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 θ.

Download

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

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 4Sketch 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 θ.

Download

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:

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

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.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 5Shaded 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.

Download

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 6Plots 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.

Download

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 7Results 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.

Download

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

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 8Example χ–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.

Download

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 10Exploration 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).

Download

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 11Basins analysed near the Gulf of Evia, Greece, UTM zone 34 N, that interact with active normal faults previously studied by Whittaker and Walker (2015).

Download

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 12The 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.

Download

Figure 13Profile χ–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.

Download

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 14Spatial 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.

Download

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.