Articles | Volume 10, issue 4
Short communication
11 Aug 2022
Short communication |  | 11 Aug 2022

Short communication: Forward and inverse analytic models relating river long profile to tectonic uplift history, assuming a nonlinear slope–erosion dependency

Yizhou Wang, Liran Goren, Dewen Zheng, and Huiping Zhang

The long profile of rivers is shaped by the tectonic history that acted on the landscape. Faster uplift produces steeper channel segments, and knickpoints form in response to changes in the tectonic uplift rates. However, when the fluvial incision depends non-linearly on the river slope, as commonly expressed with a slope exponent of n≠1, the links between tectonic uplift rates and channel profile are complicated by channel dynamics that consume and form river segments. These non-linear dynamics hinder formal attempts to associate the form of channel profiles with the tectonic uplift history. Here, we derive an analytic model that explores a subset of the emergent non-linear dynamics relating to consuming channel segments and merging knickpoints. We find a criterion for knickpoint preservation and merging, and we develop a forward analytic model that resolves knickpoints and long profile evolution before and after knickpoint merging. We further develop a linear inverse scheme to infer tectonic uplift history from river profiles when all knickpoints are preserved. Application of the inverse scheme is demonstrated over the main trunks of the Dadu River basin that drains portions of the east Tibetan Plateau. The model infers two significant changes in the relative uplift rate history since the late Miocene that are compatible with low-temperature thermochronology. The analytic derivation and associated models provide a new framework to explore the links between tectonic uplift history and river profile evolution when the erosion rate and local slopes are non-linearly related.

1 Introduction

Bedrock rivers that incise into tectonically active highlands are sensitive to changes in the tectonic conditions (Whipple and Tucker, 1999). Upon a change in the rock uplift rate with respect to a base level, the river steepness changes (Wobus et al., 2006; Whipple and Tucker, 2002), which in turn changes the local incision rate. In particular, an increase in uplift rate generates steeper slopes that facilitate faster incision, which can eventually lead to incision–uplift equilibrium. However, equilibrium is not achieved synchronously across the river long profile. Upon a change in the tectonic uplift rates, a knickpoint forms that divides the profile into reaches with different steepness and erosion rates (Rosenbloom and Anderson, 1994; Berlin and Anderson, 2007; Oskin and Burbank, 2007). Below the knickpoint, the steepness and erosion rate have already been shaped by the new tectonic conditions, and above the knickpoint, river steepness and erosion rate correspond to the previous conditions (Niemann et al., 2001; Kirby and Whipple, 2012). The erosion rate gradient across the knickpoint promotes knickpoint migration upstream, gradually changing the proportion of the channel that is equilibrated to the new tectonic conditions. For these dynamics, knickpoints are viewed as moving boundaries that separate channel reaches, recording different portions of the tectonic uplift history (e.g., Pritchard et al., 2009; Whittaker and Boulton, 2012).

Since the links between tectonic uplift history and river shape are mediated by fluvial incision, resolving these links requires a fluvial incision theory. The stream-power incision model (SPIM) is widely used to describe detachment-limited vertical incision into channel bedrock, over long timescales (commonly beyond millennial) and large length scales (Howard and Kerby, 1983; Howard, 1994; Whipple and Tucker, 1999; Lague, 2014; Venditti et al., 2019). The SPIM represents the rate of bedrock incision, E (L/T, length/time) as a power-law function of channel slope (S=z/x, L/L) and upstream drainage area (A, L2), a proxy for both discharge and channel width (Howard and Kerby, 1983):

(1) E ( x t ) = KA ( x ) m z t , x x n ,

where x (L) denotes a spatial coordinate along the channel and t (T) is time. The channel erodibility, K (L12m/T), primarily depends on the bedrock lithology and the effective rate of precipitation (Whipple and Tucker, 1999; Snyder et al., 2000). The positive exponents, m and n, control the sensitivity of incision rate to the drainage area and slope, respectively. Assigning Eq. (1) in a topography conservation equation gives rise to a partial differential equation describing the time–space evolution of the fluvial channel long profile:

(2) z ( t , x ) t = U t , x - KA m z t , x x n ,

where U (L/T) is the rate of tectonic uplift. Notably, the formulation of Eq. (1) represents many simplifications of the processes of river bedrock incision. For example, it does not explicitly account for incision thresholds, discharge variability, sediment flux incision sensitivity and dynamic changes in channel width (Lave and Avouac, 2001; Whipple and Tucker, 2002; Duvall et al., 2004; Dibiase et al., 2010). Nonetheless, Gasparini and Brandon (2011) argued that many of these processes could still be approximated by modifying the exponents, m and n.

Equation (2) is a non-linear advection equation for the elevation, where U acts as a forcing term. Consequently, Eq. (2) predicts the first-order dynamics of bedrock rivers, whereby knickpoints form in response to uplift rate changes and migrate upstream. The relative simplicity of Eq. (2) presents a unique opportunity for an analytic exploration of channel dynamics in response to changing tectonic and environmental conditions. In particular, when the analytic solution is sufficiently simple, its representation can be used as part of forward models that predict topographic evolution (e.g., Steer, 2021) and inverse models that infer the tectonic uplift history from observations of river long profiles (Fox et al., 2015; Rudge et al., 2015; Gallen and Fernández-Blanco, 2021; Goren et al., 2022).

Previous general analytic exploration of Eq. (2) (e.g., Luke, 1972; Weissel and Seidl, 1998; Prichard et al., 2009; Royden and Perron, 2013) identified that upon a change in uplift rate that induces a long-profile steepness change, portions of the solution, representing the river profile, could form that are not strictly associated with the change in uplift rate, and portions of the solution that hold tectonic information may be lost. More specifically, when U increases and n<1 or U decreases and n>1, “stretched zones” form along the river long profile that are not associated with any particular tectonic input (Royden and Perron, 2013). When U increases and n>1 or U decreases and n<1, some portions of the channel reach are consumed at knickpoints (Royden and Perron, 2013). Unlike the non-linear cases, when n=1, stretched and consumed channel reaches do not occur, and there is a 1-to-1 mapping between the tectonic uplift history and the river long profile. For this reason, so far, only analytic solutions that assume slope–incision linearity (n=1) were adapted into forward (Steer, 2021) and inverse models (for a recent review see, Goren et al., 2022) of tectonically forced fluvial landscape evolution.

While some field studies support the slope–incision linearity assumption (e.g., Wobus et al., 2006; Ferrier et al., 2013; Schwanghart and Scherler, 2020), a growing body of work shows that n could be different than unity and is mostly inferred to be >1 (Whipple et al., 2000; Harkins et al., 2007; Lague, 2014; Harel et al., 2016). From a process perspective, large values of n were suggested to stem from incision thresholds, small discharge variability and dynamic channel narrowing (Anthony and Granger, 2007; Ouimet et al., 2009; Dibiase et al., 2010; Lague, 2014; Gallen and Fernández-Blanco, 2021).

When n=1, it is well accepted that, under a well-constrained erodibility, a full tectonic uplift rate history can be retrieved from river long profiles (e.g., Goren et al., 2022, and references therein). However, when n≠1, the potential formation of stretched zones and consumption of channel segments challenge the links between river long profiles and the tectonic uplift history. On the one hand, some studies (e.g., Kirby and Whipple, 2012) proposed that, even for n≠1, knickpoint ages could be determined based on the known channel incision rates up- and down-stream of the knickpoints by using paleo-channel projection, and other studies attempted a non-linear inversion to infer uplift histories with variable values of n (Pritchard et al., 2009; Roberts and White, 2010; Paul et al., 2014). On the other hand, Royden and Perron (2013) showed that information of tectonic uplift history could be entirely lost when reaches of the channel profile are fully consumed. Therefore, the questions of to what extent the channel long profile records and preserves a full tectonic uplift rate history and if and how this history can be retrieved when n≠1 are still outstanding.

The current study addresses these questions by developing an analytic description of the evolution of channel long profiles for the cases where channel reaches may be consumed; namely U(t) is a staircase decreasing function and n<1, or U(t) is a staircase increasing function and n>1. The latter scenario is particularly applicable for tectonically active and rejuvenated landscapes. Unlike previous analytic explorations (e.g., Luke, 1972; Weissel and Seidl, 1998; Royden and Perron, 2013) that solved for the long profile as a whole, the current analysis focuses on knickpoint kinematics from a Lagrangian perspective that follows the knickpoints along their migration path. With this approach, we develop a criterion for knickpoint preservation and merging, a forward analytic model that can propagate knickpoints beyond merging, and a linear inverse model constrained by knickpoint preservation. The current study focuses on the theory and model derivation, and the operation of the inverse model is demonstrated along the Dadu River basin that drains the steep margins of the east Tibetan Plateau.

2 Theoretical background

The SPIM model, Eq. (1), predicts that for channel segments that erode at the uniform rate, the channel slope scales as a power-law function of the drainage area:

(3) z x = k s A - θ .

