Drainage reorganization induces deviations in width-area-slope scaling of valleys and channels

The width of valleys and channels affects the hydrology, ecology, and geomorphic functionality of drainage networks. Valley and channel widths are often estimated through a power-law scaling between width (W) and drainage area (A), and where lithologic variability or differential uplift rates dominate, width was suggested to scale with both 15 slope (S) and drainage area, through the relation W = kbA S. However, in fluvial systems that experience drainage reorganization, abrupt changes in drainage area distribution can result in widths that are disproportional to their drainage areas. Consequently, in such cases, the width-area-slope scaling is expected to deviate relative to drainages that did not experience reorganization. To explore the effect of reorganization on width-area-slope scaling, we studied 12 valley sections in the Negev desert, 20 Israel, categorized into undisturbed, beheaded, and reversed valleys. We found that the drainage area exponent, b, differs between valley categories, and that reversed valleys are characterized by a negative b exponent, indicating valley narrowing with increasing drainage area. A detailed study of a reversed valley reveals that unlike the negative b exponent that links drainage area to valley width, the relation between drainage area and channel width is best fitted with a positive b exponent. This difference indicates that the timescale of channel width adjustment to post25 reorganization drainage area distribution is faster than that of the valley width adjustment. We find that the difference in channel width across the divide causes a step change in unit stream power between the adjusted reserved channel and the unadjusted, beheaded channel. Gradients in width and unit stream power across the divide, lead to a widthfeedback that promotes ongoing divide migration and reorganization. The identified distinct width-area-slope scaling of reorganized valleys could assist in recognizing and constraining the 30 dynamics of landscapes influenced by drainage reorganization and likely has critical implications for the distribution of erosion rates in reorganized landscapes. https://doi.org/10.5194/esurf-2021-105 Preprint. Discussion started: 21 January 2022 c © Author(s) 2022. CC BY 4.0 License.


Introduction
The width of channels and their hosting valleys control river's dynamics and functionality with far-reaching implications across a wide range of disciplines, from flood and seismic hazards (e.g., Lóczy et al., 2009;Mashael Al, 35 2010;Morell et al., 2020;Sampson et al., 2015), to river ecosystems and habitats (e.g., Beeson et al., 2018;Brussock et al., 1985;May et al., 2013;Sweeney et al., 2004), and to hydrological modeling (e.g., Looper et al., 2012).Valley and channel width further play a central role in landscape evolution, as they respond to changes in the forcing that act on the landscape ֵ (Amos and Burbank, 2007;Fisher et al., 2013;Hancock and Anderson, 2002).The ratio between valley width, which integrates the channels, terraces, and floodplains, on the one hand, and morphometric valley 40 properties such as depth or fill thickness, on the other hand, are used to elucidate drainage evolution over geological timescales (Gibling, 2006;Schumm and Ethridge, 1994) and for inferring climate changes (Dury, 1964;Hancock and Anderson, 2002;Marcotte et al., 2021) and tectonic variations (Giaconia et al., 2012).The channel width is a key component in landscape evolution for its control on the shear stress, sediment transport capacity, and erosion rate (Whittaker et al., 2007b;Yanites et al., 2010).Particularly, many landscape evolution and hydrologic models 45 approximate the local erosion rate as a function of the channel stream power per channel width (i.e., unit stream power, Harbor, 1998;Lague, 2014;Magilligan et al., 2015;Turowski, 2018;Yanites, 2018).
The wide use of valley and channel width across various disciplines highlights the value of high-resolution width measurements, which could vary by several orders of magnitude within a single basin and across basins and landscapes (Schumm and Ethridge, 1994).However, producing high-resolution field-based width measurement is 50 challenging and time-consuming.In recent years, a growing body of work focused on developing tools for automatic width extraction based on remotely-sensed data (e.g., Fisher et al., 2013;Gilbert et al., 2016;Hilley et al., 2020;Monegaglia et al., 2018;Roux et al., 2014;Rowland et al., 2016).Although these tools enabled a significant advancement in river research and management, they commonly focus on specific types of river morphology, and require parameter calibrations, as well as human supervision (Fryirs et al., 2019;Golly and Turowski, 2017).Due to 55 these limitations, width of natural channels and valleys is more commonly estimated based on the widely recognized scaling relationships between valley and channel widths and fundamental basin properties such as discharge (or its proxy, drainage area), which could be relatively easily measured (e.g., Lavé and Avouac, 2001;Wobus et al., 2006).Furthermore, these scaling relationships are used frequently in landscape evolution models, where the width is parametrized based on the drainage area (e.g., Goren et al., 2014;Lague et al., 2014;Shobe et al., 2017;Yanites et al., https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.2013).Beyond steady conditions, several studies demonstrated that width can adjust dynamically in response to transient environmental conditions (Snyder and Kammer, 2008;Whipple et al., 2013;Yanites, 2018), and that in such landscapes, a more complex scaling involving width, area, and slope might be more applicable.

