the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
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
Huiping Zhang
The long profile of rivers is shaped by the tectonic history that acted on the landscape. Faster uplift produces steeper channel segments, and knickpoints form in response to changes in the tectonic uplift rates. However, when the fluvial incision depends nonlinearly 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 nonlinear 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 nonlinear 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 lowtemperature 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 nonlinearly related.
 Article
(3312 KB)  Fulltext XML

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 streampower incision model (SPIM) is widely used to describe detachmentlimited 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 powerlaw 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 nonlinear advection equation for the elevation, where U acts as a forcing term. Consequently, Eq. (2) predicts the firstorder 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ándezBlanco, 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 longprofile 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 nonlinear cases, when n=1, stretched and consumed channel reaches do not occur, and there is a 1to1 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ándezBlanco, 2021).
When n=1, it is well accepted that, under a wellconstrained 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 downstream of the knickpoints by using paleochannel projection, and other studies attempted a nonlinear 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 powerlaw function of the drainage area:
Notably, the powerlaw 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 baselevel 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 steadystate 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 slopebreak 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, slopebreak 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 detachmentlimited conditions might be violated. This behavior is beyond the scope of the current analysis.
To predict the retreat rate of slopebreak 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 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
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 rewritten 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 convexup 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 fastmigrating 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 convexup 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 concaveup 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 slopebreak 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 slopebreak 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 closedform analytic solution for slopebreak 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 piecewise 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 longprofile 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 steadystate 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 1D upwind firstorder finitedifference 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 timedependent tectonic uplift history U(t), and the block has a uniform erodibility. The model infers a single bestfit 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 Kindependent version for the knickpoint age and uplift rate. To derive a Kindependent 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 Kindependent uplift rate, Eq. (7) is divided by a steepness index scale factor, ${A}_{\mathrm{0}}^{m/n}$, (${L}^{\mathrm{2}m/n}$), yielding a nondimensional Kindependent 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 Kindependent but also nindependent. 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 slopebreak 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 slopebreak knickpoints. Several algorithms have been previously proposed to identify slopebreak 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ándezBlanco, 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 lowresolution numerical model that solves Eq. (2). The model is used to generate 10 river profiles with variable channel length and drainage area distribution with prechosen 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 prechosen 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 bestfit solution has two division points. Figure 5d shows the inferred history, indicating that the twodivisionpoint 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 tectonogeomorphic 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 lowtemperature 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 catchmentwide 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 tectonogeomorphic 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 nondimensional 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}Bederived 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}Bederived 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 tectonogeometric 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 presentday 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 nonunity 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 piecewise 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, welldefined 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 spaceinvariant uplift (also consistent with a baselevel fall). This assumption, which is commonly referred to as block uplift conditions, is more likely to hold over discrete, welldefined 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 processbased 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 slopebreak 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 slopebreak 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 highelevation 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ándezBlanco, 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 selfadjusting 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ándezBlanco, 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 slopebreak knickpoint retreat model under the assumption of spaceinvariant 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 concaveup (n<1, decreasing U) knickpoints. Knickpoint celerity depends on the streampower 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 lowtemperature 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 rewrite $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 nondimensional 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/esurf108332022supplement.
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 1D 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 MunozSalinas, E.: Knickpoint retreat and landscape evolution of the Amatlan de Canas halfgraben (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/esurf72112019, 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ándezBlanco, D.: A new datadriven 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 southcentral Crete, Greece, Earth Surf. Dynam., 5, 161–186, https://doi.org/10.5194/esurf51612017, 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/B9780128182345.000754, 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 streamgradient 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 streampower approach and the χ method, Earth Surf. Dynam., 4, 1–9, https://doi.org/10.5194/esurf412016, 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 detachmentlimited 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/00167606(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 longterm stream erosion models, Earth Surf. Process. Landf., 26, 1317–1332, 2001.
Oskin, M. E. and Burbank, D.: Transient landscape evolution of basementcored 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 baselevel 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/esurf912392021, 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 sedimentfluxdependent 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
 Slopebreak 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
 Slopebreak 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