# Using Disorder Metrics to Distinguish Discharge-Driven From Drainage Area-Driven Incision and Quantify Deviations in Channel Steepness

Marina Ruiz Sánchez-oro, Simon M Mudd, Boris Gailleton

► To cite this version:

Marina Ruiz Sánchez-oro, Simon M Mudd, Boris Gailleton. Using Disorder Metrics to Distinguish Discharge-Driven From Drainage Area-Driven Incision and Quantify Deviations in Channel Steepness. Journal of Geophysical Research: Earth Surface, 2024, 129 (9), pp.e2023JF007553. 10.1029/2023jf007553 . insu-05322412

HAL Id: insu-05322412

https://insu.hal.science/insu-05322412v1

Submitted on 20 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Earth Surface

## RESEARCH ARTICLE
10.1029/2023JF007553

**Key Points:**

* Discharge-driven incision can be identified a *posteriori* using channel disorder metrics in modeled landscapes but not in natural topography

* The choice of concavity index ($\theta$) can distort the normalized channel steepness index more than the choice of incision type

* Lithology clouds the rainfall signal in channel profile disorder, needing erosion rates to fully understand the incision mechanism

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Ruiz Sánchez-Oro,
marina.ruiz.so@ed.ac.uk

**Citation:**

Ruiz Sánchez-Oro, M., Mudd, S. M., & Gailleton, B. (2024). Using disorder metrics to distinguish discharge-driven from drainage area-driven incision and quantify deviations in channel steepness. *Journal of Geophysical Research: Earth Surface*, 129, e2023JF007553. https://doi.org/10.1029/2023JF007553

Received 21 NOV 2023
Accepted 7 AUG 2024

**Author Contributions:**

**Conceptualization:** Marina Ruiz Sánchez-Oro, Simon M. Mudd
**Data curation:** Marina Ruiz Sánchez-Oro
**Funding acquisition:** Simon M. Mudd
**Investigation:** Marina Ruiz Sánchez-Oro
**Methodology:** Marina Ruiz Sánchez-Oro, Simon M. Mudd, Boris Gailleton
**Project administration:** Simon M. Mudd
**Software:** Marina Ruiz Sánchez-Oro, Simon M. Mudd, Boris Gailleton
**Supervision:** Simon M. Mudd
**Validation:** Marina Ruiz Sánchez-Oro
**Visualization:** Marina Ruiz Sánchez-Oro

© 2024. The Author(s).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Using Disorder Metrics to Distinguish Discharge-Driven From Drainage Area-Driven Incision and Quantify Deviations in Channel Steepness

Marina Ruiz Sánchez-Oro¹ ID, Simon M. Mudd¹ ID, and Boris Gailleton² ID

¹School of GeoSciences, University of Edinburgh, Edinburgh, UK, ²CNRS - Geosciences Rennes, Université de Rennes, Rennes, France

**Abstract** The rate of channel incision in bedrock rivers is often described using a power law relationship that scales erosion with drainage area. However, erosion in landscapes that experience strong rainfall gradients may be better described by discharge instead of drainage area. In this study, we test if these two end member scenarios result in identifiable topographic signatures in both idealized numerical simulations and in natural landscapes. We find that in simulations using homogeneous lithology, we can differentiate a *posteriori* between drainage area and discharge-driven incision scenarios by quantifying the relative disorder of channel profiles, as measured by how well tributary profiles mimic both the main stem channel and each other. The more heterogeneous the landscape becomes, the harder it proves to identify the disorder signatures of the end member incision rules. We then apply these indicators to natural landscapes, and find, among eight test areas, no clear topographic signal that allows us to conclude a discharge or area-driven incision rule is more appropriate. We then quantify the distortion in the channel steepness index induced by changing the incision rule. Distortion in the channel steepness index can also be driven by changes to the assumed reference concavity index, and we find that distortions in the normalized channel steepness index, frequently used as a proxy for erosion rates, is more sensitive to changes in the concavity index than to changes in the assumed incision rule. This makes it a priority to optimize the concavity index even under an unknown incision mechanism.

**Plain Language Summary** Rivers erode into mountains as a result of both the sediment transported as the water flows downstream and the amount of water that the river transports. The amount of rainfall that each part of the river receives affects how much the channel cuts into the rock. In this study, we assess whether it is possible to identify rivers where rainfall has left a measurable imprint on the erosion history. Through computer simulations, we measure the river fingerprint in the landscape through calculations involving how tributary and main channel slopes compare, and how quickly a channel steepens as it travels away from the headwaters. We extract these fingerprints and search for them in the more complex natural landscapes, measuring how much they change under heavy rainfall. We find that these fingerprints are camouflaged by factors such as changes in rock types, making it a challenge to identify them without field observations.

## 1. Introduction

Physical intuition suggests that, if all other factors are equal, a steeper river will erode faster than a gentler one. This basic relationship has been proposed by geomorphologists for over a century (Gilbert, 1877). It is unusual, however, to find two channels identical in their properties with the exception of their gradient. Headwaters are, for example, frequently steeper than the downstream rivers they feed. In the early 1960s geomorphologists realized that gradient could be related to drainage area in a power law with a negative exponent (Hack, 1960; Morisawa, 1962). This basic relationship was formalized by Flint (1974):

$$S = k_s A^{-\theta} \quad (1)$$

where the concavity index, $\theta$, describes how fast the gradient of the river changes downstream and the constant $k_s$, the channel steepness index, describes the dependence of gradient normalized for drainage area. We can further set $\theta$ to a fixed reference value ($\theta_{ref}$), after which we denote the channel steepness index with $k_{sn}$. This normalization allows us to compare the relative steepness of rivers with different drainage areas (Flint, 1974).
**Writing – original draft:** Marina Ruiz Sánchez-Oro

**Writing – review & editing:** Simon M. Mudd, Boris Gailleton

Numerous studies have found that this relative steepness ($k_{sn}$) is positively correlated with measured erosion rates in upland landscapes (Adams et al., 2020; DiBiase et al., 2010; Gailleton et al., 2021; Harel et al., 2016; Harries et al., 2021; Kirby & Whipple, 2012; Peifer et al., 2021; Wobus et al., 2006). In regions without data on erosion rates we might therefore use $k_{sn}$ as a proxy for erosion rate (Kirby & Whipple, 2012). However, this is not quite as straightforward as it sounds. If we return to physical intuition, consider two channels with all properties equal apart from the amount of water they convey (quantified as, e.g., their mean annual discharge or some other statistical representation of runoff). Many authors have proposed relationships between channel incision and the physical properties of bedrock rivers, and these proposals include the influence of sediment supply (dependent on A), sediment transport capacity (dependent on hydraulic conditions and gradient), shear stress (dependent on flow depth, and thus hydraulic conditions), stream power (again, hydraulic conditions), and thresholding behavior of all the above factors (Baynes et al., 2020; Finnegan et al., 2007; Howard, 1987; Johnson & Whipple, 2010; Wobus et al., 2010). This means that computation of $k_{sn}$ based on drainage area rather than hydraulic conditions (e.g., discharge) may not represent the incision process.

An additional factor affecting channel steepness and river erosion is lithology. Recently, the role of lithology on the evolution of landscapes has been the focus of many studies, especially related to normalized channel steepness index. Gallen (2018) shows that lithological heterogeneities can be the dominant controlling mechanism in certain tectonically quiescent areas where 5-fold erodibility differences between rock types would overwhelm any present signal, such as a climatic imprint. Leonard et al. (2023a, 2023b) highlight the need to control for rock type when relating erosion rates with climatic variables such as mean annual rainfall. Under strong rainfall gradients, Be-10 erosion rates show different relations with fluvial relief (from normalized channel steepness index) depending on the underlying lithology. Thus, controlling for this when inferring erosion rates becomes particularly important when basins are subject to rainfall gradients. Marder and Gallen (2023) show further evidence for the need to control for bedrock and constrain erosion rates to detect the climatic signal through fluvial relief using Be-10 erosion rates, which can otherwise not be disentangled.

One reason erosion rate proxies have tended toward using an area-based $k_{sn}$ is because drainage area is trivial to extract from topographic data. Discharge ($Q$) records are not always easy to obtain, and gauging stations are at points rather than distributed throughout the landscape. However, various global data sets, for example, TRMM (Kummerow et al., 2000) and GPM (Skofronick-Jackson et al., 2017) have made it relatively simple to estimate and aggregate precipitation over a basin, over a variety of timescales. This means that, given a linear relationship between aggregate precipitation and discharge, based simply on a mass balance, it is now quite easy to obtain an estimate of a discharge proxy. Note that because we are not including water loss through evapotranspiration and infiltration, these values may differ from the true runoff by a factor hydrologists call the runoff ratio (Rossi et al., 2016).

These precipitation data sets enable us to calculate the channel steepness index, $k_{sn}$, based on effective discharge rather than drainage area. And indeed, a number of recent authors have taken this approach (e.g., Adams et al., 2020; Babault et al., 2018; Harries et al., 2021; Leonard & Whipple, 2021; Leonard et al., 2023b). If precipitation rates are uniform across a catchment, the drainage area-gradient relationship will have the same spatial pattern as the drainage area-discharge relationship. But rainfall can be influenced by mountains (Anders et al., 2006; Bookhagen & Burbank, 2006; Bookhagen & Strecker, 2008; Craddock et al., 2007; Gasparini & Whipple, 2014; Roe et al., 2002), meaning that the patterns of $k_{sn}$ might be different if one uses A or a more direct estimate of $Q$ that incorporates spatially varied precipitation.

One might assume that using $Q$ is always “better” than A. But there are some reasons why erosion rates, and gradients, might be more sensitive to drainage area than discharge. The main reason for this is that water does not erode the bed of rivers, sediment does. And the amount of sediment fluxing through any part of the channel in a steadily uplifting mountain range should depend on drainage area and not discharge. In addition, rivers transporting gravel will alter their geometry, for example, their width (Dunne & Jerolmack, 2020; Pfeiffer et al., 2017; Phillips & Jerolmack, 2016), to accommodate sediment supply and this could cause a damping effect on the relationship between discharge and erosion rates. So although it intuitively might make sense to always use discharge-based calculations of $k_{sn}$, we are not, at present, certain if this is better than a calculation using A.

Any proposed erosion rule, be it area or discharge driven, can be transformed into a prediction of topography. For example, the most basic erosion law incorporating gradient and area takes the form (Howard & Kerby, 1983):
$$E = KA^mS^n \qquad (2)$$

where K is an erodibility and m and n are empirical coefficients. In a variety of derivations of Equation 2 the A term includes specific power law relationships between area and both discharge and channel width (e.g., Whipple & Tucker, 1999). If one rearranges Equation 2 to isolate S and compares the result with Equation 1, one can see that, if Equation 2 is correct, $k_{sn} = (E/K)^{1/n}$ and $\theta = m/n$. Thus Equation 2 predicts a power law relationship between gradient and drainage area where erosion rates are invariant in time and space. It can also be shown that this power law relationship holds in segments of constant erosion that move upstream in transiently eroding landscapes (L. Royden & Perron, 2013). Various threshold models have been proposed as more complex versions of Equation 2, but much of their relevant behavior can be captured by altering the n exponent (Gasparini & Brandon, 2011). While this approach is not without controversy (Lague, 2014), it is at least not clearly falsified by relationships between topographic data and measured erosion rates.

Such topographic predictions suggest a basic test: if there are strong gradients in precipitation, we expect the topographic outcomes of incision rules that are driven by either A or Q to differ. In this contribution we first explore the question: if we know erosion rates are driven by either A or Q, could we tell the difference just based on topography? We use numerical experiments to answer this question. We explore the extent to which heterogeneity in uplift rates and erodibilities can cloud this signal. Finally, we explore real landscapes to see if we can find locations, based on our proposed metrics, where it is clear that an area-based or discharge-based calculation of $k_{sn}$ is more appropriate.

## 2. Methods

Our study includes three components. First, we perform a series of numerical experiments using a simple landscape evolution model with different imposed incision rules, and in addition alter other model parameters such as the rainfall gradient. The aim of these simulations is to produce landscapes under idealized and controlled conditions against which metrics for determining the most likely incision rule may be tested. We then develop metrics that allow us to test if a particular incision rule better describes observed topography. Finally, we deploy these metrics on real landscapes.

### 2.1. Numerical Simulations

To simulate landscapes where bedrock channels incise through an uplifting landscape we must make choices about the incision rule. For the purposes of simplicity, we use the basic form of the stream power incision model (Equation 2), which can emulate different incision mechanisms by altering the n exponent (Gasparini & Brandon, 2011).

Equation 2 is not influenced by discharge. For discharge-based incision, we follow other authors (e.g., Adams et al., 2020) and replace drainage area with a proxy for discharge, which we compute with the substitution $Q = A \times \text{Rainfall}$ where the rainfall is converted into runoff and accumulated downstream, yielding:

$$E = K_{lp}Q^mS^n \qquad (3)$$

In reality, the discharge (and sometimes the exponents m and n) will be modulated by other features such as evapotranspiration and infiltration rates. Other climate-independent variable such as bedrock erodibility or channel geometry are subsumed into the parameter $K_{lp}$ (Adams et al., 2020).

The model then simulates topographic evolution with a simple mass balance:

$$\frac{\partial z}{\partial t} = U - E \qquad (4)$$

where U is the uplift rate and the erosion rate, E, is solved by either Equation 2 or 3. Simulations are performed using the Fastscape framework (Bovy & Lange, 2023). Fastscape uses the methods developed by Braun and Willett (2013) that takes advantage of graph theory to efficiently solve Equation 2 (or 3) with an implicit finite different scheme.
In order to implement Equation 3, we must assign a precipitation pattern. In our experiments, we explore different precipitation gradients. Models do exist where precipitation depends on elevation and prevailing wind direction, for example, the model of Smith and Barstad (2004) which was used by Han et al. (2015). Such models involve multiple parameters, and to obtain the desired precipitation gradient one must engage in fine tuning each of these parameters. For consistency across experiments, we decided instead to increase precipitation linearly as a function of distance from the mountain front, with precipitation gradients comparable to real orographic rainfall gradients. Leonard and Whipple (2021) found that a linear rainfall gradient can be a good approximation to studying orographic precipitation patterns in natural landscapes.

