Drainage divide networks, Part 1: Identification and ordering in digital elevation models

. We propose a novel way to measure and analyse networks of drainage divides from digital elevation models. We developed an algorithm that extracts drainage divides, based on the drainage basin boundaries defined by a stream network. In contrast 10 to streams, there is no straightforward approach to order and classify divides, although it is intuitive that some divides are more important than others. A meaningful way of ordering divides is the average distance one would have to travel down on either side of a divide to reach a common stream location. However, because measuring these distances is computationally expensive and prone to edge effects, we instead sort divide segments based on their tree-like network structure, starting from endpoints at river confluences. The sorted nature of the network allows assigning distances to points along the divides, which can be 15 shown to scale with the average distance downslope to the common stream location. Furthermore, because divide segments tend to have characteristic lengths, an ordering scheme in which divide orders increase by one at junctions mimics these distances. We applied our new algorithm to the Big Tujunga catchment in the San Gabriel Mountains of southern California and studied the morphology of the drainage divide network. Our results show that topographic metrics, like the downstream flow distance to a stream, and hillslope relief attain characteristic values that depend on the drainage area threshold used to 20 derive the stream network. Portions along the divide network that have lower than average relief or are closer than average to streams are often distinctly asymmetric in shape, suggesting that these divides are unstable. Our new and automated approach thus helps to objectively extract and analyse divide networks from digital elevation models.


Introduction
Drainage divides are fundamental elements of the Earth's surface. They define the boundaries of drainage basins and thus form barriers for the transport of solutes and solids by rivers. It has long been recognized that drainage divides are not static through time, but that they are mobile and migrate laterally (e.g., Gilbert, 1877). The lateral migration of divides is a consequence of spatial gradients in surface uplift (positive or negative) as well as stream captures. These frequently accompany tectonic 30 deformation due to shearing, stretching, and rotating stream networks (Bonnet, 2009;Castelltort et al., 2012;Goren et al., 2015;Forte et al., 2015;Guerit et al., 2018), but recent studies have shown that even in tectonically inactive landscapes, drainage divides migrate over prolonged periods of time (Beeson et al., 2017). Such behaviour is consistent with the notion that small and local perturbations can trigger nonlocal responses with potentially large effects on drainage form and area (Fehr et al., 2011;O'Hara et al., 2019). At regional scales, mobile divides can lead to profound changes in drainage configurations 35 and subsequent alterations of base-level and sediment dispersal to sedimentary basins. For example, Cenozoic building of the Eastern Tibetan Plateau margin has been proposed to account for major reorganization of large East Asian river systems and associated changes in sediment delivery to marginal basins (Clark et al., 2004;Clift et al., 2006). Moreover, changes in drainage area that are associated with migrating divides affect river incision rates (Willett et al., 2014), and thus the topographic development of landscapes, which potentially confounds their interpretation in the context of climatic and tectonic changes 40 (Yang et al., 2015).
Although divide networks might be thought of as mirrors of stream networks, there exist fundamental differences between the two. Starting at channel heads, i.e., the tips of stream networks, streams always flow downhill and the upstream area 60 monotonically increases. Stream networks are therefore directed networks that have a tree-like structure, and a natural order, which has been quantified in different ways (e.g., Strahler, 1954;Shreve, 1966). Divide networks, however, are neither directed nor rooted, and they may even contain cycles. They do not obey any monotonic trends in elevation, or other topographic properties that could be easily measured. As a consequence, their ordering is less straightforward. Nevertheless, it is intuitive that some divides (e.g., a continental divide) should have a different order than others. In addition, the structure 65 of divide networks could be important in their susceptibility to drainage captures. For example, higher-order divides may record perturbations longer, as they are farther away from the base level and thus cannot adjust as quickly as lower-order divides. Furthermore, where higher-order divides are close to higher-order streams, drainage-capture events would result in profound changes in drainage area and thus a greater impact on stream discharge and power (e.g., Willett et al., 2014).
In this study, we propose to measure and analyse networks of drainage divides to address questions like: How is the geometry 70 of a divide network related to that of a stream network? Do similar scaling relationships apply? And: Can the divide network be used to infer catchment/drainage dynamics? Empirically driven answers to these questions require tools to study drainage divides, most efficiently from DEMs. We present our study in two separate papers. In the following, part 1, we present a new approach that allows identification and ordering of drainage divides in a DEM. We investigate ways of ordering drainage divide networks and analyse basic statistical as well as topographic properties with a natural example from the Big Tujunga 75 catchment in the San Gabriel Mountains in Southern California. In part 2 of this study (Scherler and Schwanghart, submitted), we present the results from numerical experiments with a landscape evolution model that we conducted to examine the response of drainage divide networks to perturbations.

Drainage divides in digital elevation models 80
Drainage divides are the boundaries between adjacent drainage basins and thus their determination is based on the definition of drainage basins. In a gridded DEM, drainage basins are generally defined through the use of flow direction algorithms. The D8 flow direction algorithm  assigns flow from each pixel in a DEM to one of its eight neighbours in the direction of the steepest descent. As a result, each pixel is associated with a distinct upstream, or uphill, drainage basin. In contrast, multiple flow direction algorithms, such as the D∞ flow direction algorithm (Tarboton, 1997), split 85 the flow from one pixel to several others, which results in some pixels contributing to more than one drainage basin (Schwanghart and Heckmann, 2012). In the following, we only consider drainage basins derived from the D8 flow direction algorithm. In gridded DEMs, flowpaths derived from this algorithm run along pixel centers (Armitage, 2019), and there exists the possibility that two flowpaths run parallel to each other in neighboring pixels. As a consequence, drainage basin boundaries, and thus divides, must be located in between DEM pixels, and have infinitesimal width (Fehr et al., 2009). Another important 90 consequence is that divides will have only two possible orientations that are parallel to pixel boundaries. Our definition of divides is different from one, in which divides are linked to the highest points (pixels) on interfluves (Haralick, 1983). In the case of multiple flow directions, for example, a meaningful position of a drainage divide would be the place within a pixel that partitions the pixel area according to the flow contributions to adjacent drainage basins.
For a given point in a channel network, its drainage basin is uniquely defined to be the upstream area of that point. The drainage 95 divide of that basin, however, does not intersect the channel itself. We thus define drainage divides as lines (or graphs) that mark the margin of drainage basins and that do not cross rivers (Figure 2). When derived from a DEM, these graphs consist of nodes and edges: nodes are located on pixel corners and edges follow pixel boundaries. A meaningful property of divide nodes and edges is that they should not coincide with nodes or edges of the drainage network. When applying the D8 flow routing algorithm to a gridded DEM, with square elevation cells, however, this requirement poses a problem, due to the different pixel 100 connectivity of divides and rivers. Whereas divide nodes can be connected to only four cardinal neighbours, river nodes can be connected to eight different neighbours. In consequence, divide nodes may exist that coincide with diagonal edges of drainage networks (Figure 2). In a gridded DEM this issue could be resolved with a D4 flow direction algorithm; or, more generally, this issue could be avoided if flow is only allowed orthogonal to cell boundaries. In our approach, we nonetheless adopt the D8 flow direction algorithm and allow for spatial congruence of streams and divides. In practice, such issues mainly 105 arise near confluences ( Figure 2).

