Short communication: Forward and inverse models relating river long profile to monotonic step-changes in tectonic rock uplift rate history: A theoretical perspective under a nonlinear slope-erosion dependency

. The long profile of rivers is widely considered as a recorded of tectonic uplift rate. Knickpoints form in response to rate changes and faster rates produce steeper channel segments. However, when the exponent relating fluvial incision to river slope, n , is not unity, the links between tectonic rates and channel profile are complicated by channel dynamics that consume and form river segments. Here, we explore non-linear cases leading to channel segment consumption and develop a Lagrangian analytic model for knickpoint migration. We derive a criterion for knickpoint preservation and merging, and 15 develop a forward analytic model that resolves knickpoint and long profile evolution before and after knickpoint merging. We further propose a linear inverse scheme to infer tectonic history from river profiles when all knickpoints are preserved. Our description provides a new framework to explore the links between tectonic uplift rates and river profile evolution when n is not unity. The knickpoint merging analysis further emphasizes a critical property of the links between tectonic and long profile evolution when 𝑛 ≠ 1 . Each tectonic history is associated with a single, well-defined river profile at any given time. However, any particular river profile could be generated by infinitely many tectonic histories. All histories except for one lead to knickpoint merging dynamics. The linear inverse model that we develop here finds the single history for which all knickpoint are preserved. While this inverse approach is highly restrictive, it finds the correct solution when only a single 300 knickpoint exists 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 history that shaped the fluvial landscape.


Introduction 20
Bedrock rivers that incise into tectonically active highlands are sensitive to changes in the tectonic conditions (Whipple and Tucker, 1999;Kirby et al., 2003). Upon a change in the rock uplift rate with respect to a base level, the river steepness changes (Wobus et al., 2006;Kirby and Whipple, 2001;Whipple and Tucker, 2002), which in turn, changes the local incision rate. Particularly, an increase in uplift rate generates steeper slopes that facilitate faster incision, overall promoting incision-uplift equilibrium. However, equilibration is not achieved synchronously across the river long profile. Upon a 25 tectonic change, a knickpoint forms that divides the profile to 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, while above the knickpoint, river steepness and erosion rate correspond to the previous conditions (Niemann et al., 2001;Kirby and Whipple, 2012). The https://doi.org/10.5194/esurf-2021-101 Preprint. Discussion started: 21 December 2021 c Author(s) 2021. CC BY 4.0 License.
Previous, general analytic exploration of equation (2) (e.g., Luke, 1972;Weissel and Seidl, 1998;Prichard et al., 2009;Royden and Perron, 2013) have identified that upon a tectonic change that induces a long-profile steepness change, portions 60 of the solution, representing the river profile, could form that are not strictly associated with the tectonic change, and, portions of the solution that hold tectonic information may be lost. More specifically, when U increases and < 1 or U decreases and > 1, 'stretched zone' along a river long profile form that are not associated with any particular tectonic input . When U increases and > 1 or U decreases and < 1, some portions of the channel reach are consumed at knickpoints . Unlike the non-linear cases, when = 1, stretched and 65 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 ( = 1) were applied as part of forward (Steer, 2021) and inverse models (for a recent review see, Goren et al., 2021) of tectonically forced fluvial landscape evolution.
However, while some field studies support the linearity assumption (e.g., Wobus et al. 2006, Ferrier et al. 201370 Schwanghart and Scherler 2020), a growing body of work show that n could be different than unity Harkins et al., 2007;Lague, 2014;Harel et al., 2016). Particularly, significant incision thresholds and relatively small discharge variability are expected to lead to > 1 (Anthony and Granger, 2007;Ouimet et al., 2009;Dibiase et al., 2011;Lague, 2014). The current study addresses this gap by developing a simple analytic description of the evolution of channel long profile for the cases where channel reaches may be consumed, namely, U(t) is a staircase decreasing function and < 1, 75 or U(t) is a staircase increasing function and > 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 in a Lagrangian perspective that follows the knickpoint along their migration path. With this approach we develop a criterion for knickpoint preservation and merging, a simple and easy implement forward analytic model, and a linear inverse model 80 constrained by knickpoint preservation. The current study focuses on the analytic derivations and their implications, intentionally leaving field applications to future studies.