Notably, the power-law scaling in Eq. (3) was originally identified based on topographic data (e.g., Morisawa, 1962; Hack, 1973; Flint, 1974) and is thus independent of any incision model. In the context of the SPIM, θ=m/n and ks=(E/K)1/n (L2m/n) are commonly referred to as the channel concavity and steepness indices, respectively (Wobus et al., 2006). An alternative perspective to Eq. (3) emerges when integrating it along the channel, while assuming constant E/K. Following such an integration, a linear relation emerges between the elevation, z, and the parameter χ (L) (Perron and Royden, 2013):


where zb is the base-level elevation, and the area scale factor A0 (L2) is introduced to maintain the χ dimensions to length. The parameter χ depends only on the drainage area distribution along the channel, and it can easily be calculated for any m/n as part of basic morphometric analysis (Perron and Royden, 2013). When setting A0=1 L2, the slope of the χz plot becomes channel steepness index, ks.

Under steady-state conditions, when dz/dt=0 and E=U, the SPIM steepness index becomes a function of the tectonic uplift rate:

(6) k s = ( U / K ) 1 / n .

When U varies in time, Eq. (6) can be used to express transient conditions, where a channel segment is eroding at a rate that corresponds to some previous uplift rate, Up (Niemann et al., 2001; Goren et al., 2014). In this case, its steepness index could be expressed as

(7) k s _ p = ( U P / K ) 1 / n .
3 Slope-break knickpoint migration

A slope-break knickpoint occurs when there is an abrupt change in the slope and steepness index along a channel long profile (Wobus et al., 2006; Haviv et al., 2010). Within the scope of the SPIM, slope-break knickpoints are commonly associated with a step change in the rate of base level lowering. When the rate increases, the slope and steepness index downstream the knickpoint are greater, and the slope break is convex upward. When the rate decreases, the slope and steepness index below the knickpoint are smaller, and the slope break would appear as a concave kink along the overall concave channel profile. In this latter case, alluviation might occur below the knickpoint, and the assumption of detachment-limited conditions might be violated. This behavior is beyond the scope of the current analysis.

To predict the retreat rate of slope-break knickpoints, we develop a model based on long profile linearization in the proximity of the knickpoint as shown in Fig. 1.

Figure 1(a) Schematics of a channel profile evolution in response to an increase in the relative uplift rate from U0 to U1 (revised from Goren et al., 2014). The blue solid line shows the steady-state channel under uplift rate U0. The black solid and gray dashed lines show the transient channel at time t and t+dt. The black dashed line shows the final steady-state channel under uplift rate U1. (b) Schematics of knickpoint retreat (revised from Wang et al., 2017). Points A and D are the knickpoint positions at time t and t+dt. Evolution of the channel profile in the time step dt is shown as the transition from zt to zt+dt. The black dashed line AG is parallel to the x axis.


Figure 1a shows the predicted channel profile evolution following a step increase in the rock uplift rate from U0 to U1 and n>1. The figure emphasizes that below and above the knickpoint the channel segments erode at rates that correspond to the new (U1) and old (U0) uplift rates, respectively, and their corresponding steepness indices are ks_1=(U1/K)1/n and ks_0=(U0/K)1/n. Figure 1b shows the linearized channel segments near the knickpoint. The river profile varies from zt to zt+dt during time step dt, accompanied by the knickpoint migrating from points A to D. Segment DG represents the vertical change in knickpoint location, and it can be expressed as

(8) DG = z t + d t x + d x - z t + d t x = z x 1 v H d t ,

where vH is the horizontal retreat velocity for the knickpoint (hereafter, knickpoint celerity). Figure 1b shows that

(9) DG = DB + BG ,

where DB is a function of the difference between the present uplift rate (U1) and previous river incision rate, U0:

(10) DB = U 1 - U 0 d t .

The segment BG is the elevation difference between points A and B:

(11) BG = z x 0 v H d t .

Combining Eqs. (8)–(11), we solve for the knickpoint celerity:

(12) v H = U 1 - U 0 z x 1 - z x 0 ,

which resembles the derivation of Whipple and Tucker (1999). Assigning Eqs. (1) and (6)–(7) into (12), vH can be re-written as

(13) v H = K ( k s _ 1 n - k s _ 0 n ) ( k s _ 1 - k s _ 0 ) A m / n = k s _ 1 n ( 1 - γ 0 _ 1 n ) k s _ 1 ( 1 - γ 0 _ 1 ) KA m / n = k s _ 1 n - 1 ( 1 - γ 0 _ 1 n ) ( 1 - γ 0 _ 1 ) KA m / n ,

where γ0_1=ks_0/ks_1. Accordingly, the fluvial response time (Whipple and Tucker, 1999) of the knickpoint, τ(xp), is expressed as

(14) τ ( x p ) = 0 x p 1 v H d x = 0 x p k s _ 1 ( 1 - γ 0 _ 1 ) k s _ 1 n ( 1 - γ 0 _ 1 n ) 1 KA m / n d x = k s _ 1 ( 1 - γ 0 _ 1 ) k s _ 1 n ( 1 - γ 0 _ 1 n ) 1 KA 0 m / n χ ( x p ) .

The response time is the time for a perturbation, e.g., a knickpoint, to propagate from the river outlet (x=0) to its present location xp. Alternatively, τ(xp) can also be thought of as the knickpoint age (Gallen and Wegmann, 2017), or the time before the present when the knickpoint was formed at the river outlet.

Equations (8)–(14) are developed for the migration of a single knickpoint based on a Lagrangian perspective, i.e., in the reference frame of the migrating knickpoint. Accordingly, Eqs. (13)–(14) predict that knickpoint celerity and response time depend only on the steepness indices immediately above and below the knickpoint and are independent of the steepness indices at lower reaches below lower, newer knickpoints. This means that as long as knickpoints do not merge, as discussed in the following section, knickpoint celerity and response time are not affected by later changes in the tectonic uplift rate and channel steepness.

Equations (13)–(14) reveal that knickpoint dynamics depends on both the slope exponent, n, and the steepness ratios, γ. Notably, although the derivations in this section are based on convex-up knickpoint (increasing U and n>1), Eqs. (12)–(14) are also valid for concave knickpoints (decreasing U and n<1; see details in Supplement Sect. S1). For n=1, vH and τ(xp) are independent of the steepness indices and their ratio. Section S2 in the Supplement compares the current derivation to previous models of knickpoint celerity (Rosenbloom and Anderson, 1994; Weissel and Seidl, 1998; Oskin and Burbank, 2007; Royden and Perron, 2013; Castillo et al., 2017).

4 Knickpoint preservation and merging

When more than a single knickpoint propagates upstream a channel profile andn≠1, the sensitivity of knickpoint celerity to ks and γ leads to potentially complex interactions between the knickpoints. Considering the case of n>1 and two knickpoints that formed by two step increases in tectonic uplift rate: kp1 formed when U0 changed to U1 and kp2 formed when U1 changed to U2 (U2>U1>U0); then the celerity of knickpoint kp2 is larger than that of kp1, the distance between them gradually decreases and the channel segment between them is consumed (see a detailed derivation in Appendix A). Consequently, depending on the knickpoints' relative celerity and the channel length, kp2 can eventually reach kp1, and the two knickpoints merge (referred to as consuming knickpoint in Royden and Perron, 2013). Here, “consumption” is reserved for channel segments that are shortened by a fast-migrating knickpoint, and “merging” is reserved for knickpoints to highlight the different dynamics of the merged knickpoint from that of two knickpoints that joined to form it. To elucidate knickpoint merging dynamics, we derive an expression for the time of knickpoint merging. Assuming that kp1 formed at time t=0 and that kp2 formed at time t=T1, Eq. (14) is used to express the χ values of the two knickpoints at any time T>T1 as

(15) χ kp 2 = T K k s _ 2 n ( 1 - γ 1 _ 2 n ) k s _ 2 ( 1 - γ 1 _ 2 ) and χ kp 1 = ( T + T 1 ) K k s _ 1 n ( 1 - γ 0 _ 1 n ) k s _ 1 ( 1 - γ 0 _ 1 ) ,

where γ1_2=ks_1/ks_2. Knickpoint merging occur at time Tm when χ(kp1)=χ(kp2). The ratio Tm/T1 is expressed as

(16) T m / T 1 = γ 1 _ 2 n 1 - γ 0 _ 1 n γ 1 _ 2 1 - γ 0 _ 1 / 1 - γ 1 _ 2 n 1 - γ 1 _ 2 - γ 1 _ 2 n 1 - γ 0 _ 1 n γ 1 _ 2 1 - γ 0 _ 1 .

Equation (16) predicts that the timing of knickpoint merging depends on the ratios of channel steepness indices but not on the steepness indices themselves. We present a detailed description of the relationship between Tm/T1, the slope exponent and the steepness ratios in Figs. 2 and 3.

