**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

^{1},

^{2},

^{3},

^{1}

**Yizhou Wang et al.**Yizhou Wang Liran Goren Dewen Zheng and Huiping Zhang

^{1},

^{2},

^{3},

^{1}

^{1}State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China^{2}Department of Earth and Environmental Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel^{3}Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, China

^{1}State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China^{2}Department of Earth and Environmental Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel^{3}Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, China

**Correspondence**: Yizhou Wang (wangyizhou2016@outlook.com)

**Correspondence**: Yizhou Wang (wangyizhou2016@outlook.com)

Received: 13 Dec 2021 – Discussion started: 21 Dec 2021 – Revised: 01 Jul 2022 – Accepted: 19 Jul 2022 – Published: 11 Aug 2022

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.

- Article
(3312 KB) -
Supplement
(653 KB) - BibTeX
- EndNote

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=\partial z/\partial x$, $L/L$) and upstream drainage area (*A*, *L*^{2}), a
proxy for both discharge and channel width (Howard and Kerby, 1983):

where *x* (*L*) denotes a spatial coordinate along the channel and *t* (*T*) is time.
The channel erodibility, *K* (${L}^{\mathrm{1}\text{\u2013}{\mathrm{2}}^{m}}/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:

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.

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:

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,
$\mathit{\theta}=m/n$ and ${k}_{\mathrm{s}}=(E/K{)}^{\mathrm{1}/n}$ (${L}^{\mathrm{2}m/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 *z*_{b} is the base-level elevation, and the area scale factor
*A*_{0} (*L*^{2}) 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 *A*_{0}=1 *L*^{2}, the slope of the *χ*−*z* plot becomes channel steepness
index, *k*_{s}.

Under steady-state conditions, when $\mathrm{d}z/\mathrm{d}t=\mathrm{0}$ and *E*=*U*, the SPIM steepness
index becomes a function of the tectonic uplift rate:

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, *U*_{p} (Niemann et al., 2001; Goren et al.,
2014). In this case, its steepness index could be expressed as

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 1a shows the predicted channel profile evolution following a step
increase in the rock uplift rate from *U*_{0} to *U*_{1} and *n*>1. The figure
emphasizes that below and above the knickpoint the channel segments erode
at rates that correspond to the new (*U*_{1}) and old (*U*_{0}) uplift rates,
respectively, and their corresponding steepness indices are
${k}_{s\mathit{\_}\mathrm{1}}=({U}_{\mathrm{1}}/K{)}^{\mathrm{1}/n}$ and ${k}_{s\mathit{\_}\mathrm{0}}=({U}_{\mathrm{0}}/K{)}^{\mathrm{1}/n}$. Figure 1b shows the linearized channel segments near
the knickpoint. The river profile varies from *z*_{t} to
*z*_{t+dt} during time step d*t*, 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

where *v*_{H} is the horizontal retreat velocity for the knickpoint
(hereafter, knickpoint celerity). Figure 1b shows that

where DB is a function of the difference between the present uplift rate
(*U*_{1}) and previous river incision rate, *U*_{0}:

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

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

which resembles the derivation of Whipple and Tucker (1999). Assigning
Eqs. (1) and (6)–(7) into (12), *v*_{H} can be re-written as

where ${\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}={k}_{\mathrm{s}\mathrm{\_}\mathrm{0}}/{k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}$. Accordingly, the fluvial response time
(Whipple and Tucker, 1999) of the knickpoint, *τ*(*x*_{p}), is
expressed as

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 *x*_{p}.
Alternatively, *τ*(*x*_{p}) 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,
*v*_{H} and *τ*(*x*_{p}) 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).

When more than a single knickpoint propagates upstream a channel profile
and *n*≠1, the sensitivity of knickpoint celerity to *k*_{s} 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: kp_{1} formed
when *U*_{0} changed to *U*_{1} and kp_{2} formed when
*U*_{1} changed to *U*_{2} (${U}_{\mathrm{2}}>{U}_{\mathrm{1}}>{U}_{\mathrm{0}}$); then the
celerity of knickpoint kp_{2} is larger than that of
kp_{1}, 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, kp_{2} can eventually
reach kp_{1}, 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 kp_{1} formed at time *t*=0 and
that kp_{2} formed at time *t*=*T*_{1}, Eq. (14) is used to express the *χ*
values of the two knickpoints at any time *T*>*T*_{1} as

where ${\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}={k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}/{k}_{\mathrm{s}\mathrm{\_}\mathrm{2}}$. Knickpoint merging occur at time
*T*_{m} when *χ*(kp_{1})=*χ*(kp_{2}).
The ratio ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ is expressed as

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 ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$, 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 ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ 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 ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$. To explore the
dependency of ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ 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 ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ 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.

For the case of concave-up consuming knickpoints (*n**<*1 and decreasing
*U*), Fig. 3a shows that the ratio ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ (when *γ*_{1_2}=*γ*_{0_1}) increases
with increasing *n*, and for a constant *n*, a higher steepness ratio leads to a
higher ${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ ratio. This means that a lower uplift rate, *U*_{2}
(with a lower steepness index below knickpoint kp_{2}), leads to a shorter
time to knickpoint merging *T*_{m}. 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
${T}_{\mathrm{m}}/{T}_{\mathrm{1}}$ is more sensitive to *γ*_{0_1} than
to *γ*_{1_2}, indicating that the preservation time
of kp_{1} is more sensitive to its own celerity than to that of the
younger knickpoint.

We note that when *n*=1, *χ*(kp_{1})>*χ*(kp_{2}) 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 *k*_{s_0} and *k*_{s_2},
respectively. Based on Eq. (13), the instantaneous merged knickpoint
celerity becomes

where ${\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{2}}={k}_{\mathrm{s}\mathrm{\_}\mathrm{0}}/{k}_{\mathrm{s}\mathrm{\_}\mathrm{2}}$. The channel reach between the two
knickpoints is fully consumed, and the channel profile holds no record of
*U*_{1}. 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
*U*_{0} to *U*_{2}.

The elevation change of slope-break knickpoint, $z\left(t,x\right)=z\left[t,x={x}_{\mathrm{p}}\left(t\right)\right]$, formed by a step increase in the uplift rate from
*U*_{0} to *U*_{1}, can be expressed as

where $\frac{\mathrm{d}x}{\mathrm{d}t}=\frac{\mathrm{d}{x}_{\mathrm{p}}\left(t\right)}{\mathrm{d}t}={v}_{\mathrm{H}}$ is the knickpoint celerity. Combining Eqs. (2), (13) and (18) yields

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

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

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\left(t{x}_{\mathrm{p}}\left(t\right)\right)={\int}_{\mathrm{0}}^{t}U\left({t}^{\prime}\right)\mathrm{d}{t}^{\prime}$.

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, kp_{1} and
kp_{2}, generated at times *t*=0 and *t*=*T*_{1}, respectively, by
two step increases in *U*, ${U}_{\mathrm{2}}>{U}_{\mathrm{1}}>{U}_{\mathrm{0}}$ and *n*>1. The time of merging,
*T*_{m}, measured with respect to *T*_{1} is constrained by
Eq. (16). For any time $t<{T}_{\mathrm{m}}+{T}_{\mathrm{1}}$, the
elevations of kp_{1} and kp_{2} are predicted by Eq. (21), and when
assigning the knickpoint ages *t*= *τ*(*x*_{p1}) and
$t-{T}_{\mathrm{1}}=\phantom{\rule{0.125em}{0ex}}\mathit{\tau}\left({x}_{\mathrm{p}\mathrm{2}}\right)$, which corresponds to the
time since the change in *U*(*t*) that generated the knickpoints. Upon merging,
for $t>{T}_{\mathrm{m}}+{T}_{\mathrm{1}}$, the elevation of the merged
knickpoint, *z*_{kp_12}, with respect to the formation
time of kp_{1} (*t*=0) can be expressed as follows.

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

where again *t*= *τ*(*x*_{p}) is the knickpoint age. After
merging, for $t>{T}_{\mathrm{m}}+{T}_{\mathrm{1}}$

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.

## 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, $\left({z}_{\mathrm{1}},{\mathit{\chi}}_{\mathrm{1}}\right),\phantom{\rule{0.125em}{0ex}}\left({z}_{\mathrm{2}},{\mathit{\chi}}_{\mathrm{2}}\right),\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\left({z}_{q-\mathrm{1}},{\mathit{\chi}}_{q-\mathrm{1}}\right)$, which are duplicated
across basins and tributaries. The basin outlets are at $\left({z}_{\mathrm{0}}=\mathrm{0},{\mathit{\chi}}_{\mathrm{0}}=\mathrm{0}\right)$, and the highest *χ* channel head is identified with $\left({z}_{q},{\mathit{\chi}}_{q}={\mathit{\chi}}_{\mathrm{max}}\right)$. The knickpoints are used to divide the
*χ*−*z* space into segments. Segment *j*, between $\left({\mathit{\chi}}_{j-\mathrm{1}},{\mathit{\chi}}_{j}\right)$, is characterized by a uniform
steepness index that shaped the river profile during time interval $\left({t}_{j-\mathrm{1}},{t}_{j}\right)$, where time *t*_{j}=*τ*_{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 *t*_{j}=*τ*_{j} based on Eq. (14) and *U*_{j} 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, ${\mathrm{KA}}_{\mathrm{0}}^{m/n}\frac{{k}_{\mathrm{s}\mathrm{\_}j}^{n}\left(\mathrm{1}-{\mathit{\gamma}}_{j}^{n}\right)}{{k}_{\mathrm{s}\mathrm{\_}j}\left(\mathrm{1}-{\mathit{\gamma}}_{j}\right)}$, ($L/T$). The scaled knickpoint age, ${t}_{j}^{*}={\mathit{\tau}}_{j}^{*}$, with dimensions of length, becomes

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, ${A}_{\mathrm{0}}^{m/n}$, (${L}^{\mathrm{2}m/n}$),
yielding a non-dimensional *K*-independent uplift rate:

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, $({U}_{j}^{*},\phantom{\rule{0.125em}{0ex}}{\mathit{\tau}}_{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, ${\mathit{\chi}}_{j}\phantom{\rule{0.125em}{0ex}}(j=\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}}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
${\mathit{\tau}}_{j}^{*}={\mathit{\chi}}_{j}$. Then, linear regression is applied
in the *χ*−*z* domain, independently for each segment. The slope of the
regression is identified as *k*_{s_j}, from which
${U}_{j}^{*}$ 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:

where *z*_{i} and $\stackrel{\mathrm{\u0303}}{{z}_{i}}$ 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. $\stackrel{\mathrm{\u0303}}{{z}_{i}}$ is obtained by integrating *U*^{*} along the
*χ*(or *t*^{*}) axis:

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 (*U*_{j}, *t*_{j}) based
on the scaled history $({U}_{j}^{*},\phantom{\rule{0.125em}{0ex}}{t}_{j}^{*})$ 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:
$\widehat{{z}_{i}}\left(\mathrm{perturbed}\right)={z}_{i-\mathrm{1}}+({z}_{i+\mathrm{1}}-{z}_{i-\mathrm{1}})\cdot \mathrm{rand}[\mathrm{0},\mathrm{1}],\phantom{\rule{0.125em}{0ex}}$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.

## 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).

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 ^{10}Be-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=\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{6}}$ m^{0.1} a^{−1}, and
the second set is *n*=2 and $K=\mathrm{4.01}\times {\mathrm{10}}^{-\mathrm{9}}$ m^{−0.8} a^{−1}. The
erodibility coefficients were inferred based on regressions through the
^{10}Be-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.

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 *n*1 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*, (${U}_{\mathrm{2}}>{U}_{\mathrm{1}}>{U}_{\mathrm{0}}$)
leading to an older kp_{1} that formed at time *t*=0 and younger
kp_{2} knickpoints that formed at time *t*=*T*_{1}.
This scenario could result in knickpoints merging and complete consumption
of the middle channel reach whose steepness index was
${k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}=({U}_{\mathrm{1}}/K{)}^{\mathrm{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 *χ*(kp_{2})>*χ*(kp_{1}) at time *t*=*T*, where *T**>**T*_{1}.
Based on Eq. (15), this condition can be expressed as

The inequality in Eq. (30), describing the relation between *T*, *T*_{1}
and *U*_{1} (through *γ*_{0_1} and *γ*_{2_1}), could be used to rule out potential lost histories
that do not obey *χ*(kp_{2})>*χ*(kp_{1}).

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.

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.

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
kp_{1}, which was formed by an uplift rate increase
from *U*_{0} to *U*_{1}, and kp_{2} formed by an
increase from *U*_{1} to *U*_{2} (${U}_{\mathrm{2}}>{U}_{\mathrm{1}}>{U}_{\mathrm{0}}$). The
celerity of the two knickpoints is expressed by Eq. (13):

Since kp_{2} is located below
kp_{1}, *A*(kp_{2}) is larger than
*A*(kp_{1}). Next, it is left to show that
$\frac{{k}_{\mathrm{s}\mathrm{\_}\mathrm{2}}^{n-\mathrm{1}}(\mathrm{1}-{\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{n})}{(\mathrm{1}-{\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}})}>\frac{{k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}^{n-\mathrm{1}}(\mathrm{1}-{\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{n})}{(\mathrm{1}-{\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}})}$. We define a variable:

Because *n*>1, we can re-write $n=\mathit{\alpha}/\mathit{\beta}$, where $\mathit{\alpha}>\mathit{\beta}>\mathrm{1}$ and *α* and *β* are both integers. Thus,

where *f*_{nume} and *f*_{deno} are the numerator and
denominator of *f*, respectively. We use the method of polynomial division:

Assigning Eq. (A4) into *f*_{nume}, we can derive

Because $({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\alpha}-\mathrm{1}-\mathit{\beta}}+({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\alpha}-\mathrm{2}-\mathit{\beta}}+\mathrm{\dots}+({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\beta}-\mathit{\beta}}$ $<\mathrm{1}+\mathrm{1}+\mathrm{\dots}+\mathrm{1}=\mathit{\alpha}-\mathit{\beta}$ and $({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{-\mathrm{1}}+({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{-\mathrm{2}}+\mathrm{\dots}+$ $({\mathit{\gamma}}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta}}{)}^{-\mathit{\beta}}>\mathrm{1}+\mathrm{1}+\mathrm{\dots}+\mathrm{1}=\mathit{\beta}$, we can derive

Again, we use polynomial division.

Assigning Eq. (A7) into *f*_{deno}, we derived

Because $({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{-\mathrm{1}}+({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{-\mathrm{2}}+\mathrm{\dots}+({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\beta}-\mathit{\alpha}}>\mathrm{1}+\mathrm{1}+\mathrm{\dots}+\mathrm{1}=\mathit{\alpha}-\mathit{\beta}$ and $({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\beta}-\mathrm{1}}+({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathit{\beta}-\mathrm{2}}+\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}}+({\mathit{\gamma}}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta}}{)}^{\mathrm{0}}<\mathrm{1}+\mathrm{1}+\mathrm{\dots}+\mathrm{1}=\mathit{\beta}$, we derived

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

Thus, *v*_{H_kp1}<*v*_{H_kp2},
kp_{2} always migrates faster than kp_{1} and given sufficient channel
length the two knickpoints will merge. The time of merging is given by
Eq. (16).

Equations (26)–(27) express the scaled uplift rate history as a series of
values (*U*^{*}, *τ*^{*}) describing the scaled age of a
knickpoint ${\mathit{\tau}}_{j}^{*}$ and the non-dimensional uplift rate,
${U}_{j}^{*}$, that operated between scaled time ${\mathit{\tau}}_{j}^{*}$
and ${\mathit{\tau}}_{j-\mathrm{1}}^{*}$, where *τ*^{*} is identified
with the *χ* axis and ${\mathit{\tau}}_{\mathrm{0}}^{*}=\mathrm{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 $\stackrel{\mathrm{\u0303}}{{z}_{j}}$ and the elevations of other
pixels $\stackrel{\mathrm{\u0303}}{{z}_{i}}$. 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 *z*_{1} to be

where ${\mathit{\gamma}}_{\mathrm{1}}={k}_{\mathrm{s}\mathrm{\_}\mathrm{2}}/{k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}$, *t*_{1} is the age of the
knickpoint, and *U*_{1} is the uplift rate that generated the knickpoint.
According to Eqs. (26)–(27), *t*_{1} and *U*_{1} are defined as

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

Then, using the definition ${k}_{\mathrm{s}\mathrm{\_}\mathrm{1}}={U}_{\mathrm{1}}^{*}\cdot {A}_{\mathrm{0}}^{-m/n}$ (Eq. 27), Eq. (B4) can be simplified to

where ${t}_{\mathrm{0}}^{*}={\mathit{\chi}}_{\mathrm{0}}=\mathrm{0}$ (basin outlet).

Then, assuming that Eq. (29) holds for knickpoint *j* with elevation
*z*_{j}, we prove the induction step for knickpoint *z*_{j+1}. Noting that
${z}_{j+\mathrm{1}}=\phantom{\rule{0.125em}{0ex}}{z}_{j}+({z}_{j+\mathrm{1}}-{z}_{j}$), we evaluate the elevation difference
between the two knickpoints *j* and *j*+1 following Eq. (21) as

Arranging Eq. (B6),

Based on Eqs. (26)–(27) and the scaling ${k}_{\mathrm{s}\mathrm{\_}j}={U}_{j}^{*}\cdot {A}_{\mathrm{0}}^{-m/n}$, we define

where ${\mathit{\gamma}}_{j}={k}_{\mathrm{s}\mathrm{\_}j+\mathrm{1}}/{k}_{\mathrm{s}\mathrm{\_}j}$. Substituting Eqs. (B8)–(B10) into (B7),

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

since ${\mathit{\gamma}}_{j}\cdot {U}_{j}^{*}=\frac{{k}_{\mathrm{s}\mathrm{\_}j+\mathrm{1}}}{{k}_{\mathrm{s}\mathrm{\_}j}}\cdot {k}_{\mathrm{s}\mathrm{\_}j}\cdot {A}_{\mathrm{0}}^{m/n}={U}_{j+\mathrm{1}}^{*}$ Eq. (B12) becomes

and

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

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

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

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 https://srtm.csi.cgiar.org/srtmdata/; SRTM Data, 2022; Farr et al., 2007).

The supplement related to this article is available online at: https://doi.org/10.5194/esurf-10-833-2022-supplement.

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.

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.

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

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, https://doi.org/10.1126/sciadv.aaz3166, 2020.

Akaike, H.: A new look at the statistical model identification, IEEE Trans. Autom. Control, 19, 716–723, https://doi.org/10.1109/TAC.1974.1100705, 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, https://doi.org/10.1029/2006JF000553, 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, https://doi.org/10.1029/2009JF001426, 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, https://doi.org/10.1029/2003JF000086, 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, https://doi.org/10.1029/2005RG000183, 2007.

Ferrier, K. L., Huppert, K. L., and Perron, J. T.: Climatic control of bedrock river incision, Nature, 496, 206–209, https://doi.org/10.1038/nature11982, 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, https://doi.org/10.1002/2014JF003196, 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, https://doi.org/10.5194/esurf-7-211-2019, 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, https://doi.org/10.1029/2020JF006060, 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, https://doi.org/10.1029/2020JF005651, 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, https://doi.org/10.5194/esurf-5-161-2017, 2017.

Gasparini, N. M. and Brandon, M. T.: A generalized power law approximation for fluvial incision of bedrock channels, J. Geophys. Res., 116, F02020, https://doi.org/10.1029/2009JF001655, 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, https://doi.org/10.1016/B978-0-12-818234-5.00075-4, 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, https://doi.org/10.1002/esp.5295, 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, https://doi.org/10.1016/j.geomorph.2016.05.035, 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, https://doi.org/10.1029/2006JF000570, 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, https://doi.org/10.1029/2008JF001187, 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, https://doi.org/10.5194/esurf-4-1-2016, 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, https://doi.org/10.1016/j.jsg.2012.07.009, 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, https://doi.org/10.1029/2019GL086882, 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, https://doi.org/10.1029/2018JF004880, 2019.

Morisawa, M. E.: Quantitative Geomorphology of Some Watersheds in the Appalachian Plateau, Geol. Soc. Am. Bull., 73, 1025–1046, https://doi.org/10.1130/0016-7606(1962)73[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, https://doi.org/10.1029/2006JF000563, 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, https://doi.org/10.1130/G30013A.1, 2009.

Paul, J. D., Roberts, G. G., and White, N.: The African landscape through space and time, Tectonics, 33, 898–935, https://doi.org/10.1002/2013TC003479, 2014.

Perron, J. T. and Royden, L.: An integral approach to Bedrock River profile analysis, Earth Surf. Process. Landf., 38, 570–576, https://doi.org/10.1002/esp.3302, 2013.

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

Roberts, G. G. and White, D.: Estimating uplift rate histories from river profiles using African examples, J. Geophys. Res., 115, B02406, https://doi.org/10.1029/2009JB006692, 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, https://doi.org/10.1002/jgrf.20031, 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, https://doi.org/10.1002/2017JF004195, 2017.

Schwanghart, W. and Scherler, D.: Divide mobility controls knickpoint migration on the Roan Plateau (Colorado, USA), Geology, 48, 698–702, https://doi.org/10.1130/G47054.1, 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, https://doi.org/10.1002/2016JF004138, 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, https://srtm.csi.cgiar.org/srtmdata/ (last access: 25 March 2022), 2022.

Steer, P.: Short communication: Analytical models for 2D landscape evolution, Earth Surf. Dynam., 9, 1239–1250, https://doi.org/10.5194/esurf-9-1239-2021, 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, https://doi.org/10.1002/2014GL062383, 2015.

Venditti, J. G., Li, G., Deal, E., Dingle, E., and Church, M.: Struggles with stream power: Connecting theory across scales, Geomorphology, 366, 106817, https://doi.org/10.1016/j.geomorph.2019.07.004, 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, https://doi.org/10.1029/2021GL092587, 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, https://doi.org/10.1029/2000JB000044, 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, https://doi.org/10.1029/2011JF002157, 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, https://doi.org/10.1130/2006.2398(04), 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, https://doi.org/10.1016/j.epsl.2019.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, https://doi.org/10.1002/2014TC003671, 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, https://doi.org/10.1016/j.epsl.2017.02.025, 2017.

- Abstract
- Introduction
- Theoretical background
- Slope-break knickpoint migration
- Knickpoint preservation and merging
- A forward analytic model for knickpoint and channel long profile evolution
- An inverse model to infer tectonic uplift rate history
- Discussion
- Conclusions
- Appendix A: A mathematical demonstration of knickpoint merging
- Appendix B: Calculating the predicted elevations based on the scaled relative uplift rate history
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement

- Abstract
- Introduction
- Theoretical background
- Slope-break knickpoint migration
- Knickpoint preservation and merging
- A forward analytic model for knickpoint and channel long profile evolution
- An inverse model to infer tectonic uplift rate history
- Discussion
- Conclusions
- Appendix A: A mathematical demonstration of knickpoint merging
- Appendix B: Calculating the predicted elevations based on the scaled relative uplift rate history
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement