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

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

$$\begin{array}{}\text{(1)}& E\left(xt\right)=\mathrm{KA}(x{)}^m{\left[{\displaystyle \frac{\partial z\left(t,x\right)}{\partial x}}\right]}^n,\end{array}$$

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

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial z(t,x)}{\partial t}}=U\left(t,x\right)-{\mathrm{KA}}^m{\left[{\displaystyle \frac{\partial z\left(t,x\right)}{\partial x}}\right]}^n,\end{array}$$

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:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial z}{\partial x}}=k_{\mathrm{s}}{A}^{-\mathit{\theta }}.\end{array}$$

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

$$\begin{array}{}\text{(4)}& {\displaystyle }z\left(x\right)=z_{\mathrm{b}}+{\left({\displaystyle \frac{E}{{\mathrm{KA}}_{\mathrm{0}}^m}}\right)}^{\frac{\mathrm{1}}{n}}\mathit{\chi }\left(x\right),\text{(5)}& {\displaystyle }\mathit{\chi }\left(x\right)=\underset{x_{\mathrm{b}}}{\overset{x}{\int }}{\left({\displaystyle \frac{A_{\mathrm{0}}}{A\left(x{}^{\prime }\right)}}\right)}^{m/n}\mathrm{d}{x}^{\prime },\end{array}$$

where z_b is the base-level elevation, and the area scale factor A₀ (L²) 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₀=1 L², 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:

$$\begin{array}{}\text{(6)}& k_{\mathrm{s}}=(U/K{)}^{\mathrm{1}/n}.\end{array}$$

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

$$\begin{array}{}\text{(7)}& k_{\mathrm{s}\mathrm{\_}\mathrm{p}}=(U_{\mathrm{P}}/K{)}^{\mathrm{1}/n}.\end{array}$$

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 U₀ to U₁ (revised from Goren et al., 2014). The blue solid line shows the steady-state channel under uplift rate U₀. 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 U₁. (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 z_t to $z_{t{+}_{\mathrm{d}t}}$. 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 U₀ to U₁ and n>1. The figure emphasizes that below and above the knickpoint the channel segments erode at rates that correspond to the new (U₁) and old (U₀) 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 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

$$\begin{array}{}\text{(8)}& \mathrm{DG}=z_{t+\mathrm{d}t}\left(x+\mathrm{d}x\right)-z_{t+\mathrm{d}t}\left(x\right)={\left({\displaystyle \frac{\partial z}{\partial x}}\right)}_{\mathrm{1}}\cdot v_{\mathrm{H}}\cdot \mathrm{d}t,\end{array}$$

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

$$\begin{array}{}\text{(9)}& \mathrm{DG}=\mathrm{DB}+\mathrm{BG},\end{array}$$

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

$$\begin{array}{}\text{(10)}& \mathrm{DB}=\left(U_{\mathrm{1}}-U_{\mathrm{0}}\right)\cdot \mathrm{d}t.\end{array}$$

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

$$\begin{array}{}\text{(11)}& \mathrm{BG}={\left({\displaystyle \frac{\partial z}{\partial x}}\right)}_{\mathrm{0}}\cdot v_{\mathrm{H}}\cdot \mathrm{d}t.\end{array}$$

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

$$\begin{array}{}\text{(12)}& v_{\mathrm{H}}={\displaystyle \frac{\left(U_{\mathrm{1}}-U_{\mathrm{0}}\right)}{{\left(\frac{\partial z}{\partial x}\right)}_{\mathrm{1}}-{\left(\frac{\partial z}{\partial x}\right)}_{\mathrm{0}}}},\end{array}$$

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

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{v}_{\mathrm{H}}& ={\displaystyle \frac{K(k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n-k_{\mathrm{s}\mathrm{\_}\mathrm{0}}^n)}{(k_{\mathrm{s}\mathrm{\_}\mathrm{1}}-k_{\mathrm{s}\mathrm{\_}\mathrm{0}})}}{A}^{m/n}={\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}}{\mathrm{KA}}^{m/n}\\ & ={\displaystyle \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}})}}{\mathrm{KA}}^{m/n},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathit{\tau }\left(x_{\mathrm{p}}\right)& =\underset{\mathrm{0}}{\overset{x_{\mathrm{p}}}{\int }}{\displaystyle \frac{\mathrm{1}}{v_{\mathrm{H}}}}\mathrm{d}x=\underset{\mathrm{0}}{\overset{x_{\mathrm{p}}}{\int }}{\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}^{m/n}}}\mathrm{d}x\\ & ={\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}_{\mathrm{0}}^{m/n}}}\cdot \mathit{\chi }\left(x_{\mathrm{p}}\right).\end{array}\end{array}$$

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

4 Knickpoint preservation and merging

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₁ formed when U₀ changed to U₁ and kp₂ formed when U₁ changed to U₂ ($U_{\mathrm{2}}>U_{\mathrm{1}}>U_{\mathrm{0}}$); then the celerity of knickpoint kp₂ is larger than that of kp₁, 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₂ can eventually reach kp₁, 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₁ formed at time t=0 and that kp₂ formed at time t=T₁, Eq. (14) is used to express the χ values of the two knickpoints at any time T>T₁ as

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}& \mathit{\chi }\left({\mathrm{kp}}_{\mathrm{2}}\right)=T\cdot K{\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{2}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n)}{k_{\mathrm{s}\mathrm{\_}\mathrm{2}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}})}}\mathrm{and}\\ & \mathit{\chi }\left({\mathrm{kp}}_{\mathrm{1}}\right)=(T+T_{\mathrm{1}})K{\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}},\end{array}\end{array}$$

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₁)=χ(kp₂). The ratio $T_{\mathrm{m}}/T_{\mathrm{1}}$ is expressed as

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{T}_{\mathrm{m}}/T_{\mathrm{1}}& ={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n\right)}{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}\right)}}/\left({\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}\right)}}\right.\\ & \left.-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n\right)}{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}\right)}}\right).\end{array}\end{array}$$

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.

Figure 2 The 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, kp₁ (upper) and kp₂ (lower), generated by two step increases in tectonic uplift rates, T₁ corresponds to the duration between the formation of kp₁ and the formation of kp₂, and T_m corresponds to the time from the emergence of kp₂ to its merging with kp₁. γ1_2 is the ratio of steepness indices above and below kp₂, and γ0_1 is the ratio of steepness indices above and below kp₁.

Figure 3 The 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 $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₂ (with a lower steepness index below knickpoint kp₂), 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₁ is more sensitive to its own celerity than to that of the younger knickpoint.

We note that when n=1, χ(kp₁)>χ(kp₂) 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

$$\begin{array}{}\text{(17)}& v_{\mathrm{H}\mathrm{\_}\mathrm{after}\mathrm{\_}\mathrm{merger}}={\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{2}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}}^n)}{k_{\mathrm{s}\mathrm{\_}\mathrm{2}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}})}}\mathrm{KA}(x_{\mathrm{p}}{)}^{m/n},\end{array}$$

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₁. 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₀ to U₂.

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

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₀ to U₁, can be expressed as

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}={\displaystyle \frac{\partial z}{\partial t}}+{\displaystyle \frac{\partial z}{\partial x}}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}},\end{array}$$

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

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}z(t,x_{\mathrm{p}}(t\left)\right)}{\mathrm{d}t}}=U\left(t\right)-{\mathrm{KA}}^m{\left(k_{\mathrm{s}\mathrm{\_}\mathrm{1}}{A}^{-\frac{m}{n}}\right)}^n\\ & +k_{\mathrm{s}\mathrm{\_}\mathrm{1}}{A}^{-\frac{m}{n}}{\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n\right)}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}\right)}}{\mathrm{KA}}^{\frac{m}{n}}=U\left(t\right)-U_{\mathrm{1}}\\ & +U_{\mathrm{1}}{\displaystyle \frac{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(20)}& z\left(t,x_{\mathrm{p}}\left(t\right)\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\left[U\left(t{}^{\prime }\right)-U_{\mathrm{1}}+U_{\mathrm{1}}{\displaystyle \frac{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}}\right]\mathrm{d}t{}^{\prime }.\end{array}$$

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

$$\begin{array}{}\text{(21)}& z\left(t,x_{\mathrm{p}}\left(t\right)\right)=\underset{\mathrm{0}}{\overset{t}{\int }}U\left(t{}^{\prime }\right)\mathrm{d}t{}^{\prime }+\left[{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}\right)}}-\mathrm{1}\right]\cdot U_{\mathrm{1}}\cdot t.\end{array}$$

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₁ and kp₂, generated at times t=0 and t=T₁, 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₁ is constrained by Eq. (16). For any time $t<T_{\mathrm{m}}+T_{\mathrm{1}}$, the elevations of kp₁ and kp₂ are predicted by Eq. (21), and when assigning the knickpoint ages t= τ(x_p1) and $t-T_{\mathrm{1}}=\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₁ (t=0) can be expressed as follows.

$$\begin{array}{}\text{(22)}& {\displaystyle }\begin{array}{rl}{z}_{\mathrm{kp}\mathrm{\_}\mathrm{12}}& =z_{\mathrm{1}}\left(T_{\mathrm{m}}+T_{\mathrm{1}}\right)+z_{\mathrm{12}}\mathrm{or}{z}_{\mathrm{kp}\mathrm{\_}\mathrm{12}}=z_{\mathrm{2}}(T_{\mathrm{m}}+T_{\mathrm{1}})\\ & +z_{\mathrm{12}},\end{array}\text{(23)}& {\displaystyle }\mathrm{where}\left\{\begin{array}{l}{z}_{\mathrm{1}}\left(t=T_{\mathrm{m}}+T_{\mathrm{1}},x_{\mathrm{p}\mathrm{1}}\right)=\underset{\mathrm{0}}{\overset{T_{\mathrm{m}}+T_{\mathrm{1}}}{\int }}U\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }+\\ \left[\frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}\right)}-\mathrm{1}\right]\cdot U_{\mathrm{1}}\cdot \left(T_{\mathrm{m}}+T_{\mathrm{1}}\right)\\ z_{\mathrm{2}}\left(t=T_{\mathrm{m}}+T_{\mathrm{1}},x_{\mathrm{p}\mathrm{2}}\right)=\underset{T_{\mathrm{1}}}{\overset{T_{\mathrm{m}}+T_{\mathrm{1}}}{\int }}U\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }+\\ \left[\frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}\right)}-\mathrm{1}\right]\cdot U_{\mathrm{2}}\cdot T_{\mathrm{m}}\\ z_{\mathrm{12}}(t>T_{\mathrm{m}}+T_{\mathrm{1}},x_{\mathrm{p}\mathrm{12}})=\underset{T_{\mathrm{m}}+T_{\mathrm{1}}}{\overset{t}{\int }}U\left(t{}^{\prime }\right)\mathrm{d}{t}^{\prime }+\\ \left[\frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}}\right)}-\mathrm{1}\right]\cdot U_{\mathrm{2}}\cdot (t-\left(T_{\mathrm{m}}+T_{\mathrm{1}}\right))\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(24)}& x_{\mathrm{p}}\left(t\right)={\mathit{\chi }}^{-\mathrm{1}}\left[{\mathrm{KA}}_{\mathrm{0}}^{m/n}t{\displaystyle \frac{k_{s\mathit{\_}\mathrm{1}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)}{k_{s\mathit{\_}\mathrm{1}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}}\right],\end{array}$$

where again t= τ(x_p) is the knickpoint age. After merging, for $t>T_{\mathrm{m}}+T_{\mathrm{1}}$

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{x}_{\mathrm{p}}\left(t\right)& ={\mathit{\chi }}^{-\mathrm{1}}\left\{\left[{\mathrm{KA}}_{\mathrm{0}}^{m/n}{T}_{\mathrm{m}}{\displaystyle \frac{k_{s\mathit{\_}\mathrm{2}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n)}{k_{s\mathit{\_}\mathrm{2}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}})}}\right]\right.\\ & \left.+\left[{\mathrm{KA}}_{\mathrm{0}}^{m/n}(t-T_{\mathrm{m}}-T_{\mathrm{1}}){\displaystyle \frac{k_{s\mathit{\_}\mathrm{2}}^n(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}}^n)}{k_{s\mathit{\_}\mathrm{2}}(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{2}})}}\right]\right\}.\end{array}\end{array}$$

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⁻¹. Then, we set two step increases in the rates to be 0.5 mm a⁻¹ at 0.8 Ma and 1.0 mm a⁻¹ 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 4 Comparison 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=k_{\mathrm{a}}(L-x{)}^h$. The model parameters are n=3, $K=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{11}}$ m^−1.7 a⁻¹, m=1.35, k_a=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 U₀=0.10 mm a⁻¹ prior to 0.8 Ma, U₁=0.5 mm a⁻¹ between 0.5–0.8 Ma and U₂=1.0 mm a⁻¹ 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 x−z (a) and χ−z (b) domains. Panel (c) depicts the river χ−z long profiles offset in elevation, demonstrating knickpoint merging dynamics. Knickpoint kp₁ formed at 0.8 Ma, as a response to the increase in the uplift rate from U₀ to U₁. Knickpoint kp₂ formed at 0.5 Ma, due to uplift rate increase from U₁ to U₂. At ∼0.2 Ma, kp₁ merged with kp₂

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, $\left(z_{\mathrm{1}},{\mathit{\chi }}_{\mathrm{1}}\right),\left(z_{\mathrm{2}},{\mathit{\chi }}_{\mathrm{2}}\right),\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 }}_{\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

$$\begin{array}{}\text{(26)}& t_j^{*}=t_j\cdot {\mathrm{KA}}_{\mathrm{0}}^{m/n}\cdot {\displaystyle \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)}}={\mathit{\chi }}_j.\end{array}$$

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:

$$\begin{array}{}\text{(27)}& U_j^{*}=A_{\mathrm{0}}^{-m/n}{k}_{\mathrm{s}\mathrm{\_}j}=A_{\mathrm{0}}^{-m/n}\cdot (U_j/K{)}^{\mathrm{1}/n}.\end{array}$$

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^{*},{\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(j=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }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:

$$\begin{array}{}\text{(28)}& \mathrm{misfit}={\displaystyle \frac{\mathrm{1}}{N/M}}\sqrt{{\sum }_{i=\mathrm{1}}^N(z_i-\stackrel{\mathrm{̃}}{z_i}{)}^{\mathrm{2}}},\end{array}$$

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

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{z_i}& =\underset{\mathrm{0}}{\overset{t_i^{*}={\mathit{\chi }}_i}{\int }}{U}^{*}\left(t^{{*}^{\prime }}\right)\mathrm{d}{t}^{{*}^{\prime }}={\sum }_{a=\mathrm{1}}^j{U}_a^{*}(t_a^{*}-t_{a-\mathrm{1}}^{*})\\ & +U_{j+\mathrm{1}}^{*}\left(t_i^{*}-t_j^{*}\right),\end{array}\end{array}$$

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^{*},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}],$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 5 Inversion 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, k_a∼2 to 1.55 and h∼0.67 to 4.27. The stream-power parameters are n=2, $K=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ m^−0.8 a⁻¹ and m=0.9. The time (dt) and space (dx) steps are 100 years and 100 m, respectively. The applied uplift rate history is U₀=0.05 mm a⁻¹ prior to 2.5 Ma, U₁=0.2 mm a⁻¹ between 0.5–2.5 Ma and U₂=1.0 mm a⁻¹ 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 6 Application 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⁻¹ before the middle Miocene and gradually increased to ∼0.35 mm a⁻¹ 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 ¹⁰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⁻¹, and the second set is n=2 and $K=\mathrm{4.01}\times {\mathrm{10}}^{-\mathrm{9}}$ m^−0.8 a⁻¹. The erodibility coefficients were inferred based on regressions through the ¹⁰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⁻¹) and the present-day rate (∼0.375 mm a⁻¹). With n=2, the oldest relative uplift rate (no more than 0.05 mm a⁻¹) is slower by approximately a factor of 10 with respect to recent rates (∼0.5 mm a⁻¹). 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, ($U_{\mathrm{2}}>U_{\mathrm{1}}>U_{\mathrm{0}}$) leading to an older kp₁ that formed at time t=0 and younger kp₂ knickpoints that formed at time t=T₁. 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₂)>χ(kp₁) at time t=T, where T>T₁. Based on Eq. (15), this condition can be expressed as

$$\begin{array}{}\text{(30)}& T_{\mathrm{1}}/T<{\displaystyle \frac{\mathrm{1}-{\mathit{\gamma }}_{\mathrm{2}\mathit{\_}\mathrm{1}}^n}{\mathrm{1}-{\mathit{\gamma }}_{\mathrm{2}\mathit{\_}\mathrm{1}}}}\cdot {\displaystyle \frac{\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}}{\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n}}-\mathrm{1},\mathrm{with}{\mathit{\gamma }}_{\mathrm{2}\mathit{\_}\mathrm{1}}=k_{\mathrm{s}\mathrm{\_}\mathrm{2}}/k_{\mathrm{s}\mathrm{\_}\mathrm{1}}.\end{array}$$

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

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

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& v_{\mathrm{H}\mathrm{\_}\mathrm{kp}\mathrm{1}}={\displaystyle \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}})}}\mathrm{KA}({\mathrm{kp}}_{\mathrm{1}}{)}^{m/n}\mathrm{and}\\ & v_{\mathrm{H}\mathrm{\_}\mathrm{kp}\mathrm{2}}={\displaystyle \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}})}}\mathrm{KA}({\mathrm{kp}}_{\mathrm{2}}{)}^{m/n}.\end{array}\end{array}$$

Since kp₂ is located below kp₁, A(kp₂) is larger than A(kp₁). 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:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}f& ={\displaystyle \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}})}}/{\displaystyle \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}})}}\\ & ={\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}-n}}}{\displaystyle \frac{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^n)/(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}})}{(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^n)/(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}})}}.\end{array}\end{array}$$

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,