Theoretical background
The SPIM model, equation (1) predicts that for channel segments that erode at the uniform rate (in space and time), the channel slope scales as a power-law function of the drainage area (Wobus et al., 2006;Cyr et al., 2010): 85 where = / and s = ( / ) 1/ (L 2m/n ) are commonly referred to as the channel concavity and steepness indices, respectively (Wobus et al., 2006). An alternative perspective to equation (3)  while assuming constant E/K. Following such an integration, a linear relation emerges between the elevation, z, and the parameter χ (L) : 90 where zb is the base-level elevation, and the factor A0 (L 2 ) is introduced to maintain the χ dimensions to length. The parameter χ depends only on the drainage area distribution along the channel, which can easily be calculated for any m/n as part of basic morphometric analysis . When setting 0 = 1 L 2 , the slope of the χ-z plot becomes 95 channel steepness index, ks.
Under steady-state conditions, when / = 0 and = , the SPIM steepness index becomes a function of the tectonic uplift rate: When U varies in time, equation (6) can be used to express transient conditions, where a channel segment is eroding at a rate 100 that corresponds to some previous uplift rate, Up (Niemann et al., 2001;Goren et al., 2014). In this case, its steepness index could be expressed as:

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 105 profile (Haviv et al., 2010). Within the scope of the SPIM, slope-break knickpoints are commonly associated with a step change in the tectonic uplift rate. When the rate increases, the slope and steepness index below 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 110 behaviour is beyond the scope of the current analysis. Also outside of the scope of the current study, are the cases of 'stretch zones' (U increases and < 1, and, U decreases and > 1) .
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 Figure 1.  Figure 1a shows the predicted channel profile evolution following a step increase in the rock uplift rate from U0 to U1 and 115 > 1. The figure emphasizes that below and above the knickpoint, the channel segments erode at rates that correspond to the new (U1) and old (U0) uplift rates, respectively, and their corresponding steepness indices are _1 = ( 1 / ) 1/ and _0 = ( 0 / ) 1/ . Figure 1b shows the linearized channel segments near the knickpoint. The river profile varies from zt to zt+dt during time step dt, accompanied by the knickpoint migrating from point A to D. Segment DG represents the vertical change in knickpoint location, and it can be expressed as: 120 where vH is the horizontal velocity for the knickpoint retreat (hereafter, knickpoint celerity). Figure 1b shows that: Where DB is a function of the difference between present uplift rate (U1) and previous river incision rate, U0: The segment BG is the elevation difference between points A and B: Combining equations (8-11), we solve for the knickpoint celerity: which resembles the derivation of Whipple and Tucker (1999). Assigning equations (1, 6-7) into (12), vH can be re-written as: 130 where 0_1 = s_0 / s_1 . Accordingly, the fluvial response time, τ(xp) is expressed as: The response time is the time for a perturbation, e.g., a knickpoint, to propagate from the river outlet ( = 0) to its present location xp. Alternatively, τ(xp), can also be thought of as the knickpoint age (Gallen and Wegmann, 2017), or the time before 135 the present when the knickpoint formed at the river outlet.
Importantly, 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, equations (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 are not merging, as discussed in the following section, knickpoints 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 > 1), equations (12-14) are valid also for concave knickpoints (decreasing U and < 1, see details in supplementary Text S1). For = 1, vH 145 and τ(xp) are independent of the steepness indices and their ratio. Supplementary Text S2 compares the current derivation to previous models of knickpoint celerity (Rosenbloom and Anderson, 1994;Weissel and Seidl, 1998;Oskin and Burbank, 2007;Castillo et al., 2017).