Figure 2 shows the results for convex-up consuming knickpoints (n>1 and increasing U). When γ1_2=γ0_1, the ratio Tm/T1 decreases with n (Fig. 2a). This means that a higher slope exponent reduces the life expectancy of knickpoints. Figure 2a also shows that for a constant n, lower steepness index ratio leads to lower values of Tm/T1. To explore the dependency of Tm/T1 on γ1_2 and γ0_1, we fix n=2 and vary each of the steepness ratios independently (Fig. 2b, c). Comparing Fig. 2b and c, it is found that Tm/T1 is more sensitive to γ1_2 than to γ0_1, indicating that the celerity of the younger knickpoint has a greater control over the timing of knickpoint merging.

Figure 2The duration of convex knickpoint preservation as a function of slope exponent n (a), γ1_2 (b) and γ0_1 (c). In (a), the steepness ratios are equal, γ1_2=γ0_1. In (b) and (c), n=2. The analysis assumes two knickpoints, kp1 (upper) and kp2 (lower), generated by two step increases in tectonic uplift rates, T1 corresponds to the duration between the formation of kp1 and the formation of kp2, and Tm corresponds to the time from the emergence of kp2 to its merging with kp1. γ1_2 is the ratio of steepness indices above and below kp2, and γ0_1 is the ratio of steepness indices above and below kp1.


Figure 3The duration of concave knickpoint preservation as a function of the slope exponent n (a), γ1_2 (b) and γ0_1 (c), under decreasing U and n<1. In (a), the steepness ratios are equal, γ1_2=γ0_1. In panels (b) and (c), n=0.5.


For the case of concave-up consuming knickpoints (n<1 and decreasing U), Fig. 3a shows that the ratio Tm/T1 (when γ1_2=γ0_1) increases with increasing n, and for a constant n, a higher steepness ratio leads to a higher Tm/T1 ratio. This means that a lower uplift rate, U2 (with a lower steepness index below knickpoint kp2), leads to a shorter time to knickpoint merging Tm. In Fig. 3b, c, n is fixed at 0.5, and the steepness ratios change. Here as well, an inverse dependency is observed with respect to the convex slope-break knickpoints, showing that Tm/T1 is more sensitive to γ0_1 than to γ1_2, indicating that the preservation time of kp1 is more sensitive to its own celerity than to that of the younger knickpoint.

We note that when n=1, χ(kp1)>χ(kp2) always holds, indicating that within the framework of the linear SPIM, knickpoints are always preserved and merging cannot occur.

Upon knickpoint merging, only a single knickpoint propagates along the channel, and the steepness indices above and below the merged knickpoint correspond to ks_0 and ks_2, respectively. Based on Eq. (13), the instantaneous merged knickpoint celerity becomes

(17) v H _ after _ merger = k s _ 2 n ( 1 - γ 0 _ 2 n ) k s _ 2 ( 1 - γ 0 _ 2 ) KA ( x p ) m / n ,

where γ0_2=ks_0/ks_2. The channel reach between the two knickpoints is fully consumed, and the channel profile holds no record of U1. Consequently, evaluating the merged knickpoint age by using Eq. (14) and the steepness indices above and below the merged knickpoint does not yield a meaningful answer. The reason is that upon merging, the steepness indices above and below the merged knickpoint correspond to the steepness indices above the older knickpoint and below the younger knickpoint, respectively. Critically, the channel profile does not hold any evidence that knickpoints have merged, and the river profile would be indistinguishable from a case of a single step increase in uplift rate from U0 to U2.

5 A forward analytic model for knickpoint and channel long profile evolution

The elevation change of slope-break knickpoint, zt,x=zt,x=xp(t), formed by a step increase in the uplift rate from U0 to U1, can be expressed as

(18) d z d t = z t + z x d x d t ,

where dxdt=dxp(t)dt=vH is the knickpoint celerity. Combining Eqs. (2), (13) and (18) yields

(19) d z ( t , x p ( t ) ) d t = U t - KA m k s _ 1 A - m n n + k s _ 1 A - m n k s _ 1 n 1 - γ 0 _ 1 n k s _ 1 1 - γ 0 _ 1 KA m n = U t - U 1 + U 1 ( 1 - γ 0 _ 1 n ) ( 1 - γ 0 _ 1 ) .

Integrating Eq. (19) to solve for the knickpoint elevation leads to

(20) z t , x p ( t ) = 0 t U t - U 1 + U 1 ( 1 - γ 0 _ 1 n ) ( 1 - γ 0 _ 1 ) d t .

As long as knickpoints do not merge, the second and third terms of the integrand in Eq. (20) are time invariant, and the elevation of the knickpoint could be more simply expressed as

(21) z t , x p t = 0 t U t d t + 1 - γ 0 _ 1 n 1 - γ 0 _ 1 - 1 U 1 t .

Equation (21) predicts the elevation of knickpoints for all values of n as the sum of the time integral over the uplift rate history and a term that depends on the steepness index ratio. Before knickpoint merging, Eqs. (14) and (21) represent a closed-form analytic solution for slope-break knickpoint positions (χ and elevation). When only a single knickpoint propagates along the channel profile, Eqs. (14) and (21) reduce to the simpler form presented in Mitchell and Yanites (2019) (see details in Sect. S3 in the Supplement). When n=1, Eq. (21) become a function of the uplift history only (Goren et al., 2014), z(txpt)=0tUtdt.

Next, we combine Eq. (21), which is conditioned by knickpoint preservation, with Eq. (16) that predicts the duration of preservation to generate a piecewise solution for knickpoint elevation before and after knickpoint merging. We consider the case of two knickpoints, kp1 and kp2, generated at times t=0 and t=T1, respectively, by two step increases in U, U2>U1>U0 and n>1. The time of merging, Tm, measured with respect to T1 is constrained by Eq. (16). For any time t<Tm+T1, the elevations of kp1 and kp2 are predicted by Eq. (21), and when assigning the knickpoint ages t= τ(xp1) and t-T1=τxp2, which corresponds to the time since the change in U(t) that generated the knickpoints. Upon merging, for t>Tm+T1, the elevation of the merged knickpoint,  zkp_12, with respect to the formation time of kp1 (t=0) can be expressed as follows.


Before merging, the horizontal position of the knickpoints can be expressed as the inverse of Eq. (14):

(24) x p t = χ - 1 KA 0 m / n t k s _ 1 n ( 1 - γ 0 _ 1 n ) k s _ 1 ( 1 - γ 0 _ 1 ) ,

where again t= τ(xp) is the knickpoint age. After merging, for t>Tm+T1

(25) x p t = χ - 1 KA 0 m / n T m k s _ 2 n ( 1 - γ 1 _ 2 n ) k s _ 2 ( 1 - γ 1 _ 2 ) + KA 0 m / n ( t - T m - T 1 ) k s _ 2 n ( 1 - γ 0 _ 2 n ) k s _ 2 ( 1 - γ 0 _ 2 ) .

While Eqs. (22)–(25) present a simple case of two merging knickpoints, it is possible to use Eq. (16) to calculate the timing and order of multiple knickpoint merging events, including the merger of already merged knickpoints, and to develop a tailored piece-wise analytic solution for their elevation.

When deriving an analytic solution for the channel long profile as a function of time, Eqs. (21)–(23) are used for knickpoint elevation, Eqs. (24)–(25) are used for the knickpoint x positions and Eq. (14) is used for the knickpoint χ values. The channel profile between knickpoints is represented in the χz domain as a linear line connecting the knickpoints. We use our analytic forward model to illustrate long-profile and knickpoint time evolution before and after knickpoint merging (Fig. 4). In the artificial case, the long profile (Fig. 4a) and χz plot (Fig. 4b) show a river that was originally under steady-state conditions with the tectonic uplift rate of 0.10 mm a−1. Then, we set two step increases in the rates to be 0.5 mm a−1 at 0.8 Ma and 1.0 mm a−1 at 0.5 Ma. The changes in the uplift rates produce two knickpoints. The lower knickpoint generated by the higher uplift rate migrates faster than the higher one and consumes the higher knickpoint at ∼0.2 Ma, causing knickpoint merging (Fig. 4c). After ∼0.2 Ma, only one knickpoint is observed on the river long profile. To demonstrate the validity of the analytic forward model, Fig. 4 also compares the analytic solution and a 1-D upwind first-order finite-difference solver of Eq. (2) and shows a consistency.


Figure 4Comparison between the analytic forward model and a numerical model for the evolution of a river long profile, with drainage area set by Hack's law, A=ka(L-x)h. The model parameters are n=3, K=2×10-11 m−1.7 a−1, m=1.35, ka=3.5 and h=1.7. The time (dt) and space (dx) steps in the numerical model are 10 a and 10 m, respectively. The applied uplift rate history is U0=0.10 mm a−1 prior to 0.8 Ma, U1=0.5 mm a−1 between 0.5–0.8 Ma and U2=1.0 mm a−1 between 0.0–0.5 Ma. L, total river length, is 6 km (hillslope length is 1 km). The analytic (colored, solid) and numerical solutions (black, dashed) match in the xz (a) and χz (b) domains. Panel (c) depicts the river χz long profiles offset in elevation, demonstrating knickpoint merging dynamics. Knickpoint kp1 formed at 0.8 Ma, as a response to the increase in the uplift rate from U0 to U1. Knickpoint kp2 formed at 0.5 Ma, due to uplift rate increase from U1 to U2. At ∼0.2 Ma, kp1 merged with kp2


