Criteria and Tools for Determining Drainage Divide Stability 1 2

Abstract Watersheds are the fundamental organizing units in landscapes and thus the controls on drainage divide location and mobility are an essential facet of landscape evolution. Additionally, many common topographic analyses fundamentally assume that river network topology and divide locations are largely static, allowing channel profile form to be interpreted in terms of spatio-temporal patterns of rock uplift rate relative to base level, climate, or rock properties. Recently however, it has been suggested that drainage divides are more mobile than previously thought and that divide mobility, and resulting changes in drainage area, could potentially confound interpretations of river profiles. Ultimately, reliable metrics are needed to diagnose the mobility of divides as part of routine landscape analyses. One such recently proposed metric is cross-divide contrasts in χ, a proxy for steady-state channel elevation, but cross-divide contrasts in a number of topographic metrics show promise. Here we use a series of landscape evolution simulations in which we induce divide mobility under different conditions to test the utility of a suite of topographic metrics of divide mobility and for comparison with natural examples in the eastern Greater Caucasus Mountains, the Kars Volcanic Plateau, and the western San Bernadino Mountains. Specifically, we test cross-divide contrasts in mean gradient, mean local relief, channel bed elevation, and χ all measured at, or averaged upstream of, a reference drainage area. Our results highlight that cross-divide contrasts in χ only faithfully reflect current divide mobility when uplift, rock erodibility, climate, and catchment outlet elevation are uniform across both river networks on either side of the divide, otherwise a χ-anomaly only indicates a possible future divide instability. The other metrics appear to be more reliable representations of current divide motion, but in natural landscapes, only cross-divide contrasts in mean gradient and local relief appear to consistently provide useful information. Multiple divide metrics should be considered simultaneously and across-divide values of all metrics examined quantitatively as visual assessment is not sufficiently reliable in many cases. We provide a series of Matlab tools built using TopoToolbox to facilitate routine analysis.