Knickpoint preservation and merging
When more than a single knickpoint propagates upstream a channel profile and ≠ 1, the sensitivity of knickpoint celerity 150 to ks and γ leads to potentially complex interactions between the knickpoints. Considering the case of > 1 and two knickpoints that formed by two step-increase in tectonic uplift rate: kp 1 formed when U0 changed to U1 and kp 2 formed when U1 changed to U2 ( 2 > 1 > 0 ), then the celerity of knickpoint kp 2 is larger than that of kp 1 , and the distance between them gradually decreases (see detailed demonstration 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 155 consuming knickpoint in Royden and Perron, 2013). To elucidate knickpoint merging dynamics, we derive an expression for the time of knickpoint merging. Assuming that kp 1 formed at time = 0 and that kp 2 formed at time = T 1 , equation (14) is used to express the values of the two knickpoints at any time T 2 > T 1 as: , and (kp 1 ) = (T 2 + T 1 ) s_1 (1− 0_1 ) where 1_2 = s_1 / s_2 . Knickpoints merging occur at time T 2_m when (kp 1 ) = (kp 2 ). The ratio T 2_m /T 1 is expressed as: 160 Equation (16) predicts that the timing of knickpoint merging depends on the ratios of channel steepness indices but not on steepness indices themselves. We present a detailed description of the relationship between T 2_m /T 1 , slope exponent, and the steepness ratios (Figure 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 T2_m/T1 165 decreases with n ( Figure S2a). This means that a higher slope exponent reduces the life expectancy of knickpoints. Figure 2a also shows that for a constant n, lower steepness indices ratio leads to lower T2_m/T1. To explore the dependency of T2_m/T1 on γ1_2 and γ0_1, we fix n = 2 and vary each of the steepness ratios independently (Figure 2b it is found that T2_m/T1 is more sensitive to γ1_2 than to γ0_1, indicating that the celerity of the younger knickpoint has a greater control over the timing of knickpoint merging. 170 For the case of concave-up consuming knickpoints (n < 1 and decreasing U), figure 3a shows that the ratio T2_m/T1 (when γ1_2 = γ0_1) increases with increasing n, and for a constant n, a higher steepness ratio leads to a higher T2_m/T1 ratio. This means that a lower uplift rate, U2 (with a lower steepness index below knickpoint kp2) leads to a shorter time to knickpoint merging T2_m. In Figures 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 T2_m/T1 is more sensitive to γ0_1 than to γ1_2, 175 indicating that the preservation time of kp1 is more sensitive to its own celerity than to that of the younger knickpoint.
We note that when = 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 knickpoint correspond to ks_0 and ks_2, respectively. Based on equation (13), the instantaneous merged knickpoint celerity 180 becomes: where 0_2 = s_0 / s_2 . The channel reach that used to stretch between the two knickpoints is fully consumed, and the channel profile holds no record of U1. Consequently, evaluating the merged knickpoint age by using equation (14) and the steepness indices above and below the merged knickpoint does not yield a meaningful answer. The reason is that upon 185 merging, the steepness indices above and below the merged knickpoint change. Critically, the channel profile does not hold any clue for the event of knickpoint merging, and the river profile would be indistinguishable from a case of a single step increase in uplift rate from U0 to U2.

A forward analytic model for knickpoint and channel long profile evolution
The elevation change of slope-break knickpoint, ( , ) = [ , = ( )], formed by a step-increase in uplift rate from U0 190 to U1, can be expressed as: where ( ) = H is the knickpoint celerity. Combining equations (2), (13) and (18) yields: Integrating equation (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 equation (20) are time invariant, and the elevation of the knickpoint could be more simply expressed as: Equation (21)  Next, we combine equation (21), which is conditioned by knickpoint preservation, with equation (16) that predicts the duration of preservation to generate a piecewise solution for knickpoint elevation beyond merging. We consider the case of two knickpoints, kp1 and kp2, generated by two step-increase of 2 > 1 > 0 , and > 1. The time of merging, T2_m, is 205 constrained by equation (16). For any < T 2_m + T 1 , the elevations of kp1 and kp2 is predicted by equation (21), when assigning the knickpoint ages, = ( ), which corresponds to the time since the change in U(t) that generated the knickpoint. Upon merging, when > T 2_m + T 1 , the elevation of the merged knickpoint, kp_12 , with respect to the formation time of kp1 ( = 0) can be expressed as: kp_12 = 1 (T 2_m + T 1 ) + 12 or kp_12 = 2 (T 2_m + T 1 ) + 12 , (1− 0_1 ) − 1] • 1 • (T 2_m + T 1 ) 2 ( = T 2_m + T 1 , 2 ) = ∫ ( ′) (1− 1_2 ) − 1] • 2 • (T 2_m ) 12 ( > T 2_m + T 1 , 12 ) = ∫ ( ′ ) Before merging, the horizontal position of the knickpoints can be expressed as the inverse of equation (14): where again, = ( ), is the knickpoint age. After merging, for > T 2_m + T 1  (14) is used for the knickpoint 220 values. The channel profile between knickpoints is represented in the -z domain as a linear line connecting the knickpoints. To illustrate long-profile and knickpoint time evolution before and after knickpoint merging and to demonstrate the validity of the analytic forward model, Figure 4 shows the consistence between the analytic solution and a 1-D upwind first-order finite-difference solver of equation (2), for a channel that experiences two step-increases in U.

An inverse model to estimate tectonic uplift rate history 225
Here, the analytic solution for knickpoint evolution is used to derive a linear inverse model for retrieving the tectonic uplift history from river long profile. The inverse model relaxes the critical assumption of = 1 that was a precondition for previous linear inverse models (Goren et al. 2021) and allows inferring the uplift history for any value of n, under two assumptions: First, if > 1, U(t) is a monotonically increasing staircase function and if < 1, U is decreasing. Second, all the knickpoints are preserved within the time resolved by the model. The model is based on the block uplift assumption, 230 whereby a suite of basins and tributaries experience and respond to the same time-dependent tectonic history U(t). The model infers the best fit U(t) based on the long profiles of the tributaries and basins.
Changes in U through time emerge as a series of knickpoints with elevations and χ values, ( 1 , χ 1 ), ( 2 , χ 2 ), … ( −1 , χ −1 ), which are duplicated across the tributaries and basins. The basin outlets are at ( 0 = 0, χ 0 = 0) and the highest χ channel head is identified with ( , χ = χ ). The knickpoints are used to divide the χ-z space into segments. Segment j, between 235 (χ −1 , χ ), is characterized by a uniform steepness index that shaped the river profile during time interval ( −1 , ), where time tj is the age of knickpoint j. Knickpoint ages can be constrained from equation (15), and the uplift rate responsible for the formation of each knickpoint can be constrained based on the steepness index below the knickpoints by using equation (7). Consequently, a full uplift history, with discrete step-increase can be derived.
A difficulty may arise because tj in equation (15) and Uj in equation (8) depend on the erodibility, K, whose value is 240 commonly poorly constrained. Thus, following Goren et al. (2014), we present a K-independent version for the knickpoint age and uplift rate. Scaling equations (15)  Equations (26-27) produce a non-dimensional uplift rate history, ( * , * ), without any prior information on K, as long as it is spatially uniform.
We propose the following steps for the application of the inverse model. First, the data of basins and tributaries is considered in the χ-z domain, where it is divided into q segments along the χ space, χ ( = 0, 1, 2, … ). The division points are considered to be slope-break knickpoints that formed in response to a step-increase in uplift rate. The scaled age of the 250 knickpoints is solved based on equation (26) as * = χ . Second, linear regression is applied in the χ-z domain, independently for each segment. The slope of the regression is identified as s_ , from which * can be derived based on equation (27). Third, conversion from ( * , * ) to a dimensional history ( , ) by solving equations (26-27), after K and n are independently constrained (e.g. Dibiase et al., 2010;Ma et al., 2020).
The first step of dividing the χ-z domain into segments calls for some consideration. First, calculating the χ value requires 255 calibrating for the concavity, 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., Hergartena et al., 2016;Gaillton et al., 2021) that finds the m/n that minimizes the scatter in the χ-z domain. Second, segment division should ideally be based on division points that represent true slope-break knickpoints. Many algorithms have been previously proposed to identify slope-break knickpoints (e.g., Mudd et al., 2014). Here, we suggest a different approach that rely on the simplicity and efficiency of the 260 inverse model. The inversion procedure could be run many times, while choosing the division points randomly. Inversion results could be evaluated by comparing the measured profiles and the profile predicted by our forward model. 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. Here, we consider a misfit function that penalizes models with more knickpoints (more parameters) for their excess complexity: 265 where zi and ̃ are the measured and predicted elevations at pixel i, respectively. N is the total number of data along the river long profiles, and = is the number of division points, or the number of parameters.
To demonstrate the applicability of the inverse method, we use a low resolution numerical model (which suffers from numerical diffusion) to generate ten river profiles with variable channel length and drainage area distribution. These rivers 270 respond to the same uplift rate history, with two step-increases in the rate forming two knickpoints in each profile.
Knickpoints do not merge over the timeframe of U(t) application (Figure 5a and b). To artificially increase the noise in the data, the elevations are perturbed by random errors (equation 29): ̂(perturbed) = −1 + ( +1 − −1 ) * rand(1), where rand(1) is a random number between 0 and 1. Inversion in applied to the data with 1-6 division points. For each 275 number of division points, 5000 realizations of the inversion model were performed with different random position of the division points. Figure 5c shows the minimal misfit (equation 28) achieved for each number of division points, indicating that the best fit solution has division points corresponding to the two knickpoints. Figure 5d compares between the applied and inferred histories, showing that the 2 division points inversion correctly infers the applied history.

