Tectonic geomorphology at small catchment sizes – extensions of the stream-power approach and the χ method

Introduction Conclusions References Tables Figures


Introduction
The vast majority of the approaches used to derive information on tectonic processes from topography are based on the analysis of longitudinal river profiles.The fundamental relationship between channel slope S and upstream catchment size A, which is used to infer such information, dates back to a seminal empirical study of Hack (1957) and is often referred to as Flint's law (Flint, 1974).The parameters k s and θ denote steepness index and concavity index, respectively.Understanding and quantitative interpretation of Eq. (1) hinges on the stream-power approach (e.g., Howard, 1994;Whipple and Tucker, 1999;Whipple, 2004;Wobus et al., 2006), where it is assumed that the rate of fluvial erosion in a bedrock channel depends on the product A θ S. In this context, Eq. ( 1) reflects a constant erosion rate along the river as it occurs, e.g., in equilibrium with homogeneous uplift.
In the simplest version of the stream-power approach it is assumed that the erosion rate E is linearly proportional to A θ S. The more general approach implements a power-law relationship where K denotes erodibility.The arbitrary reference catchment size A 0 has been introduced as a scaling parameter Published by Copernicus Publications on behalf of the European Geosciences Union.
S. Hergarten et al.: Tectonic geomorphology at small catchment sizes in order to avoid an odd physical dimension of K. Using this scaling, K describes the erosion rate at a catchment size A 0 and a (hypothetical) channel slope of 1.Although called erodibility, K not only refers to the properties of the channel bed but also contains the effect of precipitation, as the erosion rate in principle depends on the discharge instead of the catchment size.
Physically based models of bedrock incision suggest that the concavity index θ of a steady-state bedrock river under homogeneous conditions depends not only on the constitutive laws of the erosion process but also on the crosssectional geometry of the channels (e.g., Whipple, 2004;Whipple et al., 2013;Lague, 2014).This explains some variation in θ around the value θ ≈ 0.5 originally found by Hack (1957) or around the reference value θ ref = 0.45 being widely assumed for perfect bedrock channels under homogenous steady-state conditions (Whipple et al., 2013;Lague, 2014).
A range of θ between about 0.4 and 0.7 has been found under relatively homogeneous conditions (e.g., Whipple, 2004;Whipple et al., 2013), while a wider range from less than 0.2 in steep headwater channels to more than 1 in some alluvial channels has been reported (Brummer and Montgomery, 2003;Montgomery, 2001;Sofia et al., 2015).Apparent variations in θ may also arise from spatial inhomogeneity or non-steady topography.Analyzing channel slopes at constant catchment sizes, Hergarten et al. (2010) found a strong positive correlation between surface elevation and slope in several orogens, suggesting a correlation between uplift rate and elevation.This correlation will lead to a higher apparent concavity index when following individual rivers, which may explain why the majority of the values of θ found in nature are greater than θ ref = 0.45.
Compared to the concavity index θ , less is known about the exponent n as it cannot be determined from individual equilibrium river profiles under uniform conditions.According to Eq. ( 4), the exponent n can be determined by comparing river segments being in equilibrium with different uplift rates, and the results tentatively suggest that n should not be far away from 1 (Wobus et al., 2006).
Using Eq. ( 2), the evolution of the surface elevation H (x, t) along the stream profile through time under a given uplift rate U follows the partial differential equation where the linear coordinate x follows the upstream direction of the considered river.Both U and K may vary spatially and temporally.
The simplest interpretation of Eq. ( 3) refers to steady-state topography where uplift and erosion are in local equilibrium.Under these conditions, the ratio of uplift rate and erodibility can be directly obtained from the steepness index (Eq.1) according to The most interesting applications of the stream-power erosion equation (Eq.3), however, concern nonequilibrium river profiles due to temporally changing uplift rates or due to climate-induced changes in the erodibility.If such changes are discontinuous, they result in distinct knickpoints propagating in the upstream direction.

