Estimating the disequilibrium in denudation rates due to divide migration at the scale of river basins

Basin-averaged denudation rates may locally exhibit a wide dispersion, even in areas where the topographic steady state is supposedly achieved regionally. This dispersion is often attributed to the accuracy of the data or to some degree of natural variability of local erosion rates which can be related to stochastic processes such as landsliding. Another physical explanation of this dispersion is local and transient disequilibrium between tectonic forcing and erosion at the scale of catchments. Recent studies have shown that basin divide migration can potentially induce such perturbations, and they propose metrics to assess divide mobility based on cross-divide contrasts in headwater topographic features. Here, we use a set of landscape evolution models assuming spatially uniform uplift, rock strength and rainfall to assess the effect of divide mobility on basin-wide denudation rates. We propose using basin-averaged aggressivity metrics based on cross-divide contrasts (1) in channel χ , an integral function of position in the channel network; (2) in channel local gradient; and (3) in channel height, measured at a reference drainage area. From our simulations, we show that the metric based on differences in χ is the most reliable to diagnose local disequilibrium. The other metrics are more suitable for relatively active tectonic regions such as mountain belts, where contrasts in local gradient and elevation are more important. We find that the ratio of basin denudation associated with drainage migration to uplift can reach a factor of 2, regardless of the imposed uplift rate, erodibility, diffusivity coefficient or critical hillslope gradient. A comparison with field observations in the Great Smoky Mountains (southern Appalachians, USA) underlines the difficulty of using the metric based on χ , which depends on the – poorly constrained – elevation of the outlet of the investigated catchment. Regardless of the considered metrics, we show that observed dispersion is controlled by catchment size: a smaller basin may be more sensitive to divide migration and hence to disequilibrium. Our results thus highlight the relevance of divide stability analysis from digital elevation models as a fundamental preliminary step for basin-wide denudation rate studies based on cosmogenic radionuclide concentrations.

physical exp lanation of this dispersion is local and transient disequilibriu m between tectonic forcing and erosion at the scale of catchments. Recent studies have shown that divide migration can potentially induce such perturbations and they propose metrics to assess divide mob ility based on cross -divide contrasts in headwater topographic features. Here, we use a set of landscape evolution models assuming spatially uniform uplift, rock strength and rainfall to assess the effect of divide mobility on basin-wide denudation rates. We propose the use of basin-averaged aggressivity met rics based on cross-divide 15 contrasts (1) in channel , an integral function of position in the channel network, (2) in channel local gradient and (3) in channel height, measured at a reference drainage area. Fro m our simulations, we show that the metric based on differences in is the most reliable to diagnose local disequilibriu m. The other metrics are mo re suitable for relatively active tectonic regions as mountain belts, where contrasts in local gradient and elevation are more important. We find that the ratio o f basin denudation associated with drainage migrat ion to uplift can reach a factor of t wo, regardless the imposed uplift rate, 20 erodibility, d iffusivity coefficient and critical h illslope gradient. A co mparison with field observations in the Great Smo ky Mountains (southern Appalachians, USA) underlines the difficulty of using the metric based on χ, wh ich depends on thepoorly constrained-elevation of the outlet of the investigated catchment. Regardless of the considered met rics, we show that observed dispersion is controlled by catch ment size: a smaller basin may be more s ensitive to divide mig ration and hence to disequilibriu m. Our results thus highlight the relevance of divide stability analysis fro m dig ital elevation models as a 25 fundamental preliminary step for basin-wide denudation rate studies based on cosmogenic radio nuclide concentrations.