Width scaling relations 65
The common approach for width estimation relies on the seminal work of Leopold and Maddock (1953), who used empirical data to establish a power-law relation between the channel width, W [m] and discharge, Q [m 3 /sec].Based on the well-established correlation between discharge and drainage area, A [km], (e.g., Dunne and Leopold, 1978), the scaling between channel width and drainage area is often expressed as  =   *    (1) 70 Leopold and Maddock's relation, (Eq.(1)), was initially established for alluvial rivers, where  was found to be ~0.5.
Several studies suggested that the relation between the valley width and drainage area follows a similar power-law scaling (Beeson et al., 2018;Brocard and van der Beek, 2006;Langston and Temme, 2019;Langston and Tucker, 2018;May et al., 2013;Phillips, 2011;Schanz and Montgomery, 2016;Snyder et al., 2003;Tomkin et al., 2003).In these cases, the reported range of the exponent d was significantly wider, between 0.1 (Langston and Temme, 80 2019) and 1.18 (Beeson et al., 2018).Here too, the range was attributed to differences in the properties of the valleybounding rocks (Keen-Zebert et al., 2017;Langston and Temme, 2019;Schanz and Montgomery, 2016), or, in some high relief landscapes, to the spatial distribution of deep-seated landslides (Beeson et al., 2018;May et al., 2013).
While the applicability of the simple width-area power low scaling (Eq.( 1)) was demonstrated in many settings (Montgomery and Gran, 2001;Whitbread et al., 2015), field observations show that it is not applicable across 85 all landscapes.Particularly, the scaling was demonstrated to fail when channels respond to transient conditions, such https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.
as crossing areas of enhanced uplift rate (e.g., due to local faulting or folding) (Allen et al., 2013;Finnegan et al., 2005;Lavé and Avouac, 2001;Montgomery, 2004;Whittaker et al., 2007a;Yanites et al., 2010), inducing dynamic changes to channels local slope and width (Amos and Burbank, 2007;Kirby and Ouimet, 2011;Whipple et al., 2000;Whittaker et al., 2007aWhittaker et al., , 2007b;;Yanites, 2018).Consequently, substantial research focused on understanding dynamic 90 adjustments of bedrock channel width in response to uplift rate changes (e.g., Turowski, 2018;Wobus et al., 2006;Yanites, 2018).Finnegan et al. (2005) developed a model for the case of a channel that crosses terrains with variable rock uplift rates.Adopting Manning's equation (Manning et al., 1890) and assuming a constant bankfull width-to-depth ratio along the channel, Finnegan's model predicted that the channel width not only depends on drainage area, as in 95 Eq. ( 1), but also on the channel slope, S [m/m]: Based on theoretical calculations, the exponents b and c in Finnegan's model equal 0.38 and -0.19, respectively.These values were also validated in the field (Finnegan et al., 2005) and were later applied by following studies that examined width-area-slope scaling in field settings with alternating lithologies and variable uplift rates (Kirby and Ouimet, 2011;100 Spotila et al., 2015).Whittaker et al. (2007a) andAttal et al. (2008) studied settings of transient channels adjusting to changing tectonic forcing and suggested a more negative value of -0.44 for c, reflecting the dominant role of the slope in these settings.

Drainage reorganization and width scaling 105
Another form of variability and transiency in rivers is associated with drainage reorganization, which is widely recognized as an important process affecting the evolution of fluvial systems (e.g., Bishop 1995, Willett et al. 2014, Fan et al., 2018;Prince et al., 2011).Reorganization occurs when drainage divides shift through time (Bishop, 1995;Davis, 1889), change basin geometry, and consequently induce changes in the discharge and drainage area distribution along the channels (e.g., Menier et al., 2017;Pechlivanidou et al., 2019).The scaling of Eq. ( 1) and (2) predicts that 110 the addition or reduction of drainage area will result in widening or narrowing, respectively, of the channel and valley.persist over long time scales (e.g.,Jones, 2018;Snyder and Kammer, 2008).Valley width, particularly, is likely to adjust at slower rates than channel width because it represents the channel location integrated over long periods, commonly in the scale of >1ky (Hancock and Anderson, 2002;Langston and Tucker, 2018;Schumm and Ethridge, 115 1994;Tomkin et al., 2003).Nonetheless, despite the intrinsic links between drainage reorganization and drainage area change, to the best of our knowledge, the effect of drainage reorganization on fluvial width scaling has not yet been evaluated.Recognizing and accounting for these scaling deviations could be highly consequential for predicting the hydrologic and geomorphic functionality of reorganized drainages.
We hypothesize that the scaling between drainage area and width could express the delayed response of 120 width adjustment to changes in drainage area following the reorganization of drainage networks.Accordingly, the coefficient and exponent values that relate drainage area with valley and channel width in reorganized drainages could meaningfully deviate from drainages that did not experience reorganization.
To test this hypothesis and to evaluate the effect of reorganization on the width-area-slope scaling, we analyzed and compared the geometry of reorganized and undisturbed drainages in the southern Negev desert, Israel, 125 where drainage reorganization is well established by field observations (Avni et al., 2000;Ginat et al., 2000Ginat et al., , 2002;;Harel et al., 2019).In the current analysis, we aim to (1) explore if and how the scaling between valley width, drainage area, and slope varies between reorganized and non-reorganized drainages and between drainages that experienced different modes of reorganization; (2) compare the adjustment of channel width relative to valley width following reorganization; and (3) examine implications of the width-area-slope scaling for reorganized drainages on landscape 130 evolution.