Discussion and conclusion 280
The current analysis explores river long profile evolution in response to temporal step-changes in the tectonic rock uplift rate U(t) and a non-unity slope exponent, leading to consuming channel segments  and merging knickpoints. The approach we adopt here, of resolving knickpoint kinematics in a Lagrangian frame of reference, allows us to constrain the timing of knickpoint merging and the elevation of knickpoint before and after merging. The finding that despite channel reach consumption, knickpoint celerity depends only on the channel steepness below and above the 285 knickpoint, allow 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.
Analytic solutions of long profile evolution can significantly expedite forward and inverse tectonicfluvial landscape evolution. However, so far, analytic solutions were used in such models only under linear assumption (Pritchard et al., 2009;Fox et al. 2014;Goren et al., 2014Goren et al., , 2021Rudge et al., 2015;Steer et al. 2021). The simple analytic derivation that we 290 present here can expand the domain of parameters for which analytic solutions are used in such models, by including new geomorphic scenarios with ≠ 1 (with the restriction of increasing U(t) for > 1 and decreasing U(t) for < 1). For example, inverse models that are based on Bayesian statistics (Fox et al., 2015;Gallen and Ferná ndez-Blanco, 2021), which have gain recent popularity could become significantly more efficient and accurate when the forward model is represented with an analytic solution. 295 The knickpoint merging analysis further emphasizes a critical property of the links between tectonic and long profile evolution when ≠ 1. Each tectonic history is associated with a single, well-defined river profile at any given time.
However, any particular river profile could be generated by infinitely many tectonic histories. All histories except for one lead to knickpoint merging dynamics. The linear inverse model that we develop here finds the single history for which all knickpoint are preserved. While this inverse approach is highly restrictive, it finds the correct solution when only a single 300 knickpoint exists 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 history that shaped the fluvial landscape.

Code and data availability
This study has no complex codes or data sharing issue, because all the figures can be re-produced by solving the related equations. 335

Author contribution
Both YW and LG developed these models, conducted the study, and designed structure of the manuscript. YW solved these equations, and wrote the manuscript. LG wrote the 1-D numerical codes and helped to test the analytical derivations. The manuscript was mostly revised, polished, and improved by LG. DZ and HZ provided valuable suggestions and made some revisions. 340 Gallen for constructive instructions on an earlier version of this manuscript. This work was supported by the National Science Foundation of China (41802227).