6 An inverse model to infer tectonic uplift rate history

6.1 Description of the inversion algorithm

Here, the analytic solution for knickpoint evolution is used to derive a linear inverse model for retrieving the tectonic uplift history from the river long profile. The inverse model relaxes the critical assumption of n=1 that was a precondition for previous linear inverse models (Goren et al., 2022) and allows inference of the uplift history for any value of n, under two assumptions: first, if n>1, U(t) is a monotonically increasing staircase function and if n<1, U is a monotonically decreasing function. Second, all the knickpoints are preserved within the time resolved by the model. The model is based on the block uplift assumption, whereby a suite of basins and tributaries experience and respond to the same time-dependent tectonic uplift history U(t), and the block has a uniform erodibility. The model infers a single best-fit history, U(t), based on the long profiles of the analyzed rivers and tributaries.

Changes in U through time emerge as a series of knickpoints with elevations and χ values, z1,χ1,z2,χ2,zq-1,χq-1, which are duplicated across basins and tributaries. The basin outlets are at z0=0,χ0=0, and the highest χ channel head is identified with zq,χq=χmax. The knickpoints are used to divide the χz space into segments. Segment j, between χj-1,χj, is characterized by a uniform steepness index that shaped the river profile during time interval tj-1,tj, where time tj=τj is the age of knickpoint j. The steepness indices of channel segments between the knickpoints are used for constraining knickpoint ages based on Eq. (14). The uplift rate responsible for the formation of each knickpoint is constrained based on the steepness index below the knickpoint by using Eq. (7). Consequently, a full uplift rate history, subject to the assumptions of no merged knickpoints and a staircase uplift change, can be derived.

Difficulty may arise because tj=τj based on Eq. (14) and Uj in Eq. (7) depend on the erodibility, K, whose value is commonly poorly constrained. Thus, following Goren et al. (2014), we present a K-independent version for the knickpoint age and uplift rate. To derive a K-independent knickpoint age, Eq. (14) is multiplied by an erosion rate scale factor, KA0m/nks_jn1-γjnks_j1-γj, (L/T). The scaled knickpoint age, tj*=τj*, with dimensions of length, becomes

(26) t j * = t j KA 0 m / n k s _ j n 1 - γ j n k s _ j 1 - γ j = χ j .

Namely, the scaled knickpoint age corresponds to the χ value of the knickpoint. To derive a K-independent uplift rate, Eq. (7) is divided by a steepness index scale factor, A0m/n, (L2m/n), yielding a non-dimensional K-independent uplift rate:

(27) U j * = A 0 - m / n k s _ j = A 0 - m / n ( U j / K ) 1 / n .

These specific scaling choices allow the use of Eqs. (26)–(27) to predict knickpoint elevations with natural dimensions as explained in Appendix B. Importantly, Eqs. (26)–(27), which describe the scaled uplift rate history, (Uj*,τj*), are not only K-independent but also n-independent. This means that as long as K and n are spatially uniform, a scaled uplift rate history could be inferred without prior knowledge of K and n.

We propose the following three steps for the application of the inverse model. First, the data of basins and tributaries are considered in the χz domain. Calculating the χ value requires calibrating for the concavity index, m/n. We propose a tributary and basin collapse approach (e.g., Perron and Royden, 2013; Goren et al., 2014; Shelef et al., 2018) or the disorder approach (e.g., Hergarten et al., 2016; Gailleton et al., 2021), which finds the m/n that minimizes the scatter in the χz domain.

Second, the χz domain is divided into q segments along the χ space, χj(j=0,1,2,q). The division points are considered to be slope-break knickpoints that formed in response to step changes in uplift rate. The scaled age of the knickpoints is defined based on Eq. (26) as τj*=χj. Then, linear regression is applied in the χz domain, independently for each segment. The slope of the regression is identified as ks_j, from which Uj* is defined based on Eq. (27).

Segment division should ideally be based on division points that represent true slope-break knickpoints. Several algorithms have been previously proposed to identify slope-break knickpoints (e.g., Mudd et al., 2014). Here, we suggest a different approach that relies on the simplicity and efficiency of the inverse model. We propose running the inversion procedure many times, while choosing the number and location of division points randomly. The quality of the solution with a specific number and location of division points could be evaluated based on an optimization criterion, such as a misfit. Mudd et al. (2014) used the Akaike information criterion (Akaike, 1974) to balance the goodness of fit against model complexity. Here, we consider a simpler misfit function that penalizes models with more knickpoints (more parameters) for their excess complexity:

(28) misfit = 1 N / M i = 1 N ( z i - z i ̃ ) 2 ,

where zi and zĩ are the measured and predicted elevations at pixel i, respectively. N is the total number of data along the river long profiles, and M=q is the number of division points, or the number of parameters. zĩ is obtained by integrating U* along the χ(or t*) axis:

(29) z i ̃ = 0 t i * = χ i U * ( t * ) d t * = a = 1 j U a * ( t a * - t a - 1 * ) + U j + 1 * t i * - t j * ,

where pixel i is located between knickpoints j and j+1. Appendix B derives a proof for Eq. (29).

The third step is introducing natural dimensions to the tectonic uplift history by solving Eqs. (26)–(27) for (Uj, tj) based on the scaled history (Uj*,tj*) and after constraining K and n independently. K and n could be constrained though, for example, with correlations between locally measured steepness index and erosion rates or uplift rates, following Eq. (7). Inferences of erosion and uplift rates could rely, for example, on detrital cosmogenic radionuclide concentrations (e.g., Ouimet et al., 2009; Dibiase et al., 2011; Harel et al., 2016; Hilley et al., 2019; Adams et al., 2020) or dated uplifted terraces (e.g., Whittaker and Boulton, 2012; Gallen and Fernández-Blanco, 2021).

In the following, the inversion procedure is demonstrated for both numerical data and natural data from the Dadu River basin. For the numerically based demonstration, we use a low-resolution numerical model that solves Eq. (2). The model is used to generate 10 river profiles with variable channel length and drainage area distribution with pre-chosen model parameters of n, m and K (Fig. 5). These rivers respond to the same uplift rate history, with two step increases in the uplift rate forming two knickpoints in each profile. Knickpoints do not merge over the timeframe of model application (Fig. 5a and b). In addition to the numerical diffusion inherent to the model, to artificially increase the noise in the data, the elevations are perturbed by random errors: zi^(perturbed)=zi-1+(zi+1-zi-1)rand[0,1],where rand[0,1] is a random number between 0 and 1. Inversion is applied to the data while using the known pre-chosen values of n, m and K and attempting a variable number of division points between 1–6. For each number of division points, 5000 realizations of the inversion are performed with different random positions of the division points. Figure 5c shows the minimal misfit (Eq. 28) achieved for each number of division points, indicating that the best-fit solution has two division points. Figure 5d shows the inferred history, indicating that the two-division-point inversion correctly infers the applied history.

Figure 5Inversion of numerical rivers with n=2. (a–b) River profiles and χz plots of numerically generated rivers with length that ranges between 12–21 km, ka∼2 to 1.55 and h∼0.67 to 4.27. The stream-power parameters are n=2, K=1×10-8 m−0.8 a−1 and m=0.9. The time (dt) and space (dx) steps are 100 years and 100 m, respectively. The applied uplift rate history is U0=0.05 mm a−1 prior to 2.5 Ma, U1=0.2 mm a−1 between 0.5–2.5 Ma and U2=1.0 mm a−1 between 0.0–0.5 Ma. Profiles are shown with added elevation noise. See text for details. (c) Elevation misfit, Eq. (28), as a function of the number of division points. (d) The inferred scaled (black line) and dimensional (red line) uplift history, based on two division points. Introducing correct dimensions was achieved by using the known model n and K.


6.2 Application of the inverse model to the Dadu River basin

As a second demonstration of the n≠1 inverse model, we applied it to the Dadu River basin that drains portions of the east Tibetan Plateau (Fig. 6a). The main stream (Rivers 1 and 2 in Fig. 6a) of the Dadu River basin originates from the interior of the plateau (with an elevation over 5000 m) and runs across the steep plateau margin flowing in the N–S direction. Near the city of Shimian, the main stream turns eastwards and flows into the Sichuan Basin (with elevation of ∼500 m).

Figure 6Application of the inverse model to the Dadu River basin. (a) The main trunks of the Dadu River basin (shaded light blue area) and nearby rivers. Rivers labeled 1 to 5 are used for inversion. Main faults (red and black) are based on Zhang et al. (2015, 2017). (b) The χz plots of the main trunks that are used in the inversion (thin, colored lines). The thick gray line represents the best-fit inferred χz profile with two division points. The gray arrows indicate the position of the knickpoints on the inferred χz profile. (c) Elevation misfit as a function of the number of division points. (d) The non-dimensional uplift history with the best-fit solution using two division points. (e) The inferred dimensional uplift history with n=1 (gray) and n=2 (red).