We impose rainfall gradients that will mimic natural orographic patterns found in natural landscapes. In the driest scenario, we impose a rainfall gradient across the simulated mountain range of 1 m/yr (i.e., the precipitation is 1 m/yr greater at the peak rainfall location than at the minimum rainfall location). In our wettest scenario, the increase in rainfall reaches 10 m/yr across the simulation domain, corresponding to the highest rainfall achieved in the Bhutan Himalayas (Anders et al., 2006; Grujic et al., 2018).

## 2.1.1. Homogeneous Lithologies

We run experiments on a regular grid ($\Delta x = 30$ m) over a mountain range that is 15 km by 30 km (see Table S2 in Supporting Information S1 for full parameter details). This size was selected as a compromise between the number of basins that could be formed during a simulation and computational expense. We begin each landscape with a random surface generated using the diamond square algorithm (Fournier et al., 1982; Perron & Royden, 2013) with noise ranging from 0 to 1 m. We choose this initial condition as it produces a greater variety of channel network structures than the more widely used white noise. Each simulation is run to steady state, where the change in erosion is balanced by uplift.

For a given uplift and precipitation gradient, we calculate the landscape resulting from each of the two incision laws as described by Equations 2 and 3. Initial experiments use uniform erodibility to simulate homogeneous lithology. The boundary conditions include a fixed elevation on the east and west edges (which allow flux to exit the model) and periodic boundaries at the north and south of the model domain. The resulting mountain range emerges in the north-south direction.

Within the discharge-driven model, we set up one simulation for each rainfall scenario, starting from a gradient of 0 m/yr to 10 m/yr (increasing from east to west in the simulation domain), for a total of 11 simulations. The precipitation runs explore how different rainfall amounts affect river incision mechanisms and whether larger gradients generate a stronger signal in the landscape.

Both the erodibility coefficient and the uplift rate are kept constant across all simulations at $3 \times 10^{-8} m^{1-2m}yr^{-1}$ and $1 \times 10^{-5} m yr^{-1}$ respectively, yielding landscape reliefs within the ranges of those found on Earth. Similarly, *n* is chosen as 1 and *m* as 0.45 to keep the *m/n* ratio equal to 0.45, which is the central tendency of the concavity index across a large number of global landscapes (Gailleton et al., 2021; Kirby & Whipple, 2012; Tucker & Whipple, 2002). When studying the channel steepness index across multiple basins, we set a reference value, $\theta_{ref}$ (Gailleton et al., 2021), of 0.45 to establish a comparison between basins.

## 2.1.2. Heterogeneous Lithologies

To study the role of lithology under the signature of rainfall patterns we run a set of simulations with spatially varying values of erodibility, *K*. These simulations are designed to be closer to natural landscapes, where the lithological composition is more complex.

Within this study, we choose two lithological units reflecting the properties of a harder and a softer rock, with erodibility values ranging between $10^{-7} < K < 5 \times 10^{-8} m^{1-2m}yr^{-1}$. This is in line with Forte et al. (2016), Bernard et al. (2021), and Peifer et al. (2021), who determined that hard rocks can be related by a factor from 2 to 10 in erodibility to softer rocks.

We simulate three heterogeneous lithologies scenarios: (a) dense blob (Figure 5a), (b) sparse blob (Figure 5b), and (c) striped lithologies (Figure 5c).
<caption>**Table 1**</caption>

Geographical Areas Chosen Along With the Number of Selected Basins in Each Regions, Area, and Rainfall Range Across the Basin, From the Outlet to the Headwaters of the Catchment

<table><thead><tr><td>Location</td><td>N basins</td><td>Area (km²)</td><td>Precipitation range (m/yr)</td></tr></thead><tbody><tr><td>Andes, Southern Perú</td><td>5</td><td>29,979</td><td>0.687–3.983</td></tr><tr><td>Andes, Northern Argentina</td><td>7</td><td>5,932</td><td>0.045–0.100</td></tr><tr><td>North Qinling Mts, China</td><td>14</td><td>30,832</td><td>0.734–0.938</td></tr><tr><td>Kaçkar Mts, Turkey</td><td>8</td><td>4,279</td><td>0.784–1.673</td></tr><tr><td>Colorado Front Range, USA</td><td>7</td><td>9,282</td><td>3.726–4.220</td></tr><tr><td>Alburz Mts, Iran</td><td>7</td><td>8,167</td><td>0.357–0.849</td></tr><tr><td>Massif Central, France</td><td>5</td><td>12,745</td><td>0.885–1.111</td></tr><tr><td>Pyrénées, Spain-France</td><td>5</td><td>6,632</td><td>0.851–1.024</td></tr></tbody></table>

Note. We choose a varied range of area sizes and precipitation gradients to study prevalent trends across regions. The data is extracted from the 30 m Copernicus DEM (European Space Agency, 2020) and NASA's Global Precipitation Measurement Mission (GPM) (Skofronick-Jackson et al., 2017).

1. Dense blob lithology: small round hard rock blobs ($K = 5 \times 10^{-8} m^{1-2m}yr^{-1}$) of a few square meters in diameter, generated using Perlin noise (Perlin, 1985). The rest of the simulation domain transitions smoothly to soft rock blobs of the same order of magnitude. This is an example of fast-changing lithological units, such as the study area in the Colorado Front Range.

2. Sparse blob lithology: 4 large round hard rock blobs ($K = 5 \times 10^{-8} m^{1-2m}yr^{-1}$) evenly distributed on the landscape domain. They have a diagonal dipole configuration with high K at the top left and bottom right, and low K at the opposite corners. There is a transition area of medium K rocks in the center of the domain of around 15 km. This scenario shows slowly changing rock hardness, such as the chosen region in the Argentinian Andes.

3. Striped lithology: straight vertical band of hard rock ($K = 5 \times 10^{-8} m^{1-2m}yr^{-1}$) in the center of our simulated mountain range, surrounded by soft rock ($K = 1 \times 10^{-7} m^{1-2m}yr^{-1}$) to emulate cases like the Pyrénées. The width of the hard rock band is 7 km and the two soft rock bands are 3.5 km wide.

### 2.1.3. Varying the Simulation Concavity Index, $\theta$ (i.e., m/n Ratio)

We assess the impact of running the homogeneous lithology simulation with a further two choices of m/n: 0.35 and 0.55. The rationale for these experiments is to determine if changes to the imposed m/n ratio causes distortions to $k_{sn}$ of the same magnitude as those induced by changing the incision rule. The broad effect on the landscapes with different values of m/n compared to the base case of m/n = 0.45 is that at 0.35 the landscape is smoother with lower relief and higher drainage density with sinuous tributaries, whereas at m/n = 0.55 the landscape has sharper features, with higher relief and lower drainage density, forming straighter tributaries. For detailed results of these simulations, see Text S6, Figures S23–S26, Tables S8 and S9 in Supporting Information S1.

### 2.1.4. Natural Landscapes

We have chosen eight natural landscapes across the globe for this analysis on the basis of precipitation gradients. However, the areas experiencing the largest rainfall gradients are often active orogens and/or experience intense volcanic activity, such as the Himalayas (Taylor & Yin, 2009) or Hawai'i (Moore & Clague, 1992). To introduce some controls over the natural forcings, we decided to favor places with lower lithological and tectonic complexity over those with the largest orographic rainfall gradients.

We analyze basins with a minimum drainage area of $10^7 m^2$. We incorporate all channel pixels within tributaries that have a source area greater than 1.35 $km^2$ (which corresponds to 1,500 pixels in topographic data with 30 m grid spacing). Below is a description of the lithological and tectonic setting of each of the chosen areas. Table 1 summarizes the rainfall ranges and size of the chosen study locations. The natural rainfall gradients use the mean annual precipitation calculated from NASA's Global Precipitation Measurement Mission (GPM) (Skofronick-Jackson et al., 2017) between the years 2001–2021 over the extent of each region.
Eastern Cordillera, Northwest Argentina (Figure 10): this area sits west of the Puna Plateau at 4,000 m, with surrounding peaks above 5,000 m. The main lithological units are sedimentary layers of siltstone and mudstone, with localized metamorphic units (Aceñolaza & Miller, 1982; Escayola et al., 2011; Hartmann & Moosdorf, 2012; Rahl et al., 2018). As part of the Andes, tectonics are controlled by the subduction of the Nazca plate under the South American plate, with no active faults through the study area (Styron & Pagani, 2020).

Subandes, Southwest Perú (Figure S13 in Supporting Information S1): the study area lies in the Madre de Dios basin, where the dominant lithological units are siliciclastic sediments and plutonic rocks. The area has peaks above 5,000 m (Aceñolaza & Miller, 1982; Hartmann & Moosdorf, 2012; McQuarrie et al., 2005; Safran et al., 2005; Whipple & Gasparini, 2014). Similar as in Argentina, the tectonics are controlled by the Nazca plate, with active thrust faults at the foreland of the mountain range (Styron & Pagani, 2020).

Pyrénées, Spain/France (Figure S14 in Supporting Information S1): created from the subduction of the Iberian plate under the European plate, the Pyrénées has peaks above 3,000 m. It is bound by the Aquitaine basin in the north and the Ebro basin in the south, where the sedimentary sections accumulate. Metamorphic rocks and interspersed volcanics dominate the central area of the mountain belt (Bernard et al., 2019; Hartmann & Moosdorf, 2012; Morris & Yell, 1998).

Massif Central, France (Figure S15 in Supporting Information S1): located west of the Alps, it is a region of anomalously high elevation in its context, with peaks of 1,700 m. The Cevenne fault goes through the area on the northeast and the relief is otherwise controlled by the Rhône river (Olivetti et al., 2016; Sanchis & Séranne, 2000). The area is dominated by metamorphic and plutonic rocks in the north and central regions, with sedimentary rocks elsewhere (Hartmann & Moosdorf, 2012).

Alburz Mountains, Iran (Figure S16 in Supporting Information S1): bounded by the Caspian sea in the north, this mountain range reaches elevations above 5,000 m within 100 km of the sea. Several thrust faults cross the region, with a dominant volcanic lithology mixed with smaller areas of sedimentary and some metamorphic rocks (Allen et al., 2003; Hartmann & Moosdorf, 2012; Rezaeian et al., 2012).

Qinling Mountains, China (Figure S17 in Supporting Information S1): part of the Central China orogenic belt, it is bounded in the north by the Lingbao-Lushan-Wuyang thrust fault and south by the Luonan-Luanchuan fault (Dong et al., 2015, 2022). The thrust fault separates areas of unconsolidated sediment with metamorphic and plutonic rocks, mixed with some harder sedimentary regions (Hartmann & Moosdorf, 2012). The elevation reaches 3,000 m.

Kaçkar Mountains, Turkey (Figure S18 in Supporting Information S1): part of the Eastern Pontides, they are located north of the North Anatolian Fault, draining to the eastern Black Sea (Espurt et al., 2014; Hippolyte et al., 2017). The highest peaks reach over 3,500 m. The area is dominated by plutonic and volcanic rocks (Hartmann & Moosdorf, 2012).

Rocky Mountains, CO, USA (Figure S19 in Supporting Information S1): the study area sits in the Southern Rocky Mountains, with an average height of 3,000 m, between the Great Plains on the east and the Colorado Plateau (≈ 2,000 m) on the west. The area is traversed by Laramide thrust faults and Quaternary faults (Karlstrom et al., 2012, 2022). The study area is mostly composed of metamorphic and plutonic rocks, with sedimentary rocks on the Great Plains (Hartmann & Moosdorf, 2012).

## 2.2. Metrics to Quantify Topographic Outcomes

Our models simulate river incision, and so we use river profiles to explore topographic outcomes of simulations. Because gradients should scale by either A or Q (depending on the incision rule) we use a coordinate transformation, first proposed by (L. H. Royden et al., 2000), that integrates either A or Q along the river profile:

$$ \chi_A = \int_{x_0}^{x} \left( \frac{A_0}{A(x)} \right)^{\theta} dx \qquad (5) $$

$$ \chi_Q = \int_{x_0}^{x} \left( \frac{Q_0}{Q(x)} \right)^{\theta} dx \qquad (6) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of how $\chi$ profiles and $k_{sn}$ change based on the choice of incision scenario and concavity index values. We present $\chi$ profiles resulting from Q-driven simulations with a precipitation gradient of 10 m/yr across the model domain (increasing east to west), under initial $m/n = 0.35$. The basin shown is located on the drier side of the domain, draining to the east. The expected true solution is shown in (a), with perfect collinearity under $\chi_Q$ and $\theta_{best} = 0.35$, matching the simulation parameters. Panel (b) shows a slightly more collinear profile, in this case only having the incorrect incision. Panel (c) shows distortions using both the incorrect incision rule and the wrong $\theta$. Panel (d) captures the correct incision (discharge) but uses $\theta = 0.45$, where we see the largest increase in disorder and changes in $k_{sn}$.</figcaption></figure>

This transformation has various useful features. The gradient in $\chi_A$-elevation space is equivalent to $k_s$ if $A_0 = 1 \text{ m}^2$ (e.g., Mudd et al., 2014; Perron & Royden, 2013; L. Royden & Perron, 2013), and the gradient in $\chi_Q$-elevation space is equivalent to the metric $k_{s-q}$ where gradient is scaled by Q instead of A in the form:

$$S = k_{s-q} Q^{-\theta} \qquad (7)$$