The χ transformation and its limitation
Recently, the so-called χ plot (or χ transformation) has introduced the perhaps most important methodic progress in evaluating and interpreting longitudinal river profiles since the seminal work of Howard (1994).It transforms the upstream coordinate x into a new coordinate χ in such a way that the inherent curvature of equilibrium profiles due to the reduction in catchment size in the upstream direction vanishes.The catchment size A can be eliminated from Eq. ( 3) if the transformation satisfies the condition which can be achieved by where x 0 is an arbitrary reference point.As the channel slope is S = ∂H ∂x , the erosion rate (Eq.2) can be written in the form Thus, the local erosion rates is directly related to the slope of the river profile in the H vs. χ representation, and Eq.(3) simplifies to The solutions of this equation and their potential for unraveling the uplift and erosion history have been discussed by Royden and Perron (2013), and a formal inversion procedure for the linear case (n = 1) has been presented by Goren et al. (2014).
The most striking property of the χ transformation is immediately recognized in Eq. ( 9): if U and K are spatially homogeneous, all upstream paths starting from x 0 are described by the same differential equation, so that the H vs. χ curves of all tributaries must collapse with the H vs. χ curve of the main stream.Conversely, spatial inhomogeneity results in a deviation of the curves belonging to different branches that increases in the upstream direction.Thus, a narrow bunch of H vs. χ curves with a nonlinear overall shape is the fingerprint of temporal variations under spatially homogeneous conditions, while a wide but overall straight bunch points towards spatial heterogeneity under steady-state conditions.This simple interpretation, however, only holds as long as the drainage pattern has not changed in the past, since changes in catchment sizes also result in deviations between different branches (Willett et al., 2014;Yang et al., 2015).
Since a clear distinction requires the consideration of a large number of tributaries, the inherent limitation of the stream-power approach to the fluvial regime also limits the χ method.As addressed in several studies, Flint's law (Eq.1), and thus the stream-power erosion equation (Eq.2) with a constant concavity index θ , breaks down at small catchment sizes where lower limits between about 0.1 and 5 km 2 have been reported (Montgomery and Foufoula-Georgiou, 1993;Stock and Dietrich, 2003;Wobus et al., 2006).
The transition from a fluvial regime at large catchment sizes to a regime dominated by hillslope processes is explored by an example from Taiwan in Fig. 1.Based on the recently released SRTM1 digital elevation model (DEM) with a mesh width of 1 arcsecond, flow directions (D8 algorithm; O' Callaghan and Mark, 1984), catchment sizes, and channel slopes were computed for the entire island after filling all local depressions.For comparison, the same analysis was performed on the older SRTM3 DEM with a mesh width of 3 arcseconds.The mean slope (black markers) follows Eq. ( 1) well above a catchment size of a few square kilometers with a steepness index θ ref = 0.45.Clear deviations from this behavior are visible at catchment sizes below about 2 km 2 in the Taiwan data set.These deviations are even more pronounced in the finer SRTM1 data set than in the SRTM3 data set, suggesting that they indeed arise from a limitation of the stream-power approach and not from the potentially inadequate representation of the drainage pattern on coarse DEMs (Stock and Dietrich, 2003).
On the other hand, the number of nodes with a catchment size of A or larger roughly decreases like A −0.5 (Maritan et al., 1996).For a DEM with a mesh width of 1 arcsecond, this means that only some 2 % of all DEM nodes have a catchment size A ≥ 2 km 2 ; thus, about 98 % of all nodes cannot be used in the χ method here.Therefore, an extension of the χ method towards smaller catchment sizes is desirable for deriving maps of χ in order to separate temporal variations from effects of spatial heterogeneity.Beyond increasing the pure data density, it is also helpful with regard to the contest of catchments brought into discussion by Willett et al. (2014).Growing or shrinking catchments are characterized by the curvature in the H vs. χ plot becoming more and more significant close to the migrating drainage divide.Therefore, analyzing the migration of drainage divides quantitatively requires a χ transform free of any bias at small catchment sizes induced by the limited applicability of the stream-power law.

Extending the χ method to small catchment sizes
In the following we present two extensions of the basic relationship between channel slope and catchment size (Eq. 1) towards small catchment sizes and their implementation in the χ method.In most applications of the χ method, the concavity index θ is considered an adjustable parameter and used to improve either the straightness of the H vs. χ plot or the collinearity with tributaries.In the following, the approach with adjustable concavity index θ is denoted χ θ , while χ represents the version with the reference value θ ref = 0.45.However, the curvature of the data in Fig. 1 already suggests that the adjustment of θ may only introduce a limited improvement at small catchment sizes compared to the reference value θ ref .
The approaches presented in the following are intended to be as simple as possible.First, we aim at a representation by a uniform equation without distinguishing different regimes, although the domain below (concerning catchment size) the region where Flint's law holds is sometimes described as the debris flow regime (Stock and Dietrich, 2003).Second, it should involve as few parameters as possible in order to limit the numerical effort of parameter estimation.
The data shown in Fig. 1 suggest that the erosion rate still depends on the catchment size at least for A ≥ 0.01 km 2 , but this dependence is weaker than predicted by the streampower law.A simple modification of Eq. ( 1) consists in adding a constant value a to the catchment size, i.e., to assume The respective modification of the χ transformation is obtained by replacing A in Eq. ( 6) with A + a: This extension can be either considered a one-parameter approach where a is an adjustable parameter, while θ = θ ref is pre-defined, or a two-parameter approach with both a and θ being free parameters.For consistency, the latter is denoted χ θa in the following.
As an alternative approach, a constant value can be added to the term A θ .In order to avoid odd physical dimensions, this term is written in the form b θ , where b has the dimension of an area.With this extension, Eq. ( 1) turns into and the erosion rate (Eq.2) becomes In the linear case (n = 1), this extension can be interpreted as an erosion rate consisting of two additive components being both proportional to the channel slope.One of them depends on the catchment size according to the stream-power law, while the second one is independent of the catchment size and may correspond, for example, to hillslope erosion.Equation ( 12) is essentially the same as the empirical relationship suggested by Stock and Dietrich (2003) for debris flow valleys.Here, s 0 is the hypothetic slope at the valley head (A = 0), and a 2 is the counterpart of the concavity index θ.
The parameters a 1 and s 0 are related to those from Eq. ( 12) by a 1 = b −θ and s 0 = k s b −θ .In this sense the difference between the approaches only concerns the considered regime and the definition of the parameters.We use a parameter b characterizing the catchment size where fluvial erosion and the sum of surface processes independent of the catchment size contribute equally to total erosion, while Stock and Dietrich (2003) used a more abstract parameter a 1 .
The respective modified χ transformation reads Similar to the first approach, χ b (x) refers to the oneparameter version with θ = θ ref in the following, while χ θb denotes the two-parameter version with adjustable parameters b and θ .
As shown by the red and green lines in Fig. 1, both extensions with θ = θ ref do not capture the behavior of the mean slope at small catchment sizes perfectly.While the first version (χ a ) should be better at catchment sizes moderately below the range where the original stream-power approach is valid, the second version (χ b ) should be preferable if the entire range shown in Fig. 1 is considered.
Each of the approaches contains one or two adjustable parameters (a, b and/or θ ) where the optimum value differs from catchment to catchment.As already pointed out by Perron and Royden (2013) for the one-parameter version χ θ with variable θ , determining the respective optimum parameter value is nontrivial.
In the simplest situation, a steady-state topography under homogeneous uplift and erodibility, the H vs. χ plot should be a straight line.Here, the R 2 value (coefficient of determination) of a linear fit and Pearson's correlation coefficient provide equivalent objective functions to be minimized.However, this may lead to systematic bias for transient topographies.An extension based on fitting piecewise linear functions was recently suggested by Mudd et al. (2014), but this algorithm may become numerically expensive, in particular if applied to a large number of catchments or if two adjustable parameters are involved.
Including small catchment sizes in the analysis even facilitates the determination of the adjustable parameters since the collinearity of a large bunch of lines in the H vs. χ plot can be tested.Therefore, a criterion that measures how well the data follow a monotonic relationship between H and χ without being too sensitive to the shape of this relationship (such as R 2 and Pearson's correlation coefficient preferring linear relations) should be used.Spearman's rank correlation coefficient is the most widely used criterion in this context.Here, both the H values and the χ values are sorted independently.Then, an H rank and a χ rank are assigned to each data point, and the correlation coefficient of the two ranks is computed.However, this approach suffers from the χ rank being a discontinuous function of the χ values and thus of the adjustable parameters.As a consequence, the rank correlation coefficient is a piecewise constant function of the parameters with a huge number of discontinuities, which makes its numerical maximization at least theoretically problematic.This problem could be avoided by considering the correlation between the χ values themselves and the H rank as the elevations are fixed values.However, this would introduce a bias towards a certain overall relationship depending on the hypsographic curve of the catchment, so that there is no advantage to the R 2 value or Pearson's correlation coefficient (preferring a linear relationship).
Due to the problems with the rank correlation coefficients discussed above, we suggest an alternative criterion for assessing the collinearity of all rivers in the H vs. χ plot.In a first step, all pairs of χ i and H i are sorted in order of increas-Earth Surf.Dynam., 4, 1-9, 2016 ing H , and the sum is computed.This sum becomes minimal if H and χ are related monotonically and increases with each pair of subsequent points where χ decreases.However, S linearly scales with the absolute χ values; therefore, minimizing S would introduce a bias towards parameters leading to small overall χ values.Thus, S must be rescaled appropriately.As the χ values start from zero, the lowest possible value of S is the maximum χ value, χ max , occurring for the given parameter set.Thus, the straightforward rescaled objective function to be minimized is denoted χ disorder in the following.A perfect monotonic relationship between H and χ is characterized by D = 0. Some attention should be paid to pairs of identical elevation values occurring frequently in integer-valued DEMs.
Here we suggest to assume that all χ values belonging to the same elevation are always in ascending order, so that they do not increase D artificially.

Results and discussion
In the following we compare the different approaches χ a , χ b , χ θ , χ θa , and χ θb using the recently released SRTM1 DEM with a mesh width of 1 arcsecond.Taiwan was selected as a region with high tectonic activity where glaciation only affected rather small regions around the highest mountains (Ono et al., 2005).Therefore, Taiwan should be an almost perfect example of a fluvial landscape.
In order to get a sufficient number of catchments of similar sizes where each catchment contains a significant portion in the fluvial regime, a procedure to automatically delineate catchments with a size A ≈ 100 km 2 was implemented.In a first step, all sites with catchment sizes A < 100 km 2 where the catchment size of the respective flow target is greater than 100 km 2 , and where the site itself makes the largest contribution to its flow target, are determined.These points or, more precisely, their flow targets are considered the base points (x 0 ) of the respective catchments.The drainage pattern is then followed in the upstream direction down to a catchment size of 0.01 km 2 , and the different methods are applied to each of the catchments.All DEM nodes without valid elevation data or where the surface elevation had to be increased when filling local depressions were disregarded.
The topography of Taiwan yields 89 catchments meeting these criteria, with each of them containing between 6464 and 27 732 valid SRTM1 DEM nodes.These catchments are shown in Fig. 2, while Fig. 3 displays the resulting cumulative distribution of the χ disorder of the five approachesχ a , χ b , χ θ , χ θa , and χ θb -compared to the reference χ.As expected, the three approaches involving an adjustable parameter are much better than χ , and χ b is the best among these www.earth-surf-dynam.net/4/1/2016/Earth Surf.Dynam., 4, 1-9, 2016 approaches if the entire range 0.01 km 2 ≤ A ≤ 100 km 2 is considered.For 46 out of the 89 catchments, χ b yields the best approximation (lowest χ disorder) among the three oneparameter approaches, while χ a was found to be best for 17 catchments, and χ θ for 26 catchments.As an immediate consequence of the additional parameters, both two-parameter approaches χ θa and χ θb yield a further improvement.The benefit is, however, weaker than that of the one-parameter approaches towards the version χ involving no adjustable parameters.The difference between χ θa and χ θ b seems to be negligible.
Figures 4 and 5 show two examples of catchments and their representation in the H vs. χ plot.The first one (Fig. 4) is the catchment with the lowest χ disorder in all approaches except χ (ranging from D = 199 for χ θb to D = 230 for χ θ ) located in the mountain range.The χ disorder of the second example (Fig. 5) is more than 2 times higher than in the first catchment (ranging from D = 474 for χ θa to D = 502 for χ θ ).These values are close to the two-thirds quantile of the respective distribution, which means that about two-thirds of the 89 considered catchments have a lower χ disorder than this example.
Taking into account the width of the H vs. χ bunches, none of the two catchments shows a significant overall curvature in the H vs. χ plot except for the lowermost region of the first example (Fig. 4).However, this small part of the catchment is located at the edge of the mountain belt and even in an anthropogenically disturbed region.This finding suggests that spatial heterogeneity has a stronger effect on the H vs. χ plot than potential temporal changes in uplift rate in the past.This heterogeneity may be due to spatial variations in uplift rate, precipitation or, perhaps most likely at the catchment scale, the resistance of the rocks to erosion.
The relevance of spatial heterogeneity to the H vs. χ plot is obviously related to the topology of the drainage network.The example shown in Fig. 4 is an elongated catchment consisting of one main river and several smaller tributaries.The drainage network of the other example (Fig. 5) is characterized by confluences of rivers of similar sizes, and thus subcatchments with comparable χ values may occur at quite large spatial distances.For such a topology, spatial heterogeneity will likely generate diverging segments in the overall H vs. χ plot.In this example, it is readily recognized that a large part of the heterogeneity even arises from a small region when plotting χ (here χ θb ) and elevation in a map (Fig. 6).In a region east of the center of the map, the contour lines of χ are at significantly higher elevations than elsewhere in the domain.This behavior corresponds to the single very steep tributary visible in Fig. 5. Therefore, a limited region east of the center of the map seems to be characterized by a lower erodibility than the rest of the domain.
If only a smaller range of catchment sizes is considered, the differences between the methods partly disappear.Figure 7 shows the same analysis applied down to catchment sizes of 1 km 2 instead of 0.01 km 2 .Here the approach χ θ outcompetes the other one-parameter concepts by providing a better approximation in 78 out of 89 catchments and even comes close to the two-parameter approaches.As the extensions involving the parameters a and b were designed to capture the behavior at small catchment sizes, this result is not surprising.
It is also immediately recognized that restricting the lower limit of catchment size reduces the absolute values of the χ disorder.The reduction mainly arises from the number of valid DEM nodes decreasing by more than 1 order of magnitude to a range from 670 to 2174.However, the χ disorder is diminished by less than 1 order of magnitude here.This nonlinear scaling is presumably related to the scale dependence of spatial heterogeneity.Nonlinear scaling properties can also be expected with respect to the total catchment size and to the resolution of the DEM, but investigating this in detail would go beyond the scope of this study.Thus the χ disorder is found to be a good criterion for comparing different extensions of the χ method and for determining the respective parameter values, but it cannot be used for comparing catchments of different sizes and data obtained from different DEMs.
As a second example we consider the European Alps as an orogen that was heavily affected by glacial erosion in the past.For simplicity, we define the region as the domain inside the 600 m elevation contour line as previously done by Hergarten et al. (2010).Although the properties having an influence on the absolute values of the χ disorder (upper and lower limit of catchment size and DEM resolution) are the same as for the Taiwan example, the values of the χ disorder www.earth-surf-dynam.net/4/1/2016/Earth Surf.Dynam., 4, 1-9, 2016 ≤ A ≤ 100 km 2 .Each curve describes the relative number of the catchments with a χ disorder lower than or equal to the value D on the x axis.
shown in Fig. 8 are significantly higher than those obtained for Taiwan.This increase is presumably related to glaciation causing strong local deviations from fluvial topography.The observed differences between the different approaches, however, persist or become even more pronounced.Here, the method χ b yields the best results among the one-parameter approaches for more than 75 % of all catchments (280 out of 371).
Beyond the goodness of the fit expressed by the χ disorder, the best-fit parameter values may also be taken into account when comparing the different approaches.In this context, the concavity index θ is the most important parameter as it has already been addressed in numerous studies on larger catchment sizes.Figure 9 compares the statistical distributions of the best-fit θ values for the three methods involving θ as an adjustable parameter in the Taiwan example.If θ is the only parameter (χ θ ), the best-fit θ values tend to be below the widely used reference value θ ref = 0.45.This effect becomes more pronounced if points with small catchments size are included.The median for 0.01 km 2 ≤ A ≤ 100 km 2 is θ = 0.33, and 82 out of 89 catchments have θ < 0.45.As the deviations from Flint's law at small catchment sizes can only be compensated for by smaller θ values here, the significant bias towards smaller θ values found for χ θ is not surprising.In return, the two-parameter approach χ θa exhibits a tendency towards values θ > θ ref , reflected in median of θ = 0.56.The other two-parameter approach, χ θb , yields best-fit θ values with a median of θ = 0.47 close to the reference value θ ref = 0.45.While χ θa and χ θb are evenly matched with respect to the χ disorder on average, χ θa obviously needs artificially increased θ values for achieving the best fit.The approach χ θb turns out to be more robust against this bias, although some tendency towards larger θ values occurs if the catchment sizes are restricted to a narrower Each curve describes the relative number of the catchments with an estimated concavity index lower than or equal to the value θ on the x axis.Solid lines refer to fits over the entire range 0.01 km 2 ≤ A ≤ 100 km 2 , while dashed lines correspond to fits for A ≥ 1 km 2 only.range (here, 1 km 2 ≤ A ≤ 100 km 2 ).Under this aspect, the approach χ θb should be superior to χ θa if a wide range of catchment sizes is taken into account.

Conclusions
We have presented and investigated several concepts of extending Flint's law and the χ method towards small catchment sizes.Including points with small catchment sizes into the analysis of stream profiles strongly increases the data density and thus allows for a better distinction between effects of temporal changes in uplift rate or climate and spatial heterogeneity.
Among the approaches considered in this study, an extension of Flints's law similar to an equation originally suggested for debris flow channels (Stock and Dietrich, 2003) turned out to be the most suitable concept if a wide range of catchment sizes is included.The respective definition of the extended χ transform (Eq.15) can be implemented either as a two-parameter approach where both θ and b are adjustable parameters or as a one-parameter approach where b is variable and θ = θ ref with θ ref = 0.45 or any other fixed reference value.
Minimizing the χ disorder defined in Eq. ( 18) provides a simple way to determine the best values of the adjustable parameters.The χ disorder refers to the collinearity of tributaries and does not require any further assumptions such as spatial homogeneity or an uplift rate being constant over distinct time intervals and should thus be applicable in a wide context.

Figure 1 .
Figure 1.Relationship between mean channel slope and catchment size for the topography of Taiwan.Channel slopes and catchment sizes were derived from the SRTM1 and SRTM3 DEMs, and mean slopes were obtained from logarithmic bins with a factor of √ 2 (black markers).

Figure 2 .Figure 3 .
Figure 2. Map of the 89 considered catchments in Taiwan with catchment sizes A ≈ 100 km 2 .The two catchments bordered in magenta and yellow are considered in detail in Figs.4-6.

Figure 4 .
Figure 4.The mountainous catchment in Taiwan with the lowest χ disorder.(a) Topography and drainage pattern for A ≥ 0.01 km 2 .The largest river is drawn in light blue.(b) H vs. χ plots of the main river.χ 0 refers to θ = 0, so that χ 0 = x, and the plot describes the original river profile.(c) H vs. χ plots for the entire drainage network.The plots are shifted horizontally in order to avoid overlapping curves.The black lines show the part of the drainage network with A ≥ 1 km 2 .

Figure 5 .
Figure 5.A catchment in Taiwan with a rather high χ disorder.(a) Topography and drainage pattern for A ≥ 0.01 km 2 .The largest river is drawn in light blue.(b) H vs. χ plots of the main river.χ 0 refers to θ = 0, so that χ 0 = x, and the plot describes the original river profile.(c) H vs. χ plots for the entire drainage network.The plots are shifted horizontally in order to avoid overlapping curves.The black lines show the part of the drainage network with A ≥ 1 km 2 .

Figure 6 .Figure 7 .
Figure 6.Map of elevation (coded by color) and χ θb values (contour lines) of the catchment considered in Fig. 5.The contour line interval is 0.5 km, and the lines of χ θb = 2.5 km and χ θb = 5 km are emphasized in white.

Figure 8 .
Figure 8. Cumulative distribution of the χ disorder for the 371 considered catchments in the European Alps for 0.01 km 2≤ A ≤ 100 km 2 .Each curve describes the relative number of the catchments with a χ disorder lower than or equal to the value D on the x axis.

Figure 9 .
Figure 9. Cumulative distribution of the concavity index θ for the 89 considered catchments in Taiwan.Each curve describes the relative number of the catchments with an estimated concavity index lower than or equal to the value θ on the x axis.Solid lines refer to fits over the entire range 0.01 km 2 ≤ A ≤ 100 km 2 , while dashed lines correspond to fits for A ≥ 1 km 2 only.