Two main tectono-geomorphic events were suggested to control the late Cenozoic erosional history of the Dadu River basin. First, a regional cooling event dated to the late Miocene was inferred based on synchronous rapid exhumation from Shimian and upstream as recorded by low-temperature thermochronology and was attributed to be a response to the regional tectonic uplift that initiated at about 9–12 Ma (e.g., Tian et al., 2015; Zhang et al., 2017; Yang et al., 2019). Second, a major capture event of the upper Dadu River that used to drain through the Anninghe and was redirected to the Yangtze River near Shimian (Clark et al., 2004) was dated to the early Pleistocene by using provenance analysis and thermal modeling (Yang et al., 2019) together with inversion of detrital apatite fission track (AFT) ages in the modern Anninghe River basin (Wang et al., 2021).

Ma et al. (2020) performed a linear inversion on all the streams of the Dadu River basin while assuming n=1 and equal segment length in the χ domain. According to their inversion results, the uplift rates were ∼0.05 mm a−1 before the middle Miocene and gradually increased to ∼0.35 mm a−1 from 12–15 Ma until the present. However, a correlation between catchment-wide denudation rates and steepness indices by Ouimet et al. (2009) indicates that the slope exponent in the region is likely >1. This means that the true history probably deviates from that inferred with the assumption of n=1 (Goren et al., 2014). We explore the long profile of the Dadu River basin through the inversion procedure proposed here with both n=1 and n>1 with the goal of identifying changes in the basin relative uplift rate history and exploring their relation to the previously inferred tectono-geomorphic history. Here, relative uplift rate refers to the uplift rate experienced by the inverted rivers relative to their local base level (Goren et al., 2014) at Shimian.

We inverted five long main trunks of the Dadu River basin, which all drain to the same base level and generate a uniform trend in the χ-elevation domain, when the concavity index m/n=0.45 (Fig. 6b). The inversion was repeated with 1–10 division points. For each number of division points, 5000 realizations of the inversion were performed with different random positions of the division points. For each inverse model run, the elevation of the modeled rivers was calculated using Eq. (29), and the elevation misfit on all data points was calculated following Eq. (28). Figure 6c shows the elevation misfit as a function of the number of division points. The minimum misfit corresponds to two division points. The non-dimensional uplift history that corresponds to the minimal misfit is presented in Fig. 6d.

To introduce natural dimensions to the uplift rate history, the slope exponent, n, and erodibility coefficient, K, need to be constrained. For that, we rely on the correlation between 10Be-derived erosion rates at tributary basins and steepness indices reported by Ouimet et al. (2009), which could be consistent with a slope exponent ranging between n=1–4 (n=2 yields the best correlation coefficient). We introduce natural dimensions to the scaled uplift rate history with two sets of parameters. The first set is n=1 and K=1.25×10-6 m0.1 a−1, and the second set is n=2 and K=4.01×10-9 m−0.8 a−1. The erodibility coefficients were inferred based on regressions through the 10Be-derived erosion rates – steepness index data of Ouimet et al. (2009), while fixing the value of n. With both n=1 and n=2, the inferred histories predict significant increases in the relative uplift rate at ∼8 and 1–2 Ma (Fig. 6e), consistent with the timing of the tectono-geometric events seen in low temperatures (Tian et al., 2015; Zhang et al., 2015; Yang et al., 2019). In particular, knickpoint ages are expected not to be older than the inferred signal age with thermochronology. While the different sets of n and K predict that the tectonic uplift rate changes at approximately the same times, the inferred relative uplift rate values are different. With n=1, the relative uplift rate increases by no more than a factor of 4 between the oldest inferred rate at the late Miocene (∼0.1 mm a−1) and the present-day rate (∼0.375 mm a−1). With n=2, the oldest relative uplift rate (no more than 0.05 mm a−1) is slower by approximately a factor of 10 with respect to recent rates (∼0.5 mm a−1). The greater change in relative uplift rate and the faster recent relative uplift rate with n=2 are both more consistent with Ouimet et al. (2009) inferred distribution of erosion rate between the higher and the lower reaches of the Dadu River basin.

7 Discussion

The analysis presented here explores river long profile evolution in response to temporal step changes in the tectonic rock uplift rate U(t) with a non-unity slope exponent, which can lead to consuming channel segments (Royden and Perron, 2013) and merging knickpoints. The approach we adopt, of resolving knickpoint kinematics in a Lagrangian frame of reference, allows us to constrain the timing of knickpoint merging and the elevation and position of knickpoints before and after merging. The finding that despite channel reach consumption knickpoint celerity depends only on the channel steepness below and above the knickpoint allows us to develop a piece-wise analytic solution that represents the evolution of knickpoints and channel long profile through time, before and after knickpoint merging.

The analysis of merging knickpoints further emphasizes a critical property of the links between tectonic and long profile evolution when n≠1. Each tectonic uplift history is associated with a single, well-defined river profile at any given time. Therefore, the forward model that we develop here could be used without any restrictions. The inverse inference, however, has a different property, whereby any particular river long profile could be generated by many tectonic uplift histories (as demonstrated by the evolution depicted in Fig. 4). If a tectonic uplift history has occurred without merging of knickpoints, our method can reconstruct this history. However, a tectonic history that results in knickpoints merging cannot be recovered using our linear inversion method. More specifically, when our inverse approach is applied to a river long profile, the outcome will be the one history for which all knickpoints are preserved, although this inferred outcome might not be the real history that shaped the profile. While this inverse approach is highly restrictive, it finds the correct solution when only a single knickpoint group existed in the data. We further suggest that when a small number of knickpoint groups is identified in the data, the solution of this simple inverse model could still be highly informative as a preliminary guess for the tectonic uplift rate history that shaped the fluvial landscape.

7.1 Assumptions underlying the analytic derivation and models

A basic assumption underlying the analytic derivation and particularly the forward and inverse models is that the channel system experiences space-invariant uplift (also consistent with a base-level fall). This assumption, which is commonly referred to as block uplift conditions, is more likely to hold over discrete, well-defined tectonic domains with relatively little internal complexity rather than over large length scales (Goren et al., 2022). However, larger domains could also experience spatial uniformity in the uplift history. One way to test for this uniformity is to explore the χz space of the rivers and tributaries. If they all collapse along a single trend, then they likely represent channels responding to block uplift conditions. Figure 6b demonstrates this for the Dadu River basin tributaries. Despite the length scale of hundreds of kilometers for the Dadu basin, the five tributaries that we analyzed for the relative uplift rate history collapse on a single trend in the χz domain, which we interpret to support the block uplift conditions for these tributaries. Generally, when applying the inverse model over a branching network of channels, then the inferred uplift rate history smoothens local variabilities that may exist in the true uplift rate signal. The inverse solution may then be regarded as an “average” from which local histories slightly deviate.

The analytic derivation lacks a process-based perspective of knickpoint migration and instead relies on a simplified stream power parameterization of knickpoint dynamics. Consequently, a second major assumption, with specific impact on the inverse model, is that the natural knickpoints analyzed for changes in the tectonic uplift history are indeed slope-break knickpoints, which were formed following a change in the tectonic uplift rate. Knickpoints may also form by autogenic processes (e.g., Scheingross and Lamb, 2017) or due to spatial changes in the uplift rate (Wobus et al., 2006), rock erodibility (Kirby and Whipple, 2012) or local hydrologic conditions (Hamawi et al., 2022). However, when analyzing a branching channel network, it is relatively easy to distinguish between migrating slope-break knickpoints which were formed due to a regional uplift rate change and locally controlled knickpoints. The migrating knickpoints share approximately similar χ and elevation values across tributaries and basins (under a block uplift assumption), while the latter do not (e.g., Hamawi et al., 2022).

The current derivation focuses on particular combinations of tectonic uplift histories and slope exponent with either increasing U and n>1 or decreasing U and n<1. While these combinations may appear restrictive, the former combination likely describes many (if not the majority) of the dynamic high-elevation landscapes that are dissected by bedrock rivers. Recent studies have found that such landscapes represent rejuvenated response to recent faster uplift rate (e.g., Whittaker and Boulton, 2012; Harkins et al., 2007; Ouimet et al., 2009; Gallen and Fernández-Blanco, 2021). These landscapes are further characterized by convex upward knickpoints, pointing at n≥1. This is in a general agreement with the recent global compilation by Harel et al. (2016), who argued that n>1 characterizes most fluvial drainages.

7.2 Future work

When U increases and n1 or U decreases and n>1, stretched zones form along the river long profile that contain no information about tectonic uplift history. Instead, they represent self-adjusting fluvial dynamics. Royden and Perron (2013) derived an analytic solution for the channel profile along a stretched zone. Future studies that combine solutions for stretched zones together with the Lagrangian perspective developed here for consuming channel reaches and merging knickpoints could promote the derivation of efficient forward and possibly inverse models that allow for a general uplift rate history.

With n≠1, fluvial dynamics could lead to consuming channel reaches and eventually merging knickpoints. While the inverse model cannot resolve merging knickpoint dynamics, the forward model resolves knickpoint evolution through and beyond merging. This means that the forward model can be used to test any tectonic scenario, including those that lead to knickpoint merging, and identify those scenarios that are consistent with the remaining knickpoints and steepness indices observed in any particular fluvial landscape.