when $Q_0 = 1$ unity in the units of Q used to calculate both $\chi_Q$ and $k_{s-q}$ (e.g., Adams et al., 2020; D'Arcy & Whittaker, 2014; Harries et al., 2021; Leonard & Whipple, 2021; Leonard et al., 2023b). These steepness metrics can be used both in steady state and transient landscapes (L. Royden & Perron, 2013). Another advantage of using the $\chi$ transformation is that tributaries to the main channel at the same elevation yield the same $\chi$ value regardless of their drainage area, assuming they have the same uplift and erosion rates. Given a landscape in steady state, the main channel and its tributaries should follow the same linear relationship on a $\chi-z$ plot, assuming the same erosion and uplift rates (Perron & Royden, 2013) and the optimal concavity index value.

We exploit this latter feature in our efforts to discriminate, topographically, between incision rules. Regardless of the incision rule, selection of the incorrect value of $\theta$ will result in tributaries that are not collinear, introducing distortions in channel steepness index (Gailleton et al., 2021; Goren et al., 2014; Harries et al., 2021; Hergarten & Robl, 2022; Mudd et al., 2018; Perron & Royden, 2013). Computing $\chi_A$ or $\chi_Q$ also affects the channel steepness index values and the spread of the data, which can lead to different patterns of the channel steepness index and spreading the data in $\chi$ space (Figure 1). Because $k_{sn}$ values are used to infer relative erosion rates across tectonically active regions, distortions to the interpretations of topographic patterns (e.g., Gailleton et al., 2021).
We quantify the spread of the $\chi-z$ profiles using a disorder metric, first proposed by Goren et al. (2014) and further developed by Hergarten et al. (2016) and Mudd et al. (2018). One begins by ranking every point in the channel network by increasing elevation, and then checks to see if the associated $\chi$ coordinates are similarly ranked (or not):

$$R = \sum_{i=1}^{N} |\chi_{s,i+1} - \chi_{s,i}|, \qquad (8)$$

where the subscript $s,i$ represents the $i^{th}$ $\chi$ coordinate that has been sorted by its elevation ($\chi_{s,i}$). This sum, $R$, is minimal if elevation and $\chi$ are related monotonically. However, it scales with the absolute values of $\chi$, which are sensitive to the concavity index (see Equations 5 and 6), so following Hergarten et al. (2016) we scale the disorder metric, $D$, by the maximum value of $\chi$ in the tributary network ($\chi_{max}$):

$$D(\theta) = \frac{1}{\chi_{\max}(\theta)} \left( \sum_{i=1}^{N} |\chi_{s,i+1}(\theta) - \chi_{s,i}(\theta)| - \chi_{\max}(\theta) \right). \qquad (9)$$

We use the method of Mudd et al. (2018) to constrain uncertainty of this metric by creating subset networks formed from the trunk stream and every possible combination of three tributaries in a particular basin. This creates a population of $D$ values for a given basin from which a median and interquartile range may be reported.

We normalize the disorder values across all tributary combinations, obtaining:

$$D^*(\theta) = \frac{D(\theta)}{D_{\max}(\theta)} \qquad (10)$$

where $D_\theta$ is the disorder for each tributary combination and $D_{\max}(\theta)$ is the maximum disorder over all combinations.

In our analysis, we aim to decipher whether we can identify a *posteriori* the signal from a rainfall gradient in a landscape through topographic metrics. To quantify this, we focus on the effect of rainfall in the $\chi-z$ profiles and in the disorder metrics. To emulate what the analysis would look like if we did not know the incision rule, we calculate the $\chi$ profiles in two ways for each simulation scenario—regardless of what the actual imposed incision rule is for a given numerical experiment. We calculate these metrics for each basin in the simulation draining to the edge and reaching the main drainage divide (Figure 2). For each of the basins simulated in each model scenario, we calculate the following:

1. $\chi_A$: assumes a drainage area-driven incision (Equation 5).

2. $\chi_Q$: assumes a discharge-driven incision (Equation 6).

For each of the $\chi$ cases, we pick the optimal concavity, $\theta_{best}$, by minimizing the disorder metric, $D(\theta)$ (Equation 9). This leads to the most collinear river profile configuration (Gailleton et al., 2021; Mudd et al., 2018).

## 2.3. Statistical Analysis

For each incision scenario, we calculate minimum normalized disorder values ($D^*(\theta)$) corresponding to each of the basins under each of the four incision and $\chi$ scenarios (Table 2). Each basin has a disorder parabola from which its $\theta_{best}$ is derived. The point of minimum disorder happens at $\theta_{best}$. We calculate the disorder parabola assuming both A-driven and Q-driven incision, regardless of the “true” model incision, obtaining values for $D_A^*$ and $D_Q^*$ at their corresponding $\theta_{best}$. We calculate the absolute error between the minimum disorder values:

$$\Delta D^* = D_{\chi_A}^* - D_{\chi_Q}^* \qquad (11)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** River networks generated with the numerical model under two discharge-driven incision scenarios, with rainfall ranges across the simulation domain between (a) 1–2 m/yr and (b) 1–11 m/yr in the east-west direction. The smaller rainfall ranges in (a) lead to higher relief, steeper channels and a symmetric drainage divide. The higher rainfall in (b) leads to lower relief and a more sinuous and asymmetric drainage divide.</figcaption></figure>

This allows us to quantify whether $\Delta D^*$ can be used to determine if A and Q-driven basins are statistically distinguishable.

Since we are dealing with non-parametric distributions, we take the median of $\Delta D^*$ to quantify whether the A-based incision models may be distinguished from the Q-based incision models. We represent the distributions with kernel density estimates (KDE) (Cox, 2007; Silverman, 1998). If a percentage smaller than 5% is shared between the two distributions we consider them to be distinguishable from each other with 95% confidence.

## 2.4. Measuring the Effects of Rainfall and $\theta$ on $k_{sn}$

In our numerical landscapes, we impose an erosion rule and then quantify the extent to which profiles are disordered when applying the appropriate or incorrect $\chi$ transformation (i.e., using either A or Q). In real landscapes, however, we can only infer which of the two $\chi$ transformations is correct based upon relative disorder, and indeed neither may be correct.

Another way of understanding how much uncertainty is introduced by not knowing the appropriate $\chi$ transformation is to quantify the degree to which $k_{sn}$ varies under different $\chi$ transformations. We do this using distortion metrics, introduced by Gailleton et al. (2021). These metrics were developed for analyses of basins where the concavity index, $\theta$, is uncertain. If $\theta$ in the $\chi$ transformation is uncertain, then changes to the $k_{sn}$ values will depend on drainage area (see Equation 5), and so this metric first measures the ratio of $k_{sn}$ values in the upstream part of the catchment to those in downstream parts of the catchment. We can then compare the ratio calculated with one value of $\theta$ against another, reported as a ratio, so this metric is a ratio of ratios (see Text S1 in Supporting Information S1 for full explanation and derivation). We also follow this approach for A and Q-based $\chi$ transformations because in the presence of a rainfall gradient upstream catchments will be affected differently from downstream catchments.

To give a concrete example, if we apply this metric to a simulated landscape where the “correct” $\chi$ transformation is used and the landscape is at steady state, the ratio of $k_{sn}$ between the headwaters and downstream will be 1. If, using the “wrong” $\chi$ transformation the headwaters are 20% steeper, the second ratio is 1.2 and the ratio of ratios (our distortion metric) is also 1.2. The reason we do this instead of simply comparing $k_{sn}$ values is because $k_{sn}$ will vary in space due to differences in erosion and lithology, so we are trying to separate the area-driven (or discharge driven) change in the differences in $k_{sn}$ driven by choices of $\chi$ transformation. A lengthy discussion of the rationale

<table><caption>**Table 2**
Summary of the Four Scenarios Considered for Each Basin</caption><thead><tr><th>A-driven incision</th><th>Q-driven incision</th></tr></thead><tbody><tr><td>$\chi_A, \chi_Q$</td><td>$\chi_A, \chi_Q$</td></tr></tbody></table>

Note. Under each incision scenario (A-driven, Q-driven incision) we calculate $\chi$ with ($\chi_Q$) and without rainfall ($\chi_A$).
and merits of this approach can be found in Gailleton et al. (2021). For this study, we quantify distortion depending on concavity index (Case i), incision scenario (Case ii) or both (Case iii).

# 3. Results

In this section we present the numerical experiments and explore the effects of rainfall gradients on topographic metrics such as channel steepness and concavity index. We cover both homogeneous and heterogeneous lithologies.

## 3.1. Homogeneous Lithology

### 3.1.1. Drainage Area or Discharge-Driven Incision?

Figure 3 shows the distributions of $\Delta D^*$ for the A-driven and Q-driven simulations. We have performed a KDE fitting assuming both a Q-driven incision rule (Equation 3) and a A-driven incision rule (Equation 2). $\chi_Q$ changes in each row of the figure for both the A and Q-driven simulations to incorporate the increase in rainfall gradients up to 10 m/yr across the simulation domain. This illustrates how $\Delta D^*$ reacts to assuming a different incision rule from the true model incision. Each kernel estimate is obtained for the collection of basins in the simulated domain. We obtain the median (solid line), 95th and 5th (dashed lines) percentiles for each distribution, repeating this procedure for each of the rainfall gradients. We consider $\Delta D^* = 0$ as the point where calculating $\chi_Q$ or $\chi_A$ would have no effect on the minimum disorder.

Basins where the $\chi$ and the model incision scenario match have lower $D^*$ values. This is represented by the A-driven distributions lying on the negative side of the x-axis, meaning that the disorder value calculated using $\chi_A$ is lower than the disorder calculated using $\chi_Q$. That is, unsurprisingly, using the version of $\chi$ that corresponds with the imposed incision rule resulting in less disordered channel profiles. The opposite is true of the Q-driven basins, which lie on the positive x-axis range, where the disorder is lower when calculating $\chi_Q$ than $\chi_A$. This is true of all three indicated metrics: median, 95th and 5th percentiles. None of the distribution tails overlap inside of the 5–95th percentile ranges. For every rainfall distribution illustrated, the drainage area and the discharge driven incision are statistically far enough apart to be considered distinct distributions.

As rainfall gradients increase, the medians of the distributions diverge. For smaller rainfall gradients, the distribution medians appear closer together but still outside the 95% of each other, which we consider an indication they are statistically distinct. The most significant differences arise from the discharge-driven scenarios, with the changes in the drainage area $\Delta D^*$ evolving slower with increased rainfall.

For every rainfall distribution illustrated under homogeneous lithology, calculating $\chi_Q$ or $\chi_A$ leads to statistically distinct $\Delta D^*$ distributions for the A and the Q-driven scenarios. Disorder can thus be used as a tool to recognize the dominating incision rule in numerical simulations: when $\Delta D^* < 0$, the incision is drainage area-driven, whereas $\Delta D^* > 0$ implies that discharge is the main incision mechanism.

### 3.1.2. $k_{sn}$ Distortion

We quantify $k_{sn}$ distortion based on the cases outlined in Section 2.4. We remind the reader of the three distortion cases we consider in this study (a) change in $\theta$, (b) change in incision rule, and (c) change in both. Similar patterns in changes in $k_{sn}$ can arise from either of the three cases.

Figure 4 illustrates each of the distortion cases along a series of rainfall gradients for the Q-driven simulations. Figures 4a and 4b correspond to changes in $\theta$ (Cases $i_Q$ and $i_A$). Figure 4c represents a change in incision rule (Case ii) and Figure 4d reflects a combined change in incision rule and $\theta$ values (Case iii). In all plots we indicate the no-distortion scenario at $y = 0$ with a solid black line. Each of the points in the plots correspond to a basin, colored by its draining direction (eastward or westward draining). In the following description, we use $\theta_{best}$ to refer to each basin's optimal $\theta$.

Figure 4a reflects data from $k_{sn-q}$ calculations, with the rainfall range used for the discharge-driven simulations on the x-axis. The $k_{sn-q}$ distortion remains close to 1, indicating that using $\theta = 0.45$ as opposed to $\theta_{best}$ has minimal effects in the basin-averaged $k_{sn}$ values. This can be explained by referring to the model set up. In the discharge-driven model with $m/n = 0.45$, with no other external factors to disturb equilibrium, we obtain steady state
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of $\Delta D^*$ kernel density estimates distributions (with medians and percentiles) fitted through all basins in A (red) and Q-driven (blue) incision models under homogeneous lithology with initial $m/n = 0.45$. The number of basins in each row corresponds to those shown in Figure 2. In each row, the Q-driven incision model is subject to a different rainfall gradient, starting from 1 m/yr up to 10 m/yr across the simulation domain. The larger the rainfall gradient, the more separated the distribution medians become, with the Q-driven distributions becoming flatter. Regardless of the rainfall gradient, the distributions medians (with their 5–95th percentiles) lie on opposite sides of the $\Delta D^* = 0$ threshold, indicating that the two models are always distinguishable with 95% confidence.</figcaption></figure>

channel profiles with $\theta_{best} = 0.45$. This value is obtained from disorder minimization including rainfall in the calculations. Figure 4b reflects the result of using a A-driven incision rule (not including a rainfall range) in the disorder minimization procedure when calculating $\theta_{best}$. Starting with the base case of constant rainfall of 1 m/yr across the simulation domain (with a null rainfall gradient), distortion gradually increases with rainfall ranges.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** *k*<sub>sn</sub> distortion ratios for the Q-driven incision scenarios under homogeneous lithology and initial *m*/*n* = 0.45. A *k*<sub>sn</sub> distortion ratio of 1 (solid black line) indicates no change in *k*<sub>sn</sub> under a given concavity and incision rule assumption. Each point represents a basin, colored by its draining direction in model coordinates: eastward or westward draining. Since the rainfall gradient increases linearly toward the west, east-draining basins experience less rainfall than west-draining basins. Panel (a) indicates that no distortion in *k*<sub>sn</sub> occurs when the concavity index and the incision case assumptions match those explicitly defined in the model. Panels (b–d) show the possible distortion scenarios that one might encounter under different assumptions from those given in the model. In general, the larger the rainfall gradient the greater the *k*<sub>sn</sub> distortion ratios become. West-draining (wetter) basins show smaller levels of distortion. Panel (b) highlights the effects of optimizing concavity index under an incorrect incision scenario, (c) assumes concavity index is kept at 0.45 but the incision scenario changes and (d) comprises the effects of *θ* optimization under different assumptions of incision scenarios, where we see the largest *k*<sub>sn</sub> distortions of up to 34%.</figcaption></figure>

The distortion values reach 23% where θ<sub>best</sub> > 0.45 (above 1) and 15% when θ<sub>best</sub> < 0.45 (below 1). As the rainfall range increases, θ<sub>best</sub> values diverge from 0.45. Figure 4c shows the distortion in case ii, which results from a difference in the assumed incision rule at a fixed concavity index (θ = 0.45). The distortion pattern follows a similar path to Figure 4b: as rainfall ranges increase, so does the k<sub>sn</sub> distortion. In both cases, distortion reaches 11%, with differences originating from the magnitude and the direction of the rainfall gradient. Figure 4d depicts distortion case iii: the effects of both a change in θ and a change in the incision rule. The percentage of k<sub>sn</sub> distortion is larger than in the other three scenarios: a 34% increase at the highest point against 23% and 11% in Figures 4b and 4c respectively. This arises from the basins having different θ<sub>best</sub> depending on the incision rule used to calculate the k<sub>sn</sub> distortion. The effect in the distortion is additive, meaning that compared to the correct case for both incision and θ (k<sub>sn-q</sub> (θ = θ<sub>best</sub>)), optimizing θ for the wrong incision case would lead to the greatest distortion out of the four cases considered. In Figures 4b–4d, we see a separation between the eastward and the westward draining slopes which originates from direction of the rainfall gradient. Westward facing slopes have higher rainfall amounts at the bottom of the basin, whereas eastward facing slopes experience high rainfall at the top of the basins. This creates a change in the distortion value with respect to the null distortion case at y = 0.

Distortion in the channel steepness index can be caused both by incorrectly assuming concavity indices and incorrectly assuming incision rules. That is, to calculate k<sub>sn</sub> one must set a value of θ, and one must choose whether or not to incorporate a proxy for discharge, and either of these assumptions may not be the best reflection of reality in any given landscape. Optimizing θ will decrease the distortion in k<sub>sn</sub> values, for cases when the incision mechanism is unknown (Figure 4c). If we identify an incision mechanism, then distortion can also be
<figure><img src="image_7.png" /><figcaption>**Figure 5.** River networks and basin outlines generated with the numerical model under a Q-driven incision scenario, with rainfall ranges across the simulation domain between 1 and 11 m/yr in the east-west direction. Each panel depicts the simulations evolved with a different lithology type: (a) Dense blob, (b) Sparse blob, and (c) Striped lithology. Panel (a) shows a highly irregular drainage divide, with variations in $k_{sn}$ within basins corresponding to the blob locations. Panel (b) shows a more regular drainage divide. $k_{sn}$ is more homogeneous within each basin, with small variability between the basins that contain the larger blobs. Panel (c) shows a clear fringe of high channel steepness in the center of the simulation, where the band of harder rock lies. The differences in $k_{sn}$ at either side of this band are due to rainfall differences. We also show the three types of lithologies used to evolve the landscapes above: dense blobs (d), sparse blobs (e) and striped (f).</figcaption></figure>

decreased by using either $k_{sn}$ (Figure 4b) or $k_{sn-q}$ (Figure 4a), depending on the case. However, we find that optimizing $\theta$ for an incorrect incision case leads to the highest distortion values (Figure 4d).

## 3.2. Heterogeneous Lithology

We have shown that it is possible to differentiate a signal from an A-driven incision scenario from a Q-driven scenario, under homogeneous lithology. However, many natural landscapes contain a range of lithologies with different erodibilities. We complement our analysis with data from simulations with heterogeneous lithologies, as described in Section 2.1.2. Figures 5a–5c shows how the different lithologies affect basin formation in the
<figure><img src="image_8.png" /></figure>

**Figure 6.** Comparison of $\Delta D^*$ kernel density estimates distributions medians and 5–95th percentiles fitted through all basins in A (red) and Q-driven (blue) incision models. We show the results under three simulated heterogeneous lithologies under an initial $m/n = 0.45$. In each row, the Q-driven incision model is subject to a different rainfall gradient, starting from 1 m/yr up to 10 m/yr across the simulation domain. More heterogeneous lithologies, such as the dense blob and the striped lithologies lead to less distinguishability between A and Q-driven incision cases: their distribution medians lie each other's 5–95th percentile bounds. The sparse blob simulations (the least heterogeneous) show a similar behavior to the homogeneous lithology case where A and Q-driven incision case basins are located at either side of the $\Delta D^* = 0$ divide, albeit with some overlap in the 5–95th percentiles.

simulation domain. We also include in Figures 5d–5f the simulated underlying lithologies. We highlight the sudden differences in $k_{sn}$ in (a) where the higher erodibility blobs are located. (b) Portrays a more gradual change in lithologies, leading so less changes in channel steepness within each basin. The striped lithology case (c) shows step changes in $k_{sn}$, aligning with the two-fold change in rock erodibility between the central and the side bands. We note as well the differences in the divide location and sinuosity between the cases, with the dense blob lithology (a) leading to the largest asymmetry between the west and the east draining basins.

In Figure 6 we explore whether $\Delta D^*$ shows significant changes in A-driven distributions compared to Q-driven cases for the three modeled heterogeneous lithology scenarios.

In contrast with the homogeneous lithology case, the dense blob A and Q-driven basins are not statistically distinct. Both distribution medians fall within each other's 95% percentiles, regardless of the rainfall range. This is also the case for the striped lithology basins. The more heterogeneous the erodibility pattern is, the harder it becomes to distinguish between the original incision rules. Hence, it becomes more challenging to differentiate A-driven from Q-driven incision rules using the $\Delta D^*$ method. The sparse blob lithology model has a clearer differentiation in $\Delta D^*$ for each incision mechanism, with the $\Delta D^* = 0$ line separating the A and the Q-driven basins. However, the confidence is lower than for the homogeneous lithology case, since the 5–95th percentiles intercept, yielding a lower confidence threshold.

We also quantify how heterogeneous lithologies impact the distortion of $k_{sn}$ values. We start by looking at the result of adding lithology to the discharge-driven model scenario and compare the results to Figure 4. We highlight that in the heterogeneous lithology cases, we do not expect $\theta$ to be 0.45 at the base rainfall case (1 m/yr constant across the domain). We have introduced a perturbation, leading the system to diverge from $\theta = 0.45$, which we obtained in an unperturbed, steady state scenario. This is reflected in a distortion between 10% and 20% at the base rainfall scenario in all distortion cases. Irrespective of the range in rainfall, this distortion remains within the same bounds. We see consistently higher distortion values for all rainfall ranges in Figure 7b, as a result of using the incorrect incision scenario when optimizing $\theta$. In Figure 7c we see the distortion associated with case ii (caused by changes in incision rule), when keeping $\theta = 0.45$. As the rainfall range increases, the distortion shifts
<figure><img src="image_9.png" /><figcaption>**Figure 7.** $k_{sn}$ distortion ratios for the *Q*-driven incision scenarios under a heterogeneous lithology consisting of densely packed blobs of different erodibilities and initial $m/n = 0.45$. A $k_{sn}$ distortion ratio of 1 (solid black line) indicates no change in $k_{sn}$ under a given concavity and incision rule assumption. Each point indicates a basin, colored by its draining direction in model coordinates: eastward or westward draining. Since the rainfall gradient increases linearly toward the east, west-draining basins experience less rainfall than west-draining basins. Panel (a) shows distortion of up to 20%, in both eastward and westward draining basins, assuming the correct incision scenario (*Q*) and optimizing $\theta$. Panels (b–d) show the possible distortion scenarios that one might encounter under different assumptions from those given in the model. Panel (b) highlights the effects of optimizing concavity index under an incorrect incision scenario, (c) assumes concavity index is kept at 0.45 but the incision scenario changes and (d) comprises the effects of $\theta$ optimization under different assumptions of incision scenarios, where we see the largest $k_{sn}$ distortions of up to 54%. The rainfall effects are less visible than in the homogeneous lithology case, with no clear trend associated with increases in rainfall in plots (a)–(c). Panel (d) shows a clear divide between east and west-draining basins, and an increase in distortion as the rainfall gradient increases. In this lithological scenario, westward (wetter) basins display larger $k_{sn}$ distortions.</figcaption></figure>

toward being larger, up to 44% in east-draining basins. When we optimize $\theta$ and compare between the two incision scenarios in Figure 7d, the distortion in both the west and the east side is increasing at a similar rate with increasing rainfall ranges. This is the only distortion case where we see a systematic increase in distortion with rainfall ranges. Only in this case do we see a clear differentiation between east and west draining basins. In panels (a)–(c) we observe that basins are not strictly separated by draining location as a result of the direction of the rainfall gradient. Instead, there is a mixture of basins crossing the null distortion line implying that changes in erodibility are masking the clear basin differentiation we see in Figure 4.

Heterogeneities in lithology add to the distortions in channel steepness index caused by rainfall ranges. The signal from the rainfall is reflected as a systematic increase in distortion in *A*-driven cases. In *Q*-driven scenarios the systematic increase of distortion with rainfall range is masked by the lithological heterogeneities when comparing between $\theta$ values. Optimizing $\theta$ will remove the distortion that we incorporate by using $\theta = 0.45$, which is more prominent in heterogeneous than homogeneous lithologies. Making a decision on the type of incision, will also reduce the distortion. However, this can end up with an increased distortion if the incision type assumed for calculating $\chi$ is incorrect and $\theta$ is subsequently optimized using channel information obtained from that incorrect incision type.
<figure><img src="image_10.png" /></figure>

**Figure 8.** Comparison of median $\Delta D^*$ values for all basins in each of the mountain ranges, calculated using real rainfall gradients. We highlight the median and the 5-95th percentile bounds with a horizontal solid line. Areas where the $\Delta D^*$ is greater than 0 identify as being predominately forced by rainfall. Conversely, if $\Delta D^* < 0$ other forcings aside from rainfall (such as tectonics or lithology) are dominating the signal. Based on the $\Delta D^*$ values, some cases show a slight preference toward drainage $\chi_A$ (N. Andes, S. Andes and Kaçkar) whereas the rest favor a discharge-driven description, using $\chi_Q$. However, the 5-95th percentiles always cross the $\Delta D^* = 0$ point, which separates the A and the Q-driven scenarios, meaning that the results are largely basin dependent.

### 3.3. Natural Landscape Study Cases

We have quantified how changes in rainfall rates across a landscape can lead to differences in $k_{sn}$ and $\theta_{best}$ under various incision assumptions. Our numerical experiments suggest that it is possible to statistically distinguish A-driven and Q-driven incision landscapes in a lithologically homogeneous model using disorder metrics. Experiments with heterogeneous lithologies capture how the rainfall signal is obscured by additional forcings, namely spatial variation in erodibility, a common feature in natural landscapes. Having isolated these signals in experiments, we now apply this knowledge to real landscapes, where we only have a *posteriori* topographic states from which to infer the landscape incision rule. If the incision of natural landscapes is driven by *Q*, we expect $D^*(\theta_{best,Q}) < D^*(\theta_{best,A})$ and $\Delta D^* > 0$. For each basin, we repeat the calculations outlined in Section 2.2, and plot the median values of normalized distortion ($\Delta D^*$) and their distributions in Figure 8. This metric quantifies which of the two assumed incision scenarios results in lower distortion of the channel network. If all basins are positive, or all are negative, we view this as an indication that one of the two assumed incision signals better describes the observed topographic data.

We see that in Figure 8 most natural landscapes have basins that straddle $\Delta D^* = 0$, meaning that we cannot determine the most likely incision scenario with any confidence. Only the areas in the Rockies and the Northern Andes show statistically significant imprints on $\Delta D^*$, favoring respectively discharge and drainage area as the incision mechanism yielding the least disorder. The remaining sites show only a slight trend toward an incision scenario (discharge - Pyrénées, Qinling, Massif Central, Alburz- and drainage area—Southern Andes, Kaçkar-), with the distributions centered largely around $\Delta D^* = 0$ and spreading at either side of the incision case division.

If we are not able to distinguish the correct incision scenario using the disorder metrics, how uncertain will our interpretations of channel steepness become? In Figure 9 we investigate distortions in the channel steepness index, colored by the predominant lithology in each basin, following the classification from the GLIM data set (Hartmann & Moosdorf, 2012). Whilst some areas share the same lithology in all basins (Kaçkar, Alburz—volcanic rocks), other regions show larger heterogeneity, with basins of sedimentary, metamorphic and volcanic composition (Qinling).

We highlight the difference in scale of the results shown here. In the modeling experiments, some scenarios show a distortion of no more than 10% under low rainfall (e.g., Q-driven homogeneous lithology (Figure 4a), A-driven homogeneous lithology (Figure S11b in Supporting Information S1), and A-driven dense blob lithology (Figures S2a–S2c in Supporting Information S1)). In real landscapes, we observe larger distortion values, in some instances reaching 50% in Figures 9a and 9d, and above 75% distortion in Figures 9b and 9c. In natural basins, using $\theta = 0.45$ becomes a large source of distortion to the channel steepness index regardless of the assumed incision scenario. This is a consequence of the $\theta_{best}$ values diverging from 0.45, which corresponds to the $m/n$ value of the models. Compared to these large distortion values incurred by using $\theta = 0.45$, the distortion values arising from choosing A-driven or Q-driven incision scenarios fall considerably. We see most distortion falls within the 25% bounds, regardless of the mountain range, although with different amounts of spread around the null distortion line at $y = 1$. This shows that in natural landscapes, optimizing the concavity index has a larger impact on $k_{sn}$ values than the choice of incision scenario, at least where there are no constraints on erosion rates (e.g., Leonard et al., 2023a). Regarding the impact of lithology on $k_{sn}$, we observe that in some basins, such as the only siliciclastic sedimentary basin in the Massif Central, the distortion is larger than its neighboring metamorphic basins. In other cases, such as the Rockies, the metamorphic basin has distortion comparable to the volcanic basins in cases *ii* and *iii*, suggesting that the choice of incision rule introduces little distortion. However, the choice of $\theta$ places the metamorphic basin on the opposite side of the null distortion line, implying different a concavity trend from its neighbors. Other mountain ranges, such as Alburz or Qinling, both show large amounts of distortion, despite Alburz basins composition being homogeneous and Qinling having large heterogeneities.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** *k*<sub>sn</sub> distortion ratios for the selected natural landscapes, where each basin is colored by their dominating lithology. A *k*<sub>sn</sub> distortion ratio of 1 (solid black line) indicates no change in *k*<sub>sn</sub> under a given concavity and incision rule assumption. The four cases (a)–(d) illustrate how one might introduce distortions in *k*<sub>sn</sub> when we do not know what the dominating incision scenario is, as it often happens. We show that the *k*<sub>sn</sub> distortion introduced by choosing θ = 0.45 instead of θ<sub>best</sub>, reaches 50% under Q-driven incision (a) and 81% under A-driven incision (b). Fixing θ to 0.45 under different incision scenarios can lead to distortions of up to 92% (c). Panel (d) depicts the distortion introduced by optimizing θ under both incision regimes, with distortions of up to 51%, the lowest out of all cases.</figcaption></figure>

3.3.1. Argentinian Andes Study Site

We choose the study area in the Northern Argentinian Andes to illustrate how the topographic metrics react to A
and Q incision rule assumptions. In depth topographic metrics results for the remaining areas are in the Supporting
Information S1.

The study area is situated in the Salta region of Northern Argentina, near the border with Bolivia. The lithology of the region is shown in Figure 10a, where metamorphic and sedimentary rocks dominate the landscape. The basins cross a range of rock hardnesses, with sedimentary units dominating in the southern basins and northern basins dominated by metamorphics. Over a distance of 30 km we see rainfall gradients of 0.05 m/yr (Figure 10b). This area has a very strong topographic gradient, with relief ranging from 1,500 m to over 5,000 m. Figure 10c shows how the rainfall swath profile is linearly decreasing with elevation.

To better understand the behavior of ΔD* and k<sub>sn</sub> distortion ratios, which are at the core of our analysis, we have plotted maps of localized k<sub>sn</sub> differences (Figure 11) and χ − z profiles (Figure 12). Figure 11 includes the same distortion cases as in Figure 9, but in this case we skip the step where we average over upper and the lower reaches of each basin to obtain a single distortion ratio per basin. Instead, we use k<sub>sn</sub> local differences to visualize which points in the basins are most affected by concavity optimization and changes in incision rule. We see from Figures 11a and 11b that optimizing concavity index affects k<sub>sn</sub> differences both at a basin wide scale and localized points,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Summary of the topographic and climatic characteristics for the study area in the Southern Andes, Argentina. The base maps are in coordinates of UTM20S. Panel (a) shows the lithological composition from the GLIM database (Hartmann & Moosdorf, 2012), following their nomenclature. Most basins present metamorphic, sedimentary rocks or a mix of both. Panel (b) includes the mean annual precipitation from NASA GPM (Skofronick-Jackson et al., 2017), averaged between 2001 and 2021. We highlight the basins selected for this study and the channel networks. Panel (c) shows 30 km swath profiles for elevation and rainfall from E-W over the black dashed rectangle in (b). The elevation data is from the Copernicus 30 m DEM (European Space Agency, 2020). Due to differences in resolution between elevation (30 m) and rainfall (10 km), the bin spacing has been accordingly scaled. The elevation increases up to around 5,000 m, linearly over the first 30 km, whilst the rainfall decreases linearly throughout.</figcaption></figure>

in some cases coinciding with changes in lithology. Most basins show values between 0.1 and 0.2 throughout, with jumps up to 0.5 (note a value of 0 would mean that $k_{sn}^*$ is the same for each term of the difference). Case *ii* shows differences of less than 0.1 in most basins except from some very localized points shared also by panels (A)-(B). This is the scenario that originates the least changes in $k_{sn}^*$; changing the incision rule tends to cause very localized changes in $k_{sn}^*$ but does not originate large changes in the overall $k_{sn}^*$ values. Case *iii* shows some basins remaining practically unchanged, as in case *ii* but others showing differences similar to those in cases $i_A$ and $i_Q$. We see that some spots are more susceptible to changes in $k_{sn}^*$, regardless of the cause (either concavity index changes or incision rule). However, changes that happen throughout the basin seem to be linked to changes in concavity index rather than different incision rule assumptions.

We can also investigate the channel-wise behavior through the $\chi-z$ profiles and $k_{sn}^*$ values in Figure 12. Figures 12a and 12b show the profiles for $\theta_{best}$. These have generally higher normalized steepness values than the profiles calculated using $\theta = 0.45$, both for the *A* and the *Q* assumed incision cases. The highest $k_{sn}^*$ points correspond to breaks in the profiles, or knickpoints. Some of those are placed in the smaller tributaries to the main stem. Both top panels, calculated with $\theta_{best}$, show more collapsed profiles than the bottom panels, which use $\theta = 0.45$. Figures 12a and 12c, which assume a *Q*-driven incision rule, show more channel collinearity than Figures 12b and 12c, which assume a *A*-driven incision mechanism. For this particular basin, $\chi_Q$ calculated with $\theta_{best}$ has the most collapsed profile (Figure 12a).

# 4. Discussion

## 4.1. Channel Steepness Index, Concavity Index, and Erosion Rates

Our numerical experiments show that in lithologically homogeneous landscapes it is possible to distinguish between *A* and *Q*-driven landscapes from topographic metrics alone, using the proposed disorder method. Distortions in normalized channel steepness index can go above 30% in our simulations, with variations depending on the lithology and the intensity of the rainfall gradient. Homogeneous lithology leads to $k_{sn}$ distortion that increases monotonically with increasing rainfall ranges (Figures 4b–4d), whereas heterogeneous lithologies muddle rainfall-related signals (see Figures 7a–7c). This means that lithology needs to be accounted for before one can test for a rainfall signal.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Maps showing the four cases of *k*<sub>sn</sub> distortion for selected river basins in a region of the Southern Andes, Argentina. The river network is selected based on the criteria for calculating the normalized *k*<sub>sn</sub> distortion outlined in S1. We show the extracted basins over the lithological map of the region, where sedimentary (brown) and metamorphic (green) rocks dominate. Panels (a) and (b) explore the effect of optimizing the basin concavity index. Panel (c) and (d) highlight changes due to concavity index and assumed incision rule. The differences in *k*<sub>sn</sub> depend strongly on the basin, with the largest changes occurring at localized points, sometimes corresponding to changes in lithology. Panel (b) has generally higher *k*<sub>sn</sub> difference values than *k*<sub>sn-q</sub> (a) and (c and d). Panel (c) shows the overall least difference out of all four cases, with the high valued points at similar locations as in (a) and (b).</figcaption></figure>

Natural landscapes show more climatic and lithological complexity than our numerical models, leading to further difficulties in identifying the least disordered incision scenario (Figure 8). By choosing a small number of basins across a variety of rainfall gradients and lithological compositions, we were hoping to identify a disorder pattern that would clearly favor either a *A* or *Q*-based incision rule following the classification system derived from our simulations. Even though this method proved non-conclusive, were are now able to provide an estimate of how much distortion in normalized channel steepness is introduced assuming different incision rules or concavity indices (Figure 9). For $\theta = \theta_{best}$, the variability between $k_{sn}$ and $k_{sn-q}$ is smaller than for $\theta = 0.45$. This means that as long as the concavity index has been constrained, the channel steepness index variations between *A* and *Q*-scenarios will be smaller than if $\theta = 0.45$ was chosen. We note that the mountain ranges where the distributions diverge for *A* and *Q* cases often correspond to the ones where the disorder is high, demonstrating once more that using unconstrained values of $\theta$ to extract geomorphometrics can amplify process-linked biases.

Aside from the distortion caused by rainfall gradients and different rock types, normalized channel steepness can also be affected by the value of the concavity index. In simulations depicted in Figure 4 we have imposed a *Q*-driven incision rule, with $\theta = 0.45$. We find that the distortions associated with choosing $\theta = 0.45$ (as opposed to $\theta_{best}$) are similar in magnitude to distortions resulting from the choice of incision rule (Figures 4b and 4c). In a real landscape we cannot define the “true” incision rule or the “true” value of $\theta$ because of the complex interplay of climate, lithology and tectonics. To do this, we would need to resort to field measurements of erosion or uplift rates, topography alone would not be enough. However, we can determine which incision scenario and which
<figure><img src="image_14.png" /><figcaption>**Figure 12.** $\chi$ profiles for a basin in the Argentinian Andes, under different concavity index and incision scenarios which use the mean annual precipitation (m/yr) from the area. The profiles are colored by $k_{sn}^*$ values. Disorder decreases when we optimize the concavity index ($\theta_{best}$) (a) and (b)), increasing the collinearity of the channels. $k_{sn}^*$ is higher for localized areas in (a) and (b). In panels (c) and (d) the $\chi$ profiles are more fanned, which indicates higher disorder. Compared to the effects that changes in concavity index have in the profiles, the choice of incision rule introduces secondary differences in the channels.</figcaption></figure>

concavity index provides the most accurate representation of the geomorphic processes affecting measured erosion rates, through the incision rules and $\theta$ that we test in this contribution.

Determining the incision rule that is more consistent with observed topography will help in interpreting fluvial landscapes in areas with large rainfall gradients. The distortion that we have found in our study can then be translated into uncertainties in erosion rates. Adams et al. (2020) found that Be-10 erosion rates in the Himalayas and Bhutan are better constrained when using $k_{sn-q}$ instead of $k_{sn}$. Recent work by Leonard et al. (2023a, 2023b) corroborates this finding for the north-central Andes, a region with high rainfall gradients, where only $k_{sn-q}$ is able to correctly match erosion rates of different lithologies.

Our study highlights the need to obtain field erosion rates to verify normalized channel steepness values. We have quantified the distortion that different incision rule and concavity index assumptions can introduce in normalized channel steepness index. We find that changes in the concavity index can lead to changes in normalized channel steepness index which are often larger than those caused by spatially varying rainfall, suggesting a potential cumulative effect of the different biases. With complex lithology and rainfall patterns, it is challenging to use our proposed disorder method to assess the suitability of an incision rule, given no additional erosion rate information.

Studies of topography have linked $\theta$ values to different rainfall (Harries et al., 2021; Zaprowski et al., 2005) and lithological regimes (Duvall et al., 2004). Other theoretical studies have shown that even under the assumption that the stream power family of incision rules are correct (e.g., Lague, 2014), several other processes such as the nature of the river bed (Howard & Kerby, 1983; Tucker & Whipple, 2002), the sediment availability (Wickert & Schildgen, 2019), and the tectonic setting (Kirby & Whipple, 2001) can also affect the value of $\theta$ in different reaches of the rivers in a way that would become hard to disentangle from spatial changes in rainfall. Han et al. (2015) explores changes that spatially varied rainfall can have on the concavity index through modeling
experiments. In natural landscapes, changes around their chosen value of 0.5 vary comparatively less due to rainfall gradients than changes caused by factors such as the type of channel bed (Whipple & Tucker, 2002), sediment amounts (Gasparini et al., 2007) or uplift rate changes (Kirby & Whipple, 2001).

In our model simulations, we have assumed a spatially and temporally constant uplift. Since variations in uplift can affect $k_{sn}$ values, we decided to keep uplift constant and focus only on heterogeneities in rainfall and lithology. In natural landscapes, uplift can vary and some studies have shown that varying uplift clouds both climatic and lithologic signals (Godard et al., 2014; Seybold et al., 2021). Work by Leonard et al. (2023a) and (Leonard et al., 2023b) decouples the climate effects in tectonically active regions, such as the Andes, providing evidence that rainfall gradients can indeed be detected despite the presence of strong uplift fields. We formulate our study assuming that the disorder metric is robust against tectonics (Mudd et al., 2018). In this contribution we have quantified the distortion introduced by the choice of incision rule and concavity index. Any changes in $k_{sn}$, $k_{sn-q}$ or $\theta$ linked to the presence of varying uplift is intrinsically included in the reported distortion values.

## 4.2. Lithology

Our simulations show that heterogeneous lithology patterns have a strong effect in $k_{sn}$ and disorder metrics, clouding the signal from a rainfall gradient which is otherwise visible under homogeneous lithology using the proposed disorder method. Many studies focus on the effects that lithology has on the evolution of landscapes, and how we can quantify it through topographic metrics such as normalized channel steepness index. This is particularly important in ancient orogens where changes in topography will likely be due to lithological or climatic changes (Clementucci et al., 2022). We discuss some of the most relevant work for our study below.

Originally suggested by Mills (2003), Gallen (2018) and Gallen and Thigpen (2018) quantify how differences in rock type map to changes in $k_{sn}$. Similarly, Leonard et al. (2023a) assess whether using $k_{sn}$ or $k_{sn-q}$ leads to better groupings of lithological units under strong rainfall gradients. In our simulated landscapes with heterogeneous lithologies, we set erodibilities between $5 - 10 \times 10^{-8}$, a two fold maximum increase which is in line with erodibility changes in areas such as the Pyrénées (Bernard et al., 2019). Given that we are not able to clearly discern the rainfall signal in these simulations, in particular under the dense blob and the striped lithology scenarios (Figure 6) we must pose the question of under what specific condition are we able to detect these climatic signals using the proposed disorder method.

The disorder method assesses whether rainfall gradients leave traces in the river profiles which cause measurable differences in $\chi$ profiles using a constrained concavity index. Under homogeneous lithology, it is possible to use disorder to distinguish between the *Q* and the A-driven incision cases. However, when heterogeneous lithology is introduced, this climatic signal is no longer clearly visible. There are two possible reasons why this is the case. First, the changes in erodibility might be too large for any climatic signal to show. This has been recorded in the studies by Gallen (2018), in the Appalachians, where 5 fold erodibility changes cloud any climatic or uplift signal. Clementucci et al. (2022) show that in relict landscapes, the prevailing lithology has a stronger control over erodibility than climate. Second, the disorder method, and consequently $\chi$ profiles, where lithology is not explicitly controlled for, are not able to resolve the rainfall gradients.

In our simulations, we have accounted for the rainfall gradients by explicitly incorporating *Q* in the calculations of $\chi$ (Equation 6) and slope (Equation 7). Lithology, on the other hand, is implicitly incorporated in the parameters describing fluvial morphology such as $k_{sn}$ or $\theta$. This means that unless lithology can be accounted for by selecting areas where the bedrock is largely uniform (such as done in the study by Marder and Gallen (2023)) or, alternatively, rock type is explicitly incorporated in the $\chi$ formulation, the visibility of the climate signal is largely dependent on the intensity of lithological heterogeneities. Alternative methods, such as an inverse approach, that various studies have used are covered in the limitations sections.

## 4.3. Rainfall and Climate

In this work, we use a linear function of rainfall with topography to simulate the orographic effect. This captures two precipitation patterns: increasing toward and away from the drainage divide. Similar linear patterns are used across the literature as a simplification of the complex orographic rainfall patterns seen in nature (Anders et al., 2006; Bookhagen & Strecker, 2008; Leonard & Whipple, 2021). In our simulations, we have seen how the overall relief decreases as precipitation rates increase. This has been quantified in natural landscapes by Marder
and Gallen (2023). However, the relief on each side of the divide varies depending on whether rainfall increases or decreases toward the outlet. The divide migrates slightly toward the drier side which presents steeper channels and is bounded by higher relief. In natural landscapes, the relationship between relief and rainfall at a local scale is complex and heavily influenced by local processes (Champagnac et al., 2012; Montgomery et al., 2001).

We observe different levels of channel steepness index distortion for the wet and the dry side of the simulated mountain ranges. This allows us to use the distortion ratios to identify the direction of the rainfall gradient, either increasing or decreasing with elevation. This difference between the wet and the dry side of natural landscapes has been studied in depth in natural settings such as Hawai'i (Ferrier et al., 2013).

In the modeling framework, we have assumed that the rainfall pattern remains constant throughout time, but this is not necessarily the case in natural landscapes. Even small climate changes can lead to changes in the rainfall pattern and discharge amounts at different parts of the catchment, altering the local erosion rates and displacing the landscape from an equilibrium state (Leonard & Whipple, 2021). Changes in atmospheric circulation patterns or temperature changes (Bradley, 2015; Herman et al., 2013; Ward & Galewsky, 2014) are some examples of cases when other climate variable can affect erosion rates. Glaciations are also important when considering a large portion of the Earth's landscapes, and they also show a relation with the erosion rates and the relief of the landscape.

## 4.4. Limitations

One of the limitations of our study lies in the precipitation treatment, both in the modeling study and in the natural landscapes. While it is possible to reproduce a realistic rainfall pattern using the Fastscape module adapted from Smith and Barstad (2004), it is computationally expensive and requires knowledge of the wind patterns of the region of interest, which change at different atmospheric layers. We assume that the orographic rainfall can be approximated in the modeling framework by a linear rainfall trend that does not change through the simulation, regardless of the relief (Roe et al., 2003). In real orographic rainfall scenarios, there would be a positive feedback between the rainfall and the topography, which our model does not capture.

In our natural sites, we also assume that the rainfall pattern, derived from the average rainfall from 20 years (2001–2021), is representative of the precipitation at that site throughout its history, over geologic timescales. Climate patterns have changed throughout the centuries driven by changes in atmospheric conditions, solar irradiance, and biosphere and ocean changes (Bradley, 2015). The regions of study lie within the mid-latitudes, where the climate has remained largely unchanged over long timescales (Roe et al., 2002). However, due to using rainfall data from the 21st century, recent changes in rainfall patterns due to human induced climate change are present in data from the past years compared to data prior to industrial revolution (Trenberth, 2011).

The question of whether mean annual precipitation should be used to describe the climate of a region is also a highly debated topic, with studies suggesting that it is the storms and extreme events which contribute the most to mean annual precipitation and do the most erosive work (Deal et al., 2017, 2018; DiBiase & Whipple, 2011; Rasmussen et al., 2016; Sorensen & Yanites, 2019). Other studies prefer the use of mean annual precipitation (Adams et al., 2020; Anders et al., 2006; Armitage et al., 2011; D'Arcy & Whittaker, 2014; Gasparini & Whipple, 2014; Leonard & Whipple, 2021; Rossi et al., 2016), especially when capturing the erosion work over a longer climatic trend or incorporating it in long-term landscape evolution models.

In our study, we assume that aggregate precipitation is related linearly with discharge. We do not consider additional factors that might affect this relation such as infiltration or evapotranspiration, which would alter the runoff ratio. Several studies have analyzed this relationship, finding that runoff ratios are smallest in arid regions (Chang et al., 2014; Rossi et al., 2016), with the largest variability in discharge values. DiBiase and Whipple (2011) and Lague et al. (2005) explore how runoff variability affect channel steepness introducing a complex model accounting for stochastic high-discharge events. These approaches explore in-depth how runoff can affect channel morphology and play a role in estimating correct erosion rates. Since the simple theorized relationships we use are not able to clearly differentiate between Q and A-driven incision scenarios, we consider that introducing these complex hydrologic relations would cloud the parameter space and introduce unnecessary uncertainties in our analysis.

In the modeling framework, many processes have been simplified. We have already mentioned the rainfall and discharge patterns. The representation of lithological units, the exclusion of sediment supply and the homogeneous uplift, have all been choices made to isolate the climatic signal as much as possible. Future work would
include modeling to control varying uplift and isolate its signal from climatic and lithological forcings. We are also assuming that the detachment limited stream power model forms a good basis for how rainfall interacts with uplift and erosion, which many studies support (e.g., Adams et al., 2020; Gasparini & Whipple, 2014; Harries et al., 2021; Leonard & Whipple, 2021; Leonard et al., 2023b), while acknowledging it still does not fully explain all geomorphic processes at the landscape scale.

We explore the use of normalized channel steepness and disorder metrics to study the effects of rainfall gradients on landscapes. We have framed this study as a minimization problem involving the disorder metric and the concavity index. Previous studies propose alternative frameworks for minimization problems, such as the inversion methods initially developed by Pritchard et al. (2009). These methods reconstruct the erosion history and relief of a catchment from an initial assumption of uplift and fluvial incision based on the current day topographic shape of the landscape. The work by Roberts and White (2010), Fox et al. (2014), and Rudge et al. (2015) describe in detail the inverse method and its applicability. Inversion approaches have successfully reconstructed the erosional and uplift history of lithologically heterogeneous landscapes, assuming that lithological forcings play a role in rock uplift that is otherwise not possible to disentangle using channel steepness (Fisher et al., 2022; Gallen, 2018, following the approach by Goren et al. (2014)). Similar to our study, these methods attempt to constrain a certain element of the landscape (e.g., lithology or tectonics), through imposing of a series of spatio-temporal variations in the forcing in question. In a similar fashion, it would be an interesting avenue to isolate the effects of rainfall on a natural landscape using inverse methodologies, to reconstruct the climate regime that would lead to the topography we now observe.

## 5. Conclusions

In this study, we explore if topographic metrics can be used to determine whether channel incision is most closely related to drainage area or discharge (or some proxy thereof), given a rainfall pattern. Many past papers quantify channel steepness calculated based on drainage area as an indicator of river incision rates (e.g., Harel et al., 2016; Kirby & Whipple, 2012), but now that precipitation records are more readily available (Skofronick-Jackson et al., 2017), we can question whether adding empirical rainfall gradients to the equation will yield different topographic outcomes in river channels.

In a simple numerical model with homogeneous lithology, normalized disorder metrics ($D^*$ and $\Delta D^*$) (Gailleton et al., 2021; Goren et al., 2014; Mudd et al., 2018) can be used to distinguish between A- and Q-driven incision, as they scale monotonically with the imposed rainfall gradient. However, when adding heterogeneous lithology, the incision signal resulting from rainfall gradients is obscured. The signature of rainfall gradients is also masked in the disorder metrics by the choice of concavity index. We have quantified the effects that optimizing the concavity index $\theta$ can have in the channels, concluding that using the standard value of $\theta = 0.45$ amplifies the distortion caused by rainfall effects.

In natural landscapes, we cannot establish a general topographic rule to distinguish whether A or Q-based metrics are better suited to describe landscapes subject to rainfall gradients, assuming that no erosion or uplift information is known about the region. Based on the differences in normalized disorder values of catchments, some favor A-driven incision over Q-driven incision and vice versa. Heterogeneities in bedrock erodibility cloud the rainfall signal, and these can only be constrained by field measurements of erosion rates. Given that we are not able to separate those cases topographically, we quantify how much distortion we would introduce in channel steepness index if we failed to identify the incision mechanism that would lead to the correct erosion rates. Our results suggest that in most basins we would see maximum changes in normalized channel steepness index of up to 25%, which does not constitute enough to drastically change the interpretation of erosion rates across the landscape. We compare this to distortions in $k_{sn}$ of 50% obtained from using $\theta = 0.45$ instead of optimizing the concavity index. We suggest readers to use $\theta_{best}$ as an efficient method to reduce distortions already introduced by an unknown incision mechanism.

## Data Availability Statement

Analyses have been run using open source software (lsdtopotools v0.9, Mudd et al., 2023; lsdtopytools). Precipitation data was retrieved using the package gpmprecipitation_tools. Visualization scripts, data and model workflows are available online (Ruiz Sánchez-Oro, 2024).
Acknowledgments

Funding for MRSO was provided by the Natural Environment Research Council (NERC NE/S007407/1). This research has been supported by the H2020 European Research Council (Grant 803721). We thank Sean F. Gallen, an anonymous reviewer, the Associate Editor and the Editor, Amy East, for their thoughtful and thorough comments that have greatly improved the paper. Perceptually uniform scientific color maps (Crameri, 2023) are used in this study to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies (Crameri et al., 2020).

References

Aceñolaza, F. G., & Miller, H. (1982). Early Paleozoic orogeny in southern South America. *Precambrian Research*, **17**(2), 133-146. https://doi.org/10.1016/0301-9268(82)90052-3

Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, A. M., & Hodges, K. V. (2020). Climate controls on erosion in tectonically active landscapes. *Science Advances*, **6**(42), 3166-3182. https://doi.org/10.1126/sciadv.aaz3166

Allen, M. B., Ghassemi, M. R., Shahrabi, M., & Qorashi, M. (2003). Accommodation of late Cenozoic oblique shortening in the Alborz range, northern Iran. *Journal of Structural Geology*, **25**(5), 659-672. https://doi.org/10.1016/S0191-8141(02)00064-0

Anders, A. M., Roe, G. H., Hallet, B., Montgomery, D. R., Finnegan, N. J., & Putkonen, J. (2006). Spatial patterns of precipitation and topography in the Himalaya. *Special Paper of the Geological Society of America*, **398**, 39-53. https://doi.org/10.1130/2006.2398(03)

Armitage, J. J., Duller, R. A., Whittaker, A. C., & Allen, P. A. (2011). Transformation of tectonic and climatic signals from source to sedimentary archive. *Nature Geoscience*, **4**(4), 231-235. https://doi.org/10.1038/ngeo1087

Babault, J., Viaplana-Muzas, M., Legrand, X., Van Den Driessche, J., González-Quijano, M., & Mudd, S. M. (2018). Source-to-sink constraints on tectonic and sedimentary evolution of the western Central Range and Cenderawasih Bay (Indonesia). *Journal of Asian Earth Sciences*, **156**, 265-287. https://doi.org/10.1016/j.jseas.2018.02.004

Baynes, E. R., Lague, D., Steer, P., Bonnet, S., & Illien, L. (2020). Sediment flux-driven channel geometry adjustment of bedrock and mixed gravel-bedrock rivers. *Earth Surface Processes and Landforms*, **45**(14), 3714-3731. https://doi.org/10.1002/esp.4996

Bernard, T., Sinclair, H. D., Gailleton, B., & Fox, M. (2021). Formation of longitudinal river valleys and the fixing of drainage divides in response to exhumation of crystalline basement. *Geophysical Research Letters*, **48**(8), e2020GL092210. https://doi.org/10.1029/2020GL092210

Bernard, T., Sinclair, H. D., Gailleton, B., Mudd, S. M., & Ford, M. (2019). Lithological control on the post-orogenic topography and erosion history of the Pyrenees. *Earth and Planetary Science Letters*, **518**, 53-66. https://doi.org/10.1016/j.epsl.2019.04.034

Bookhagen, B., & Burbank, D. W. (2006). Topography, relief, and TRMM-derived rainfall variations along the Himalaya. *Geophysical Research Letters*, **33**(8), L08405. https://doi.org/10.1029/2006GL026037

Bookhagen, B., & Strecker, M. R. (2008). Orographic barriers, high-resolution TRMM rainfall, and relief variations along the eastern Andes. *Geophysical Research Letters*, **35**(6), L06403. https://doi.org/10.1029/2007GL032011

Bovy, B., & Lange, R. (2023). fastscape-lem/fastscape: Release v0.1.0. Zenodo. https://doi.org/10.5281/ZENODO.8375653

Bradley, R. S. (2015). Climate and climatic variation. *Paleoclimatology*, 13-54. https://doi.org/10.1016/B978-0-12-386913-5.00002-8

Braun, J., & Willett, S. D. (2013). A very efficient *O(n)*, implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. *Geomorphology*, **180-181**, 170-179. https://doi.org/10.1016/j.geomorph.2012.10.008

Champagnac, J.-D., Molnar, P., Sue, C., & Herman, F. (2012). Tectonics, climate, and mountain topography. *Journal of Geophysical Research*, **117**(B2), B02403. https://doi.org/10.1029/2011JB008348

Chang, H., Johnson, G., Hinkley, T., & Jung, I.-W. (2014). Spatial analysis of annual runoff ratios and their variability across the contiguous U.S. *Journal of Hydrology*, **511**, 387-402. https://doi.org/10.1016/j.jhydrol.2014.01.066

Clementucci, R., Ballato, P., Siame, L. L., Faccenna, C., Yaaqoub, A., Essaifi, A., et al. (2022). Lithological control on topographic relief evolution in a slow tectonic setting (Anti-Atlas, Morocco). *Earth and Planetary Science Letters*, **596**, 117788. https://doi.org/10.1016/j.epsl.2022.117788

Cox, N. J. (2007). Kernel estimation as a basic tool for geomorphological data analysis. *Earth Surface Processes and Landforms*, **32**(12), 1902-1912. https://doi.org/10.1002/esp.1518

Craddock, W. H., Burbank, D. W., Bookhagen, B., & Gabet, E. J. (2007). Bedrock channel geometry along an orographic rainfall gradient in the upper Marsyandi River valley in central Nepal. *Journal of Geophysical Research*, **112**(3), F03007. https://doi.org/10.1029/2006JF000589

Crameri, F. (2023). Scientific colour maps. Zenodo. https://doi.org/10.5281/zenodo.8409685

Crameri, F., Shephard, G. E., & Heron, P. J. (2020). The misuse of colour in science communication. *Nature Communications*, **11**(1), 5444. https://doi.org/10.1038/s41467-020-19160-7

D'Arcy, M., & Whittaker, A. C. (2014). Geomorphic constraints on landscape sensitivity to climate in tectonically active areas. *Geomorphology*, **204**, 366-381. https://doi.org/10.1016/j geomorph.2013.08.019

Deal, E., Braun, J., & Botter, G. (2018). Understanding the role of rainfall and hydrology in determining fluvial erosion efficiency. *Journal of Geophysical Research: Earth Surface*, **123**(4), 744-778. https://doi.org/10.1002/2017JF004393

Deal, E., Favre, A.-C., & Braun, J. (2017). Rainfall variability in the Himalayan orogen and its relevance to erosion processes. *Water Resources Research*, **53**(5), 4004-4021. https://doi.org/10.1002/2016WR020030

DiBiase, R. A., & Whipple, K. X. (2011). The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate. *Journal of Geophysical Research*, **116**(F4), 4036. https://doi.org/10.1029/2011JF002095

DiBiase, R. A., Whipple, K. X., Heimsath, A. M., & Ouimet, W. B. (2010). Landscape form and millennial erosion rates in the San Gabriel Mountains, CA. *Earth and Planetary Science Letters*, **289**(1-2), 134-144. https://doi.org/10.1016/j.epsl.2009.10.036

Dong, Y., Shi, X., Sun, S., Sun, J., Hui, B., He, D., et al. (2022). Co-evolution of the Cenozoic tectonics, geomorphology, environment and ecosystem in the Qinling Mountains and adjacent areas, Central China. *Geosystems and Geoenvironment*, **1**(2), 100032. https://doi.org/10.1016/j.geoeco.2022.100032

Dong, Y., Zhang, X., Liu, X., Li, W., Chen, Q., Zhang, G., et al. (2015). Propagation tectonics and multiple accretionary processes of the Qinling Orogen. *Journal of Asian Earth Sciences*, **104**, 84-98. https://doi.org/10.1016/j.jseas.2014.10.007

Dunne, K. B. J., & Jerolmack, D. J. (2020). What sets river width? *Science Advances*, **6**(41), eabc1505. https://doi.org/10.1126/sciadv.abc1505

Duvall, A., Kirby, E., Burbank, D., Duvall, C., Kirby, E., & Burbank, D. (2004). Tectonic and lithologic controls on bedrock channel profiles and processes in coastal California. *Journal of Geophysical Research*, **109**(F3), 3002. https://doi.org/10.1029/2003JF000086

Escayola, M. P., van Staal, C. R., & Davis, W. J. (2011). The age and tectonic setting of the Puncoviscana Formation in northwestern Argentina: An accretionary complex related to Early Cambrian closure of the Puncoviscana Ocean and accretion of the Arequipa-Antofalla block. *Journal of South American Earth Sciences*, **32**(4), 438-459. https://doi.org/10.1016/j.jsames.2011.04.013

Espurt, N., Hippolyte, J.-C., Kaymakci, N., & Sangu, E. (2014). Lithospheric structural control on inversion of the southern margin of the Black Sea Basin, Central Pontides, Turkey. *Lithosphere*, **6**(1), 26-34. https://doi.org/10.1130/L316.1

European Space Agency, S. (2020). Copernicus DEM GLO-30 Dataset. European Space Agency. https://doi.org/10.5270/ESA-c5d3d65

Ferrier, K. L., Huppert, K. L., & Perron, J. T. (2013). Climatic control of bedrock river incision. *Nature*, **496**(7444), 206-209. https://doi.org/10.1038/nature11982

Finnegan, N. J., Sklar, L. S., & Fuller, T. K. (2007). Interplay of sediment supply, river incision, and channel morphology revealed by the transient evolution of an experimental bedrock channel. *Journal of Geophysical Research*, **112**(F3), F03S11. https://doi.org/10.1029/2006JF000569
Fisher, J. A., Pazzaglia, F. J., Anastasio, D. J., & Gallen, S. F. (2022). Linear inversion of fluvial topography in the Northern Apennines: Comparison of base-level fall to crustal shortening. *Tectonics*, 41(11), e2022TC007379. https://doi.org/10.1029/2022TC007379

Flint, J. J. (1974). Stream gradient as a function of order, magnitude, and discharge. *Water Resources Research*, 10(5), 969–973. https://doi.org/10.1029/WR010i005p00969

Forte, A. M., Yanites, B. J., & Whipple, K. X. (2016). Complexities of landscape evolution during incision through layered stratigraphy with contrasts in rock strength. *Earth Surface Processes and Landforms*, 41(12), 1736–1757. https://doi.org/10.1002/esp.3947

Fournier, A., Fussell, D., & Carpenter, L. (1982). Computer rendering of stochastic models. *Communications of the ACM*, 25(6), 371–384. https://doi.org/10.1145/358523.358553

Fox, M., Goren, L., May, D. A., & Willett, S. D. (2014). Inversion of fluvial channels for paleorock uplift rates in Taiwan. *Journal of Geophysical Research: Earth Surface*, 119(9), 1853–1875. https://doi.org/10.1002/2014JF003196

Gailleton, B., Mudd, S. M., Clubb, F. J., Grieve, S. W., & Hurst, M. D. (2021). Impact of changing concavity indices on channel steepness and divide migration metrics. *Journal of Geophysical Research: Earth Surface*, 126(10), e2020JF006060. https://doi.org/10.1029/2020JF006060

Gallen, S. F. (2018). Lithologic controls on landscape dynamics and aquatic species evolution in post-orogenic mountains. *Earth and Planetary Science Letters*, 493, 150–160. https://doi.org/10.1016/J.EPSL.2018.04.029

Gallen, S. F., & Thigpen, J. R. (2018). Lithologic controls on focused erosion and intraplate earthquakes in the eastern Tennessee seismic zone. *Geophysical Research Letters*, 45(18), 9569–9578. https://doi.org/10.1029/2018GL079157

Gasparini, N. M., & Brandon, M. T. (2011). A generalized power law approximation for fluvial incision of bedrock channels. *Journal of Geophysical Research*, 116(2), F02020. https://doi.org/10.1029/2009JF001655

Gasparini, N. M., & Whipple, K. X. (2014). Diagnosing climatic and tectonic controls on topography: Eastern flank of the northern Bolivian Andes. *Lithosphere*, 6(4), 230–250. https://doi.org/10.1130/L322.1

Gasparini, N. M., Whipple, K. X., & Bras, R. L. (2007). Predictions of steady state and transient landscape morphology using sediment-flux-dependent river incision models. *Journal of Geophysical Research*, 112(F3), F03S09. https://doi.org/10.1029/2006JF000567

Gilbert, G. (1877). *Geology of the Henry Mountains* (Technical Report). Government Printing Office.

Godard, V., Bourlès, D. L., Spinabella, F., Burbank, D. W., Bookhagen, B., Fisher, G. B., et al. (2014). Dominance of tectonics over climate in Himalayan denudation. *Geology*, 42(3), 243–246. https://doi.org/10.1130/G35342.1

Goren, L., Fox, M., & Willett, S. D. (2014). Tectonics from fluvial topography using formal linear inversion: Theory and applications to the Inyo Mountains, California. *Journal of Geophysical Research: Earth Surface*, 119(8), 1651–1681. https://doi.org/10.1002/2014JF003079

Grujic, D., Govin, G., Barrier, L., Bookhagen, B., Coutand, I., Cowan, B., et al. (2018). Formation of a rain shadow: O and H Stable isotope records in authigenic clays from the Siwalik Group in Eastern Bhutan. *Geochemistry, Geophysics, Geosystems*, 19(9), 3430–3447. https://doi.org/10.1029/2017GC007254

Hack, J. T. (1960). Interpretation of erosional topography in humid temperate regions. *American Journal of Science*, 258(A), 80–97.

Han, J., Gasparini, N. M., & Johnson, J. P. (2015). Measuring the imprint of orographic rainfall gradients on the morphology of steady-state numerical fluvial landscapes. *Earth Surface Processes and Landforms*, 40(10), 1334–1350. https://doi.org/10.1002/esp.3723

Harel, M. A., Mudd, S. M., & Attal, M. (2016). Global analysis of the stream power law parameters based on worldwide ¹⁰Be denudation rates. *Geomorphology*, 268, 184–196. https://doi.org/10.1016/j geomorph.2016.05.035

Harries, R. M., Gailleton, B., Kirstein, L. A., Attal, M., Whittaker, A. C., & Mudd, S. M. (2021). Impact of climate on landscape form, sediment transfer and the sedimentary record. *Earth Surface Processes and Landforms*, 46(5), 990–1006. https://doi.org/10.1002/esp.5075

Hartmann, J., & Moosdorf, N. (2012). The new global lithological map database glim: A representation of rock properties at the earth surface. *Geochemistry, Geophysics, Geosystems*, 13(12), Q12004. https://doi.org/10.1029/2012GC004370

Hergarten, S., & Robl, J. (2022). The linear feedback precipitation model (LFPM 1.0) – A simple and efficient model for orographic precipitation in the context of landform evolution modeling. *Geoscientific Model Development*, 15(5), 2063–2084. https://doi.org/10.5194/gmd-15-2063-2022

Hergarten, S., Robl, J., & Stüwe, K. (2016). Tectonic geomorphology at small catchment sizes-extensions of the stream-power approach and the χ method. *Earth Surface Dynamics*, 4, 1–9. https://doi.org/10.5194/esurf-4-1-2016

Herman, F., Seward, D., Valla, P. G., Carter, A., Kohn, B., Willett, S. D., & Ehlers, T. A. (2013). Worldwide acceleration of mountain erosion under a cooling climate. *Nature*, 504(7480), 423–426. https://doi.org/10.1038/nature12877

Hippolyte, J.-C., Müller, C., Sangu, E., & Kaymakci, N. (2017). Stratigraphic comparisons along the Pontides (Turkey) based on new nanno-plankton age determinations in the Eastern Pontides: Geodynamic implications. *Geological Society, London, Special Publications*, 428(1), 323–358. https://doi.org/10.1144/SP428.9

Howard, A. D. (1987). Modelling fluvial systems: Rock-, gravel- and sand-bed channels. In K. Richards (Ed.), *River channels* (pp. 69–94).

Howard, A. D., & Kerby, G. (1983). Channel changes in badlands. *GSA Bulletin*, 94(6), 739–752. https://doi.org/10.1130/0016-7606(1983)94(739:CCIB)2.0.CO;2

Johnson, J. P. L., & Whipple, K. X. (2010). Evaluating the controls of shear stress, sediment supply, alluvial cover, and channel morphology on experimental bedrock incision rate. *Journal of Geophysical Research*, 115(F2), F02018. https://doi.org/10.1029/2009JF001335

Karlstrom, K. E., Coblentz, D., Dueker, K., Ouimet, W., Kirby, E., Van Wijk, J., et al. (2012). Mantle-driven dynamic uplift of the Rocky Mountains and Colorado Plateau and its surface response: Toward a unified hypothesis. *Lithosphere*, 4(1), 3–22. https://doi.org/10.1130/L150.1

Karlstrom, K. E., Wilgus, J., Thacker, J. O., Schmandt, B., Coblentz, D., & Albonico, M. (2022). Tectonics of the Colorado Plateau and its margins. *Annual Review of Earth and Planetary Sciences*, 50(1), 295–322. https://doi.org/10.1146/annurev-earth-032320-111432

Kirby, E., & Whipple, K. (2001). Quantifying differential rock-uplift rates via stream profile analysis. *Geology*, 29(5), 415–418. https://doi.org/10.1130/0091-7613(2001)029(0415:QDRURV)2.0.CO;2

Kirby, E., & Whipple, K. X. (2012). Expression of active tectonics in erosional landscapes. *Journal of Structural Geology*, 44, 54–75. https://doi.org/10.1016/J.JSG.2012.07.009

Kummerow, C., Simpson, J., Thiele, O., Barnes, W., Chang, A. T., Stocker, E., et al. (2000). The status of the tropical rainfall measuring mission (TRMM) after two years in orbit. *Journal of Applied Meteorology*, 39(12 PART 1), 1965–1982. https://doi.org/10.1175/1520-0450(2001)040(1965:tsottr)2.0.CO;2

Lague, D. (2014). The stream power river incision model: Evidence, theory and beyond. *Earth Surface Processes and Landforms*, 39(1), 38–61. https://doi.org/10.1002/esp.3462

Lague, D., Hovius, N., & Davy, P. (2005). Discharge, discharge variability, and the bedrock channel profile. *Journal of Geophysical Research*, 110(F4), F04006. https://doi.org/10.1029/2004JF000259
Leonard, J. S., & Whipple, K. X. (2021). Influence of spatial rainfall gradients on river longitudinal profiles and the topographic expression of spatially and temporally variable climates in mountain landscapes. *Journal of Geophysical Research: Earth Surface, 126*(12), e2021JF006183.
https://doi.org/10.1029/2021JF006183

Leonard, J. S., Whipple, K. X., & Heimsath, A. M. (2023a). Controls on topography and erosion of the north-central Andes. *Geology, 52*(2), 153–158.
https://doi.org/10.1130/G51618.1

Leonard, J. S., Whipple, K. X., & Heimsath, A. M. (2023b). Isolating climatic, tectonic, and lithologic controls on mountain landscape evolution. *Science Advances, 9*(3), eadd8915.
https://doi.org/10.1126/sciadv.add8915

Marder, E., & Gallen, S. F. (2023). Climate control on the relationship between erosion rate and fluvial topography. *Geology, 51*(5), 424–427.
https://doi.org/10.1130/G50832.1

McQuarrie, N., Horton, B. K., Zandt, G., Beck, S., & DeCelles, P. G. (2005). Lithospheric evolution of the Andean fold–thrust belt, Bolivia, and the origin of the Central Andean plateau. *Tectonophysics, 399*(1), 15–37.
https://doi.org/10.1016/j.tecto.2004.12.013

Mills, H. H. (2003). Inferring erosional resistance of bedrock units in the east Tennessee mountains from digital elevation data. *Geomorphology, 55*(1), 263–281.
https://doi.org/10.1016/S0169-555X(03)00144-2

Montgomery, D. R., Balco, G., & Willett, S. D. (2001). Climate, tectonics, and the morphology of the Andes. *Geology, 29*(7), 579–582.
https://doi.org/10.1130/0091-7613(2001)029(0579:CTATMO)2.0:CO;2

Moore, J. G., & Clague, D. A. (1992). Volcano growth and evolution of the island of Hawaii. *GSA Bulletin, 104*(11), 1471–1484.
https://doi.org/10.1130/0016-7606(1992)104(1471:VGAEOT)2.3:CO;2

Morisawa, M. E. (1962). Quantitative geomorphology of some watersheds in the Appalachian Plateau. *Geological Society of America Bulletin, 73*(9), 1025–1046.
https://doi.org/10.1130/0016-7606(1962)73[1025:QGOSWI]2.0:CO;2

Morris, S., & Yell. (1998). Exhumation of the Pyrenean orogen: Implications for sediment discharge. *Basin Research, 10*(1), 69–85.
https://doi.org/10.1046/j.1365-2117.1998.00053.x

Mudd, S. M., Attal, M., Milodowski, D. T., Grieve, S. W., & Valters, D. A. (2014). A statistical framework to quantify spatial variation in channel gradients using the integral method of channel profile analysis. *Journal of Geophysical Research: Earth Surface, 119*(2), 138–152.
https://doi.org/10.1002/2013JF002981

Mudd, S. M., Clubb, F. J., Gailleton, B., & Hurst, M. D. (2018). How concave are river channels? *Earth Surface Dynamics, 6*(2), 505–523.
https://doi.org/10.5194/ESURF-6-505-2018

Mudd, S. M., Clubb, F. J., Grieve, S. W. D., Milodowski, D. T., Gailleton, B., Hurst, M. D., et al. (2023). LSDtopotools/LSDTopoTools2: LSDTopoTools2 v0.9 (Version v0.9) [Computer software]. Zenodo. https://doi.org/10.5281/ZENODO.8076231

Olivetti, V., Godard, V., & Bellier, O. (2016). Cenozoic rejuvenation events of Massif Central topography (France): Insights from cosmogenic denudation rates and river profiles. *Earth and Planetary Science Letters, 444*, 179–191.
https://doi.org/10.1016/j.epsl.2016.03.049

Peifer, D., Persano, C., Hurst, M. D., Bishop, P., & Fabel, D. (2021). Growing topography due to contrasting rock types in a tectonically dead landscape. *Earth Surface Dynamics, 9*(2), 167–181.
https://doi.org/10.5194/ESURF-9-167-2021

Perlin, K. (1985). Image synthesizer. *Computer Graphics, 19*(3), 287–296.
https://doi.org/10.1145/325165.325247

Perron, J. T., & Royden, L. (2013). An integral approach to bedrock river profile analysis. *Earth Surface Processes and Landforms, 38*(6), 570–576.
https://doi.org/10.1002/esp.3302

Pfeiffer, A. M., Finnegan, N. J., & Willenbring, J. K. (2017). Sediment supply controls equilibrium channel geometry in gravel rivers. *Proceedings of the National Academy of Sciences, 114*(13), 3346–3351.
https://doi.org/10.1073/pnas.1612907114

Phillips, C. B., & Jerolmack, D. J. (2016). Self-organization of river channels as a critical filter on climate signals. *Science, 352*(6286), 694–697.
https://doi.org/10.1126/science.aad3348

Pritchard, D., Roberts, G. G., White, N. J., & Richardson, C. N. (2009). Uplift histories from river profiles. *Geophysical Research Letters, 36*(24), L24301.
https://doi.org/10.1029/2009GL040928

Rahl, J. M., Harbor, D. J., Galli, C. I., & O'Sullivan, P. (2018). Foreland basin record of uplift and exhumation of the Eastern Cordillera, northwest Argentina. *Tectonics, 37*(11), 4173–4193.
https://doi.org/10.1029/2017TC004955

Rasmussen, K. L., Chaplin, M. M., Zuluaga, M. D., & Houze, R. A. (2016). Contribution of extreme convective storms to rainfall in South America. *Journal of Hydrometeorology, 17*(1), 353–367.
https://doi.org/10.1175/JHM-D-15-0067.1

Rezaeian, M., Carter, A., Hovius, N., & Allen, M. B. (2012). Cenozoic exhumation history of the Alborz Mountains, Iran: New constraints from low-temperature chronometry. *Tectonics, 31*(2), TC2004.
https://doi.org/10.1029/2011TC002974

Roberts, G. G., & White, N. (2010). Estimating uplift rate histories from river profiles using African examples. *Journal of Geophysical Research, 115*(B2), B02406.
https://doi.org/10.1029/2009JB006692

Roe, G. H., Montgomery, D. R., & Hallet, B. (2002). Effects of orographic precipitation variations on the concavity of steady-state river profiles. *Geology, 30*(2), 143–146.
https://doi.org/10.1130/0091-7613(2002)030(0143:EOOPVO)2.0:CO;2

Roe, G. H., Montgomery, D. R., & Hallet, B. (2003). Orographic precipitation and the relief of mountain ranges. *Journal of Geophysical Research, 108*(B6), 2315.
https://doi.org/10.1029/2001JB001521

Rossi, M. W., Whipple, K. X., & Vivoni, E. R. (2016). Precipitation and evapotranspiration controls on daily runoff variability in the contiguous United States and Puerto Rico. *Journal of Geophysical Research: Earth Surface, 121*(1), 128–145.
https://doi.org/10.1002/2015JF003446

Royden, L., & Perron, J. T. (2013). Solutions of the stream power equation and application to the evolution of river longitudinal profiles. *Journal of Geophysical Research: Earth Surface, 118*(2), 497–518.
https://doi.org/10.1002/jgrf.20031

Royden, L. H., Clark, M. K., & Whipple, K. X. (2000). Evolution of river elevation profiles by bedrock incision: Analytical solutions for transient river profiles related to changing uplift and precipitation rates. *Eos, Transactions American Geophysical Union, 81*(48), Fall Meeting Supplement, Abstract T62F-09.
https://doi.org/10.1002/eost.v81.16

Rudge, J. F., Roberts, G. G., White, N. J., & Richardson, C. N. (2015). Uplift histories of Africa and Australia from linear inverse modeling of drainage inventories. *Journal of Geophysical Research: Earth Surface, 120*(5), 894–914.
https://doi.org/10.1002/2014JF003297

Ruiz Sánchez-Oro, M. (2024). Marinaruizso/jgr: paper. Jgrearthsurfacepaper_v0.1 [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.10909989

Safran, E. B., Bierman, P. R., Aalto, R., Dunne, T., Whipple, K. X., & Coffee, M. (2005). Erosion rates driven by channel network incision in the Bolivian Andes. *Earth Surface Processes and Landforms, 30*(8), 1007–1024.
https://doi.org/10.1002/esp.1259

Sanchis, E., & Séranne, M. (2000). Structural style and tectonic evolution of a polyphase extensional basin of the Gulf of Lion passive margin: The Tertiary Alès basin, southern France. *Tectonophysics, 322*(3), 219–242.
https://doi.org/10.1016/S0040-1951(00)00097-4

Seybold, H., Berghuijs, W. R., Prancevic, J. P., & Kirchner, J. W. (2021). Global dominance of tectonics over climate in shaping river longitudinal profiles. *Nature Geoscience, 14*(7), 503–507.
https://doi.org/10.1038/s41561-021-00720-5

Silverman, B. W. (1998). *Density estimation for statistics and data analysis*. Routledge. https://doi.org/10.1201/9781315140919
Skofronick-Jackson, G., Petersen, W. A., Berg, W., Kidd, C., Stocker, E. F., Kirschbaum, D. B., et al. (2017). The Global Precipitation Measurement (GPM) mission for science and society. *Bulletin of the American Meteorological Society*, 98(8), 1679–1695. https://doi.org/10.1175/BAMS-D-15-00306.1

Smith, R. B., & Barstad, I. (2004). A linear theory of orographic precipitation. *Journal of the Atmospheric Sciences*, 61(12), 1377–1391. https://doi.org/10.1175/1520-0469(2004)061(1377:ALTOOP)2.0:CO;2

Sorensen, C. S., & Yanites, B. J. (2019). Latitudinal trends in modern fluvial erosional efficiency along the Andes. *Geomorphology*, 329, 170–183. https://doi.org/10.1016/j geomorph.2018.12.030

Styron, R., & Pagani, M. (2020). The GEM global active faults database. *Earthquake Spectra*, 36(1_suppl), 160–180. https://doi.org/10.1177/8755293020944182

Taylor, M., & Yin, A. (2009). Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism. *Geosphere*, 5(3), 199–214. https://doi.org/10.1130/GES00217.1

Trenberth, K. (2011). Changes in precipitation with climate change. *Climate Research*, 47(1), 123–138. https://doi.org/10.3354/cr00953

Tucker, G. E., & Whipple, K. X. (2002). Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison. *Journal of Geophysical Research*, 107(B9), ETG 1-1–ETG 1-16. https://doi.org/10.1029/2001JB000162

Ward, D. J., & Galewsky, J. (2014). Exploring landscape sensitivity to the Pacific Trade Wind Inversion on the subsiding island of Hawaii. *Journal of Geophysical Research: Earth Surface*, 119(9), 2048–2069. https://doi.org/10.1002/2014JF003155

Whipple, K. X., & Gasparini, N. M. (2014). Tectonic control of topography, rainfall patterns, and erosion during rapid post–12 Ma uplift of the Bolivian Andes. *Lithosphere*, 6(4), 251–268. https://doi.org/10.1130/L325.1

Whipple, K. X., & Tucker, G. E. (1999). Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. *Journal of Geophysical Research*, 104(B8), 17661–17674. https://doi.org/10.1029/1999JB900120

Whipple, K. X., & Tucker, G. E. (2002). Implications of sediment-flux-dependent river incision models for landscape evolution. *Journal of Geophysical Research*, 107(B2), ETG 3-1–ETG 3-20. https://doi.org/10.1029/2000JB000044

Wickert, A. D., & Schildgen, T. F. (2019). Long-profile evolution of transport-limited gravel-bed rivers. *Earth Surface Dynamics*, 7(1), 17–43. https://doi.org/10.5194/ESURF-7-17-2019

Wobus, C. W., Tucker, G. E., & Anderson, R. S. (2006). Self-formed bedrock channels. *Geophysical Research Letters*, 33(18), 18408. https://doi.org/10.1029/2006GL027182

Wobus, C. W., Tucker, G. E., & Anderson, R. S. (2010). Does climate change create distinctive patterns of landscape incision? *Journal of Geophysical Research*, 115(4), F04008. https://doi.org/10.1029/2009JF001562

Zaprowski, B. J., Pazzaglia, F. J., & Evenson, E. B. (2005). Climatic influences on profile concavity and river incision. *Journal of Geophysical Research*, 110(F3), 3004. https://doi.org/10.1029/2004JF000138