Introduction
Topographic steady state, in which average topography is constant over time, is one of the key concepts of modern geomorphology (e.g. Gilbert, 1877;Hack 1960;Montgomery, 2001). Though simple, this paradig m provides a useful framework to study landscape evolution related to tectonic and/or climatic forcing (e.g. Willett et al., 2001;Reinhart and 2 Ellis, 2015), to spatial variations in rock strength (Perne et al., 2017) or to the geomet ry of active crustal structures (Lave and Avouac, 2001;Stolar et al., 2007;Scherler et al., 2014;Le Rou x-Mallouf et al., 2015). To define topographic steady state, the temporal and spatial scales of the processes involved are essential parameters. Co mpared to large scale geodynamic processes operating over 1-100 Myr timescales, river incision and sediment transport are rapid processes driving landscapes to stable forms over this long t imescale, whereas rapid climatic fluctuations during the Quaternary may prevent the 35 occurrence of steady-state conditions in modern landscapes (Whipple, 2001).
The timescale of div ide migrat ion has received increasing attention in the recent years. Although rivers exh ibit a rap id adjustment to tectonic or climatic changes to maintain their profiles, Whipple et al. (2017) show that divides continue to migrate over t ime periods of 10 6 -10 7 years as response to the same changes . This suggests that long-term transience might be pervasive in the planar structure of landscapes, even in the absence of new variat ions in landscape characteristics or forcings 40 (e.g. tectonic or climate) (Hasbargen and Paola, 2000;Hasbargen and Paola, 2003;Pelletier, 2004, Dah lquist et al., 2018. In addition to the influence of spatial variability o f rock uplift rate, rock strength or rainfall (e.g. Reiner et al., 2003;Godard et al., 2006;Miller et al., 2013), th is long timescale could also explain the persistence of spatial variations in denudation r ates observed in tectonically inactive orogens which achieved regional-scale topographic steady state (Willett et al., 2014).
As an examp le, in the Great Smoky Mountains in the southern Appalachians, uplift and erosion rates integrated over varying 45 time periods fro m 10s kyr to 100 Myr give a similar average magn itude of ca. 0.03 mm.yr -1 (Mat mon et al., 2003a, b andPortenga andBierman, 2011). These results suggest a regional quasi-topographic steady state over the last ~180 Myr, maintained by the isostatic response of the thickened crust since the end of the Appalachian orogeny (Mat mon et al., 2003a, b). Beyond this average value, indiv idual basin-wide denudation rates exh ibit a strong dispersion (up to a factor of two, Fig.   1), which is not related to spatial variation in rainfall o r in erodib ility of the substrate (Matmon et al., 2003b). In a recent 50 study, Willett et al. (2014) assess divide mobility fro m the contrast in the channel head topographic metric , taken here as a proxy for steady-state river profile elevation (Perron and Royden, 2012;Royden and Perron, 2013), and propose an explanation in which a significant part of the observed dispersion in denudation rates could be due to drainage div ide migration associated with contrasting erosion rates across divides.
More recently, to characterize d ivide migrations Forte and Whipple (2018) introduced other metrics, referred as " Gilbert 55 metrics" (Gilbert, 1877), based on the cross -divide contrast in channel head local gradient and elevation. This last study indeed focused on cross -divide contrasts in headwater basin shape. Here, we propose extending these approaches by modeling divide migrat ion and by developing new metrics to assess divide stability at the scale of the entire watershed, which are an expansion of the aggressivity metric initially suggested by Willett et al. (2014). We use these metrics to assess the effect of persistent divide mob ility on basin-averaged erosion rates at a timescale of 10 4 yr. We use numerical landscape 60 evolution models, taking into account both hillslope diffusion and fluvial incision. For the sake of simp licity and to avoid the influence of other factors such as topography, lithology, climate or vegetation, we restrict our analysis to synthetic orogen s with spatially uniform uplift, rock strength and rainfall. After a brief presentation of our landscape evolution model (LEM), we describe the methods developed to assess basin-wide denudation rates and aggressivity metrics, such as average cross -divide contrasts in channel , gradient and elevation. Next, we investigate transient time and location of morphologic 65 adjustments to divide migrations. We explore the relevance and co mp lementarity o f tested relat ive stability metrics between neighboring basins. We then investigate the impact of uplift rate, erodib ility and hillslope processes on the dynamics of divide migrat ion and associated denudation rates. Finally, we apply our approach to the basin -wide denudation rates dataset of Mat mon (2003a,b) in the case of the Great Smo ky Mountains and propose new criteria to guide future sampling strategies to assess basin-wide denudation rates from river sands. 70

Landscape Evolution Model (LEM)
We use TTLEM (TopoToolbox Landscape Evolution Model) (Campforts et al., 2017), a landscape evolution model based on the Matlab function library TopoToolbox 2 (Schwanghart and Scherler, 2014). This LEM uses a finite volu me method (Camp forts and Govers, 2015) to solve the following equation of mass conservation for rock/regolith subject to uplift and 75 denudation: where / is the variation of elevation with time, ( / ) is the change of elevation due to tectonic horizontal 80 advection, is the rock uplift rate, / is the density ratio between the bedrock and the regolith. We use a linear formulation of hillslope diffusion (Culling, 1963) limited by a critical slope :

85
where is the flu x of soil-regolith material. When slopes values exceed , they are readjusted to the critical value by using a modified version of the excess topography algorithm (Blöthe et al., 2015). The diffusivity gives the rate of soil-regolith material creep. Its magnitude ranges fro m 10 -3 to 10 -1 m².yr -1 in natural settings and varies with soil thickness, lithology and vegetation (Roering et al., 1999;Jungers et al., 2009;West et al., 2013;Richardson et al., 2019). Hillslope diffusion is implemented in TTLEM using an implicit scheme, which is unconditionally stable at large t ime steps (Pellet ier, 2008). Non-90 linear diffusion formulat ion (Perron, 2011) is also implemented in TTLEM. However, we favored the use of a linear diffusion with a critical slope, which is mo re convenient for the time step used in our simu lations (=5000 yr) and the set of parameters considered (see section 2.2). Due to the relatively coarse spatial resolution of our models (= 90m), any of these diffusion formu lations generate negligible topographic differences on the direct vicinity of crest lines (Roering et al., 1999, Campforts et al., 2017 and do not affect our results (see Fig. S1). Fluvial incision is calculated with a stream power law: 95 where Is the erodibility coefficient reflecting climate, hydraulic roughness, sediment load and lithology. Its value ranges between 10 -16 and 10 0 m (1-2m) .yr −1 (Kirby and Whipple, 2001;Harel et al., 2016). is the upstream area. is the along 100 stream distance from the outlet of the river. and are two parameters that are usually reported as a / ratio ranging between 0.35 and 0.8. The river incision law is implemented in TTLEM using an explicit scheme based on a higher-order flu x-limit ing finite volu me method that is total variation d iminishing (TVD -FVM) [see Campforts and Govers (2015) and Campforts et al. (2017) for further details]. Its main advantage is to eliminate numerical d iffusion, which is present in most other schemes solving differential equations of river incision. This last point has a significant impact on the accuracy of 105 basin-wide simulated denudation rates, making TTLEM a well-suited LEM for the purpose of this study.

Geometry and meshing
Since the computation is performed using a discretized land surface, smaller mesh sizes lead to detailed topography but lengthen the computation time and memo ry requirements. Hereinafter, we consider a reference square landscape model o f 50 110 km side with a grid resolution of 90 m, which is a good compromise between computation time (3-5 hours on a PC workstation) and the total amount of basins that can be studied (>1000). Our results are not affected when the grid resolution is 30 m nor when the model size is 100x100 km (see Fig. S2).

Boundary conditions
In order to isolate the effect of d ivide migrat ions on the variability of basin-wide denudation rates, we explore simple models 115 with constant and spatially uniform uplift and precipitation rates and we assume no horizontal advection ( / ) = 0. We use a Dirichlet boundary condition: simu lation edges are not affected by uplift on a one pixel band to represent a stable bas e level for rivers. The model presents no initial topography, except for gauss ian noise ranging between 0 and 50 m so as to initiate a random fluvial network.
In order to better constrain the variab ility of our results under similar co nditions, we run for each model five simulat ions using the same parameters, but with different initial random top ographies.

Timescale 130
The total duration of simu lations is 10 Myr. The imp licit scheme used to simulate linear hillslope processes provides stable solutions regardless of the time step. In contrast, the explicit scheme used to model fluvial incision requires a time step that satisfies the Courant-Friedrich-Lewy criterion. Hereinafter, we choose a time step of 5000 yr for h illslope diffusion. Our results are not affected when using a smaller (i.e. 1000 yr) (see Fig. S1). Incision computation is nested in this time step and uses another time step that is automatically determined to assure model stability (Campforts et al., 2017). 135

Basin-wide denudation rates
We derive basins from the synthetic DEMs (Digital Elevation Models) using an accumulation map co mputed with a single flow d irection algorith m implemented in TopoToolbox (Sch wanghart and Scherler, 2014). Next, we calculate for each basin the variation in average elevation over a time interval of 10 kyr. The drainage network migrates during the simu lation, so we 140 only survey the basins that keep the same outlet location during this time interval. Furthermo re, due to divide mobility, the geometry of watersheds can also change. Hence, we measure the average difference in elevation inside the basin perimeter after 10 kyr. Here we only assess the surface uplift (England and Molnar, 1990). To appro ximate the denudation rates for each basin, we sum the surface uplift with the rock uplift rate and divide the result by the time interval. By considering the relatively small period over wh ich we integrate denudation (10 kyr), we then assume that these 145 approximations have a neglig ible impact on the results. If the basin is in a topographic steady state, is equal to zero and is equal to the background uplift rate. Thus, a positive (negative) value of traduce a deficit (an excess) of denudation.
Calculated that way, is sensitive to divide migration but also to transient features like knickpoints that migrate along the river network. In our simulat ions, knickpoints may develop due to (1) the dissection of the initial flat surface or (2) discrete drainage captures (see Sec. 3.1). We use the knickpointfinder algorithm imp lemented in TopoToolbox (Schwanghart and 150 Scherler, 2014) to identify the affected basins.

From cross-divide metrics to basin averaged aggressivity metrics
Most recent studies have focused on the relationship between drainage divide mobility and headwater across-divide contrast in either , gradient, elevation or local relief values (e.g. Whipple et al., 2017;Fo rte and Whipple, 2018). Here, in line with Willett et al. (2014, see Supp. Mat. therein) we focus on the specific in fluence of divide migration on denudation rates at the 155 scale of the entire stream basin. Our approach aims to integrate cross -divide contrasts in drainage network properties along the entire basin perimeter. We then obtain basin-averaged aggressivity metrics that determine if a watershed is either growing or shrinking (Willett et al. 2014).
First, we assess , local topographic gradient and height of the drainage network at a reference drainage area ( Fig.   2). Ideally, must be equal to the area at which channelization occurs (Forte and Whipple, 2018). Ho wever, it is 160 challenging to locate the accurate position of channel heads (Clubb et al., 2016). Hence, we use a constant value of set to 1 km². The parameter , is an integral function of position along the channel network (Perron and Royden, 2012) described by the equation: where ( ) is the upstream drainage area at the location , 0 is an arbitrary scaling area set to 1 km². The over ratio refers here to the reference concavity of an equilibrated river profile. Its value is set to 0.5 in accordance with the model parameters. For each independent drainage network, we integrate fro m the outlet , located at the model boundary (< 1 m high), to the channel heads. Local gradient is determined for each DEM pixel fro m its eight-connected neighbors. Height is 170 simply extracted from the DEM.
Then, we calcu late the difference of metrics ( , and ) across the segments of divide shared by two reference basins.
Finally, the aggressivity metric is obtained by averaging these across -divide differences along the perimeter of each samp led basin (Fig. 2). This way, the sign of the aggressivity metric in a basin corresponds to the differen ce of the averaged value of considered metric d ifference ( , and ) in this basin with respect to its neighbours. This method has the advantage of 175 weightening ind ividual d ivide segments by the number of p ixels they contain, and then to provide a robust assessment of the basin aggressivity. Aggressivity metrics based on , and are hereafter referred to as , and , respectively.
However, due to topology issues, some parts of the perimeter o f the sampled basin s may be not shared by two reference basins (Fig. 2). We quantify this incompleteness by assessing the ratio of documented pixels over the total amount of pixel along the basin perimeter. We refer to this ratio as the "confidence index" CI , assuming that an higher CI is associated with a 180 more robust basin aggresivity assessment.

Evolution of reference model
A detailed analysis of the DEM suggests that during the initial phase, the flat init ial surface (Fig.3a) is progressively uplifted to form a p lateau. At the same time the edges of this plateau are gradually regressively eroded by drainage networks that 185 spread from the base level toward the center of the model (Figs. 3b and c). This transient landscape is completely dissected after 2 Myr. Fro m this time and until the end of the simu lation, landscape changes are mainly due to co mpetition between watersheds, resulting in continuous divide migrations with decreasing intensity as th e model is moving toward a total topographic equilibrium ( Fig. 3d to f; Supplementary Video n°1).
To define the time period of regional steady state, we measure the average elevation, the maximu m elevation and the average 190 denudation rate over the entire model for each time step (Fig. 4a). We identify two d istinct stages during the evolution of our reference simulat ion. During the first million years, due to long wavelength topographic building, the calculated landscapes are far fro m steady state. This leads to a major increase of the mean elevation fro m 0 m to ca. 75 m. In a second stage, this trend reverses and the mean elevation decreases asymptotically toward ca. 60 m until the end of the simulation.
The evolution of the maximu m elevation follows the same pattern but can be affected by temporal changes in the location 195 and altitude of highest peaks. The maximu m elevation increases between 0 and ca. 250 m over the first 3Myr (Fig. 4a) then decreases progressively to remain at ca. 200 m during the rest of the simulation.
We compute the average denudation rate from the rock uplift rate and from average elevation change over the entire model between two time steps: where( / ) is the average surface uplift over the entire model on a time-step , is the imposed uniform uplift rate (0.1 mm.yr -1 ) and is the average "real" denudation rate. During the first 0.25 Myr, the mean denudation rate falls abruptly fro m ca. 0.6 mm.y r -1 to nearly 0 mm.yr -1 as a consequence of diffusion over the initial flat topography. After that 205 time and until the first 1 Myr, the mean denudation rate increases but remains lo wer than the uplift rate, leading to the increase in average elevation over this time period. In the follo wing 1 Myr, exceeds the uplift rate to reach up to 0.104 mm.yr -1 before it gently decreases to 0.1 mm.yr -1 until the end of the simulat ion. This shows that topography tends to -but never reaches -a strict steady state over the simulation time. Abrupt changes in after ca. 2.5, 3.5, 4,5 and 9.5 Myr (red circles in Figure 4b) are related to major local captures in the drainage network, wh ich can be observed during the model 210 evolution (red circles in Fig. 3e and f and Supplementary Video n°1).
Based on these results, we will consider that a regional topographic steady-state is reached between 1.5 and 2 Ma, when the plateau relict topography is totally eroded and begins to decrease (Figs. 3 and 4). This time is consistent with the time required to reach topographic steady state proposed from models with constant uplift rate and no horizontal advection (Willett et al., 2001). 215

Basin-wide denudation rates variability
We calculate basin-wide denudation rates upstream of each stable drainage network confluence after 2.5 Myr, 5 Myr and 10 Myr of simu lation (Figs. 5a, b and c, respectively). Regard less of the duration, we observe a significant variability in the calculated denudation rates depending on basin size. As exposed by Forte & Whipple (2018), the erosion rate contrasts across divides is spatially limited to areas very near the div ides . Thus, the variability is maximu m for small basins (ca. 1 220 km²) and decreases with increasing basin area. In our approach, small basins are nested in larger ones. Hence, these results can be related to the averaging of denudation rates along the drainage network, in agreement with the measurements of Matmon et al. (2003b). Th is variab ility also decreases with t ime (Figs. 5a-c). For basins with an excess of denudation relative to the uplift rate , the / ratio can reach up to 2.5 after 2.5 Myr but only 2 after 5 Myr and 1.7 after 10 Myr.
Basins with a denudation excess that stand out of the general trend at 10 Ma (Fig. 5c) are associated with a capture event 225 visible in figure 4b. For basins with a deficit of denudation, the evolution of the ratio is less obvious. It can be lower than 0.5 after 2.5 Myr, but increases slightly to 0.6 until 10 Myr. These results reflect a significant spatial variability of the d ifference between basin-wide denudation rates and uplift rate. To assess more accurately the temporal evolution of this variability, we calculate every 0.5 Myr for three distinct categories of basin sizes: 1-2 km 2 , 10-20 km 2 and 100-200 km 2 . We then estimate the mean absolute deviation (MAD) fro m the uplift rate by considering separately basins with a denudation in 230 excess or in deficit of uplift rate (Fig. 5d). Until 1.5 Ma, basins are located on the plateau where denudation rate is null. This The spatial variability of the denudation rates is neither homogeneous nor randomly distributed (Fig. 6a). The location of drainage basins with denudation rates far fro m the equilibriu m value of 0.1 mm.yr -1 co incides with migrating d rainage divides (Fig. 3d) and with cross-divide contrasts in channel , gradient and height (Figs. 6b-d). Fo llo wing Willett et al. 240 (2014) and Forte and Whipple (2018), the divide migrations predicted by these contrasts are consistent with the direction of divide mobility obtained fro m our model. One may note that the higher the contrast in these parameters across the divide, the higher the deviation of the denudation rate from the uplift rate, and therefore fro m topographic equilibriu m. None of the sampled basins in this data-set contains a knickpo int. Thus, these results based on simulat ions assuming uniform and constant properties as well as constant boundary conditions confirm that the dispersion observed in denudation rates is 245 primarily controlled by div ide mig ration. Basins that expand (shrink) show higher (lower) denudation rates compared to uplift rate, and are hereafter referred to as aggressors (victims), following the terminology adopted by Willett et al. (2014). Willett et al. (2014) showed that the basin-averaged cross-divide contrast in , could be used to deduce an aggressivity metric for basins. We extend this basin-scale approach to the Gilbert's metrics recently proposed by Forte and Whipple 250 (2018) including cross-divide contrast in headwater gradient and elevation.

Deviation of denudation rates from the uplift rate, and basin aggressivity
We here assess the relationship between the / ratio and these aggressivity metrics. First, to exclude variability related to both basin area and time, we focus on a single class of basins with a size of 2-4 km² gathered fro m five co mputed reference models after a simulat ion duration of 2.5 Myr. Denudation rates may be affected by knickpoints, wh ich are a source of transient perturbation at the scale of the catchment. Therefore, in order to focus only on perturbations associated with 255 drainage divide dynamics, basins that contain knickpoints are ignored. In agreement with cross -divide metrics tested by Forte and Whipple (2018), graphs in figure 7 must be divided into four quadrants. Aggressor (victim) basins have negative (positive) and values and conversely positive (negative) value (Fig. 2). Theoretically, aggressor (victim) basins have higher (lower) denudation rates than the underlying uplift rate. This result is verified for ca. 81 %, 52 % and 81 % of basins for , and , respectively. For this limited dataset, the evolution between / and both and 260 may be defined by a linear relat ionship (Fig. 7b). Co mpared to other metrics, is less sensitive to drainage migration and shows a more scattered distribution.
In natural settings, the stage of evolution of landscapes cannot be easily defined and t he total amount of basins with a specific size may be limited. The large dataset fro m our modeling can provide further insights by gathering the results obtained every 0.5 Myr for seven classes of basin areas expanding geometrically with a mu ltiply ing fact or of 2 fro m 1-2 to 265 64-128 km² (Figs. 7b ). Basins that contain knickpoints are discarded fro m the analysis. When all classes of drainage areas are co mbined together, we still obtain a clear relat ionship between aggressivity metrics and / , with 77 %, 56 % and 78 % of basins lying in aggressors or victims quadrants for , and , respectively (Fig 7b). Our results highlight the major control of basin size on the dispersion / . Part of the variability intrinsic fo r each class of basin area may in turn be explained by heterogeneities in aggressivity between different parts of a basin. Figure 7c shows that this dispersion is related 270 to the standard deviation of aggressivity metrics, , and . In other words, basins where d ifferent div ide segments migrate at different rates or in d ifferent direct ions are more scattered. The lower the confidence index, the more scattered the results (Fig. 7d). Thus, some d ispersion may co me fro m appro ximations due to undocumented divide segments performed when averaging metrics differences between reference basins (Fig. 2). One may note two different trends for victim and aggressor basins. Aggressors show a more scattered distribution for and metrics. When compared to 275 victims, these basins have hillslopes closer to the critical value (Fig. S3). Hence, the dispersion may be explained by the non-linear relationship existing between denudation rates and basin slope (Montgomery and Brandon, 200 2;Binnie et al., 2015).

Sensitivity tests 280
The reference model involves various parameters related to uplift, fluv ial incision and hillslope denudation. A systematic analysis of trade-off between all parameters is out of the scope of this manuscript. In th is section, we assess the sensitivity of the results to both tectonic and erosion processes, by studying the specific impact of uplift , erodibility , diffusivity and critical hillslope grad ient taken separately. Vary ing these parameters may change the simulat ion time required to erode the plateau associated with the initial boundary conditions . In this section, to reduce sensitivity dependence on these initial 285 conditions, we only consider results obtained between 5 Ma and 10 Ma.

Sensitivity to uplift rate
We test rock uplift rates of 0.01 mm.yr -1 , 0.1 mm.y r -1 (hereafter called reference model) and 1 mm.yr -1 to cover the range of a large variety of geodynamic settings (Champagnac et al., 2012). It is well known that a river responds to a fall in base level (due to changes in rock uplift rate or other forcing) by cutting downward into its bed, deepening and widening its active 290 channel. In our simu lations, changes in uplift rate lead to variations in the density of the drainage network. Co mpared to the reference model, an uplift rate of 1 mm.yr -1 (0.01 mm.yr -1 ) results in a decrease (increase) of drainage density. These results are consistent with previous studies that show an inverse relat ionship between drainage density and erosion rates in equilibriu m topography when using a threshold slope for diffusion processes (Tucker and Bras, 1998;Clubb et al., 2016). An increase in uplift rate favors river entrenchment leading to increase the range of and (Fig. 8). Hence, these two 295 Gilbert's metrics appear to be well suited to diagnose local disequilibriu m for higher uplift rates (i.e. ≥ 1 mm.yr -1 ).

Conversely increase uplift rate induces a lower range of values for
. This last observation is explained by the decrease of drainage density, and associated stream length.
Maximu m variability of / reaches a factor regardless the assumed uplift rate between 1 mm.yr -1 and 0.01 mm.yr -1 . The observed small differences suggest that limited uplift rate promote diffusive processes (see Sect. 4.1.3). 300

Influence of erodibility
Fluvial erosion is proportional to the erodibility coefficient that may reflect, among others, rock strength and climate. We let this parameter vary between 5.10 -6 m (1-2m) .yr -1 and 5.10 -5 m (1-2m) .yr -1 . As expected fro m (Eq .1 and 3), we find that erodibility and uplift rates have opposite effects . Lo wer (higher) values of erodibility lead to h igher (lower) average topography. Thus, an increase (decrease) in erodib ility decreases (increases) the range of all aggressivity metrics (Fig. 9).
Lower values of erodibility also increase the range of the / ratio. Models with higher (lower) erod ibility reach a quasitopographic steady state earlier (at a later stage). Hence, differences in the variability of / may be related to different stages of evolution for each models over the period we consider (5 to 10 Ma) (Fig. 5d).

Influence of hillslope processes
Hillslope denudation is proportional to the diffusivity coefficient and depends on the critical slope (Eq . 2). To test the 310 effect of hillslope processes, we let vary between 10 -3 m 2 .yr -1 and 10 -1 m 2 .yr -1 . Co mpared to the reference model, we find no differences in the case of a lower diffusivity (i.e. 0.001 m².yr -1 ) (Fig. 10c). In contrast, for models with higher diffusivity coefficient (i.e. 0.1 m².yr -1 ), this parameter has a significant effect on both the range of / and the aggressivity metric ( Fig. 10a). This result is consistent with the observations described in section 4.1.1 . It derives from a stronger impact of diffusive processes, which decrease local slopes in the vicinity of divides. In our modeling, the local slopes remain lower 315 than the fixed crit ical value. Then, assuming a crit ical slope between 20° and 40°, we find that does not affect significantly the relationship between the / ratio and the studied metrics (Fig.11).
Altogether, these sensitivity tests demonstrate the robustness of our findings. Regardless of the tested parameter values, we observe a relationship between aggressivity met rics and deviation of denuda tion rates from uplift rates. Thus, aggressivity 320 metrics are, to the first-order, reliable metrics to assess the effect of divide mobility on basin-wide denudation rates inferred fro m simulat ions. In the fo llo wing section, we apply this approach to field o bservations and discuss the consequences for sampling and interpretation.

Implications for the interpretation of basin-wide denudation rates
Over the last decades, measurements of cosmogenic radionuclides (CRN) concentration in alluvial sediments (see Granger et 325 al., 2013 and references therein), of suspended sediments (Gabet et al., 2008) and of detrital thermochronology (Huntington and Hodges, 2006) have become co mmon practices to assess basin -wide denudation rates. However, their interpretation remains debated, even in settings where topographic steady state is supposedly achieved regionally.

Application to the Great Smoky Mountains
As previously mentioned (Matmon et al., 2003a, b), while the Great Smoky Mountains in the southern Appalachians are 330 expected to be in a quasi-topographic steady state, basin-wide denudation rates show a strong dispersion up to a factor of two in co mparison to the estimated uplift rate (ca. 0.03 mm.yr -1 , see Fig. 1). We use the data associated with 40 basins orig inally sampled by Matmon et al., (2003a, b) and for which denudation rates were re-calcu lated by Portenga and Bierman (2011).
Following our method, we calcu late the three basin -averaged aggressivity metrics , and associated with these 40 catchments ( Fig. 12; See also Fig. S3). The calculat ion of requires to define the elevation of the catchment outlets 12 and the / ratio (Eq. 4). As underlined by Forte and Whipple (2018), the choice of the "correct" outle t elevation is nontrivial in natural settings. We first consider a local base level given by the Tennessee River. To test the relevance of this choice, we also use a base level located at a fixed arbitrary elevation = 400 m. We assume the same / ratio value of 0.45 as used by Willett et al. (2014) for the Great Smo ky Mountains. For all calculated metrics, the majority (ca. 58% for and ca. 66% for ) of the basins are located in the expected quadrants (see Fig. 7). However, more attention must be 340 given to the results based on . For this metric, ca. 58% of the analy zed basins lie in the expected quadrant when we consider the Tennessee river as the local base level versus ca. 68% fo r = 400 m (Fig. 12b). Although the overall results are similar, we show that the choice of a different base level leads to significant variations in for individual basins. This highlights the main weakness of the metric, wh ich is highly sensitive to the choice of the proper base level .
Nevertheless, our results confirm the findings of Willett et al. (2014), suggesting that a significant part of the data varia nce 345 observed in the Matmon et al. (2003a, b ) can be exp lained by divide migrat ion (Fig. 12), raising this possible explanation for the variability of most natural data-sets. One may note that the Southern Appalachians exh ibit migrat ing kn ickpoints that can locally affect denudation rates (Gallen et al., 2011;Gallen et al., 2013). This last point can also explain part of the observed variability in this dataset but this specific impact is beyond the scope of the present study.

350
Based on both our simulations and this field dataset, we propose to favor the us e of and . A mong the tested metrics, appears the least sensitive to disequilibrium, excepted in active mountain belts with rock uplift U ≥ 1 mm.yr -1 .

Assessment of topographic disequilibrium
Topographic steady-state is a very convenient assumption and concept to deduce the uplift pattern in mountains ranges from denudation rates, and thus to obtain significant information on the geomet ry of active structures and on orogen dynamics 355 Avouac, 2001, Godard et al., 2014 ;Scherler et al., 2014;Le Rou x-Mallouf et al., 2015). However, this assumption is seldom verified at the scale of sampled watersheds.
On the basis of our modeling, we show that the competit ion between low -order basins has a significant impact on basin-wide denudation rates. The proposed approach provides a new tool to assess the potential deviat ion fro m topographic steady state based on aggressivity metrics and drainage area, which can both be inferred fro m a simple DEM : the closer to zero the 360 aggressivity metrics and the lower the standard deviation of cross -divide met rics, the more representative of uplift rate the measured denudation rates .

Improvement of sampling strategy
Basin-wide denudation rates obtained fro m CRN concentration measurements, suspended sediments or detrital thermochronology depend on many parameters including lithology, ice cover, rainfall, landslide activity or tectonic uplift 365 (Vance et al., 2003;Bierman and Nichols, 2004;Witt mann et al., 2007;Yan ites et al., 2009;Norton et al., 2010;Godard et al., 2012;Whipp and Elhers, 2019). Hence, to unravel the influence of tectonics from other processes, a specific sampling strategy is usually reco mmended: (1) to sample catchments with homogeneous lit hologies to limit the effect of spatial variations in the abundance of target minerals in bedrock formations; (2) to select catchments with no ice cover (past or present) because the input of glacier -derived sediments can significantly co mplicate the interpretation of CRN 370 concentrations; (3) to choose areas with spatially uniform rain fall d istribution; and (4) to consider watersheds where the relative contribution of landslides to long-term landscape evolution is lo w. Un fortunately, these different criteria imp ly selection of watersheds with variable sizes. The first three criteria favor the sampling of small catch ments, whereas the last one requires basins large enough to be less affected by landslides.
Our approach suggests the need to pre-assess targeted basins for their potential div ide mob ility before samp ling for CRN 375 concentration measurements . If the objective is to quantify the background uplift rate, one should sample basins that satisfy the conditions we previously described in the current section and also display an aggressivity close to zero and with the smallest associated standard deviation. Conversely, to quantify the specific denudation rate associated with the migration of drainage divides, small aggressor or victim basins should be favored.
Based on our simulat ions, a relationship between the maximu m of erosion variab ility (0.5 and 99.5 percentiles, respectively) 380 due to divide mobility [( − )⁄ ] and the catchment size can be derived (Fig. 13). Ou r results suggest a logarithm dependence between these two parameters, regardless of the assumed , , and : [( − )⁄ ] = 1 log ( ) + 2 for 1 km² < < 100 km² , 385 with 1 and 2 two parameters that depend on balance between erosion processes, uplift rate and state of evolution of the landscape.

Conclusions
Calculations fro m a Landscape Evolution Model assuming spatially uniform uplift, rock strength and rainfall confirm that the concept of topographic steady state is relevant at the scale of entire mou ntain belts, but represents an oversimplification 390 at the scale of individual watersheds. Our simu lations underline the role of divide mobility on deviations from equilibriu m, which can lead to significant differences between tectonic uplift rate and basin -wide denudation rates even if an overall topographic steady state is achieved at large scale.
To better assess these deviations, we propose new basin-averaged aggressivity metrics -, and -based on the approach of Willett et al. (2014) and Forte and Whipple (2018). They include mean cross -divide contrasts in channel , local 395 gradient and height. Fro m our calcu lations, is the most reliable aggressivity metric to assess local disequilibriu m, but is highly depend on the chosen base level, wh ich remains hard to constrain. Gilbert's metrics and are mo re suitable for relat ively h igh uplift rate (i.e. ≥ 1 mm.yr-1). A ltogether, our met rics reveal that deviation of denudation rates fro m uplift rate related to div ide migrat ions depends on both basin aggressivity and basin area. Th is last parameter has a key control on the dispersion in / , which can reach a factor o f two, regardless of the imposed uplift rate (here 0.01-1 mm.yr -1 ), 400 erodibility (here 5.10 -6 -5.10 -5 m (1-2m) .yr -1 ), diffusivity (here 10 -3 -10 -1 m 2 .yr -1 ) and hillslope grad ient (here 20°-40°). By comparing our results to CRN measurements fro m the Great Smo ky Mountains (Matmon, 2003a, b), we show that this approach can be used to improve field sampling strategies and provides a new tool to derive a minimal uncertainty in basinwide denudation rates due to topographic disequilibrium.
For the sake of simp licity our models involve spatially ho mogenous and time invariant parameters. Additional simulat ions 405 are now needed to test this approach in more co mplex settings, including spatial and temporal variability in climate and tectonic forcing or parameters like stream power equation exponents and .

Author contribution
R.C. and M.F. in itiated this study. T.S.S. performed the simu lations and topographic analyses. All authors contributed to the writing of the paper. 410  For channel gradient, divides migrate toward the drainages that present lower values. In our study, channel , local gradient and height are measured at the outlet (indicated by red circles) of basins for a reference area (basins bounded with thin black lines). Aggressivity metrics are then calculated for a given basin (represented in grey) by averaging along its perimeter the individual across divide differences in metrics between reference basins. The proportion of perimeter which is not shared by two reference basins is measured to give a confidence index of the calculated aggressivity.