To elucidate this idea, we revisit the simple case discussed in Sect. 4, with n>1 and two step increases in U, (U2>U1>U0) leading to an older kp1 that formed at time t=0 and younger kp2 knickpoints that formed at time t=T1. This scenario could result in knickpoints merging and complete consumption of the middle channel reach whose steepness index was ks_1=(U1/K)1/n. Despite this complete consumption, some constraints could be placed on the “lost” uplift rate based on the conjecture that the two knickpoints merged and therefore χ(kp2)>χ(kp1) at time t=T, where T>T1. Based on Eq. (15), this condition can be expressed as

(30) T 1 / T < 1 - γ 2 _ 1 n 1 - γ 2 _ 1 1 - γ 0 _ 1 1 - γ 0 _ 1 n - 1 , with γ 2 _ 1 = k s _ 2 / k s _ 1 .

The inequality in Eq. (30), describing the relation between T, T1 and U1 (through γ0_1 and γ2_1), could be used to rule out potential lost histories that do not obey χ(kp2)>χ(kp1).

More generally, analytic solutions of river long profile evolution can significantly expedite forward and inverse tectonic–fluvial landscape evolution models. However, so far, analytic solutions were used in such models only under the n=1 assumption (Pritchard et al., 2009; Fox et al., 2014; Goren et al., 2014, 2022; Rudge et al., 2015; Steer et al., 2021). The simple analytic derivation that we present here can expand the domain of parameters for which analytic solutions are used in such models, by including new geomorphic scenarios with n≠1. For example, inverse models that are based on Bayesian statistics (e.g., Fox et al., 2015; Gallen and Fernández-Blanco, 2021), which have gained recent popularity, could become significantly more efficient and accurate when the forward model is represented with an analytic solution. This presents a great opportunity for future studies to combine our newly derived forward model as part of a Bayesian inversion of river long profiles.

8 Conclusions

We develop an analytic slope-break knickpoint retreat model under the assumption of space-invariant uplift rate. The model is based on a Lagrangian frame of reference and can deal with both convex- (n>1, monotonic step increase in U) and concave-up (n<1, decreasing U) knickpoints. Knickpoint celerity depends on the stream-power model slope exponent, n, and the ratio of channel steepness indices above and below the knickpoint. Consequently, for the conditions we study here, the celerity of newer knickpoints is greater than that of the older knickpoints that propagate along the same channel, and knickpoints could merge. We derive a mathematical formulation to determine the preservation duration of knickpoints before merging. We further derive an analytical forward model that solves for the evolution of the channel profile before and after knickpoint merging. Finally, assuming that all the knickpoints are preserved, we develop a linear inverse model to retrieve the tectonic uplift history from the river long profiles. The forward and inverse models are novel in their ability to treat cases in which n≠1. Applying the inverse model with n=2 to the Dadu River basin, the model inferred a relative tectonic uplift rate history that is consistent with the exhumation history recorded by low-temperature thermochronology. The analytic derivation presented here could be readily incorporated in forward and inverse tectonic–fluvial landscape evolution models achieving accurate and efficient solutions.

Appendix A: A mathematical demonstration of knickpoint merging

In this section, we show that two knickpoints formed with n>1 and step increases in U must eventually merge. The two knickpoints are denoted by kp1, which was formed by an uplift rate increase from U0 to U1, and kp2 formed by an increase from U1 to U2 (U2>U1>U0). The celerity of the two knickpoints is expressed by Eq. (13):

(A1) v H _ kp 1 = k s _ 1 n - 1 ( 1 - γ 0 _ 1 n ) ( 1 - γ 0 _ 1 ) KA ( kp 1 ) m / n and v H _ kp 2 = k s _ 2 n - 1 ( 1 - γ 1 _ 2 n ) ( 1 - γ 1 _ 2 ) KA ( kp 2 ) m / n .

Since kp2 is located below kp1, A(kp2) is larger than A(kp1). Next, it is left to show that ks_2n-1(1-γ1_2n)(1-γ1_2)>ks_1n-1(1-γ0_1n)(1-γ0_1). We define a variable:

(A2) f = k s _ 1 n - 1 ( 1 - γ 0 _ 1 n ) ( 1 - γ 0 _ 1 ) / k s _ 2 n - 1 ( 1 - γ 1 _ 2 n ) ( 1 - γ 1 _ 2 ) = 1 γ 1 _ 2 1 - n ( 1 - γ 0 _ 1 n ) / ( 1 - γ 0 _ 1 ) ( 1 - γ 1 _ 2 n ) / ( 1 - γ 1 _ 2 ) .

Because n>1, we can re-write n=α/β, where α>β>1 and α and β are both integers. Thus,

(A3) f = 1 γ 1 _ 2 1 - α / β ( 1 - ( γ 0 _ 1 1 / β ) α ) / ( 1 - ( γ 0 _ 1 1 / β ) β ) ( 1 - ( γ 1 _ 2 1 / β ) α ) / ( 1 - ( γ 1 _ 2 1 / β ) β ) = f nume f deno ,

where fnume and fdeno are the numerator and denominator of f, respectively. We use the method of polynomial division:

(A4) 1 - ( γ 0 _ 1 1 / β ) α = ( 1 - γ 0 _ 1 1 / β ) ( ( γ 0 _ 1 1 / β ) α - 1 + + ( γ 0 _ 1 1 / β ) β + + ( γ 0 _ 1 1 / β ) 0 ) 1 - ( γ 0 _ 1 1 / β ) β = ( 1 - γ 0 _ 1 1 / β ) ( ( γ 0 _ 1 1 / β ) β - 1 + ( γ 0 _ 1 1 / β ) β - 2 + + ( γ 0 _ 1 1 / β ) 0 ) .

Assigning Eq. (A4) into fnume, we can derive

(A5) f nume = ( γ 0 _ 1 1 / β ) α - 1 + + ( γ 0 _ 1 1 / β ) β + + ( γ 0 _ 1 1 / β ) 0 ( γ 0 _ 1 1 / β ) β - 1 + ( γ 0 _ 1 1 / β ) β - 2 + + ( γ 0 _ 1 1 / β ) 0 = ( γ 0 _ 1 1 / β ) α - 1 + ( γ 0 _ 1 1 / β ) α - 2 + + ( γ 0 _ 1 1 / β ) β ( γ 0 _ 1 1 / β ) β - 1 + ( γ 0 _ 1 1 / β ) β - 2 + + ( γ 0 _ 1 1 / β ) 0 + 1 = ( γ 0 _ 1 1 / β ) α - 1 - β + ( γ 0 _ 1 1 / β ) α - 2 - β + + ( γ 0 _ 1 1 / β ) β - β ( γ 0 _ 1 1 / β ) - 1 + ( γ 0 _ 1 1 / β ) - 2 + + ( γ 0 _ 1 1 / β ) - β + 1 .

Because (γ0_11/β)α-1-β+(γ0_11/β)α-2-β++(γ0_11/β)β-β <1+1++1=α-β and (γ0_11/β)-1+(γ0_11/β)-2++ (γ0_11/β)-β>1+1++1=β, we can derive

(A6) f nume < α - β β + 1 .

Again, we use polynomial division.

(A7) 1 - ( γ 1 _ 2 1 / β ) α = ( 1 - γ 1 _ 2 1 / β ) ( ( γ 1 _ 2 1 / β ) α - 1 + ( γ 1 _ 2 1 / β ) α - 2 + + ( γ 1 _ 2 1 / β ) 0 ) 1 - ( γ 1 _ 2 1 / β ) β = ( 1 - γ 1 _ 2 1 / β ) ( ( γ 1 _ 2 1 / β ) β - 1 + ( γ 1 _ 2 1 / β ) β - 2 + + ( γ 1 _ 2 1 / β ) 0 )

Assigning Eq. (A7) into fdeno, we derived


Because (γ1_21/β)-1+(γ1_21/β)-2++(γ1_21/β)β-α>1+1++1=α-β and (γ1_21/β)β-1+(γ1_21/β)β-2++(γ1_21/β)0<1+1++1=β, we derived

(A10) f deno > α - β β + 1 .

Assigning Eqs. (A6) and (A10) into (A3), we can derive

(A11) f = f nume f deno < ( α - β β + 1 ) / ( α - β β + 1 ) = 1 .

Thus, vH_kp1<vH_kp2, kp2 always migrates faster than kp1 and given sufficient channel length the two knickpoints will merge. The time of merging is given by Eq. (16).

Appendix B: Calculating the predicted elevations based on the scaled relative uplift rate history

Equations (26)–(27) express the scaled uplift rate history as a series of values (U*, τ*) describing the scaled age of a knickpoint τj* and the non-dimensional uplift rate, Uj*, that operated between scaled time τj* and τj-1*, where τ* is identified with the χ axis and τ0*=0 corresponds to the outlet. In this appendix, we prove Eq. (29) and show how the scaled uplift rate history could be used to calculate the forward model and predict knickpoint elevations zj̃ and the elevations of other pixels zĩ. These predicted elevations are used in the misfit calculation, Eq. (28), to evaluate the inversion results. Equation (29) is proved by induction.

First, we prove the base case with a single knickpoint. For this case, Eq. (21) predicts the knickpoint elevation z1 to be

(B1) z 1 = 0 t 1 U 1 d t + 1 - γ 1 n 1 - γ 1 - 1 U 1 t 1 = 1 - γ 1 n 1 - γ 1 U 1 t 1 ,

where γ1=ks_2/ks_1, t1 is the age of the knickpoint, and U1 is the uplift rate that generated the knickpoint. According to Eqs. (26)–(27), t1 and U1 are defined as


Substituting Eqs. (B2)–(B3) into (B1), we get

(B4) z 1 = 1 - γ 1 n 1 - γ 1 ( U 1 * A 0 m / n ) n K t 1 * 1 KA 0 m / n k s _ 1 1 - γ 1 k s _ 1 n 1 - γ 1 n = ( U 1 * A 0 m / n ) n t 1 * 1 A 0 m / n k s _ 1 k s _ 1 n .

Then, using the definition ks_1=U1*A0-m/n (Eq. 27), Eq. (B4) can be simplified to

(B5) z 1 = U 1 * t 1 * = U 1 * ( t 1 * - t 0 * ) ,

where t0*=χ0=0 (basin outlet).

Then, assuming that Eq. (29) holds for knickpoint j with elevation zj, we prove the induction step for knickpoint zj+1. Noting that zj+1=zj+(zj+1-zj), we evaluate the elevation difference between the two knickpoints j and j+1 following Eq. (21) as

(B6) z j + 1 - z j = t j t j + 1 U j + 1 d t + 1 - γ j + 1 n 1 - γ j + 1 - 1 U j + 1 t j + 1 - 1 - γ j n 1 - γ j - 1 U j t j .

Arranging Eq. (B6),

(B7) z j + 1 - z j = - U j + 1 t j + 1 - γ j + 1 n 1 - γ j + 1 U j + 1 t j + 1 - 1 - γ j n 1 - γ j U j t j + U j t j .

Based on Eqs. (26)–(27) and the scaling ks_j=Uj*A0-m/n, we define


where γj=ks_j+1/ks_j. Substituting Eqs. (B8)–(B10) into (B7),

(B11) z j + 1 - z j = - t j * U j * 1 - γ j 1 - γ j n γ j n + t j + 1 * U j + 1 * - t j * U j * + t j * U j * 1 - γ j 1 - γ j n .

Rearranging Eq. (B11), the knickpoint elevation difference can be written as

(B12) z j + 1 - z j = t j + 1 * U j + 1 * + t j * U j * 1 - γ j 1 - γ j n - 1 - γ j 1 - γ j n γ j n - 1 = t j + 1 * U j + 1 * - γ j t j * U j *

since γjUj*=ks_j+1ks_jks_jA0m/n=Uj+1* Eq. (B12) becomes

(B13) z j + 1 - z j = t j + 1 * U j + 1 * - t j * U j + 1 * = U j + 1 * ( t j + 1 * - t j * )


(B14) z j + 1 = z j + U j + 1 * ( t j + 1 * - t j * ) .

Therefore, based on the induction base (Eq. B5) and step (Eq. B14), knickpoint elevations can be expressed based on the scaled uplift rate history as

(B15) z j = 0 t j * U * ( t * ) d t * = a = 1 j U a * ( t a * - t a - 1 * ) .

For any pixel, i, between the knickpoint j and j+1, its elevation can be predicted based on the scaled uplift rate history as

(B16) z i ̃ = z j + k s _ j + 1 A 0 - m / n ( χ i - χ j ) = 0 t j * U * ( t * ) d t * + U j + 1 * ( χ i - χ j ) = 0 t i * = χ i U * ( t * ) d t * ,

showing that Eq. (29) holds for any pixel in the landscape.

Code and data availability

This study has no complex codes or data sharing issues. All the figures can be reproduced by solving the related equations. The DEM (digital elevation model) data used for river profile inversion (Fig. 6) are from the 90 m Shuttle Radar Topography Mission (SRTM) DEM (downloaded from; SRTM Data, 2022; Farr et al., 2007).


The supplement related to this article is available online at:

Author contributions

YW derived the theory and developed the forward and inverse analytic models with input from LG. YW designed the application of the inverse model to the Dadu River basin with input from LG. YW and LG wrote the manuscript. LG wrote the 1-D numerical code. DZ and HZ provided valuable suggestions and made some revisions.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank George E. Hilley for revising the earlier version of this paper many times with great patience and stimulating our inspiration in using the method of characteristics to solve the equation. We thank Fiona Clubb, Eitan Shelef and Sean F. Gallen for constructive instructions on an earlier version of this paper. Constructive suggestions from Philippe Steer and the anonymous reviewer and Simon Mudd (associate editor) helped to improve our paper.

Financial support

This research has been supported by the National Science Foundation of China (grant no. 41802227) and by the Israel Science Foundation (grant no. 526/19).

Review statement

This paper was edited by Simon Mudd and reviewed by Philippe Steer and one anonymous referee.


Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, A. M., and Hodges, K. V.: Climate controls on erosion in tectonically active landscapes, Sci. Adv., 6, eaaz3166,, 2020. 

Akaike, H.: A new look at the statistical model identification, IEEE Trans. Autom. Control, 19, 716–723,, 1974. 

Anthony, D. M. and Granger, D. E.: An empirical stream power formulation for knickpoint retreat in Appalachian Plateau fluviokarst, J. Hydrol., 343, 117–126, 2007. 

Berlin, M. M. and Anderson, R. S.: Modeling of knickpoint retreat on the Roan Plateau, western Colorado, J. Geophys. Res., 112, F03S06,, 2007. 

Bookhagen, B. and Burbank, D. W.: Toward a complete Himalayan hydrological budget: spatiotemporal distribution of snowmelt and rainfall and their impact on river discharge, J. Geophys. Res., 115, F03019,, 2010. 

Castillo, M., Ferrari, L., and Munoz-Salinas, E.: Knickpoint retreat and landscape evolution of the Amatlan de Canas half-graben (northern sector of Jalisco Block, western Mexico), J. S. A. Earth Sci., 77, 108–122, 2017. 

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

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

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The shuttle radar topography mission, Rev. Geophys., 45, RG2004,, 2007. 

Ferrier, K. L., Huppert, K. L., and Perron, J. T.: Climatic control of bedrock river incision, Nature, 496, 206–209,, 2013. 

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

Fox, M., Goren, L., May, D. A., and Willett, S. D.: Inversion of fluvial channels for paleorock uplift rates in Taiwan, J. Geophys. Res.-Earth, 119, 1853–1875,, 2014. 

Fox, M., Bodin, T., and Shuster, D. L.: Abrupt changes in the rate of Andean Plateau uplift from reversible jump Markov Chain Monte Carlo inversion of river profiles, Geomorphology, 238, 1–14, 2015. 

Gailleton, B., Mudd, S. M., Clubb, F. J., Peifer, D., and Hurst, M. D.: A segmentation approach for the reproducible extraction and quantification of knickpoints from river long profiles, Earth Surf. Dynam., 7, 211–230,, 2019. 

Gailleton, B., Mudd, S. M., Clubb, F. J., Grieve, S. W. D., and Hurst, M. D.: Impact of changing concavity indices on channel steepness and divide migration metrics, J. Geophys. Res.-Earth, 126, e2020JF006060,, 2021. 

Gallen, S. F. and Fernández-Blanco, D.: A new data-driven Bayesian inversion of fluvial topography clarifies the tectonic history of the corinth rift and reveals a channel steepness threshold, J. Geophys. Res.-Earth, 126, e2020JF005651,, 2021. 

Gallen, S. F. and Wegmann, K. W.: River profile response to normal fault growth and linkage: an example from the Hellenic forearc of south-central Crete, Greece, Earth Surf. Dynam., 5, 161–186,, 2017. 

Gasparini, N. M. and Brandon, M. T.: A generalized power law approximation for fluvial incision of bedrock channels, J. Geophys. Res., 116, F02020,, 2011. 

Goren, L., Fox, M., and Willett, S. D.: Tectonics from fluvial topography using formal linear inversion: theory and applications to the Inyo Mountains, California, J. Geophys. Res.-Earth, 119, 1651–1681, 2014. 

Goren, L., Fox, M., and Willett, S. D.: Linear Inversion of Fluvial Long Profiles to Infer Tectonic Uplift Histories, in: Treatise on Geomorphology (Second Edition), edited by: Shroder, J. F., Academic Press, 225–248,, 2022. 

Hack, J. T.: Studies of Longitudinal Stream Profiles in Virginia and Maryland, U.S. Geological Survey Professional Paper, 294, 45–97, 1957. 

Hack, J. T.: Stream profile analysis and stream-gradient index, J. Res. US Geol. Surv., 1, 421–429, 1973. 