Drainage divide networks
Analogous to streams, drainage divides are typically organized into tree-like networks (Figure 1), although cycles that correspond to internally-drained basins may exist. Because of the directed flow of water, stream networks can be regarded as directed graphs that start at channel heads (the leaves of the tree) and end at an outlet or a river mouth (the root of the tree). 110 Flow directions in stream networks can be easily derived from node elevations (e.g.,  and the hierarchy of streams can be related to their upstream area, for example. In contrast, drainage divides have no inherent direction, and there exists no terrain property, like elevation, that could be used to assign a direction to them. A meaningful metric for ordering divides may be the average branch length , i.e., the average distance Λ [m] one would have to travel down on either side of a divide to reach a common stream location ( Figure 3). However, measuring this distance 115 requires that the common stream location and the entire path leading to it be contained in the DEM. Because this may not be true for a significant part of the divide network in a DEM and because measuring this distance is computationally expensive , it is not very practical.
Instead, we suggest that directions can be derived from the tree-like structure of drainage divide networks. Analogous to a parcel of water that travels down a river from its source to its mouth, we propose to start at the leaves of the tree, which we 120 call the endpoints of the divide network (Figure 2), and incrementally move down the branches (Figure 4). Note that the term "move down" does not refer to elevation, but to the hierarchy of the divide network. Where two or more drainage divides meet, they form a junction. We call individual parts of drainage divides that link endpoints and junctions, junctions and junctions, or endpoints and endpoints, the drainage-divide segments, and we refer to the ends of divide segments as segment termini to avoid confusion with endpoints. At junctions with more than one unsorted divide segment, the sorting process pauses because it is 125 not obvious in which direction the sorting shall continue. However, in the absence of cycles (internally drained basins), each junction will reach a point in the sorting loop when there is only one unsorted divide segment left so that the sorting can continue ( Figure 4). This condition assures that the divide segments are correctly sorted in a tree-like manner, but it fails when encountering a cycle. As we will show later, the average branch length Λ scales linearly with the maximum distance from an endpoint along the sorted divide network, as well as the maximum number of divide segments (or junctions), both of which 130 are more easily computed. We thus propose to order the nodes and edges of the divide network, by their maximum distance from divide endpoints, measured either in map units or in the number of divide segments. From now on, we call the distance measured in map units along the directed divide network the divide distance (d d ).

Divide algorithm 135
We implemented the above-described way of extracting and ordering drainage divides from a DEM in the TopoToolbox v2 (Schwanghart and Scherler, 2014), a MATLAB-based software for topographic analysis. Figure 5 shows the workflow of our approach, which consists of the following steps: (1) For a given DEM, we first define a stream network, based on the D8 flow direction algorithm and a threshold drainage area (Figure 5a,b). The lower the threshold, the more detail will the stream and divide networks be. 140 (2) We extract drainage divides, based on drainage basin boundaries that we obtained for drainage areas at tributary junctions and drainage outlets ( Figure 5c). Initially, each drainage basin boundary is composed of one divide segment that connects two endpoints, and junctions do not yet exist. These divide segments do not cross any rivers but their nodes may coincide with stream edges (Figure 2). We remove redundant divide segments in the collection of divides, which arise from nested and adjoining drainage basins. As a result, we are left with a 145 set of unique divide segments, which, however, may be continuous across junctions or terminate where they should be continuous ( Figure 6).
(3) We next organize the collection of divide segments into a drainage divide network (Figure 5d). This is the core of the algorithm, in which we identify endpoints and junctions, merge broken divide segments, and split divide segments at junctions ( Figure 6). Our algorithm distinguishes between junctions, endpoints, and broken divide 150 segments, by computing for each node of the divide network the number of edges linked to it, the number of segment termini linked to it, and the existence and direction of a diagonal flow direction. For example, most nodes with two edges and two segment termini correspond to a broken segment and need to be merged, unless they coincide with a stream and merging them would make the resulting divide cross that stream ( Figure 6). See Appendix 7.1 for more details. 6 (4) Finally, we sort the drainage divide segments within the network (Figure 5e). The algorithm iteratively identifies segments that are connected to endpoints and removes them from the list of unsorted divide segments until no divide segments are left ( Figure 4). This step assigns a direction to each divide segment and transforms the divide network into a directed acyclic graph. For the sorted divide network, we then compute the divide distance, i.e., the maximum distance from an endpoint along the sorted divide network (Figure 5f). 160 After the sorting, we also assign orders to divide segments, based on the ordering of stream networks, first introduced by . We adopted both Strahler's (1954) and Shreve's (1966) rules of stream ordering, and added a third rule that we call Topo. All ordering schemes start with a value of one at endpoints and progressively update divide orders at junctions based on the following rules:

Strahler:
= max(min{ } + 1; ω i ) Eq. 1 Shreve: where i are the divide orders of the joining divide segments and k is the divide order of the following divide segment. In 165 the Strahler ordering scheme, the order increases by one if the joining divide segments have the same order, otherwise it remains at their maximum order. In the Shreve ordering scheme, the resulting divide order is the sum of those of the joining divide segments; and in the Topo ordering scheme, divide orders increase by one at each junction. Junctions typically link three different divide segments, but up to four can occur ( Figure 6). Based on the Strahler ordering scheme, the bifurcation ratio , can be derived from : 170 is the number of divide segments of order .
As aforementioned, our divide algorithm currently does not handle internally drained basins. Whereas the divides of internally drained basins are easy to identify, they are not easily sorted in a meaningful manner. In fact, the distance to a common stream location (Λ) is undefined for a divide of an internally drained basin. At the moment, the sorting procedure ( Figure 4) stops at such divides, because the divide segments cannot be assigned a direction. In consequence, also parts of the divide network that 175 potentially lie beyond an internally-drained basin, and for which the distance to a common stream location is defined, cannot be reached anymore. While we are working on a solution to this issue, our algorithm is currently best applied to acyclic drainage divide networks.

Topographic data and analysis
We investigated basic characteristics of drainage divide networks using a 30-m resolution DEM from the 1-arc second Shuttle 180 7 and landscape rejuvenation as the river incises into a relict pre-uplift landscape (DiBiase et al. 2015). We pre-processed the DEM by carving through local sinks (Soille et al. 2003) to avoid artificial internally drained basins, and we obtained a stream network based on a minimum upstream area of 0.1 km 2 . We note that this threshold is likely well within the zone of debris 185 flow instead of fluvial incision (Stock and Dietrich, 2003); but for the purpose of our analysis, this is irrelevant.
We analysed the divide network, its planform geometry, and its relation to topography. Planform geometry is studied using statistical analysis of the number and length of divide segments of different orders. Topographic analyses are based on metrics that we determined for the entire DEM and that we subsequently associated, or mapped, to divide edges and entire divide segments. As topographic metrics, we focus on hillslope relief (HR) and horizontal flow distance to the stream network (FD). 190 HR was defined to be the elevation difference between a point on the divide and the point on the river that it flows to. To quantify the morphologic asymmetry of a divide, we propose to use the across-divide difference in hillslope relief (∆ ), normalized by the across-divide sum in hillslope relief (∑ ), and call its absolute value the divide asymmetry index (DAI): The DAI ranges between 0 for entirely symmetric divides and 1 for the most asymmetric divides. Note that this index is based only on values of hillslope relief (HR). Theoretically, a divide with equal amounts of HR on either side of a divide, but contrasts 195 in flow distance (FD) and thus slope angle, would yield a DAI of zero. However, due to the definition of streams by a minimum drainage area, this hardly ever occurs. In addition, such cases can be identified by cross-divide differences in FD.

Basic divide statistics
We applied our divide algorithm to the Big Tujunga catchment, and the resulting divide network for different ordering schemes 200 is shown in Figure 7. Because the Shreve and Topo ordering schemes yield larger ranges in divide orders, their visualization allows for greater differentiation compared to the Strahler ordering scheme. Differences in the visual appearance of the divide network due to the ordering scheme are also apparent at the root node, i.e., the junction that is encountered last in the ordering process (black arrows in Figure 7). In the Topo ordering scheme, divide orders increase by one during each sorting cycle, so that the last divide segments will have orders that are different by not more than one. In contrast, the ordering rules of the 205 Strahler and Shreve scheme (see Eq. 1 and Eq. 2) may yield unequal orders during the sorting, so that the divide orders of the last divide segments may be different by more than one. In the Big Tujunga catchment, the basin area and thus the number of divide segments is larger north of the Big Tujunga River, compared to south of it. In consequence, both Strahler and Shreve divide orders increase more rapidly along the northern perimeter compared to the southern, and the junction encountered last during the sorting process (at the root of the tree) opposes divide segments with orders of 7 and 6 for Strahler, and 1463 and 210 772 for Shreve, in the north and south, respectively (Figure 7). For the Strahler ordering scheme, the frequency distribution of divide segments decreases exponentially with divide order (ω), which is consistent with Horton's law of stream numbers , and corresponds to a bifurcation ratio of = 3.89 ± 0.30 (standard error). The bifurcation ratio of the associated stream network is 5.39 ± 0.87.
We computed divide segment lengths for different drainage area thresholds (A min ) and show the associated empirical 215 distribution functions in Figure 8a. Divide segment lengths are not normally distributed, but can be reasonably fitted with a gamma distribution. However, the fitted gamma distributions predict systematically higher probabilities for shorter divide segments and lower probabilities for longer divide segments compared to the actual data. For the different drainage area thresholds that we tested, the shape parameter of the fitted gamma distribution (k) attains values that range between 0.87 and 1.75. In general, the average length of all divide segments increases with the drainage area threshold used for deriving the 220 stream network, simply because both the stream and the divide network extend to finer branches. For a drainage area threshold of 0.1 km 2 , the average length across all divide orders is 442 ± 323 m (±1σ), compared to an expected value (= k) of ~442 from the fitted gamma distribution. The average length for different divide orders tends to be slightly lower at ω < ~40 (based on the Topo ordering scheme), compared to ω ≥ ~40 (Figure 8b).
We quantified the average branch length, i.e., the average distance one would have to travel down on either side of a divide to 225 reach a common stream location (Λ), for 100 randomly chosen divide edges per divide order in the Topo ordering scheme.
Although the maximum order for which Λ can be determined (because of the size of our DEM) is limited to ω ≤ 55, results demonstrate that Λ (in km) increases linearly with 0.36 × divide order (ω), and 1.11 × divide distance (d d ) (Figure 9). The linear scaling of these two relationships is a consequence of the similarity of segment lengths for different divide orders ( Figure   8b). Whereas the Topo ordering scheme can be approximated by dividing the divide distance by the expected divide segment 230 length, this does not hold true for the Strahler and Shreve ordering schemes, which yield relationships between Λ and divide order that are non-linear (not shown).

Drainage divide network morphology of the Big Tujunga catchment
We next studied the morphology of the drainage divide network from the Big Tujunga catchment. Because the divide morphology consists of parts that lie within the catchment as well as parts that lie outside of it, we analysed the entire drainage 235 divide network from the DEM as shown in Figure 5. Although the drainage divide network is truncated along the DEM edges, the following analysis is insensitive to this issue. Figure 10 shows the drainage divide morphology of the Big Tujunga catchment, based on a stream network that was derived from a drainage area threshold of 1 km². Topographic metrics are shown for each divide edge (cf., Figure 1). Whereas across-divide mean flow distance ( ̅̅̅̅ ) varies between 0 and ~3000 m, mean hillslope relief ( ̅̅̅̅ ) varies between 0 and ~800 m. It is notable that the biggest range in values occurs at divide distances 240 <10 km, whereas divides at higher distances appear to hover around characteristic values that are controlled by the drainage area threshold used to derive the stream network ( Figure 10e). It should be noted that divide edges at low divide distances are much more abundant compared to those at higher distances, or equivalently, at higher divide orders (Figure 7). At increasingly higher drainage area thresholds however, the frequency of divides at low order and distance decreases more rapidly compared to the frequency of high orders and distances. The average (±1σ) ̅̅̅̅ and ̅̅̅̅ values for a drainage area threshold of 1 km 2 are 9 1325 ± 350 m and 341 ± 129 m, respectively. The empirically determined average ̅̅̅̅ for all tested drainage area thresholds are consistent with Hack's law (Hack, 1957), which relates the length L of the longest stream in a catchment to its drainage area A, according to = ℎ . Values of k a and h are ~1.6 and ~0.5, respectivelysimilar to values observed elsewhere (Hack, 1957). Combining ̅̅̅̅ and ̅̅̅̅ yields average slope values that vary between ~20° and ~8°, for drainage area thresholds between 0.1 km 2 and 10 km 2 . The mean slope value of the entire DEM is ~21.5°, which suggests that the lower drainage area 250 threshold better confines the divides to hillslopes as compared to lower-sloping channels.
Based on the observation of characteristic values of ̅̅̅̅ and ̅̅̅̅ , we sought to identify parts of the divide network that have anomalously low relief or are anomalously close to a stream. Instead of mean values, we turned towards across-divide minimum values of flow distance ( ) and hillslope relief ( ), as these would be more sensitive to deviations on either side of a divide. In addition, we compared these values with the divide asymmetry index (DAI), as we expected that 255 anomalous divides may also be topographically asymmetric (e.g., Whipple et al., 2017). Figure 11 shows how HR min , FD min , and DAI vary with distance along the divide network of the Big Tujunga catchment. Notable deviations from average values (cf., Figure 10) occur at ~22-25 km, ~ 41-45 km, and >54 km divide distance, and are typically associated with asymmetric divides ( Figure 11). Highly asymmetric divides are furthermore found at low divide distances (<5 km) and typically coincide with low values of , whereas could be high or low. The systematic decrease in and , concurrent with an increase 260 in the DAI at higher divide distances prompted us to enquire the geographic position of these divides and how they compare to the surrounding landscape. We thus imposed thresholds to identify anomalously low (HR min <200 m) and asymmetric divides (DAI>0.5), which are less than 1000 m from a stream (FD min <1000m) ( Figure 12). Beheaded streams as well as sharp-crested and shortened hillslopes identified in high-resolution satellite imagery (Figure 12b-e) support the impression that these divides are mobile and migrating in the direction of lower and sometimes shorter . Most of these divides can be seen to border 265 regions of contrasting local relief (Figure 12a), and many cluster along the eastern edge of the catchment. The predicted migration direction indicated in Figure 12a is derived from the orientation of the divide segments and their mean DAI magnitude. If correct, most of the divide migration along the southern and eastern edge of the catchment, from higher to lower relief, would result in area loss for the Big Tujunga catchment.

Extraction and ordering of drainage divide networks
Our new approach allows routinely extracting drainage divides from any DEM without internally drained basins. We have shown that the maximum divide distance d d , calculated as the maximum distance along the (directed) divide network from an endpoint, is a meaningful metric for ordering drainage divide networks, as it scales linearly with the average branch length Λ, i.e., the average distance one would have to travel down on either side of a divide to reach a common stream location ( Figure  275 9). In contrast to the average branch length, however, the divide distance is more easily and rapidly calculated and is less prone to edge effects that inhibit the ordering of divides (cf., . However, whenever a drainage basin intersects the edge of a DEM, it's truncation will likely produce a spurious drainage divide. Furthermore, calculated divide distances are most likely lower than they would be for a larger DEM, similar to the reduction of upstream area along a stream network. Truncated drainage basin boundaries should therefore be avoided, or discarded from analyses that rely on correct 280 divide distances. The proposed sorting procedure ( Figure 4) recovers the tree-like structure of the divide network and allows the derivation of divide orders, analogous to the well-known stream orders. Because divide segments have similar mean lengths across all divide orders (Figure 8), divide orders derived with the Topo ordering scheme can serve a similar purpose as divide distance. Shreve (1969) studied link lengths in stream networks and concluded that their distribution is better described with a gamma 285 distribution compared to an exponential or log-normal distribution. Results from the Big Tujunga catchment support this conclusion with respect to divide segment lengths, although systematic deviations can be observed ( Figure 8a). It needs to be tested with more observations whether these deviations are inherent to drainage divide networks in general, or whether they could hold clues about the dynamic state of a landscape.
An advantage of characterizing the divide network by distance instead of orders is that the divide distance is invariant with 290 respect to the chosen drainage area threshold, whereas divide orders are not, because they depend on the total number of divide segments and junctions. Further differences are apparent at the root node, which may oppose divide segments with orders that differ by more than one (Figure 7). In the case of the Big Tujunga catchment, Strahler orders are not that different across the root node, but in a different landscape that could well be the case. This issue is more prevalent in the case of Shreve ordering, but it is avoided with the Topo ordering scheme. Furthermore, the non-uniform distribution of divide segment lengths ( Figure  295 8) influences how similar or dissimilar divide distances of the meeting divide segments are at the root node. If the average divide segment length of trees that meet at the root node are different, divide distances will make a jump, even if divide orders are similar. In the Big Tujunga catchment, the divide distance jump at the root node is 5400 m.
Divide orders derived with the Strahler ordering scheme, can be used to investigate how the divide network conforms to  laws of network composition. In the Big Tujunga catchment, for example, the bifurcation ratio of the divide 300 network (R b ~3.9) is lower than that of the stream network (R b ~5.4). This may in part be due to the fact that we analysed only a part of the divide network; divide segments that originate from the main catchment boundary and extend outwards are not included in the statistics. Including those in the calculation yield R b ~4.6 for the divide network, still lower than the bifurcation ratio of the stream network. Nevertheless, these bifurcation ratios are similar to published bifurcation ratios of different natural stream networks (e.g., Tarboton et al., 1988), supporting the expected similarity of the stream and divide network topology. 305 However, more observations from different landscapes are needed to assess systematic differences and commonalities between divide and stream networks.

Drainage divide mobility in the Big Tujunga catchment
Based on the observation of characteristic values of minimum hillslope relief (300-500 m) and minimum flow distance (1000-1800 m), we identified drainage divides in the Big Tujunga catchment that are anomalously low, close to a stream, and 310 asymmetric ( Figure 11, Figure 12). These geometric properties suggest the existence of windgaps, hillslope undercutting by rivers, and spatial anomalies in erosion rates, which are diagnostic for past or ongoing mobility of drainage divides. Anomalous drainage divides are particularly frequent along the eastern edge of the catchment, where an area of low hillslope angles and local relief (Figure 12), the so-called Chilao Flats, is bordering a steep catchment to the south and east of it. This high-elevation low-relief area is thought to represent a relict peneplain surface that has been uplifted during growth of the San Gabriel 315 Mountains, and is currently destroyed by headward incision of rivers (Spotila et al., 2002;DiBiase et al., 2015). Cosmogenic 10 Be-derived erosion rates confirm lower erosion rates in the Chilao Flats area compared to the surrounding steeper catchments (DiBiase et al., 2010), which ought to drive divide migration and drainage area loss in the headwaters of the Big Tujunga catchment, consistent with our results.
We identified another stretch of anomalous divides along the southern margin of the Big Tujunga catchment (Figure 12a, d), 320 part of which is coincident with the trace of the San Gabriel Fault, which follows the orientation of the valley . Reduced relief in a ~1-km wide zone around this fault is also observed farther to the east, along the West Fork of the San Gabriel River , suggesting weaker rocks closer to the fault (e.g., . Other anomalous divides in this area, as well as along the northern margin of the Big Tujunga catchment, show signs of mobility by one-side shortened hillslopes and beheaded valleys (Figure 12b, d). We thus suggest that most, if not all of the anomalous 325 divides we identified based on hillslope relief, flow distance, and divide asymmetry, are in fact unstable and migrating with time. Because most of the peripheral divides indicate drainage area loss of the Big Tujunga catchment, these area changes ought to result in changes in stream power (Willett et al., 2014), which complicate the interpretation of stream profile knickpoints in a tectonic framework (DiBiase et al., 2015).

6
Conclusions 330 In this study, we presented an approach to objectively extract and analyse drainage divides from DEMs. We argued that divides can be ordered in a meaningful way based on the average distance one would have to travel down on either side of a divide to reach a common stream location, and we have shown that this distance can be well approximated by the maximum alongdivide distance from endpoints of the divide network, which we termed the divide distance. We have also shown that the treelike structure of divide networks lends itself for topological analysis similar to stream networks, and we introduced an ordering 335 scheme (Topo), in which divide orders increase by one at divide junctions. Because divide segments tend to have characteristic lengths, the Topo ordering scheme mimics the divide distance. Topographic analysis of the drainage divide network of the Big Tujunga catchment yielded characteristic values of flow distance and hillslope relief that can be shown to depend on the drainage area threshold, with which the stream network was derived. Based on these characteristic values and a minimum divide distance of ~5 km, below which we observed large scatter, we identified divides that have anomalously low hillslope 340 relief, are close to rivers, and are asymmetric in shape. We interpret these divides to be mobile and indicating beheaded valleys or future capture events.

Classification of divide network nodes
Once the drainage divides are defined, based on the outline of drainage basins, and redundant divide segments are removed, 345 they compose a network D = (V,E) which is defined by a set of vertices V (or nodes) and a set of edges E, each of which is associated with two distinct vertices. However, D may contain some divide segments that do not end at junctions or that terminate at nodes that are neither junctions nor endpoints. To create the divide network, we have to identify divide endpoints and junctions, as well as divide segments that need to be merged or parted. We achieve this by computing for each node x i the number of edges (1 to 4) and divide-segment termini (0 to 4) that exist in D, and identifying whether the node coincides with 350 a stream edge (0/1) ( Table A1). Based on these criteria, we classify nodes to be endpoints (EP), junctions (J), or broken segments (BS). In the case of nodes with three edges, three segment termini, and the presence of a stream edge, we also check which of these edges, if connected, would cross a stream, to distinguish the EP from the BS. After this classification, we are able to merge broken segments, split segments at junctions, and thus update D, which now contains all the divide segments that compose the drainage divide network. 355

Code availability
The divide algorithm developed in this study has been implemented in the TopoToolbox v2 (Schwanghart and Scherler, 2014).
The codes will be made available with the next TopoToolbox release, and shall be accessible through https://github.com/wschwanghart/topotoolbox.

Competing interests 360
The authors declare that they have no conflict of interest.         (±1σ) values of mean flow distance (red) and mean hillslope relief (blue) for different drainage area thresholds. Average values were determined from all divide edges at a divide distance >5 km, to minimize the influence of divides that are close to streams simply due to their proximity to confluences.  Divide migration is a time-dependent process that is difficult to quantify. While the effects of regional-scale drainage captures may be preserved within sedimentary archives (e.g., Clift et al., 2006), this is unlikely for smaller scale drainage captures or 25 gradual divide migration. In such cases, most studies rely on topographic indicators. Mobile divides are typically inferred from post-drainage capture evidence: distorted drainage structures, low divides (wind gaps), or high tributary junction angles (e.g., Clark et al., 2004) (Figure 1). However, divide mobility may also be expressed in the topography without major drainage captures or flow reversals, but as a result from more gradual migration of divides. Willett et al. (2014) inferred drainage divide mobility from across-divide differences of χ values, a proxy for steady-state river channel elevation (Perron and Royden, 2013). 30 They argue that changes in drainage area within mountain ranges, e.g., due to tectonic strain of the crust (Yang et al., 2015), may commonly lead to relative differences in incision rate and the formation of low-relief landscapes that are bordered by migrating divides. Whipple et al. (2017a,b) however argued that the time scale of such changes is too short to profoundly affect mountainous landscapes. Instead, they argue that transient low-relief landscapes, such as those in Southeastern Tibet, are more likely to be formed by regional changes in rock uplift rate and upstream propagation of knickpoints between the adjusted and 35 unadjusted parts of landscapes. They also cast doubt on the ease of comparing across-divide differences in drainage network geometry (i.e., χ values) where the common base level is far and where opposing rivers may incise into areas of different rock types, different rock uplift rates, or different climates. Whipple et al. (2017a,b) instead proposed that the shape of drainage divides themselves hold clues about their mobility (Figure 1). Amongst the topographic parameters that they tested are acrossdivide differences in channel elevation at a reference drainage area, mean headwater hillslope gradient, and mean headwater 40 local relief.
In summary, several different metrics have been proposed that may allow quantification of divide mobility in both natural and modelled landscapes. Forte et al. (2018) compared the performance of these metrics with a landscape evolution model in which they induced divide mobility, and concluded that across-divide differences in relief or gradient better depict divide motion than χ. In their analysis however, they focused on divide motion that is perpendicular to the regional drainage direction and 45 averaged divide migration rates as well as topographic metrics across the entire width of the model domain, so that each time step is associated with single values for divide migration rate, erosion rate difference, and the tested topographic metrics. In part 1 of this study (Scherler and Schwanghart, submitted), we presented a new approach for the automatic identification and ordering of drainage divide networks in a DEM. Here, we present experiments with a numerical landscape evolution model that we conducted to investigate how drainage divide networks respond to different perturbations, including fault activity and 50 differences in erodibility. In contrast to previous studies that examined the response of drainage divides to perturbations, we studied the entire drainage divide network in an objective manner and examined how different portions of the divide network respond to perturbations. In addition, we tested the utility of a new metric that quantifies the connectivity of drainage divide junctions.

Landscape evolution model
We studied the response of divide networks to stream captures and divide migration, using the TopoToolbox Landscape Evolution Model (TTLEM; Campforts et al., 2017). In our experiments, we modelled the topographic evolution of a 20 km × 20 km square block (50-m node spacing), subject to uniform rock uplift, stream power-based fluvial incision (e.g., Howard and Kerby, 1983;Whipple and Tucker, 1999), as well as hillslope diffusion (e.g., Culling, 1963): 60 where z is elevation (L), t is time (T), U is rock uplift rate (L T -1 ), A is upstream area (L 2 ), S is local channel slope (L L -1 ), K r is a parameter of the efficiency of river incision (T -1 ) (K r = 1×10 -5 yr -1 ), and m and n are dimensionless constants with values of 0.5 and 1, respectively. The last term on the right-hand side depicts elevation change due to the divergence in diffusive hillslope transport qs (L 3 L -1 T -1 ), which we consider to be a linear function of hillslope gradient: = − ∇z, where D is the diffusivity (L 2 T -1 ) of soil creep (D = 2×10 -3 m 2 yr -1 ). All four edges of the block were fixed in elevation (z = 0 m), which 65 forced rivers to flow outwards. The uplift rate (U = 1 mm yr -1 ) was constant in all models. Our choice of parameter values was guided by the study of Whipple et al. (2017b), who tested a wide range of rock uplift and erosional efficiency parameters and found almost no difference of divide mobility in models with and without hillslope diffusion and for n values of 1 and 2.
We started from a flat surface with imposed random noise and ran the experiment for 30 Myr, until the topography reached a steady state. The result of this model, which we termed 'Initialize', provided the starting point for four other models that we 70 ran for 10 Myr (Figure 2). The model 'Reference' included no further changes. In the model 'Rotating', we included a circular (10-km diameter) left-lateral strike slip fault that was active throughout the experiment. Strike-slip faults are well known for enforcing drainage captures and thus divide mobility (e.g., Castelltort et al., 2012;Duval and Tucker, 2015). Although the rotating block has, to our knowledge, no real-world equivalent, this model setup represents a convenient way of simulating extended periods of strike-slip faulting, as the fault does not intersect the model boundary (cf., Braun and Sambridge, 1997). 75 The fault slip rate was fixed at 4 mm/yr, which corresponds to an angular velocity of 8×10 -7 rad/yr, resulting in ~460° of total rotation during the model run. We note that the rotating movement requires interpolation and thus leads to numerical diffusion of elevations within the rotating disc. However, the resulting change in total volume by interpolation is <0.03% of the volume uplifted during the same time and therefore small. The model 'Inclined' included 1-km thick and 5-km spaced layers of 50-% reduced erosional efficiency of rivers (K r ), dipping 30° towards northwest. The 'Inclined' model is representative for a 80 landscape, in which rivers incise into tilted sedimentary rocks of non-uniform rock strength, similar to what has been studied by Forte and Whipple (2018). During the experiment, the combination of surface lowering and inclination resulted in the resistant layers to regularly sweep from southeast to northwest across the simulated landscape. The model 'Spheres' included 4 represent incision of a region that is characterized by country rocks with more-resistant magmatic intrusions. The expected 85 behaviour of this model is similar to the landscape response to localized perturbations studied by .

Topographic analysis
We analysed the modelled topography and the associated drainage divide network. For each modelled topography and at each time step (dt = 40,000 years), we first computed flow directions and flow accumulation, and subsequently identified the stream network using a drainage area threshold of 0.2 km 2 . We next derived the drainage divide network on the basis of the stream 90 network and using the algorithm proposed in Scherler and Schwanghart (submitted). We calculated divide distances as well as divide orders based on the Topo ordering scheme (cf., Scherler and Schwanghart, submitted). As topographic metrics, we included elevation (z), as well as hillslope relief (HR) and horizontal flow distance to the stream network (FD). HR was measured as the elevation difference between a point on the divide and the nearest river location as measured by the distance along local flow directions. We also computed χ values on either side of a divide, using a reference area A 0 of 1 m 2 , a reference 95 concavity θ ref of 0.45, and setting the base level x b to 0 at the edge of the model domain (e.g., Perron and Royden, 2013): For each divide edge, we computed these topographic metrics as well as the erosion rate (Eq. 1) for the two neighbouring pixels that belong to adjacent drainage basins and denoted the across-divide minimum, maximum, sum, difference, and average in any one metric , as , , ∑ , ∆ , and ̅ , respectively. Erosion rates were based on the erosion rate of the first downslope stream pixel to reduce the impact of local noise along hillslopes. Topographic metrics of entire divide segments are 100 based on those of the divide edges that it is composed of. For quantifying across-divide differences in topographic metrics as well as erosion rates, irrespective of the actual values, we used normalized indices of the form ∆ / ∑ . One such index that we frequently used in our study is the divide asymmetry index ( = |∆ /Σ |), which is the absolute value of the normalized hillslope relief difference, and which ranges between 0 (symmetric) and 1 (most asymmetric).
The above described across-divide differences in topographic metrics essentially aim to quantify divide mobility. In contrast, 105 Spotila (2012) studied the stability of divides and argued that divide junctions and pyramidal peaks are more stable than solitary linear divides and might therefore act as anchor points for drainage divide networks. He proposed that divide junctions are more difficult to erode than linear divides, due to their greater volume of topography per unit area, their greater mechanical stability, and their reduction of confluent flows (Spotila, 2012). He also suggested that the stability of divide junctions is related to the number of joining drainage divides. Because divide junctions obtained from our algorithm cannot connect more than 110 four divides segments (cf., Scherler and Schwanghart, submitted)and most often connect three segmentswe introduce a new metric to quantify divide junction connectivity, C J : Eq. 3.
We define C J to correspond to the sum of the ratios of the Euclidean distance, d, and the divide distance, d d , of all divide edges, n, within a specified maximum divide distance, d d,max , times the ratio of the cell size, dx, and d d,max . The dimensionless quantity C J is sensitive to the number of divides within a given divide distance from a junction, weighted by their orientation towards 115 the junction (Figure 3). The value d d,max reflects the divide distance, over which differences in junction connectivity are measured. For junctions that connect a constant number of straight and infinitely long divide segments, C J is not sensitive to the value of d d,max . However, for actual junctions, C J is typically sensitive to the value of d d,max , because as d d,max grows, increasingly more junctions are at a distance d d <d d,max of a specific junction and thus the number of divide segments grows with d d (Figure 3a). In general, C J will be sensitive to the position of a junction within the drainage divide network, if the 120 junction's maximum divide distance from an endpoint is smaller than d d,max . In other cases, C J will provide a measure of how connected a junction is within a network, or, in other words, how prominent the junction is compared to other junctions in the network. In this study, we used a d d,max value of 5 km.

3
Results of numerical experiments

General behaviour 125
The simulated landscapes along with their drainage divide networks at the end of the numerical experiments are shown in which forces them to steepen and to induce surface uplift upstream, as opposed to their appearance at drainage divides, which increases the height of the divide, but which does not induce drainage divide migration at larger scale (Figure 5d).

Network topology
We first analysed the response of the entire drainage network topology to the perturbations by quantifying the aggregated 140 length of divide segments as a function of their order (Figure 7). The first few Myr of the Initialize model are characterized by large changes in divide lengths and orders. Initially, the divide network extends to orders as high as 100, but rapidly contracts as the drainage network becomes dendritic. After about 5 Myr, the highest orders are down to 60. Subsequent changes result in some scatter of the divide lengths, but not in the range of divide orders. Compared to the Reference model, in which the divide network structure does not change anymore, the Rotating, Inclined, and Spheres models exhibit changes in the divide 145 network, mostly at divide orders greater than ~20. This observation is related to the fact that low-order divides are distributed across the entire model domain and their number is accordingly high. Any of the perturbations we imposed only affects some of these divides, and thus the impact on their average length is rather small. In contrast, high-order divides are constrained to the highest parts of the modelled land surface and their numbers are much lower. The imposed perturbations typically affect a greater portion of them and hence the scatter in divide lengths is wider. In the Rotating and Spheres models, we also observed 150 that maximum divide orders occasionally extend to higher values, but these changes are rather small. We note that the above observations also prevail when considering divide distance instead of divide order, because the two are linearly related (Scherler and Schwanghart, submitted), with divide distance ~430 m × divide order. The factor on divide order corresponds to the mean length of the divide segments.

Topographic metrics 155
We next studied how the above-described disturbances affect drainage-divide metrics during the simulations (Figure 8). For all models, we computed the averages of topographic parameters, measured at drainage divides of specific divide distance intervals (Figure 8a-d). As in the analysis of divide segment lengths by order, it should be kept in mind that the numbers of divide segments, or their aggregated lengths (Figure 7), are much higher for low orders and distances as compared to higher ones. For reference, a divide order of 20 corresponds to a divide distance of approximately 9 km. In the Initialize model, all of 160 the studied metrics attain a constant value that remains unchanged in the Reference model (Figure 8), and which may or may not depend on the divide distance. For example, the mean elevation and junction connectivity (C J ) clearly increase with divide distance, whereas the flow distance exhibits only minor dependence on divide distance, and hillslope relief appears unrelated to divide distance. The dependency of some metrics on divide distance is partly explained by the model setup. Although divides with low distances occur also at higher elevation, the bulk of them is near the model edge, close to zero elevation. In contrast, 165 divides at high distances are exclusively found near the centre of the model, where elevations are also high. Similarly, the junction connectivity (C J ) is high in the model centre, where the divides are far from most of the endpoints, which are more abundant near the edges of the model.
It is also worth noting that none of the normalized across-divide differences in the topographic metrics attain zero values in the Reference model. This means that even at topographic steady state, there exist residual across-divide differences in hillslope 170 relief, flow distance and χ. In the case of χ, these also depend on the divide distance, and are greater closer to the model edge, where divide distances are low. In the perturbed models, we observed fluctuations in all topographic metrics, although of different magnitudes. For example, comparison of across-divide differences in erosion rate with differences in hillslope relief, flow distance and χ (Figure 8e-h) shows that the normalized difference in hillslope relief (i.e., the divide asymmetry index) is sensitive to drainage divide mobility in all perturbed models, whereas across-divide differences in χ and flow distance are (C J ) metric attains values in the perturbed models that are temporally averaged quite similar to the constant values in the Reference model. In many cases, the deviations from the Reference model are greater the higher up in the divide network, i.e., for higher divide distances. This pattern is particularly well visible in the Inclined model, where the amplitudes of the oscillations in all of the parameters increase with divide distance. 180

Minimum hillslope relief and flow distance
Motivated by the observation of constant values in hillslope relief and flow distance in the Reference model, as well as in actual landscapes (Scherler and Schwanghart, submitted), and by our expectation that small values in either one would be found where one catchment loses area to another (Figure 1), we next compared how minimum hillslope relief (HR min ), minimum flow distance (FD min ), and the divide asymmetry index (DAI) vary with divide distance in the Reference model and 185 the three models with landscape perturbances (Figure 9). In contrast to the average values in Figure 8, we provide these metrics for all divide edges during the last 1 Myr of the model runs. Note also that we plotted data points in an order which brings high DAI values to the front, to better assess where asymmetric divides are located, but that in all four models, relatively high DAI points may plot on top of low DAI points. In the Reference model, both minimum hillslope relief and minimum flow distance reach relatively steady values (HR min ~250-350m; FD min ~400-600m) at a divide distance of ~5 km. At lower divide distances, 190 both HR min and FD min approach zerosimply because divides are defined to start at the stream networkand these divides can become increasingly more asymmetric. It is notable however, that some of the highest HR min and FD min values are also observed at low divide distances of approximately 1-2 km. In the three other models, the transition between quite variable divides at short distances and more steady 'background' values at higher distances appears to be preserved, but we observe generally more variability. For example, in the Rotating model, we observe divides with significantly lower HR min and FD min 195 values at higher distances. These divides are particularly prominent at a distance of ~10 km and ~15-22 km and correspond to the position of the strike-slip fault. Where HR min and FD min are low, DAI values are relatively high (divides are highly asymmetric), although there also exist divides that have high DAI but regular HR min and FD min values. In the Inclined and Spheres models, the HR min and FD min values are never as low as in the Rotating model at divide distances >5 km, which reflects the lack of drainage captures. Instead, we observe frequent excursions to both higher and lower HR min and FD min values, either 200 across all divide distances (Inclined) or at specific locations (Spheres). Deviations from the average values in the Reference model are greatest in the Spheres model compared to all other perturbed models and always correspond to peaks that grow where strong spheres are exhumed. In the three disturbance models, DAI values are generally higher than in the Reference modelalthough divides with high HR min and FD min values and at great divide distances mostly appear to have somewhat lower DAI values. 205

Junction connectivity
The spatial pattern of divide junction connectivity (C J ) values at the end of the landscape evolution experiments (Figure 10) partly follows the pattern of divide distances (Figure 4). Junctions with higher C J values tend to occur at higher elevation and at greater divide distance (d d ). In the Reference model, the highest C J values occupy the centre as well as the centres of the four quadrants of the model domain, resembling the five of a six-sided dice. In the Rotating model, these clusters of high C J 210 values are maintained, but their connection is disrupted by the strike-slip fault, which induces low C J values. The centrally located divide junctions occupy a similar range in C J values, but are all shifted to higher elevations (cf., Figure 4). Similar, although lower, offsets to higher elevations occur in the Inclined and Spheres models, where junctions coincide with rocks of reduced erodibility. In the Spheres model, the basic structure of C J values is similar to that of the Reference model, but the highest C J values are steered towards the less erodible spheres, where they also attain C J values that are distinctly higher than 215 in any of the other models (Figure 10b). In general, divide junctions with combinations of elevation and C J values that are outside the range of values observed in the Reference model are found in the most disturbed parts of the landscape ( Figure   10c). In summary, the perturbed models appear to induce mostly changes in junction elevation, whereas changes in junction connectivity (C J ) are seemingly constrained to the Spheres model.

Quantifying drainage divide mobility
The analysis of stream networks has become a standard tool for inferring tectonic forcing and landscape history (e.g., Wobus Kirby and Whipple, 2012;Demoulin, 2012;Schwanghart and Scherler, 2014). The divide network holds the potential for recording similar, but also other aspects of landscape history (e.g., Willett et al., 2014). The question is, which divide metrics are useful to analyse and what they tell us about landscape history? Our Rotating model induced relatively 225 sudden drainage captures ( Figure 6). Because such events are associated with the dissection of drainage divides, reliable indicators are values of hillslope relief (HR) and flow distance (FD) that are much lower compared to the values that divides (> ~5 km divide distance) strive for (cf., Figure 9). More gradual divide migration, however, likely lacks such simple diagnostic criteria, and in those cases, across-divide differences in topographic metrics may be more suitable indicators of divide mobility.
The most commonly used metric to infer drainage divide mobility is the across-divide difference in χ (Willett et al., 2014). 230 Although the utility of this metric has recently received some critique (Whipple et al., 2017b;Forte and Whipple, 2018), it became a popular tool for studying drainage divides. Whipple et al. (2017b) and Forte and Whipple (2018) instead advocated the use of other topographic metrics, including mean gradient, mean local relief, and channel bed elevation, measured at or upstream of a reference drainage area. We note that these latter metrics are typically highly correlated and eventually very similar to the hillslope relief and DAI metrics that we included in this study. 235 Figure 11 shows how normalized across-divide differences in χ, hillslope relief (HR), as well as flow distance (FD) compare to normalized across-divide differences in erosion rate (ER), evaluated for each divide edge from the last 1 Myr of the landscape evolution experiments Reference, Rotating, Inclined, and Spheres. We find that across-divide differences in hillslope relief (HR; Figure 11a) are most sensitive to the disturbances included in the models, whereas across-divide differences in χ are similarly sensitive to disturbances in the Rotating model, but less so in the Inclined and Spheres models (Figure 11b; Figure 9 8). Across-divide differences in flow distance (FD) are the least sensitive to disturbances in the models and show the largest scatter when compared with erosion rates (Figure 11c). However, there also exists substantial scatter in the relationship between across-divide differences in hillslope relief and erosion rate, which partly depends on the divide distance. In general, we observe that the scatter is higher for divide distances <~5 km (dark blue in Figure 11), which corresponds to the value below which we observe large variability in divide morphology, even in the Reference model (cf., Figure 9). To quantify the 245 correlation of the normalized across-divide differences in topographic metrics with normalized across-divide differences in erosion rate, we fitted a linear model to all drainage divide edges from the entire model runs, categorized into 1-km divide distance bins and show the resulting coefficients of determination (R 2 ) in Figure 12. As already suspected from Figure 8 and Figure 11, the R 2 values differ in between models and metrics, and also depend on divide distance. In general, we observe that all metrics perform poorly at divide distances < ~5 km, and that across-divide differences in flow distance perform poorly even 250 at higher distances. The highest R 2 values are linked to across-divide differences in hillslope relief, whereas across-divide differences in χ attain similar R 2 values in the Rotating model, but some of the lowest R 2 values of all metrics in the Inclined and Spheres model. This difference may be explained by the fact that in the latter two models, we introduced spatial variability in the erosional efficiency of rivers (K r ) that we did not account for in our across-divide comparison of χ, as would be required (Willett et al., 2014). In natural landscapes, however, these values are rarely well known. 255 We speculate that the influence of divide distance on topographic metric-erosion rate relationships may also account for the differences in scatter observed by Sassolas-Serrayet et al. (2019) in landscape evolution experiments similar to our Initialize model between larger and smaller basin areas. But even when excluding divides of low order or low divide distance, we still observe considerable scatter in the topographic metric-erosion rate relationships, which, at the very least, demand caution when interpreting divide morphology in terms of mobility. In this regard, Figure 5 and studying the videos of the landscape evolution 260 experiments (see Supplementary Material) is insightful: where drainage divides are migrating, one typically observes a range of across-divide topographic metric values that vary considerably during the migration. In other words, despite a continuous divide migration at large scale, there often exists small-scale variability in divide morphology that may in part be related to across-divide differences in topographic metrics lagging behind across-divide differences in erosion rate.
As a final note, we emphasize that the above observations are from our numerical experiments, which depict an idealized 265 world. It is clear that complexities present in nature, such as anisotropic and variable rock properties, hydro-climatic gradients, or biological influences on erosion processes and rates can lead to landscape patterns which bias any of the above topographic metrics, and need to be taken into account when inferring divide dynamics from divide metrics in natural landscapes.

Divide network dynamics
Stream networks tend to attain configurations that are in equilibrium with the geological and climatic environment, given an 270 initial condition (e.g., Rinaldo et al., 2014). Because drainage divides are defined by adjacent drainage basins, the geometry of divide networks should attain a similar equilibrium, which expresses itself both in the geometry of divides and in the topology of divide networks. Our numerical experiments have shown that during the initial establishment of a stream network, on a relatively flat surface, both stream and divide networks are far from their steady-state configuration and characterized by networks that extend to high orders (Figure 7), and divide distances. During the subsequent extension and shrinkage of 275 individual streams towards their steady-state configuration, the divide network contracts and primarily high-order divide segments shorten and become fewer, whereas divides of low orders maintain their frequencies (Figure 7).
In general, divide segments of high order, i.e., at great distance from endpoints, appear to be the most responsive to landscape disturbances (Figure 8). In the case of the Rotating model, this is in part expected; because the inner rotating part of the landscape contains the highest order divide segments (Figure 4b). In the cases of the Inclined and Spheres models, it may be 280 related to the increased probability of recording a disturbance, because the adjoining basins cover a larger area compared to lower-order divides. In other words, if drainage captures happen somewhere within a drainage basin, this will most likely influence divides further upstream. Over a distance of less than ~5 km from divide network endpoints, the divide segments transition from low interfluves at river junctions to high topographic ridges, as seen in the Reference model (Figure 9). In the other models, most of the investigated morphometric parameters are quite variable over the same distance and can be seen to 285 rapidly adjust to disturbances such as drainage captures or migrating divides. Such behaviour is consistent with the observation that the timescale of a rivers' response to changes in drainage area increases with the distance from the divide to the outlet of a river (Whipple et al., 2017b). To reliably distinguish the morphologic effects of real disturbances from 'noise' close to the river, a minimum divide distance of perhaps ~5 km, as in our analysis of the Big Tujunga divide network (Scherler and Schwanghart, submitted), appears appropriate. This minimum divide distance could be lower or higher, depending on factors 290 like drainage density and average hillslope relief, for example.
Our new junction connectivity index (C J ), complements existing topographic metrics in assessing divide network dynamics.
For example, the junction connectivity in our Rotating model is low along the fault (Figure 10), consistence with the absence of stable divides. In the Spheres model, however, the appearance of more resistant rocks at the surface often resulted in migration of divides towards the spheres (Figure 5c, Supplementary Material). In this case, parts of the drainage divide network 295 were mobile, not stable, but they moved towards particularly stable portions in the landscape. Therefore, the junction connectivity index (C J ) may also be interpreted as an attractor or centrality index (Phillips et al., 2015) that quantifies how strong a drainage divide network has been pulled towards and anchored at a certain junction (Spotila, 2012).

Conclusions
Based on landscape evolution model experiments in which we forced divides to migrate, we found that stable drainage divides 300 strive to attain a constant hillslope relief as well as flow distance from the nearest stream, provided a sufficiently large divide distance to avoid confounding influences near the edges of the divide network. In our experiments this distance is ~5 km from endpoints. Simple indicators of mobile divides are thus anomalously low hillslope relief or flow distance values, which could signal beheaded valleys or future capture events. Overall, drainage divides located high up in the network, i.e., at great distance from endpoints, are more vulnerable than divides closer to endpoints of the network and are more likely to record disturbance      to junctions in the perturbed models that lie higher than the junctions in the Reference models. Figure 11: Relationship between across-divide differences in the topographic metrics hillslope relief (HR), chi (χ), and flow distance (FD), with across-divide differences in erosion rates (ER) for all divide edges in the last 1 Myr of the landscape evolution experiments Reference, Rotating, Inclined, and Spheres. Marker colour denotes the divide distance (dd). Note that the markers are sorted by their divide distance to 495 show the highest divide distance above data points with lower divide distance, and that data points with lower divide distance occur beneath those with higher divide distance. Figure 12: Coefficient of determination for linear regressions between normalized across-divide differences in different topographic metrics (HR = hillslope relief; FD = flow distance; χ = chi) and normalized across-divide differences in erosion rates, as a function of drainage divide distance in the three perturbation models.