$$\begin{array}{}\text{(A3)}& f={\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}-\mathit{\alpha }/\mathit{\beta }}}}{\displaystyle \frac{(\mathrm{1}-({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }})/(\mathrm{1}-\left({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}\right)}{(\mathrm{1}-({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }})/(\mathrm{1}-\left({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}\right)}}={\displaystyle \frac{f_{\mathrm{nume}}}{f_{\mathrm{deno}}}},\end{array}$$

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

$$\begin{array}{}\text{(A4)}& \left\{\begin{array}{l}\mathrm{1}-({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }}=(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}\left)\right(({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{1}}+\mathrm{\dots }\\ +({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}})\\ \mathrm{1}-({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}=(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}\left)\right(({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{1}}+\\ ({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{2}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}})\end{array}\right..\end{array}$$

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

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& f_{\mathrm{nume}}={\displaystyle \frac{({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{1}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}{({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{1}}+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{2}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}\\ & ={\displaystyle \frac{({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{1}}+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{2}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }}}{({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{1}}+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathrm{2}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{1}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}+\mathrm{1}\\ & ={\displaystyle \frac{({\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 }}}{({\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}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A6)}& f_{\mathrm{nume}}<{\displaystyle \frac{\mathit{\alpha }-\mathit{\beta }}{\mathit{\beta }}}+\mathrm{1}.\end{array}$$

Again, we use polynomial division.

$$\begin{array}{}\text{(A7)}& \left\{\begin{array}{l}\mathrm{1}-({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }}=(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}\left)\right(({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{1}}+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{2}}\\ +\mathrm{\dots }+\left({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}\right)\\ \mathrm{1}-({\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 }}\left)\right(({\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}}\\ +\mathrm{\dots }+\left({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}\right)\end{array}\right.\end{array}$$

Assigning Eq. (A7) into f_deno, we derived

$$\begin{array}{}\text{(A8)}& {\displaystyle }\begin{array}{rl}& f_{\mathrm{deno}}={\displaystyle \frac{({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{1}}+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\alpha }-\mathrm{2}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}{({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}\cdot ({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathit{\alpha }}\\ & ={\displaystyle \frac{({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathit{\beta }-\mathit{\alpha }}}{({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}\mathrm{or}{f}_{\mathrm{deno}}& ={\displaystyle \frac{\begin{array}{c}({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}+\\ ({\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 }}\end{array}}{({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}\\ & ={\displaystyle \frac{({\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 }}}{({\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}}+\mathrm{\dots }+({\mathit{\gamma }}_{\mathrm{1}\mathit{\_}\mathrm{2}}^{\mathrm{1}/\mathit{\beta }}{)}^{\mathrm{0}}}}+\mathrm{1}.\end{array}\end{array}$$

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}}+\mathrm{\dots }+({\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

$$\begin{array}{}\text{(A10)}& f_{\mathrm{deno}}>{\displaystyle \frac{\mathit{\alpha }-\mathit{\beta }}{\mathit{\beta }}}+\mathrm{1}.\end{array}$$

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

$$\begin{array}{}\text{(A11)}& f={\displaystyle \frac{f_{\mathrm{nume}}}{f_{\mathrm{deno}}}}<({\displaystyle \frac{\mathit{\alpha }-\mathit{\beta }}{\mathit{\beta }}}+\mathrm{1})/({\displaystyle \frac{\mathit{\alpha }-\mathit{\beta }}{\mathit{\beta }}}+\mathrm{1})=\mathrm{1}.\end{array}$$

Thus, v_{H_kp1}<v_{H_kp2}, kp₂ always migrates faster than kp₁ 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 ${\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{̃}}{z_j}$ and the elevations of other pixels $\stackrel{\mathrm{̃}}{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₁ to be

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{z}_{\mathrm{1}}& =\underset{\mathrm{0}}{\overset{t_{\mathrm{1}}}{\int }}{U}_{\mathrm{1}}\mathrm{d}t+\left[{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}\right)}}-\mathrm{1}\right]\cdot U_{\mathrm{1}}\cdot t_{\mathrm{1}}\\ & ={\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}\right)}}\cdot U_{\mathrm{1}}\cdot t_{\mathrm{1}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B2)}& {\displaystyle }{t}_{\mathrm{1}}=t_{\mathrm{1}}^{*}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}_{\mathrm{0}}^{m/n}}}\cdot {\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}\right)}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}^n\right)}},\text{(B3)}& {\displaystyle }{U}_{\mathrm{1}}=(U_{\mathrm{1}}^{*}\cdot A_{\mathrm{0}}^{m/n}{)}^n\cdot K.\end{array}$$

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

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}{z}_{\mathrm{1}}& ={\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}\right)}}\cdot (U_{\mathrm{1}}^{*}\cdot A_{\mathrm{0}}^{m/n}{)}^n\cdot K\cdot t_{\mathrm{1}}^{*}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}_{\mathrm{0}}^{m/n}}}\\ & \cdot {\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}\right)}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}^n\right)}}=(U_{\mathrm{1}}^{*}\cdot A_{\mathrm{0}}^{m/n}{)}^n\cdot t_{\mathrm{1}}^{*}\cdot {\displaystyle \frac{\mathrm{1}}{A_{\mathrm{0}}^{m/n}}}\cdot {\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}}{k_{\mathrm{s}\mathrm{\_}\mathrm{1}}^n}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B5)}& z_{\mathrm{1}}=U_{\mathrm{1}}^{*}\cdot t_{\mathrm{1}}^{*}=U_{\mathrm{1}}^{*}\cdot (t_{\mathrm{1}}^{*}-t_{\mathrm{0}}^{*}),\end{array}$$

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}}=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

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{z}_{j+\mathrm{1}}-z_j& =\underset{t_j}{\overset{t_{j+\mathrm{1}}}{\int }}{U}_{j+\mathrm{1}}\mathrm{d}t+\left[{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}\right)}}-\mathrm{1}\right]\\ & \cdot U_{j+\mathrm{1}}\cdot t_{j+\mathrm{1}}-\left[{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}}-\mathrm{1}\right]\cdot U_j\cdot t_j.\end{array}\end{array}$$