We compare the Gilbert metrics to differences in the quantity 'c' across a divide, 117 (Willett et al., 2014). The derivation of and underlying rationale for the calculation of c is 118 discussed in detail in several recent publications (e.g., Harkins  (1) 124 125 where A is upstream drainage area, A 0 is a reference scaling area, q ref is a reference concavity 126 (Wobus et al., 2006), and x' is a dummy variable. A plot of channel elevation vs c for a stream 127 that is equilibrated to a spatially constant uplift rate and erosional efficiency should be a 128 straight line and under these circumstances c can be considered a proxy for steady-state 129 channel elevation. If A 0 is set to unity then the slope of the chi-z plot will equal the normalized 130 channel steepness (Wobus et al., 2006) but is dimensionless. As described by Willett  Similarly, in the eastern Greater Caucasus, c-anomalies suggest that the main divide between 162 northern and southern drainages is unstable, but is predicted to move either south using true 163 base-level ( Figure 2A) or north using constant elevation or the bedrock-alluvial transition 164 ( Figure 2B & C). This highlights that care must be exercised when choosing outlet elevations for 165 c analysis, but also that there may be non-unique answers depending on different, but still 166 reasonable, choices of outlet elevation. 167 168

Proposed Methodology for Use of Divide Metrics 169
While there are some potential problems with the use of c-maps, they are appealing as 170 a data exploration tool as they allow for quick assessment of the relative stability of a drainage 171 divide and associated drainage network. Here we develop similar maps using the three 172 "Gilbert" metrics described above. Mapping the elevation metric only requires coloring a 173 drainage network by channel elevation. The elevation metric is interpreted the same as c in c-174 maps: divides are expected to move from low to high values in the presence of an anomaly in 175 channel elevations ( Figure 1C & D). 176 For the local relief and gradient metrics, we are primarily concerned with average values 177 of these properties near the divide, so a simple strategy of coloring stream networks by 178 upstream running averages of either local relief or gradient is sufficient. These two metrics are 179 more direct proxies for erosion and as such, divides are expected to move from high to low 180 values ( Figure 1C & D). For all four metrics, we are only concerned with the values at the 181 channel heads, which are approximated by choosing a reference drainage area at which to 182 evaluate the values, which we refer to as 'stream endpoints', so for all metrics (including c) a 183 full map of values along streams are not necessary, but provide for useful visuals. 184 Visual comparisons of contrasts in colors across a divide are useful for identifying 185 potentially interesting patterns. However, the perception that a particular divide is unstable can 186 be influenced by visual bias or choices of color scales. To interrogate this further one must 187 assess the actual across-divide differences in the quantities of interest. Additionally, sometimes 188 a single drainage divide may be heterogeneous so it is useful to segment a divide and analyze 189 the stability of these sections individually. We visualize individual divide sections as histograms 190 of values at all of the stream endpoints on either side of a divide ( Figure 1D). In practice, this is 191 useful to assess the degree of overlap or separation between values on either side of a divide. 192 Along with the histograms, we calculate a mean, standard error of the mean, 95% bootstrap 193 confidence interval, and standard deviation for the population of values on either side of a 194 divide. In this study, we primarily use the conservative criteria that a divide is potentially stable 195 according to a given metric if the mean of one side of the divide is within one standard 196 deviation of the mean of the other side. These 'delta' values and their associated uncertainties 197 can then be standardized so that positive and negative delta values of the different metrics 198 indicate similar divide migration direction, providing an easy visual assessment of divide 199 rate and/or erosional efficiency be eliminated, and (2)  For analysis purposes, we segment this drainage divide into 8 sections based on visual 307 inspection of the four metrics and choose break points between portions of the divide that 308 appear to display transitions in at least one of the criteria. Results for all metrics and associated 309 river profiles for the eight divide segments are available in the supplement (Supplemental 310 Figures 1-16). In practice, while the elevation metric was useful in the model results (e.g. Figure  311 4 & 6), the results of the elevation metric are typically equivocal in natural settings we have 312 examined due to large standard deviations (Figure 7). The elevation metric, and indeed all of 313 the metrics, tend to indicate more divide mobility if the standard error of the mean is used to 314 estimate uncertainty. Unsurprisingly, the mean upstream slope and mean upstream relief 315 metrics are very similar, though the similarity of these metrics will depend on the chosen relief 316 radius (e.g., DiBiase et al., 2010). Thus, in comparing divide metrics along the length of the 317 divide for this and subsequent examples, we focus our discussion on c and relief. 318 With the exception of two segments (GC7 & GC8, Figure 7C), c always predicts 319 northward movement of the divide (using the 550m outlet elevation) whereas the relief metric 320 suggests the divide is stable within uncertainty (using the standard deviation) except for two 321 segments (GC3 and GC5, Figure 7C). The means of all metrics (except for GC7) agree in the 322 direction of divide motion and applying a less restrictive uncertainty (e.g. standard error) shows 323 more agreement between all metrics. As discussed earlier, the choice of outlet elevation for c 324 in the eastern Greater Caucasus (e.g. Figure 2) significantly influences predicted divide 325 behavior, with c suggesting southward motion of the divide if 'true base level' is used for the 326 outlet elevation ( Figure 7D). There are no quantitative estimates of erosion rates on either side 327 of the divide so we do not have a way to evaluate the 'right' answer in this setting, but 2016) directly across the divide, suggesting that a change in erosional efficiency is unlikely as a 337 driver. Thus, the simplest interpretation of these results is similar to that posited by Forte et al. 338 (2015,2014), that this indicates the presence of an uplift rate gradient that is 'holding' the 339 divide in place and that the c metric is sensitive to this and indicating the expected reaction of 340 the divide if or when this uplift rate gradient dissipates. 341 342