Geologic and Geomorphic setting
We explore the scaling between the valley and channel width, drainage area, and slope along ephemeral valleys that 135 incise into the southeastern Negev Highlands, Israel.The eastern boundary of the study area is the ~400-600m cliff of the Arava Valley (Fig. 1).The Arava Valley is a rift-like structure stretching between the Dead Sea and the Gulf of Aqaba and is the southern segment of the transtensional Dead Sea plate boundary, between the Arabia plate and Sinai sub-plate (Garfunkel, 1981;Garfunkel et al., 2014, Fig 1A).
The main drainage divide in the study area mostly separates east-flowing basins that drain across the cliff 140 toward the Arava Valley, from west-flowing, low relief channels that flow on the Negev highlands (Fig. 1b-c).The lithology exposed along the highland valleys consists primarily of Cretaceous limestone and dolomite strata (Ginat, 1991).The climate is hyper-arid with an average annual precipitation of ~<50 mm (Bitan and Rubin, 1991), typically generating one to few flashflood events per year.These climatic conditions had generally persisted through most of the Pleistocene (Amit et al., 2006(Amit et al., , 2011)), except for relatively short episodes of transition to wetter conditions (Ginat 145 et al., 2018;Vaks et al., 2013).
The eastern Negev desert has been experiencing ongoing fluvial reorganization since the late Miocene (Avni et al., 2000(Avni et al., , 2012)).Before the development of the Arava Valley, rivers that originated in the Jordanian highlands, east of the Arava Valley, flowed westward, crossing the plate boundary along the Negev highlands towards the Mediterranean (Garfunkel and Horowitz, 1966;Zilberman and Calvo, 2013).Since the Miocene, tectonic activity 150 along the Dead Sea plate boundary formed the Arava Valley that gradually became a prominent base level.
Consequently, during the Plio-Pleistocene, several large-scale capture events had redirected major drainage systems in the Negev toward the central Arava valley (Avni et al., 2000;Ginat et al., 2000Ginat et al., , 2002;;Guralnik et al., 2010).Field observations and a regional χ analysis (a morphometric parameter used to approximate the stability of drainage divides, (Willett et al., 2014)) suggest that the regional divide between the Arava Valley and the Negev highlands is still 155 actively migrating westward (Harel et al., 2019).

Categories and characteristics of valleys in the study area
To explore the effects of drainage reorganization on the scaling relations between width, drainage area, and slope, we analyzed 12 valley sections associated with different drainage reorganization categories.All sections are located 170 adjacent to the Negev-Arava drainage divide, resulting in relatively small drainage areas of 0.2-14.2km 2 .The valleys are incised into bedrock, generating a relief of several tens of meters between the valley bottoms and the highlands'  1b-c), the morphology of the valley section (Fig. 2), and additional supporting field observations: 1. Undisturbed valleys (n=4) are westward flowing sections that do not show evidence of local drainage 175 reorganization, suggesting that they were not directly affected by the formation of the Arava base level.In these valleys, the low-order (sensu Strahler) segments are characterized by a V-shaped morphology (Figs.2A and 2C) and relatively rapid downstream widening.Higher-order valleys widen downstream slower and typically have a trapezoid cross-section with a sediment-filled flat valley bed, whose relief is tens of cm (Figs.2b-d).Here, the definition of valley and channel could be ambiguous because the low local relief results in braided and dynamic flow pathways, 180 and because the entire valley is flooded during large rainstorm events (Fig. 2d).
2. Beheaded valleys (n=3) are west-flowing sections whose headwaters were beheaded.Beheading is recognized by a windgap, i.e., a valley confined drainage divide that occurs along the cliff, or that is shared with a reversed valley (described below), indicating a paleo-valley that was formed by a channel that used to drain a larger area.Close to the windgap, the beheaded valleys are characterized by a U-shaped cross-section (e.g., Figs 1b, 2e-f, 185 3a, and 6a), likely controlled by the colluvial aprons concave profile.West and downstream from the windgap, these valleys become indistinguishable from the undisturbed valleys with the trapezoid-shaped cross-section.Valley's beheading is associated with either the receding cliff and its coinciding divide, or with localized divide migration within the valley as part of a reversal process on the opposing side of the windgap (e.g., Figs. 2e and 3a, Bishop, 1995;Harel et al., 2019, Shelef andGoren, 2021). 190 3. Reversed valleys (n=5) host east-flowing channels that had reversed their drainage direction (Bishop, 1995) from west to east.These valleys are bounded between a windgap and a knickpoint generated where the channel flows across the cliff (e.g., Figs 2e, 3a and 6a).Harel et al. ( 2019) identified these sections as reversed drainages based on the presence of barbed tributaries and west grading terraces that record the antecedent valley gradient, opposite to the present-day channel's drainage direction.The reversed valley sections share windgaps with beheaded valleys, 195 indicating that they were part of the antecedent west-flowing drainage (e.g., Fan et al., 2018).In the reversed valleys, the valley morphology near the windgap is typically U-shaped, and downstream the channel incises into the alluvialcolluvial valley fill, creating cut terraces and forming a V or box-shaped channel cross-section within the broader   c-c', d-d,' and e-e