Arranging Eq. (B6),

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{z}_{j+\mathrm{1}}-z_j& =-U_{j+\mathrm{1}}\cdot t_j+{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}\right)}}\cdot U_{j+\mathrm{1}}\cdot t_{j+\mathrm{1}}\\ & -{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}}\cdot U_j\cdot t_j+U_j\cdot t_j.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B8)}& {\displaystyle }\begin{array}{rl}{U}_j\cdot t_j& =t_j^{*}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}_{\mathrm{0}}^{m/n}}}\cdot {\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}j}\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{k_{\mathrm{s}\mathrm{\_}j}^n\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\cdot (U_j^{*}\cdot A_{\mathrm{0}}^{m/n}{)}^n\\ & \cdot K=t_j^{*}\cdot U_j^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}},\end{array}\text{(B9)}& {\displaystyle }{U}_{j+\mathrm{1}}\cdot t_{j+\mathrm{1}}=t_{j+\mathrm{1}}^{*}\cdot U_{j+\mathrm{1}}^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_{j+\mathrm{1}}^n\right)}},\text{(B10)}& {\displaystyle }\begin{array}{rl}{U}_{j+\mathrm{1}}\cdot t_j& =t_j^{*}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathrm{KA}}_{\mathrm{0}}^{m/n}}}\cdot {\displaystyle \frac{k_{\mathrm{s}\mathrm{\_}j}\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{k_{\mathrm{s}\mathrm{\_}j}^n\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\cdot (U_{j+\mathrm{1}}^{*}\\ & \cdot A_{\mathrm{0}}^{m/n}{)}^n\cdot K=t_j^{*}\cdot U_j^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{k_{\mathrm{s}\mathrm{\_}j}^n\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\\ & \cdot k_{\mathrm{s}\mathrm{\_}j+\mathrm{1}}^n=t_j^{*}\cdot U_j^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\cdot {\mathit{\gamma }}_j^n,\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{z}_{j+\mathrm{1}}-z_j& =-t_j^{*}\cdot U_j^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\cdot {\mathit{\gamma }}_j^n+t_{j+\mathrm{1}}^{*}\cdot U_{j+\mathrm{1}}^{*}\\ & -t_j^{*}\cdot U_j^{*}+t_j^{*}\cdot U_j^{*}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B12)}& \begin{array}{rl}{z}_{j+\mathrm{1}}-z_j& =t_{j+\mathrm{1}}^{*}\cdot U_{j+\mathrm{1}}^{*}+t_j^{*}\cdot U_j^{*}\\ & \cdot \left[{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}-{\displaystyle \frac{\left(\mathrm{1}-{\mathit{\gamma }}_j\right)}{\left(\mathrm{1}-{\mathit{\gamma }}_j^n\right)}}\cdot {\mathit{\gamma }}_j^n-\mathrm{1}\right]\\ & =t_{j+\mathrm{1}}^{*}\cdot U_{j+\mathrm{1}}^{*}-{\mathit{\gamma }}_j\cdot t_j^{*}\cdot U_j^{*}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B13)}& z_{j+\mathrm{1}}-z_j=t_{j+\mathrm{1}}^{*}\cdot U_{j+\mathrm{1}}^{*}-t_j^{*}\cdot U_{j+\mathrm{1}}^{*}=U_{j+\mathrm{1}}^{*}\cdot (t_{j+\mathrm{1}}^{*}-t_j^{*})\end{array}$$

and

$$\begin{array}{}\text{(B14)}& z_{j+\mathrm{1}}=z_j+U_{j+\mathrm{1}}^{*}\cdot (t_{j+\mathrm{1}}^{*}-t_j^{*}).\end{array}$$

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

$$\begin{array}{}\text{(B15)}& z_j=\underset{\mathrm{0}}{\overset{t_j^{*}}{\int }}{U}^{*}\left(t^{{*}^{\prime }}\right)\mathrm{d}{t}^{{*}^{\prime }}={\sum }_{a=\mathrm{1}}^j{U}_a^{*}(t_a^{*}-t_{a-\mathrm{1}}^{*}).\end{array}$$

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

$$\begin{array}{}\text{(B16)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{z_i}& =z_j+k_{\mathrm{s}\mathrm{\_}j+\mathrm{1}}\cdot A_{\mathrm{0}}^{-m/n}\cdot ({\mathit{\chi }}_i-{\mathit{\chi }}_j)=\underset{\mathrm{0}}{\overset{t_{j^{*}}}{\int }}{U}^{*}\left(t^{{*}^{\prime }}\right)\mathrm{d}{t}^{{*}^{\prime }}\\ & +U_{j+\mathrm{1}}^{*}\cdot ({\mathit{\chi }}_i-{\mathit{\chi }}_j)=\underset{\mathrm{0}}{\overset{t_i^{*}={\mathit{\chi }}_i}{\int }}{U}^{*}\left(t^{{*}^{\prime }}\right)\mathrm{d}{t}^{{*}^{\prime }},\end{array}\end{array}$$

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