Kars Volcanic Plateau 343
The Kars Volcanic Plateau is also part of the Arabia-Eurasia collision zone, but the 344 tectonics and local geology are decidedly different than that of the Greater Caucasus to the 345 north. This portion of the collision zone has relatively low rates of active internal deformation 346 (Reilinger et al., 2006), which occur primarily on normal and strike slip faults with some 347 portions of the deformation related to local volcanic features (e.g., Dhont and Chorowicz, 2006;348 Koçyiğit et al., 2001). The Kars Plateau is part of the broader East Anatolian Plateau which lacks 349 mantle lithosphere (Zor, 2008) after a slab detachment or delamination event at ~7-8 Ma 350 (Keskin, 2003;Şengör et al., 2003). The average ~2 km high, roughly dome shaped plateau (e.g., For all the analyzed divides, the relief metric suggests they are stable using the standard 369 deviation criterion and close to stable using the standard error whereas c consistently suggests 370 that divides should move toward the center of the Kars Plateau ( Figure 9D). Using the model 371 results as a basis for interpretation suggests this is likely a case in which a contrast in either 372 erosional efficiency or rock uplift perturbs the c metric (e.g., Figure 5  each other and what is predicted from the erosion rate map ( Figure 10D). If we accept the 415 erosion rate map as accurate and that contrasts in erosion rate across a divide from this map 416 are unequivocal evidence of current or future divide motion, then despite agreement overall, 417 there are examples of both c and relief failing to correctly identify divide instability given 418 uncertainty in cross-divide differences ( Figure 10D). In detail, SB2 and SB8 are cases where c 419 agrees with erosion rates (but not relief) and SB6 is a case where relief agrees with erosion 420 rates (but not c). There are possible interpretations of these deviations, but importantly, these 421 are all cases where (1) a relatively small number of values are used to determine potential 422 divide motion and (2) the determination of divide stability or mobility is dictated by how much 423 overlap or separation in means and standard deviations are required to deem a divide stable or 424 mobile, respectively. This highlights the utility of viewing divide metrics in forms like the 425 histograms used here for evaluating confidence in a given determination and also suggests that 426 there is likely a minimum segmentation length of divides below which the data is simply too 427 noisy to make a clear determination ( Figure 10D stable portion of the divide appears to be between internal plateau streams and streams 432 draining into Big Bear Lake (SB5 & SB6, Figure 10). 433 Using the model results as a means to interpret the divide metrics would suggest that all 434 other divides are currently moving and that any spatial differences in erosional efficiency or 435 uplift rate are absent or sufficiently small such that c is still a viable metric in this setting. This is 436 consistent with known constraints from this region, specifically uniform uplift, simple bedrock 437 geology, and unique relationships between erosion rate and mean channel steepness and 438 erosion rate and mean hillslope gradient (e.g., Binnie et al. where as because c-values only require drainage area measurements, these should be 468 relatively insensitive to data resolution as long as flow routing algorithms are reasonably 469 accurate. In contrast, the divide-scale of the Gilbert metrics make them entirely insensitive to 470 any of the choice of outlet elevation issues that can potentially plague c-maps (e.g. Figure 2). It 471 is also worth noting that none of the metrics are useful for explicitly illuminating past divide 472 motion. All metrics in certain scenarios may be useful in this regard to the extent that current 473 divide motion implies some prior history of divide motion, but because none of these metrics 474 contain any time information, this assumption is hard to validate without independent 475 evidence of past divide motion. 476 In addition to considering multiple metrics, more detailed analyses of differences in 477 values across divides are necessary to fully assess divide stability. In many cases, visual 478 differences in maps of either c or the Gilbert metrics seem to suggest a robust 'anomaly' across 479 a divide, but the histogram of values or the uncertainty on delta values actually show significant 480 amounts of overlap in values, e.g. divide GC2 which in map view seems to highlight an across 481 divide difference in local relief ( Figure 7A), but in detail has relatively similar values in local 482 relief near channel heads ( Figure 7D). A lingering issue is what constitutes suitable amounts of 483 overlap in values across a divide to suggest that said divide is stable or unstable. We do not 484 have any basis for suggesting that the criteria we primarily use (i.e. neither mean value is within 485 one standard deviation of the other for a stable divide) is correct. Comparing predictions using 486 the standard deviation and standard error highlights the importance of the stability criteria, as 487 for example in the Greater Caucasus examples, using the standard deviation with the Gilbert 488 metrics suggested mostly stable divides where as using the standard error suggests more mobile divides. Generally, because standard deviations are larger than bootstrap confidence 490 intervals which are in turn larger than standard errors, using standard deviations bias results 491 towards stable divides (more possibility of overlap) and standard errors bias results towards 492 mobile divides (less possibility of overlap) with bootstrap confidence intervals representing a 493 middle ground. Choosing any estimation of uncertainty is reasonable, but we emphasize that at 494 minimum workers should specify what criteria they are using to judge relative stability or 495 mobility. 496 The software tools provided along with this work allow for relatively easy analysis of 497 drainage divide stability and hopefully will aid the addition of this analysis to routine 498 characterization of landscapes. However, this should always be done in concert with traditional 499 landscape analyses. As described above, the presence of a c-anomaly along with absence of a 500 Gilbert-anomaly at a divide indicates a spatial gradient in uplift rate, erosional efficiency, or 501 both may exist in one or both sets of the catchments that define the divide, but it doesn't 502 provide any information as to the nature of these gradients or their location. For this, maps of 503 streams colored by normalized channel steepness or examining traditional longitudinal or c-504 transformed river profiles would provide more information. Thus, we primarily view these types 505 of metrics as cursory data analysis tools to illuminate areas that necessitate deeper 506 investigation. 507 Finally, fully testing the accuracy of different metrics of divide stability fundamentally 508 requires comparing them to areas for which we have some constraints on erosion rates on 509 either side of divides and thus direct information on the degree of divide mobility. Special 510 attention should be paid to areas where the Gilbert and c metrics disagree, as understanding 511 erosion rate contrasts (or the lack of contrasts) in these settings have the greatest potential to 512 provide more general information on the utility of these metrics in different situations and thus 513 contribute to determining the most reliable topographic expression of divide mobility. 514 The most up to date version of the tools described in this paper are available on github 518 (http://github.com/amforte/DivideTools). This work and development of these tools was 519 supported by EAR-1450970 awarded to AMF and KXW. The implementation of the method to 520 control outlet elevation and check of stream completeness was adapted from methods 521 suggested by Wolfgang Schwanghart via his blog (https://topotoolbox.wordpress.com/). We 522 thank Nicole Gasparini for helpful discussions regarding implementation of the tools. 523

S1. Description of codes included in Github repository
S1.1. DivideStability -For a given area, this routine produces shapefiles (or alternatively rasters convertible to shapefiles within a GIS program) that includes a stream network with values for all four metrics (c, channel elevation, upstream mean relief, and upstream mean gradient) so that the user can produce maps of stream networks colored by these quantities (e.g. Figure 1). To aid in sensible color scaling, values for the 'Gilbert metrics' are normalized to vary between 0 and 1. Because the calculation of c is sensitive to the choice of outlet elevation (e.g. Figure 2), careful control of outlet elevation is essential for meaningful interpretations of c-anomalies.
For this reason, this function allows the user to modify a stream network to remove portions of streams below a minimum elevation. This function also checks to make sure that stream outlets along the edge of the DEM meet the criteria defined and that all drainage basins are 'complete', i.e. that the summation of drainage area is accurate and is not influenced by tributaries that are cut off. Either of these cases could result in artificial c-anomalies, but generally should have no effect on the Gilbert metrics.

S1.2. ChiGrid -The
DivideStability code calculates c along the stream network, but we find it useful to be able to visualize c and an additional metric simultaneously. This code calculates c at every pixel in the DEM so that colored stream networks can overlay this c raster. Similar to DivideStability, a minimum elevation can be specified for calculating the c raster and a check is performed to ensure included drainage basins are complete. S1.3. AcrossDivide -Tool uses the output of DivideStability and allows users to select sections of a divide of interest to perform detailed analysis of divide sections (e.g. Figure 1D). This function provides users multiple ways of defining a divide of interest, but all of them generally function on the idea that the user (or the function itself in the case of automated detection schemes) defines divides of interest by selecting the drainages that define this divide. End results are plots of the distribution and means of the values at the reference drainage area of all four metrics on either side of the divide of interest (e.g. Figure 1) along with a prediction from each metric independently regarding whether and in which direction the divide should move. The prediction of a divide stability or mobility is made on the basis of a user selected assessment of uncertainty and whether the uncertainty of the distributions overlap with the means of the opposing side of the divide. If there is overlap, the divide is considered stable, and if there is no overlap, the divide is considered mobile. This is not meant as an absolute criteria, simply a quick first order assessment. The user can choose to use the standard deviation of the population (default), the standard error on the mean, or the 95% bootstrap confidence interval determined from a 1000 iteration resampling scheme. This function also produces a list of channel head coordinates and their respective values for the four metrics that define the divide of interest.

S1.4. PlotDivideProfiles -To understand the predicted behavior of a divide, it is often necessary to consider the longitudinal profiles of the rivers in question. This functions plots c-elevation
and distance-elevation plots for the streams used to define the divide. Various plotting options exist to allow the user to plot only specific channels and to color drainages by either gradient or relief to compare predictions of individual metrics. S1.5. AlongDividePlot -If the user has defined multiple divide segments, this allows them to produce a plot similar to what is shown in the text (e.g. Figures 7C, 9D, or 10D). In detail, this function will produce three plots for each divide (made up of multiple segments), (1) a plot of divide segment means with uncertainties, (