Methods 215
We study the effect of reorganization on valley width scaling by exploring the relations between valley width, drainage area, and slope following Eq.( 2).While Eq. ( 2) was originally developed to estimate channel width, we adopt it here for valley width due to its generality.In Sect.5.2, we devote a detailed discussion to the slope effect on valley width, which distinguishes the predictions of Eqs. ( 1) and (2).To compute the coefficient kb, and exponents b and c (Eq. ( 2)) that best fit the valley sections' geometry in 220 the study area, we applied a multivariate regression over the valley width, drainage area, and slope along the studied sections.The topography data were derived from TanDEM-X DEM (Wessel, 2016) with 0.4 arcsec resolution (~11.6m in the field area).
Unlike the upstream drainage area and channel slope, which are derived by simple morphometric analysis over the DEM, defining and extracting the valley width based on DEM is not straightforward and requires a tailored 225 procedure (Clubb et al., 2017;Golly and Turowski, 2017;Hilley et al., 2020;Roux et al., 2014;Rowland et al., 2016;Sechu et al., 2021).Particularly, the location and orientation of valley width measurements require caution because the width is not necessarily well-defined in proximity to side tributaries and valley bends (e.g., Beeson et al., 2018).
To overcome these challenges, we developed a semi-automatic approach for optimal measurements of valley width.230

Valley width measurements
Valley width was measured by applying two main operations.First, a polygon representing the valley bottom of drainage basins is extracted, and second, valley width is measured over the valley bottom polygon at optimal points (Fig. 3a).The first step is achieved by applying the ArcGIS plugin VBET-'Valley bottom extractor tool' (Gilbert et al., 2016).In VBET, valley boundaries are identified based on user-defined slope thresholds, representing the 235 transition from the valley bottom to the hillslope.VBET parameters used for the current analysis are described in Table S1 in the Supplement.Importantly, these parameters were fitted to each basin separately by an iterative process of visually comparing the valley bottom polygons against 0.5m/pixel aerial orthophotos and fine-tuning the parameters to achieve the best visual fit.In six basins, this procedure was not sufficient to achieve a satisfying fit between the VBET polygon and the orthophoto, mainly due to local DEM inaccuracies.In these cases, the polygons were manually 240 edited to correct mismatches.Manual editing was based on the orthophotos, available topographic data, and field observations.Except in one case where the area difference between the original and edited polygons was 10%, the differences in the other five basins were <5% (Table S1 in the Supplement).
Width measurements over the VBET polygons were achieved by applying an ArcGIS-based algorithm that identifies optimal locations for measuring valley width.The algorithm identifies points along the valley talweg that 245 are sufficiently far from bends and confluences, such that they don't affect the valley width.The algorithm operates