Hamawi, M., Goren, L., Mushkin, A., and Levi, T.: Rectangular drainage pattern evolution controlled by pipe cave collapse along clastic dikes, the Dead Sea Basin, Israel, Earth Surf. Proc. Land., 47, 936–954,, 2022 

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

Harkins, N., Kirby, E., Heimsath, A., Robinson, R., and Reiser, U.: Transient fluvial incision in the headwaters of the Yellow River, northeastern Tibet, China, J. Geophys. Res.-Earth, 112, F03S04,, 2007. 

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

Hergarten, S., Robl, J., and Stüwe, K.: Tectonic geomorphology at small catchment sizes – extensions of the stream-power approach and the χ method, Earth Surf. Dynam., 4, 1–9,, 2016. 

Hilley, G. E., Porder, S., Aron, F., Baden, C. W., Johnstone, S. A., Liu, F., Sare, R., Steelquist, A., and Young, H. H.: Earth's topographic relief potentially limited by an upper bound on channel steepness, Nat. Geosci., 12, 828–832, 2019. 

Howard, A.: A detachment-limited model of drainage basin evolution, Water Resour. Res., 30, 2261–2285, 1994. 

Howard, A. D. and Kerby, G.: Channel changes in badlands, Geol. Soc. Am. Bull., 94, 739–752, 1983. 

Kirby, E. and Whipple, K. X.: Expression of active tectonics in erosional landscapes, J. Struct. Geol., 44, 54–75,, 2012. 

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Process. Landf., 39, 38–61, 2014. 

Lavé, J. and Avouac, J. P.: Fluvial incision and tectonic uplift across the Himalayas of central Nepal, J. Geophys. Res.-Sol. Ea., 106, 26561–26591, 2001. 

Luke, J. C.: Mathematical models for landform evolution, J. Geophys. Res., 77, 2460–2464, 1972. 

Ma, Z., Zhang, H., Wang, Y., Tao, Y., and Li, X.: Inversion of Dadu River bedrock channels for the late Cenozoic uplift history of the eastern Tibetan Plateau, Geophys. Res. Lett., 47, e2019GL086882,, 2020. 

Mitchell, N. A. and Yanites, B. J.: Spatially Variable Increase in Rock Uplift in the Northern U.S. Cordillera Recorded in the Distribution of River Knickpoints and Incision Depths, J. Geophys. Res.-Earth, 124, 1238–1260,, 2019. 

Morisawa, M. E.: Quantitative Geomorphology of Some Watersheds in the Appalachian Plateau, Geol. Soc. Am. Bull., 73, 1025–1046,[1025:QGOSWI]2.0.CO;2, 1962. 

Mudd, S. M., Attal, M., Milodowski, D. T., Grieve, S. W., and Valters, D. A.: A statistical framework to quantify spatial variation in channel gradients using the integral method of channel profile analysis, J. Geophys. Res.-Earth, 119, 138–152, 2014. 

Niemann, J. D., Gasparini, N. M., Tucker, G. E., and Bras, R. L.: A quantitative evaluation of Playfair's law and its use in testing long-term stream erosion models, Earth Surf. Process. Landf., 26, 1317–1332, 2001. 

Oskin, M. E. and Burbank, D.: Transient landscape evolution of basement-cored uplifts: example of the Kyrgyz Range, Tian Shan, J. Geophys. Res., 112, F03S03,, 2007. 

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

Paul, J. D., Roberts, G. G., and White, N.: The African landscape through space and time, Tectonics, 33, 898–935,, 2014. 

Perron, J. T. and Royden, L.: An integral approach to Bedrock River profile analysis, Earth Surf. Process. Landf., 38, 570–576,, 2013. 

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

Roberts, G. G. and White, D.: Estimating uplift rate histories from river profiles using African examples, J. Geophys. Res., 115, B02406,, 2010. 

Rosenbloom, N. A. and Anderson, R. S.: Hillslope and channel evolution in a marine terraced landscape, Santa Cruz, California, J. Geophys. Res., 99, 14013–14029, 1994. 

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

Rudge, J. F., Roberts, G. G., White, N. J., and Richardson, C. N.: Uplift histories of Africa and Australia from linear inverse modeling of drainage inventories, J. Geophys. Res.-Earth, 120, 894–914, 2015. 

Scheingross, J. S. and Lamb, M. P.: A mechanistic model of waterfall plunge pool erosion into bedrock, J. Geophys. Res.-Earth, 122, 2079–2104,, 2017. 

Schwanghart, W. and Scherler, D.: Divide mobility controls knickpoint migration on the Roan Plateau (Colorado, USA), Geology, 48, 698–702,, 2020. 

Shelef, E., Haviv, I., and Goren, L.: A potential link between waterfall recession rate and bedrock channel concavity, J. Geophys. Res.-Earth, 123, 905–923,, 2018. 

Snyder, N. P., Whipple, K. X., Tucker, G. E., and Merritts, D.: Landscape response to tectonic forcing: DEM analysis of stream profiles in the Mendocino Triple Junction region, northern California, Geol. Soc. Am. Bull., 112, 1250–1263, 2000. 

SRTM Data: 90 m Shuttle Radar Topography Mission (SRTM) DEM, (last access: 25 March 2022), 2022. 

Steer, P.: Short communication: Analytical models for 2D landscape evolution, Earth Surf. Dynam., 9, 1239–1250,, 2021. 

Tian, Y., Kohn, B. P., Hu, S., and Gleadow, A. J. W.: Synchronous fluvial response to surface uplift in the eastern Tibetan Plateau: Implications for crustal dynamics, Geophys. Res. Lett., 42, 29–35,, 2015. 

Venditti, J. G., Li, G., Deal, E., Dingle, E., and Church, M.: Struggles with stream power: Connecting theory across scales, Geomorphology, 366, 106817,, 2019. 

Wang, Y., Zhang, H., Zheng, D., Dassow, W. V., Zhang, Z., Yu, J., and Pang, J.: How a stationary knickpoint is sustained: new insights into the formation of the deep Yarlung Tsangpo Gorge, Geomorphology 285, 28–43, 2017.  

Wang, Y., Liu, C., Zheng, D., Zhang, H., Yu, J., Pang, J., Li, C., and Hao, Y.: Multistage exhumation in the catchment of the Anninghe River in the SE Tibetan Plateau: Insights from both detrital thermochronology and topographic analysis, Geophys. Res. Lett., 48, e2021GL092587,, 2021. 

Weissel, J. K. and Seidl, M. A.: Inland propagation of erosional escarpments and river profile evolution across the southeast Aus tralian passive continental margin, Geophys. Monogr., 107, 189–206, 1998. 

Whipple, K. X. and Tucker, G. E.: Dynamics of the stream power river incision model: Implications for height limits of mountain ranges, landscape response timescales and research needs, J. Geophys. Res., 104, 17661–17674, 1999. 

Whipple, K. X. and Tucker, G. E.: Implications of sediment-flux-dependent river incision models for landscape evolution, J. Geophys. Res.-Sol. Ea., 107, 2039,, 2002. 

Whipple, K. X., Hancock, G. S., and Anderson, R. S.: River incision into bedrock: Mechanics and relative efficacy of plucking, abrasion, and cavitation, Geol. Soc. Am. Bull., 112, 490–503, 2000. 

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

Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: procedures, promise, and pitfalls, in: Tectonics, Climate, and Landscape Evolution, edited by: Willett, S. D., Hovius, N., Brandon, M. T., and Fisher, D. M.: Geological Society of America, Special Paper 398, Penrose Conference Series, 55–74,, 2006. 

Yang, R., Suhail, H. A., Gourbet, L., Willett, S. D., Fellin, M. G., Lin, X., Gong, J. F., Wei, X. C., Maden, C., Jiao, R. H., and Chen, H. L.: Early Pleistocene drainage pattern changes in Eastern Tibet: Constraints from provenance analysis, thermochronometry, and numerical modeling, Earth Planet. Sc. Lett., 531, 115955,, 2019. 

Zhang, Y. Z., Replumaz, A., Wang, G. C., Leloup, P. H., Gautheron, C., Bernet, M., van der Beek, P., Paquette, L. J., Wang, A., Zhang, K. X., Chevalier, M. L., and Li, H. B.: Timing and rate of exhumation along the Litang fault system, implication for fault reorganization in Southeast Tibet, Tectonics, 34, 1219–1243,, 2015. 

Zhang, Y. Z., Replumaz, A., Leloup, P. H., Wang, G.-C., Bernet, M., van der Beek, P., Paquette, L. J., and Chevalier, M. L.: Cooling history of the Gongga batholith: Implications for the Xianshuihe Fault and Miocene kinematics of SE Tibet, Earth Planet. Sc. Lett., 465, 1–15,, 2017. 

Short summary
Abrupt changes in tectonic uplift rates induce sharp changes in river profile, called knickpoints. When river erosion depends non-linearly on slope, we develop an analytic model for knickpoint velocity and find the condition of knickpoint merging. Then we develop analytic models that represent the two-directional link between tectonic changes and river profile evolution. The derivation provides new understanding on the links between tectonic changes and river profile evolution.