Drainage area and slope extraction, and multivariate regression
Along the intersection pixels between width transects and the valley talweg, drainage area and slope were extracted based on the DEM.The drainage area was extracted from a flow accumulation raster, computed using a D8 flow routing algorithm (O'Callaghan and Mark, 1984).The slope was calculated using a 5-pixel running window along the 255 channel profile, where the slope that was calculated between the terminal pixels of the window was assigned to the central pixel.
Importantly, the extracted slope is that of the channel.In the undisturbed and the beheaded valley sections, field observations indicate that the channel and the valley slope are indistinguishable (e.g., Fig. 2b,2d and 2e).This is not the case for reversed valley sections, where the channels are decoupled from the valleys.The reversed channels 260 are incised into the bottom of paleo-valleys, whose grading may be opposite to that of the active channel (Harel et al., 2019).
Direct measurements of the reversed valley's slope are unlikely to produce reliable results because the preservation of the antecedent valley bed is, in most cases, poor and patchy.In Sect.3.3 and 5.2, we present an independent analysis dedicated to the reversed valleys and discuss the implications of using the channels' slope in the 265 valley analysis.
After extracting the valley width, drainage area, and slope at the intersection pixels, the best-fit values of kb, b, and c in Eq. (2) were obtained.The best-fit values were calculated by using a multivariate least-squares linear regression over W [m], A [km], and S [m/m] (following Attal et al., 2008;Spotila et al., 2015) that were logtransformed, (Fig. 3b-d).The data points were not binned.270

Detailed analysis of channel and valley width
For the reversed category, where the valley and the channel are decoupled, we examined how fitting the valley width compared to the channel width affects the predictors kb, b, and c in Eq. ( 2).Valley 12 (Fig. 1c and Table 1) is a thoroughly surveyed site investigated in a previous study (Harel et al., 2019) that was chosen for this analysis.The 285 channel analysis is based on a 15 cm/pixel DEM and a 3 cm/pixel orthophoto generated using a structure from motion (SfM) algorithm over drone-acquired aerial photos (80% overlap).Here, the channel bottom polygon was delineated manually based on the high-resolution DEM and orthophoto.Then, the width measurement algorithm described in Sect.3.1 was applied over the channel polygon.The multivariate regressions for the channel and the valley used the slope and drainage area based on the high-resolution DEM, following the procedure described in Sect.3.2.290

Validations and errors in the measurements and model
The main potential sources of width measurement errors originate from the DEM horizontal resolution, R (~11.6 m 2 /pixel), and the relative vertical accuracy (~2m, Wessel, 2016).To incorporate the uncertainty stemming from the horizontal resolution of the DEM, we assigned each width measurement with a constant error, evaluated as √2 * .295 To independently explore the effect of the vertical accuracy on the final width measurements, valley transects of seven selected sites were measured with a differential GPS (DGPS).The valley bottom was extracted from the transects by applying the same slope criteria as the VBET tool for that basin.The correspondence between the DEM and the DGPS-based width measurements is shown in Fig. 4: The differences between the DEM-based width measurements relative to the DGPS-based measurements range between 1.7-19.8m, and are scale-independent.300 Overall, the RMSE is 13.08 m, smaller than the error of √2 *  evaluated from the horizontal DEM resolution.The width measurements are available in Table S2 in the Supplement.

Width, drainage area, and slope scaling along valley sections
The best-fit coefficients and exponents of individual valley sections, their 95% confidence intervals, and the adjusted 310 R 2 for the regression are presented in Table 1 and Fig. 5. P-values of the predictors and of the multivariate regressions are provided in Table S3 in the Supplement.The multivariate regression results reveal a unique range of the drainage area exponents values, b, for each predefined valley category (Fig. 5b).The undisturbed valleys are characterized by the highest exponents, ranging from 320 0.25 to 0.55, whereas the b exponents of the beheaded valleys are lower, 0.17-0.23.Uniquely, the reversed valleys have negative b exponents, ranging from -0.24 to -0.71, indicating that the valley narrows with increasing drainage area in this category.In the undisturbed and beheaded valleys, the values of the slope exponent c are smaller by an order of 330 magnitude relative to the values of the area exponent, b, and are statistically insignificant (with one exception in valley 5), at the 5 % significance level (Table 1).On the other hand, in all the reversed valleys except for valley 10, the values of the c exponent are negative, have the same order of magnitude as the area exponent values, and are statistically significant (P-value <0.05).
The model performance was evaluated through the adjusted R 2 , which was higher than 0.6 in most valleys, 335 and overall >0.38 (Table 1, Fig. 5).The model P-value for all valleys was <0.05 (Table S3 in the Supplement).

Comparing valley width and channel width in a reversed drainage
We explored the effect of drainage area and slope on the width of the channel vs. the width of the valley using a dronederived high-resolution DEM of reversed valley 12 (Table 1, Fig. 6).Unlike the undisturbed and beheaded valley categories, the channels in the reversed settings are incised into the valley fill (Harel et al., 2019).The channel in valley 12 initiates east of the windgap and incises into the erodible valley fill, where it merges with short side 350 tributaries that drain the valley bottom.Further downstream, it merges with a barbed tributary that joins the valley from the north (Fig. 6a), and after ~160 m, bedrock is exposed at the base and the north bank of the channel.At this point, the reversed channel is incised ~20 m below the surfaces of the antecedent valley bottom.The channel then Field observations show that while the reversed valley narrows downstream (i.e., eastward), the channel width increases in the downstream direction (Fig. 6a).A multivariate regression over the channel data reveals a drastically different dependency between the channel's width, drainage area, and slope compared to the valley (Fig 6b).For the channel, the least-square multivariate regression (R 2 =0.68) yields a kb coefficient of 10.64±2.57(10 -6 m 1- 2b ), a positive and high b exponent of 0.58±0.17,and a relatively small, negative, statistically insignificant c exponent 360 of -0.07±0.18.In contrast, the computed values for the valley are 3.04±3.29(10 -6 m 1-2b ) for kb, a negative b exponent of -0.43±0.19,and a statistically significant c, -0.82±0.44.

Drainage reorganization affects width-area scaling
The multivariate regression results reveal that undisturbed, beheaded, and reversed valleys could be characterized by a distinct range of the drainage area exponent, b, in the width-area-slope scaling.Therefore, reorganization could meaningfully affect the width-area scaling of valleys.In our study area, the b exponent values of the undisturbed 375 valleys range between 0.25 and 0.55, consistent with previously published exponent values derived based on widtharea or width-area-slope scaling of channels and valleys (Kirby and Ouimet, 2016;Spotila et al., 2015;Whittaker et al., 2007a, Beeson et al., 2018;Brocard and van der Beek, 2006;Langston and Temme, 2019;May et al., 2013;Schanz and Montgomery, 2016;Snyder et al., 2003;Tomkin et al., 2003).The beheaded valleys in the study area are characterized by b exponent values of 0.17-0.23.While this range partly overlaps with previously published valley 380 width-area scaling (Beeson et al., 2018;Langston and Temme, 2019;Schanz and Montgomery, 2016) )), is consistent with the process of beheading.During beheading, a valley loses its narrowest headwater sections, and consequently, the beheaded valley is wider at smaller drainage areas compared to undisturbed valleys (e.g., Fig. 385   3c).Additionally, the drainage loss reduces the sediment transport capacity near the divide and may lead to aggradation that widens the valley.Thus, the log(W [m]) vs. log(A [km 2 ]) relations for beheaded valleys are characterized by a decreased slope compared to undisturbed valleys (e.g., Fig. 3b-c).
The process described above is expected to also affect the kb coefficient values of the beheaded valleys.
Whereas the median kb value is indeed somewhat higher for beheaded valleys compared to undisturbed ones (Fig. 5  390 and Table 1), the difference in kb is small.This is because the value of the kb coefficient reflects the valley width at a drainage area of 1 km 2 (given that the c exponent is approximately zero).In the study area, a 1 km 2 drainage area is reached only after the beheaded section is joined by several undisturbed tributaries, which likely obscure the beheading influence, making the kb coefficient of the undisturbed and beheaded valleys overall similar.
The negative value of the b exponent of the reversed valleys (between -0.24 and -0.71) represents downstream 395 valley narrowing, supporting the inferred reversal of these valley sections (Harel et al., 2019).Furthermore, the similar The ranges of the b exponents differ between the undisturbed and reorganized valleys, in agreement with our hypothesis that drainage reorganization modifies the scaling between valley width and drainage area.Based on these results, we suggest that scaling differences between adjacent drainages could serve as supportive evidence for drainage reorganization.Furthermore, the value of the b exponent could relate to specific categories of reorganization.The value of the b exponent for beheaded valleys is expected to be lower relative to that of undisturbed valleys.Negative 405 b exponents could indicate reversed drainages, where preserved antecedent valleys narrow with increased drainage area.Importantly, however, the width-area scaling could vary across climates (Hancock and Anderson, 2002) and lithologies (Brocard and van der Beek, 2006;Langston and Temme, 2019;Schanz and Montgomery, 2016), and is affected by structures (Keen-Zebert et al., 2017) and landslides (Beeson et al., 2018;May et al., 2013).Thus, invoking b exponents to support reorganization requires comparing the suspected reorganized valley sections to undisturbed 410 sections with similar environmental and structural conditions.Additionally, accurately identifying the extent of the reorganized valley section is crucial because sections that include more than a single category might show a noisy signal, if at all.

Slope-valley width relation in the study area
In steady-state drainage networks with uniform lithology, climate, and uplift rates, the slope, S, and the drainage area, 415 A, scale through the power-law S  A -θ (Flint, 1974).In such cases, slope can be substituted for area or vice versa, such that the width in Eq. ( 2) becomes a function of drainage area only (i.e., Eq. (2) reduces to the form of Eq. ( 1)) or of slope only.In contrast, when A and S do not covary, each of them can independently influence the width scaling (Finnegan et al., 2005;Whittaker et al., 2007b).In those cases, the width predicted by Eq. ( 2) is expected to be more accurate than that predicted by Eq. (1).420 The multivariate regression results reveal three lines of evidence indicating that the valley width is independent of the slope in the undisturbed and beheaded valley sections in the study area: (i) the slope exponent , c, is smaller by an order of magnitude with respect to the b exponent (Figure 5 statistically insignificant (with the exception of valley 5), and (iii) Least-square regressions following the slopeindependent Eq. ( 1) yield similar width predictors and R 2 values as the regression based on Eq. (2) (Fig. S6 in the 425 Supplement).However, a slope-area analysis (Fig. S7a in the Supplement) shows a low covariance between log(S [m]) and log (A [m]), (R 2 <0.24) in the studied valley sections, indicating that the slope-independent width predictions do not arise from such covariance.Alternatively, we suggest that the statistical insignificance of the c exponent in the undisturbed and beheaded valley sections stems from the scatter in the slope data (Fig. S7b in Supplement), which reduces the weight of the slope relative to drainage area in the regression.The scatter is particularly distinct due to the 430 exceptionally low slope values in these valleys, which could be <10 -4 [m/m].
As opposed to the undisturbed and beheaded valleys, in the reversed valley sections (except for valley 10), the inferred slope exponents, c, are statistically significant, negative, and with the same order of magnitude as the area exponent b.Furthermore, for these sections, the multivariate regressions based on Eq. ( 2) produced higher R 2 values than a slope-independent regression based on Eq. ( 1), indicating that the slope plays a role in predicting the reversed 435 valley width.The slope significance could be partly explained by the somewhat reduced scatter in the slope data along the reversed category relative to the other categories, likely due to the overall higher slopes along the reversed sections (Fig. S7c in the Supplement).However, from a geomorphological perspective, the measured channel slope is expected to be decoupled from the valley morphology, because in this category, the channel and the antecedent valley grade in opposite directions (Harel et al., 2019).We suggest that the slope significance in these valleys could be explained by 440 inspecting the width-slope correlation close to the downstream edge of the reversed sections, just above the knickpoint, where the valley narrows and channel steepens (e.g., valley transects in Fig 2e-f, and Fig. S7c in the Supplement).
Narrowing and steepening could reflect channel and valley response to flow acceleration close to the upper lip of the knickpoint (Haviv et al., 2010), forming a juvenile narrow valley on the expense of the paleo valley (Fig. 2f).The incision above the knickpoint that contributes to the correlation between channel slope and valley width at the 445 downstream edge of most reversed sections likely reflects a transient response to reorganization and the onset of valley width adjustment to the new drainage direction and drainage area distribution.
Valley 10 is an interesting exception in this context.Here, despite a ~80 m high knickpoint at the edge of the reversed section, slope increase and valley narrowing above the knickpoint are absent.The lack of incision above the knickpoint in valley 10 could imply a recent episode of knickpoint migration to its current location.Accordingly, in 450 this case, the multivariate regression infers a small and statistically insignificant slope exponent (Table 1).

Differing timescales of the valley and channel width adjustment
Channel width is recognized as a sensitive parameter that adjusts dynamically and relatively rapidly to changes in tectonic (Amos and Burbank, 2007;Attal et al., 2008;Morell et al., 2020;Yanites, 2018) and climatic (DiBiase and Whipple, 2011) conditions.Over longer timescales, the valley width, which depends on the time-integrated location 455 of the channel, is expected to adjust to the same changes in the boundary conditions (Hancock and Anderson, 2002;Langston and Temme, 2019;Tomkin et al., 2003).Drainage reorganization changes the drainage area distribution and, like tectonic and climatic variations, induces dynamic width changes.We suggest that the comparison between the valley and channel width in reversed valley 12 (Fig. 6) demonstrates a temporal snapshot where the channel width is adjusted to the new drainage area distribution inflicted by the drainage reversal, whereas the valley width is not yet 460 adjusted to the same forcing (Fig. 6b), in accordance with the longer timescale expected for valley adjustment.

Implications to landscape evolution
Delayed valley versus channel adjustment in response to reorganization (Fig. 6) and the diverging response of valleys of different reorganization categories (Fig. 5) can have critical implications for landscape evolution.We explore one such implication by inspecting the influence of width adjustment on a proxy for the unit stream power (ω= gQS/W 465 [Watt/m], : density, g: gravitational acceleration, Q: discharge), which is commonly used for evaluating fluvial erosion rate (e.g., Harbor, 1998;Magilligan et al., 2015).Given that g can be treated as a constant, and that Q is typically proportional to drainage area, the unit stream power is expected to be proportional to psp =AS/W [L] (Whittaker et al. (2007a)).We therefore use psp to compare the unit stream power between a reversed and beheaded valley sections that share the same windgap and flow in opposite directions (Fig. 7a).In the reversed section, psp is 470 comparably high because of the narrow width and the increased slope of the actively incising channel (Fig. 6a).In the windgap (Fig. 7b, black x symbols), where the psp values in the reversed section are larger by an order of magnitude 475 relative to the beheaded section.
The differences in psp across the windgap are likely associated with the formation of the reversed section, and feeds back with the reversal process.Harel et al. (2019) proposed that in this study area, channel reversal initiates and extends by gradual windgap migration along an antecedent valley.Linking the windgap migration to the dynamic valley and channel adjustments, we propose that in the extending reversed section, the migrating windgap increases 480 the drainage area above the knickpoint.This, in turn, is expected to increase the erosion rates of the reversed channel, while on the beheaded side, at the other side of the windgap, the valley is losing drainage area, and experiences reduced erosion rate and even deposition, that overall reduce the rate of channel adjustment to reorganization.These differences in stream power across the divide can cause feedback in which divide migration promotes rapid channel width adjustment on the reversed side (whose drainage area is increasing), while on the beheaded side (whose drainage 485 area is decreasing), adjustment is detained.This, in turn, generates a step-change in unit stream power across the windgap that could further push the windgap in the same direction.This process may be perturbed by capture of side tributaries by the reversed section, that cause abrupt changes in drainage area and erosion rate across the windgap (Shelef and Goren, 2021).
The step-change in psp values across the windgap emphasizes the importance of accurate channel width 490 estimates when exploring the evolution of landscapes undergoing drainage reorganization.More specifically, when the width is approximated based on Eq. (2) without accounting for the influence of reorganization on the channel width (i.e., by using the same width scaling across the windgap), the aforementioned step-change in psp is not recognized (Fig. 7b, blue dots).Consequently, a computed difference in erosion rates across the divide would be lower and the divide will be estimated as being more stable than it actually is (Fig. 7b).495

Conclusions
Our analysis of undisturbed and reorganized valley sections in the Negev desert shows that the scaling between valley width, drainage area and slope is modified by reorganization.The analysis further reveals that each reorganization category is associated with a distinct range of b exponent values, relating the valley width to the drainage area.In the undisturbed valleys, the range of b exponent values is consistent with values reported in previous studies (Attal et al., 515 2008;Kirby and Ouimet, 2011;Whittaker et al., 2007a).The b exponent values of the beheaded valleys are positive and smaller than in undisturbed valleys.In the reversed valleys, the b exponent values are negative, reflecting valley narrowing with increasing drainage area.We suggest that this deviation in scaling could be exploited, with caution, in future studies that aim to identify and categorize drainage reorganization.
Differences in width-area-slope scaling also occur between a channel and its hosting valley.In the reversed 520 valley section that we analyzed in detail, we found that the channel width is best fitted by a positive b exponent, reflecting the faster adjustment of channel width to the post-reorganization drainage area distribution relative to valley width (with a negative b exponent).The contrasting timescales of channel and valley width adjustment are consequential for landscape evolution.This is illustrated in a case study of reversed and beheaded valleys that diverge from a common windgap.The narrow width of the active channel in the reversed section compared to the wide channel 525 (that occupies the entire valley) of the beheaded section, results in a step-change in unit stream power across the windgap.Assuming that erosion rate correlates with unit stream power, this could lead to a divide migration feedback.
Erosion rate gradients across the windgap promote windgap migration toward the beheaded valley, which has a smaller unit stream power due to its wider channel.Windgap migration promotes rapid adjustment of the channel width on the extending reversed side, while on the beheaded side, adjustment is delayed, sustaining the gradient in unit stream 530 power.This feedback suggests that the differing response of channel and valley width in different reorganization categories could sustain ongoing divide migration, and may add to the slope and area feedbacks that were previously attributed to the process of divide migration (Shelef and Goren, 2021;Whipple et al., 2016;Willett et al., 2014).This https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.

Figure 1
Figure 1 (a) Regional map of the study area (black rectangle) near the Sinai-Arabia plate boundary (dashed white line), based on © Google earth aerial photo (2020).(b) A simplified sketch of valley categorization in the study area: Undisturbed valleys (Green, 'U' tag) are valleys that do not intersect with the cliff and are not affected by drainage reorganization.

Figure 2 :
Figure 2: Field photos and valley transects of valley sections in the study area.(a,b) Upstream and downstream segments of undisturbed valleys ((a) and (b) respectively).The drainage area in panel (a) is 0.08 km 2, and in panel (b) is 1.85 km 2 .Blue and red lines (a-a' and b-b,' respectively) mark the cross-section profiles shown in panel (c).(c) Transects of a-a' and ' emphasizing the difference between the U-shape transect near the windgap (e-e'), the V-shaped channel profile incised into the U-shaped valley terraces (d-d'), and the V-shaped transect above the knickpoint (c-c').
https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.iteratively, where each iteration identifies optimal width measurement points in valleys of a different order.The final output of the algorithm is a set of pixels located along the intersections between the valley transects, whose length represent the valley width, and the valley talweg (flowlines) (e.g., Fig 3a).The algorithm is described in detail in Section S1 and Figs.S1-S5 in the Supplement.250 https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.

Figure 4 :
Figure 4: Validation of the width measurements extracted from TanDEM-X against width measurements based on differential GPS valley transects.The maximal difference between the measurements is less than 20 m.
Unlike the b exponent values, the kb coefficient values and the slope exponent values, c, are non-unique for the categories (Fig 5a and 5c).The values of the kb coefficient, which represents the valley width at A=1 [km 2 ], ranges https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.from 56.26 to 108.73 (10 -6 m 1-2b ) in the undisturbed valleys, which overlaps with the range of the beheaded valleys, 127.74 -213.51 (10 -6 m 1-2b ) (Fig 5a).The kb coefficient values along the reversed valleys cluster into two groups, three valleys with low values of 3.04-13.72(10 -6 m 1-2b ), and two with higher values of 97.63 and 125.32, showing a similar range to the undisturbed and beheaded valleys.

Fig 6 :
Fig 6: Variations in valley and channel widths along a reversed section (valley 12 in Table 1 and Fig. 1c).(a) Valley (blue) and channel (orange) polygons and width transects sketched over a 0.5 m resolution orthophoto.(b) Width-area-slope scaling of valley width (blue), which narrows with drainage area, contrasting with the width of the reversed channel (red) that increases downstream.This difference is expressed by the drainage area exponent, b, which is positive for the channel and negative for the valley.
, it does not overlap with the b exponents values of the undisturbed valleys in our study area.The lower value of the b exponent of the beheaded valleys, reflecting a smaller change of width (log(W [m])) per unit change in drainage area (log(A [km 2 ] https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.magnitude of the b exponent values (i.e., regardless of sign) for the reversed and undisturbed valley sections is consistent with the view that (1) the geometry of the undisturbed valleys is similar to the antecedent valleys whose flow direction was later reversed (e.g., Figs 2e, 3d and 6); and (2) valley width did not drastically change following the reversal process.400 ), (ii) the values of the exponent c are https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.
contrast, psp is comparably low in the beheaded section, where anastomosing, low relief channels occupy almost the entire valley bed (Figs.2a and 2d and 7a), so that the valley width represents the effective channel width.The difference in psp between the opposing valley sections is expressed as a distinct step-change in the psp values across https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.

Figure 7 :
Figure 7: A proxy for unit stream power (psp=AS/W) along a profile from a reversed valley (# 12) to a beheaded valley (# 6) across a windgap.(a) An orthophoto of the reversed and beheaded valley sections that share a common windgap.The black line marks the profile route that follows the main channels and crosses the divide.The dashed white line marks the divide, and yellow lines mark measured width transects.In the reversed valley, east of the windgap, the active drainage is confined width feedback could be easily overlooked if channel width is parameterized based on standard scaling relations, without accounting for the deviated scaling of reorganized systems.535 Insights from this study point at several venues for future research, for example: what are the timescales over which the deviation in scaling persists and how do they vary with climate and lithology?Could the values of the b https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.exponent quantify the state of adjustment of channels and valleys?And what is the relation between the divide migration dynamics and rates to the valley and channel width adjustment?540 https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.https://doi.org/10.5194/esurf-2021-105Preprint.Discussion started: 21 January 2022 c Author(s) 2022.CC BY 4.0 License.

Table 1 :
Best-fit regressions for individual valley sections.315 a Predictor P-value >0.05.